source: git/libpolys/coeffs/longrat.cc @ 558f3cc

spielwiese
Last change on this file since 558f3cc was 558f3cc, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fix: make positive leading coeff at normalizing!
  • Property mode set to 100644
File size: 59.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
6*/
7
8#include "config.h"
9#include <misc/auxiliary.h>
10
11#ifdef HAVE_FACTORY
12#include <factory/factory.h>
13#endif
14
15#include <coeffs/longrat.h>
16
17
18// 64 bit version:
19//#if SIZEOF_LONG == 8
20#if 0
21#define MAX_NUM_SIZE 60
22#define POW_2_28 (1L<<60)
23#define LONG long
24#else
25#define MAX_NUM_SIZE 28
26#define POW_2_28 (1L<<28)
27#define LONG int
28#endif
29
30static inline number nlShort3(number x) // assume x->s==3
31{
32  assume(x->s==3);
33  if (mpz_cmp_ui(x->z,(long)0)==0)
34  {
35    mpz_clear(x->z);
36    FREE_RNUMBER(x);
37    return INT_TO_SR(0);
38  }
39  if (mpz_size1(x->z)<=MP_SMALL)
40  {
41    LONG ui=mpz_get_si(x->z);
42    if ((((ui<<3)>>3)==ui)
43    && (mpz_cmp_si(x->z,(long)ui)==0))
44    {
45      mpz_clear(x->z);
46      FREE_RNUMBER(x);
47      return INT_TO_SR(ui);
48    }
49  }
50  return x;
51}
52
53#ifndef LONGRAT_CC
54#define LONGRAT_CC
55
56#include <string.h>
57#include <float.h>
58
59#include <coeffs/coeffs.h>
60#include <reporter/reporter.h>
61#include <omalloc/omalloc.h>
62
63#include <coeffs/numbers.h>
64#include <coeffs/modulop.h>
65#include <coeffs/shortfl.h>
66#include <coeffs/mpr_complex.h>
67
68#ifndef BYTES_PER_MP_LIMB
69#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
70#endif
71
72//#define SR_HDL(A) ((long)(A))
73/*#define SR_INT    1L*/
74/*#define INT_TO_SR(INT)  ((number) (((long)INT << 2) + SR_INT))*/
75// #define SR_TO_INT(SR)   (((long)SR) >> 2)
76
77#define MP_SMALL 1
78//#define mpz_isNeg(A) (mpz_cmp_si(A,(long)0)<0)
79#define mpz_isNeg(A) ((A)->_mp_size<0)
80#define mpz_limb_size(A) ((A)->_mp_size)
81#define mpz_limb_d(A) ((A)->_mp_d)
82#define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
83#define MPZ_EXACTDIV(A,B,C) mpz_divexact((A),(B),(C))
84
85/// Our Type!
86static const n_coeffType ID = n_Q;
87
88void    _nlDelete_NoImm(number *a);
89
90/***************************************************************
91 *
92 * Routines which are never inlined by p_Numbers.h
93 *
94 *******************************************************************/
95#ifndef P_NUMBERS_H
96
97number nlShort3_noinline(number x) // assume x->s==3
98{
99  return nlShort3(x);
100}
101
102
103number nlOne=INT_TO_SR(1);
104
105#if (__GNU_MP_VERSION*10+__GNU_MP_VERSION_MINOR < 31)
106void mpz_mul_si (mpz_ptr r, mpz_srcptr s, long int si)
107{
108  if (si>=0)
109    mpz_mul_ui(r,s,si);
110  else
111  {
112    mpz_mul_ui(r,s,-si);
113    mpz_neg(r,r);
114  }
115}
116#endif
117
118static number nlMapP(number from, const coeffs src, const coeffs dst)
119{
120  assume( getCoeffType(dst) == ID );
121  assume( getCoeffType(src) == n_Zp );
122
123  number to;
124  to = nlInit(npInt(from,src), dst);
125  return to;
126}
127
128static number nlMapLongR(number from, const coeffs src, const coeffs dst);
129static number nlMapR(number from, const coeffs src, const coeffs dst);
130
131
132#ifdef HAVE_RINGS
133/*2
134* convert from a GMP integer
135*/
136number nlMapGMP(number from, const coeffs /*src*/, const coeffs /*dst*/)
137{
138  number z=ALLOC_RNUMBER();
139#if defined(LDEBUG)
140  z->debug=123456;
141#endif
142  mpz_init_set(z->z,(mpz_ptr) from);
143  //mpz_init_set_ui(&z->n,1);
144  z->s = 3;
145  z=nlShort3(z);
146  return z;
147}
148
149/*2
150* convert from an machine long
151*/
152number nlMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
153{
154  number z=ALLOC_RNUMBER();
155#if defined(LDEBUG)
156  z->debug=123456;
157#endif
158  mpz_init_set_ui(z->z,(unsigned long) from);
159  z->s = 3;
160  z=nlShort3(z);
161  return z;
162}
163#endif
164
165
166#ifdef LDEBUG
167BOOLEAN nlDBTest(number a, const char *f,const int l, const coeffs /*r*/)
168{
169  if (a==NULL)
170  {
171    Print("!!longrat: NULL in %s:%d\n",f,l);
172    return FALSE;
173  }
174  //if ((int)a==1) Print("!! 0x1 as number ? %s %d\n",f,l);
175  if ((((long)a)&3L)==3L)
176  {
177    Print(" !!longrat:ptr(3) in %s:%d\n",f,l);
178    return FALSE;
179  }
180  if ((((long)a)&3L)==1L)
181  {
182    if (((((LONG)(long)a)<<1)>>1)!=((LONG)(long)a))
183    {
184      Print(" !!longrat:arith:%lx in %s:%d\n",(long)a, f,l);
185      return FALSE;
186    }
187    return TRUE;
188  }
189  /* TODO: If next line is active, then computations in algebraic field
190           extensions over Q will throw a lot of assume violations although
191           everything is computed correctly and no seg fault appears.
192           Maybe the test is not appropriate in this case. */
193  omCheckIf(omCheckAddrSize(a,sizeof(*a)), return FALSE);
194  if (a->debug!=123456)
195  {
196    Print("!!longrat:debug:%d in %s:%d\n",a->debug,f,l);
197    a->debug=123456;
198    return FALSE;
199  }
200  if ((a->s<0)||(a->s>4))
201  {
202    Print("!!longrat:s=%d in %s:%d\n",a->s,f,l);
203    return FALSE;
204  }
205  /* TODO: If next line is active, then computations in algebraic field
206           extensions over Q will throw a lot of assume violations although
207           everything is computed correctly and no seg fault appears.
208           Maybe the test is not appropriate in this case. */
209  //omCheckAddrSize(a->z[0]._mp_d,a->z[0]._mp_alloc*BYTES_PER_MP_LIMB);
210  if (a->z[0]._mp_alloc==0)
211    Print("!!longrat:z->alloc=0 in %s:%d\n",f,l);
212
213  if (a->s<2)
214  {
215    /* TODO: If next line is active, then computations in algebraic field
216             extensions over Q will throw a lot of assume violations although
217             everything is computed correctly and no seg fault appears.
218             Maybe the test is not appropriate in this case. */
219    //omCheckIf(omCheckAddrSize(a->n[0]._mp_d,a->n[0]._mp_alloc*BYTES_PER_MP_LIMB), return FALSE);
220    if (a->z[0]._mp_alloc==0)
221      Print("!!longrat:n->alloc=0 in %s:%d\n",f,l);
222    if ((mpz_size1(a->n) ==1) && (mpz_cmp_si(a->n,(long)1)==0))
223    {
224      Print("!!longrat:integer as rational in %s:%d\n",f,l);
225      mpz_clear(a->n); a->s=3;
226      return FALSE;
227    }
228    else if (mpz_isNeg(a->n))
229    {
230      Print("!!longrat:div. by negative in %s:%d\n",f,l);
231      mpz_neg(a->z,a->z);
232      mpz_neg(a->n,a->n);
233      return FALSE;
234    }
235    return TRUE;
236  }
237  //if (a->s==2)
238  //{
239  //  Print("!!longrat:s=2 in %s:%d\n",f,l);
240  //  return FALSE;
241  //}
242  if (mpz_size1(a->z)>MP_SMALL) return TRUE;
243  LONG ui=(int)mpz_get_si(a->z);
244  if ((((ui<<3)>>3)==ui)
245  && (mpz_cmp_si(a->z,(long)ui)==0))
246  {
247    Print("!!longrat:im int %d in %s:%d\n",ui,f,l);
248    f=NULL;
249    return FALSE;
250  }
251  return TRUE;
252}
253#endif
254
255#ifdef HAVE_FACTORY
256CanonicalForm nlConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs /*r*/ )
257{
258  if (setChar) setCharacteristic( 0 );
259
260  CanonicalForm term;
261  if ( SR_HDL(n) & SR_INT )
262  {
263    term = SR_TO_INT(n);
264  }
265  else
266  {
267    if ( n->s == 3 )
268    {
269      mpz_t dummy;
270      long lz=mpz_get_si(n->z);
271      if (mpz_cmp_si(n->z,lz)==0) term=lz;
272      else
273      {
274        mpz_init_set( dummy,n->z );
275        term = make_cf( dummy );
276      }
277    }
278    else
279    {
280      // assume s==0 or s==1
281      mpz_t num, den;
282      On(SW_RATIONAL);
283      mpz_init_set( num, n->z );
284      mpz_init_set( den, n->n );
285      term = make_cf( num, den, ( n->s != 1 ));
286    }
287  }
288  return term;
289}
290
291number nlConvFactoryNSingN( const CanonicalForm n, const coeffs r)
292{
293  if (n.isImm())
294  {
295    long lz=n.intval();
296    int iz=(int)lz;
297    if ((long)iz==lz)
298    {
299      return nlInit(n.intval(),r);
300    }
301    else  return nlRInit(lz);
302    return nlInit(n.intval(),r);
303  }
304  else
305  {
306    number z = ALLOC_RNUMBER(); //Q!? // (number)omAllocBin(rnumber_bin);
307#if defined(LDEBUG)
308    z->debug=123456;
309#endif
310    gmp_numerator( n, z->z );
311    if ( n.den().isOne() )
312      z->s = 3;
313    else
314    {
315      gmp_denominator( n, z->n );
316      z->s = 0;
317    }
318    nlNormalize(z,r);
319    return z;
320  }
321}
322#endif
323number nlRInit (long i);
324
325static number nlMapR(number from, const coeffs src, const coeffs dst)
326{
327  assume( getCoeffType(dst) == ID );
328  assume( getCoeffType(src) == n_R );
329
330  double f=nrFloat(from);
331  if (f==0.0) return INT_TO_SR(0);
332  int f_sign=1;
333  if (f<0.0)
334  {
335    f_sign=-1;
336    f=-f;
337  }
338  int i=0;
339  mpz_t h1;
340  mpz_init_set_ui(h1,1);
341  while((FLT_RADIX*f) < DBL_MAX && i<DBL_MANT_DIG)
342  {
343    f*=FLT_RADIX;
344    mpz_mul_ui(h1,h1,FLT_RADIX);
345    i++;
346  }
347  number re=nlRInit(1);
348  mpz_set_d(re->z,f);
349  memcpy(&(re->n),&h1,sizeof(h1));
350  re->s=0; /* not normalized */
351  if(f_sign==-1) re=nlNeg(re,dst);
352  nlNormalize(re,dst);
353  return re;
354}
355
356static number nlMapLongR(number from, const coeffs src, const coeffs dst)
357{
358  assume( getCoeffType(dst) == ID );
359  assume( getCoeffType(src) == n_long_R );
360
361  gmp_float *ff=(gmp_float*)from;
362  mpf_t *f=ff->_mpfp();
363  number res;
364  mpz_ptr dest,ndest;
365  int size, i,negative;
366  int e,al,bl;
367  mp_ptr qp,dd,nn;
368
369  size = (*f)[0]._mp_size;
370  if (size == 0)
371    return INT_TO_SR(0);
372  if(size<0)
373  {
374    negative = 1;
375    size = -size;
376  }
377  else
378    negative = 0;
379
380  qp = (*f)[0]._mp_d;
381  while(qp[0]==0)
382  {
383    qp++;
384    size--;
385  }
386
387  e=(*f)[0]._mp_exp-size;
388  res = ALLOC_RNUMBER();
389#if defined(LDEBUG)
390  res->debug=123456;
391#endif
392  dest = res->z;
393
394  if (e<0)
395  {
396    al = dest->_mp_size = size;
397    if (al<2) al = 2;
398    dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al);
399    for (i=0;i<size;i++) dd[i] = qp[i];
400    bl = 1-e;
401    nn = (mp_ptr)omAlloc(sizeof(mp_limb_t)*bl);
402    nn[bl-1] = 1;
403    for (i=bl-2;i>=0;i--) nn[i] = 0;
404    ndest = res->n;
405    ndest->_mp_d = nn;
406    ndest->_mp_alloc = ndest->_mp_size = bl;
407    res->s = 0;
408  }
409  else
410  {
411    al = dest->_mp_size = size+e;
412    if (al<2) al = 2;
413    dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al);
414    for (i=0;i<size;i++) dd[i+e] = qp[i];
415    for (i=0;i<e;i++) dd[i] = 0;
416    res->s = 3;
417  }
418
419  dest->_mp_d = dd;
420  dest->_mp_alloc = al;
421  if (negative) dest->_mp_size = -dest->_mp_size;
422
423  if (res->s==0)
424    nlNormalize(res,dst);
425  else if (mpz_size1(res->z)<=MP_SMALL)
426  {
427    // res is new, res->ref is 1
428    res=nlShort3(res);
429  }
430  nlTest(res, dst);
431  return res;
432}
433
434//static number nlMapLongR(number from)
435//{
436//  gmp_float *ff=(gmp_float*)from;
437//  const mpf_t *f=ff->mpfp();
438//  int f_size=ABS((*f)[0]._mp_size);
439//  if (f_size==0)
440//    return nlInit(0);
441//  int f_sign=1;
442//  number work=ngcCopy(from);
443//  if (!ngcGreaterZero(work))
444//  {
445//    f_sign=-1;
446//    work=ngcNeg(work);
447//  }
448//  int i=0;
449//  mpz_t h1;
450//  mpz_init_set_ui(h1,1);
451//  while((FLT_RADIX*f) < DBL_MAX && i<DBL_MANT_DIG)
452//  {
453//    f*=FLT_RADIX;
454//    mpz_mul_ui(h1,h1,FLT_RADIX);
455//    i++;
456//  }
457//  number r=nlRInit(1);
458//  mpz_set_d(&(r->z),f);
459//  memcpy(&(r->n),&h1,sizeof(h1));
460//  r->s=0; /* not normalized */
461//  nlNormalize(r);
462//  return r;
463//
464//
465//  number r=nlRInit(1);
466//  int f_shift=f_size+(*f)[0]._mp_exp;
467//  if ( f_shift > 0)
468//  {
469//    r->s=0;
470//    mpz_init(&r->n);
471//    mpz_setbit(&r->n,f_shift*BYTES_PER_MP_LIMB*8);
472//    mpz_setbit(&r->z,f_size*BYTES_PER_MP_LIMB*8-1);
473//    // now r->z has enough space
474//    memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
475//    nlNormalize(r);
476//  }
477//  else
478//  {
479//    r->s=3;
480//    if (f_shift==0)
481//    {
482//      mpz_setbit(&r->z,f_size*BYTES_PER_MP_LIMB*8-1);
483//      // now r->z has enough space
484//      memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
485//    }
486//    else /* f_shift < 0 */
487//    {
488//      mpz_setbit(&r->z,(f_size-f_shift)*BYTES_PER_MP_LIMB*8-1);
489//      // now r->z has enough space
490//      memcpy(mpz_limb_d(&r->z)-f_shift,((*f)[0]._mp_d),
491//        f_size*BYTES_PER_MP_LIMB);
492//    }
493//  }
494//  if ((*f)[0]._mp_size<0);
495//    r=nlNeg(r);
496//  return r;
497//}
498
499int nlSize(number a, const coeffs)
500{
501  if (a==INT_TO_SR(0))
502     return 0; /* rational 0*/
503  if (SR_HDL(a) & SR_INT)
504     return 1; /* immidiate int */
505  int s=a->z[0]._mp_alloc;
506//  while ((s>0) &&(a->z._mp_d[s]==0L)) s--;
507//#if SIZEOF_LONG == 8
508//  if (a->z._mp_d[s] < (unsigned long)0x100000000L) s=s*2-1;
509//  else s *=2;
510//#endif
511//  s++;
512  if (a->s<2)
513  {
514    int d=a->n[0]._mp_alloc;
515//    while ((d>0) && (a->n._mp_d[d]==0L)) d--;
516//#if SIZEOF_LONG == 8
517//    if (a->n._mp_d[d] < (unsigned long)0x100000000L) d=d*2-1;
518//    else d *=2;
519//#endif
520    s+=d;
521  }
522  return s;
523}
524
525/*2
526* convert number to int
527*/
528int nlInt(number &i, const coeffs r)
529{
530  nlTest(i, r);
531  nlNormalize(i,r);
532  if (SR_HDL(i) &SR_INT) return SR_TO_INT(i);
533  if (i->s==3)
534  {
535    if(mpz_size1(i->z)>MP_SMALL) return 0;
536    int ul=(int)mpz_get_si(i->z);
537    if (mpz_cmp_si(i->z,(long)ul)!=0) return 0;
538    return ul;
539  }
540  mpz_t tmp;
541  int ul;
542  mpz_init(tmp);
543  MPZ_DIV(tmp,i->z,i->n);
544  if(mpz_size1(tmp)>MP_SMALL) ul=0;
545  else
546  {
547    ul=(int)mpz_get_si(tmp);
548    if (mpz_cmp_si(tmp,(long)ul)!=0) ul=0;
549  }
550  mpz_clear(tmp);
551  return ul;
552}
553
554/*2
555* convert number to bigint
556*/
557number nlBigInt(number &i, const coeffs r)
558{
559  nlTest(i, r);
560  nlNormalize(i,r);
561  if (SR_HDL(i) &SR_INT) return (i);
562  if (i->s==3)
563  {
564    return nlCopy(i,r);
565  }
566  number tmp=nlRInit(1);
567  MPZ_DIV(tmp->z,i->z,i->n);
568  tmp=nlShort3(tmp);
569  return tmp;
570}
571
572/*
573* 1/a
574*/
575number nlInvers(number a, const coeffs r)
576{
577  nlTest(a, r);
578  number n;
579  if (SR_HDL(a) & SR_INT)
580  {
581    if ((a==INT_TO_SR(1L)) || (a==INT_TO_SR(-1L)))
582    {
583      return a;
584    }
585    if (nlIsZero(a,r))
586    {
587      WerrorS(nDivBy0);
588      return INT_TO_SR(0);
589    }
590    n=ALLOC_RNUMBER();
591#if defined(LDEBUG)
592    n->debug=123456;
593#endif
594    n->s=1;
595    if ((long)a>0L)
596    {
597      mpz_init_set_si(n->z,(long)1);
598      mpz_init_set_si(n->n,(long)SR_TO_INT(a));
599    }
600    else
601    {
602      mpz_init_set_si(n->z,(long)-1);
603      mpz_init_set_si(n->n,(long)-SR_TO_INT(a));
604    }
605    nlTest(n, r);
606    return n;
607  }
608  n=ALLOC_RNUMBER();
609#if defined(LDEBUG)
610  n->debug=123456;
611#endif
612  {
613    n->s=a->s;
614    mpz_init_set(n->n,a->z);
615    switch (a->s)
616    {
617      case 0:
618      case 1:
619              mpz_init_set(n->z,a->n);
620              if (mpz_isNeg(n->n)) /* && n->s<2*/
621              {
622                mpz_neg(n->z,n->z);
623                mpz_neg(n->n,n->n);
624              }
625              if (mpz_cmp_si(n->n,(long)1)==0)
626              {
627                mpz_clear(n->n);
628                n->s=3;
629                n=nlShort3(n);
630              }
631              break;
632      case 3:
633              n->s=1;
634              if (mpz_isNeg(n->n)) /* && n->s<2*/
635              {
636                mpz_neg(n->n,n->n);
637                mpz_init_set_si(n->z,(long)-1);
638              }
639              else
640              {
641                mpz_init_set_si(n->z,(long)1);
642              }
643              break;
644    }
645  }
646  nlTest(n, r);
647  return n;
648}
649
650
651/*2
652* u := a / b in Z, if b | a (else undefined)
653*/
654number   nlExactDiv(number a, number b, const coeffs r)
655{
656  if (b==INT_TO_SR(0))
657  {
658    WerrorS(nDivBy0);
659    return INT_TO_SR(0);
660  }
661  if (a==INT_TO_SR(0))
662    return INT_TO_SR(0);
663  number u;
664  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
665  {
666    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
667    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
668    {
669      return nlRInit(POW_2_28);
670    }
671    long aa=SR_TO_INT(a);
672    long bb=SR_TO_INT(b);
673    return INT_TO_SR(aa/bb);
674  }
675  number bb=NULL;
676  if (SR_HDL(b) & SR_INT)
677  {
678    bb=nlRInit(SR_TO_INT(b));
679    b=bb;
680  }
681  u=ALLOC_RNUMBER();
682#if defined(LDEBUG)
683  u->debug=123456;
684#endif
685  mpz_init(u->z);
686  /* u=a/b */
687  u->s = 3;
688  MPZ_EXACTDIV(u->z,a->z,b->z);
689  if (bb!=NULL)
690  {
691    mpz_clear(bb->z);
692#if defined(LDEBUG)
693    bb->debug=654324;
694#endif
695    FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
696  }
697  u=nlShort3(u);
698  nlTest(u, r);
699  return u;
700}
701
702/*2
703* u := a / b in Z
704*/
705number nlIntDiv (number a, number b, const coeffs r)
706{
707  if (b==INT_TO_SR(0))
708  {
709    WerrorS(nDivBy0);
710    return INT_TO_SR(0);
711  }
712  if (a==INT_TO_SR(0))
713    return INT_TO_SR(0);
714  number u;
715  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
716  {
717    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
718    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
719    {
720      return nlRInit(POW_2_28);
721    }
722    long aa=SR_TO_INT(a);
723    long bb=SR_TO_INT(b);
724    //return INT_TO_SR((aa-(aa%bb))/bb);
725    return INT_TO_SR(aa/bb);
726  }
727  if (SR_HDL(a) & SR_INT)
728  {
729    /* the small int -(1<<28) divided by 2^28 is 1   */
730    if (a==INT_TO_SR(-(POW_2_28)))
731    {
732      if(mpz_cmp_si(b->z,(POW_2_28))==0)
733      {
734        return INT_TO_SR(-1);
735      }
736    }
737    /* a is a small and b is a large int: -> 0 */
738    return INT_TO_SR(0);
739  }
740  number bb=NULL;
741  if (SR_HDL(b) & SR_INT)
742  {
743    bb=nlRInit(SR_TO_INT(b));
744    b=bb;
745  }
746  u=ALLOC_RNUMBER();
747#if defined(LDEBUG)
748  u->debug=123456;
749#endif
750  assume(a->s==3);
751  assume(b->s==3);
752  mpz_init_set(u->z,a->z);
753  /* u=u/b */
754  u->s = 3;
755  MPZ_DIV(u->z,u->z,b->z);
756  if (bb!=NULL)
757  {
758    mpz_clear(bb->z);
759#if defined(LDEBUG)
760    bb->debug=654324;
761#endif
762    FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
763  }
764  u=nlShort3(u);
765  nlTest(u, r);
766  return u;
767}
768
769/*2
770* u := a mod b in Z, u>=0
771*/
772number nlIntMod (number a, number b, const coeffs r)
773{
774  if (b==INT_TO_SR(0))
775  {
776    WerrorS(nDivBy0);
777    return INT_TO_SR(0);
778  }
779  if (a==INT_TO_SR(0))
780    return INT_TO_SR(0);
781  number u;
782  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
783  {
784    //return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
785    if ((long)a>0L)
786    {
787      if ((long)b>0L)
788        return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
789      else
790        return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
791    }
792    else
793    {
794      if ((long)b>0L)
795      {
796        long i=(-SR_TO_INT(a))%SR_TO_INT(b);
797        if ( i != 0L ) i = (SR_TO_INT(b))-i;
798        return INT_TO_SR(i);
799      }
800      else
801      {
802        long i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
803        if ( i != 0L ) i = (-SR_TO_INT(b))-i;
804        return INT_TO_SR(i);
805      }
806    }
807  }
808  if (SR_HDL(a) & SR_INT)
809  {
810    /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
811    if ((long)a<0L)
812    {
813      if (mpz_isNeg(b->z))
814        return nlSub(a,b,r);
815      /*else*/
816        return nlAdd(a,b,r);
817    }
818    /*else*/
819      return a;
820
821    /* a is a small and b is a large int: -> a */
822    //return a;
823  }
824  number bb=NULL;
825  if (SR_HDL(b) & SR_INT)
826  {
827    bb=nlRInit(SR_TO_INT(b));
828    b=bb;
829  }
830  u=ALLOC_RNUMBER();
831#if defined(LDEBUG)
832  u->debug=123456;
833#endif
834  mpz_init(u->z);
835  u->s = 3;
836  mpz_mod(u->z,a->z,b->z);
837  if (bb!=NULL)
838  {
839    mpz_clear(bb->z);
840#if defined(LDEBUG)
841    bb->debug=654324;
842#endif
843    FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
844  }
845  if (mpz_isNeg(u->z))
846  {
847    if (mpz_isNeg(b->z))
848      mpz_sub(u->z,u->z,b->z);
849    else
850      mpz_add(u->z,u->z,b->z);
851  }
852  u=nlShort3(u);
853  nlTest(u, r);
854  return u;
855}
856
857/*2
858* u := a / b
859*/
860number nlDiv (number a, number b, const coeffs r)
861{
862  number u;
863  if (nlIsZero(b,r))
864  {
865    WerrorS(nDivBy0);
866    return INT_TO_SR(0);
867  }
868  u=ALLOC_RNUMBER();
869  u->s=0;
870#if defined(LDEBUG)
871  u->debug=123456;
872#endif
873// ---------- short / short ------------------------------------
874  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
875  {
876    LONG i=SR_TO_INT(a);
877    LONG j=SR_TO_INT(b);
878    if ((i==-POW_2_28) && (j== -1L))
879    {
880      FREE_RNUMBER(u);
881      return nlRInit(POW_2_28);
882    }
883    LONG r=i%j;
884    if (r==0)
885    {
886      FREE_RNUMBER(u); // omFreeBin((void *)u, rnumber_bin);
887      return INT_TO_SR(i/j);
888    }
889    mpz_init_set_si(u->z,(long)i);
890    mpz_init_set_si(u->n,(long)j);
891  }
892  else
893  {
894    mpz_init(u->z);
895// ---------- short / long ------------------------------------
896    if (SR_HDL(a) & SR_INT)
897    {
898      // short a / (z/n) -> (a*n)/z
899      if (b->s<2)
900      {
901        mpz_mul_si(u->z,b->n,SR_TO_INT(a));
902      }
903      else
904      // short a / long z -> a/z
905      {
906        mpz_set_si(u->z,SR_TO_INT(a));
907      }
908      if (mpz_cmp(u->z,b->z)==0)
909      {
910        mpz_clear(u->z);
911        FREE_RNUMBER(u);
912        return INT_TO_SR(1);
913      }
914      mpz_init_set(u->n,b->z);
915    }
916// ---------- long / short ------------------------------------
917    else if (SR_HDL(b) & SR_INT)
918    {
919      mpz_set(u->z,a->z);
920      // (z/n) / b -> z/(n*b)
921      if (a->s<2)
922      {
923        mpz_init_set(u->n,a->n);
924        if ((long)b>0L)
925          mpz_mul_ui(u->n,u->n,SR_TO_INT(b));
926        else
927        {
928          mpz_mul_ui(u->n,u->n,-SR_TO_INT(b));
929          mpz_neg(u->z,u->z);
930        }
931      }
932      else
933      // long z / short b -> z/b
934      {
935        //mpz_set(u->z,a->z);
936        mpz_init_set_si(u->n,SR_TO_INT(b));
937      }
938    }
939// ---------- long / long ------------------------------------
940    else
941    {
942      mpz_set(u->z,a->z);
943      mpz_init_set(u->n,b->z);
944      if (a->s<2) mpz_mul(u->n,u->n,a->n);
945      if (b->s<2) mpz_mul(u->z,u->z,b->n);
946    }
947  }
948  if (mpz_isNeg(u->n))
949  {
950    mpz_neg(u->z,u->z);
951    mpz_neg(u->n,u->n);
952  }
953  if (mpz_cmp_si(u->n,(long)1)==0)
954  {
955    mpz_clear(u->n);
956    u->s=3;
957    u=nlShort3(u);
958  }
959  nlTest(u, r);
960  return u;
961}
962
963/*2
964* u:= x ^ exp
965*/
966void nlPower (number x,int exp,number * u, const coeffs r)
967{
968  *u = INT_TO_SR(0); // 0^e, e!=0
969  if (!nlIsZero(x,r))
970  {
971    nlTest(x, r);
972    number aa=NULL;
973    if (SR_HDL(x) & SR_INT)
974    {
975      aa=nlRInit(SR_TO_INT(x));
976      x=aa;
977    }
978    else if (x->s==0)
979      nlNormalize(x,r);
980    *u=ALLOC_RNUMBER();
981#if defined(LDEBUG)
982    (*u)->debug=123456;
983#endif
984    mpz_init((*u)->z);
985    mpz_pow_ui((*u)->z,x->z,(unsigned long)exp);
986    if (x->s<2)
987    {
988      if (mpz_cmp_si(x->n,(long)1)==0)
989      {
990        x->s=3;
991        mpz_clear(x->n);
992      }
993      else
994      {
995        mpz_init((*u)->n);
996        mpz_pow_ui((*u)->n,x->n,(unsigned long)exp);
997      }
998    }
999    (*u)->s = x->s;
1000    if ((*u)->s==3) *u=nlShort3(*u);
1001    if (aa!=NULL)
1002    {
1003      mpz_clear(aa->z);
1004      FREE_RNUMBER(aa);
1005    }
1006  }
1007  else if (exp==0)
1008    *u = INT_TO_SR(1); // 0^0
1009#ifdef LDEBUG
1010  if (exp<0) Print("nlPower: neg. exp. %d\n",exp);
1011  nlTest(*u, r);
1012#endif
1013}
1014
1015
1016/*2
1017* za >= 0 ?
1018*/
1019BOOLEAN nlGreaterZero (number a, const coeffs r)
1020{
1021  nlTest(a, r);
1022  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>1L /* represents number(0) */;
1023  return (!mpz_isNeg(a->z));
1024}
1025
1026/*2
1027* a > b ?
1028*/
1029BOOLEAN nlGreater (number a, number b, const coeffs r)
1030{
1031  nlTest(a, r);
1032  nlTest(b, r);
1033  number re;
1034  BOOLEAN rr;
1035  re=nlSub(a,b,r);
1036  rr=(!nlIsZero(re,r)) && (nlGreaterZero(re,r));
1037  nlDelete(&re,r);
1038  return rr;
1039}
1040
1041/*2
1042* a == -1 ?
1043*/
1044BOOLEAN nlIsMOne (number a, const coeffs r)
1045{
1046#ifdef LDEBUG
1047  if (a==NULL) return FALSE;
1048  nlTest(a, r);
1049#endif
1050  //if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1L));
1051  //return FALSE;
1052  return  (a==INT_TO_SR(-1L));
1053}
1054
1055/*2
1056* result =gcd(a,b)
1057*/
1058number nlGcd(number a, number b, const coeffs r)
1059{
1060  number result;
1061  nlTest(a, r);
1062  nlTest(b, r);
1063  //nlNormalize(a);
1064  //nlNormalize(b);
1065  if ((a==INT_TO_SR(1L))||(a==INT_TO_SR(-1L))
1066  ||  (b==INT_TO_SR(1L))||(b==INT_TO_SR(-1L)))
1067    return INT_TO_SR(1L);
1068  if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
1069    return nlCopy(b,r);
1070  if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
1071    return nlCopy(a,r);
1072  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1073  {
1074    long i=SR_TO_INT(a);
1075    long j=SR_TO_INT(b);
1076    if((i==0L)||(j==0L))
1077      return INT_TO_SR(1);
1078    long l;
1079    i=ABS(i);
1080    j=ABS(j);
1081    do
1082    {
1083      l=i%j;
1084      i=j;
1085      j=l;
1086    } while (l!=0L);
1087    if (i==POW_2_28)
1088      result=nlRInit(POW_2_28);
1089    else
1090     result=INT_TO_SR(i);
1091    nlTest(result,r);
1092    return result;
1093  }
1094  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1095  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1096  if (SR_HDL(a) & SR_INT)
1097  {
1098    LONG aa=ABS(SR_TO_INT(a));
1099    unsigned long t=mpz_gcd_ui(NULL,b->z,(long)aa);
1100    if (t==POW_2_28)
1101      result=nlRInit(POW_2_28);
1102    else
1103     result=INT_TO_SR(t);
1104  }
1105  else
1106  if (SR_HDL(b) & SR_INT)
1107  {
1108    LONG bb=ABS(SR_TO_INT(b));
1109    unsigned long t=mpz_gcd_ui(NULL,a->z,(long)bb);
1110    if (t==POW_2_28)
1111      result=nlRInit(POW_2_28);
1112    else
1113     result=INT_TO_SR(t);
1114  }
1115  else
1116  {
1117    result=ALLOC_RNUMBER();
1118    mpz_init(result->z);
1119    mpz_gcd(result->z,a->z,b->z);
1120    result->s = 3;
1121  #ifdef LDEBUG
1122    result->debug=123456;
1123  #endif
1124    result=nlShort3(result);
1125  }
1126  nlTest(result, r);
1127  return result;
1128}
1129
1130number nlShort1(number x) // assume x->s==0/1
1131{
1132  assume(x->s<2);
1133  if (mpz_cmp_ui(x->z,(long)0)==0)
1134  {
1135    _nlDelete_NoImm(&x);
1136    return INT_TO_SR(0);
1137  }
1138  if (x->s<2)
1139  {
1140    if (mpz_cmp(x->z,x->n)==0)
1141    {
1142      _nlDelete_NoImm(&x);
1143      return INT_TO_SR(1);
1144    }
1145  }
1146  return x;
1147}
1148/*2
1149* simplify x
1150*/
1151void nlNormalize (number &x, const coeffs r)
1152{
1153  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1154    return;
1155  if (x->s==3)
1156  {
1157    x=nlShort3_noinline(x);
1158    nlTest(x,r);
1159    return;
1160  }
1161  else if (x->s==0)
1162  {
1163    if (mpz_cmp_si(x->n,(long)1)==0)
1164    {
1165      mpz_clear(x->n);
1166      x->s=3;
1167      x=nlShort3(x);
1168    }
1169    else
1170    {
1171      mpz_t gcd;
1172      mpz_init(gcd);
1173      mpz_gcd(gcd,x->z,x->n);
1174      x->s=1;
1175      if (mpz_cmp_si(gcd,(long)1)!=0)
1176      {
1177        mpz_t r;
1178        mpz_init(r);
1179        MPZ_EXACTDIV(r,x->z,gcd);
1180        mpz_set(x->z,r);
1181        MPZ_EXACTDIV(r,x->n,gcd);
1182        mpz_set(x->n,r);
1183        mpz_clear(r);
1184        if (mpz_cmp_si(x->n,(long)1)==0)
1185        {
1186          mpz_clear(x->n);
1187          x->s=3;
1188          x=nlShort3_noinline(x);
1189        }
1190      }
1191      mpz_clear(gcd);
1192    }
1193  }
1194  nlTest(x, r);
1195}
1196
1197/*2
1198* returns in result->z the lcm(a->z,b->n)
1199*/
1200number nlLcm(number a, number b, const coeffs r)
1201{
1202  number result;
1203  nlTest(a, r);
1204  nlTest(b, r);
1205  if ((SR_HDL(b) & SR_INT)
1206  || (b->s==3))
1207  {
1208    // b is 1/(b->n) => b->n is 1 => result is a
1209    return nlCopy(a,r);
1210  }
1211  result=ALLOC_RNUMBER();
1212#if defined(LDEBUG)
1213  result->debug=123456;
1214#endif
1215  result->s=3;
1216  mpz_t gcd;
1217  mpz_init(gcd);
1218  mpz_init(result->z);
1219  if (SR_HDL(a) & SR_INT)
1220    mpz_gcd_ui(gcd,b->n,ABS(SR_TO_INT(a)));
1221  else
1222    mpz_gcd(gcd,a->z,b->n);
1223  if (mpz_cmp_si(gcd,(long)1)!=0)
1224  {
1225    mpz_t bt;
1226    mpz_init_set(bt,b->n);
1227    MPZ_EXACTDIV(bt,bt,gcd);
1228    if (SR_HDL(a) & SR_INT)
1229      mpz_mul_si(result->z,bt,SR_TO_INT(a));
1230    else
1231      mpz_mul(result->z,bt,a->z);
1232    mpz_clear(bt);
1233  }
1234  else
1235    if (SR_HDL(a) & SR_INT)
1236      mpz_mul_si(result->z,b->n,SR_TO_INT(a));
1237    else
1238      mpz_mul(result->z,b->n,a->z);
1239  mpz_clear(gcd);
1240  result=nlShort3(result);
1241  nlTest(result, r);
1242  return result;
1243}
1244
1245// Map q \in QQ \to Zp
1246// src = Q, dst = Zp (or an extension of Zp?)
1247number nlModP(number q, const coeffs Q, const coeffs Zp)
1248{
1249  assume( getCoeffType(Q) == ID );
1250
1251  const int p = n_GetChar(Zp);
1252  assume( p > 0 );
1253
1254  const long P = p;
1255  assume( P > 0 );
1256
1257  // embedded long within q => only long numerator has to be converted
1258  // to int (modulo char.)
1259  if (SR_HDL(q) & SR_INT)
1260  {
1261    long i = SR_TO_INT(q);
1262    if (i<0L)
1263      return n_Init( static_cast<int>( P - ((-i)%P) ), Zp);
1264    else
1265      return n_Init( static_cast<int>( i % P ), Zp );
1266  }
1267
1268  const unsigned long PP = p;
1269
1270  // numerator modulo char. should fit into int
1271  number z = n_Init( static_cast<int>(mpz_fdiv_ui(q->z, PP)), Zp );
1272
1273  // denominator != 1?
1274  if (q->s!=3)
1275  {
1276    // denominator modulo char. should fit into int
1277    number n = n_Init( static_cast<int>(mpz_fdiv_ui(q->n, PP)), Zp );
1278
1279    number res = n_Div( z, n, Zp );
1280
1281    n_Delete(&z, Zp);
1282    n_Delete(&n, Zp);
1283
1284    return res;
1285  }
1286
1287  return z;
1288}
1289
1290#ifdef HAVE_RINGS
1291/*2
1292* convert number i (from Q) to GMP and warn if denom != 1
1293*/
1294void nlGMP(number &i, number n, const coeffs r)
1295{
1296  // Hier brauche ich einfach die GMP Zahl
1297  nlTest(i, r);
1298  nlNormalize(i, r);
1299  if (SR_HDL(i) & SR_INT)
1300  {
1301    mpz_set_si((mpz_ptr) n, SR_TO_INT(i));
1302    return;
1303  }
1304  if (i->s!=3)
1305  {
1306     WarnS("Omitted denominator during coefficient mapping !");
1307  }
1308  mpz_set((mpz_ptr) n, i->z);
1309}
1310#endif
1311
1312/*2
1313* acces to denominator, other 1 for integers
1314*/
1315number   nlGetDenom(number &n, const coeffs r)
1316{
1317  if (!(SR_HDL(n) & SR_INT))
1318  {
1319    if (n->s==0)
1320    {
1321      nlNormalize(n,r);
1322    }
1323    if (!(SR_HDL(n) & SR_INT))
1324    {
1325      if (n->s!=3)
1326      {
1327        number u=ALLOC_RNUMBER();
1328        u->s=3;
1329#if defined(LDEBUG)
1330        u->debug=123456;
1331#endif
1332        mpz_init_set(u->z,n->n);
1333        u=nlShort3_noinline(u);
1334        return u;
1335      }
1336    }
1337  }
1338  return INT_TO_SR(1);
1339}
1340
1341/*2
1342* acces to Nominator, nlCopy(n) for integers
1343*/
1344number   nlGetNumerator(number &n, const coeffs r)
1345{
1346  if (!(SR_HDL(n) & SR_INT))
1347  {
1348    if (n->s==0)
1349    {
1350      nlNormalize(n,r);
1351    }
1352    if (!(SR_HDL(n) & SR_INT))
1353    {
1354      number u=ALLOC_RNUMBER();
1355#if defined(LDEBUG)
1356      u->debug=123456;
1357#endif
1358      u->s=3;
1359      mpz_init_set(u->z,n->z);
1360      if (n->s!=3)
1361      {
1362        u=nlShort3_noinline(u);
1363      }
1364      return u;
1365    }
1366  }
1367  return n; // imm. int
1368}
1369
1370/***************************************************************
1371 *
1372 * routines which are needed by Inline(d) routines
1373 *
1374 *******************************************************************/
1375BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
1376{
1377  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1378//  long - short
1379  BOOLEAN bo;
1380  if (SR_HDL(b) & SR_INT)
1381  {
1382    if (a->s!=0) return FALSE;
1383    number n=b; b=a; a=n;
1384  }
1385//  short - long
1386  if (SR_HDL(a) & SR_INT)
1387  {
1388    if (b->s!=0)
1389      return FALSE;
1390    if (((long)a > 0L) && (mpz_isNeg(b->z)))
1391      return FALSE;
1392    if (((long)a < 0L) && (!mpz_isNeg(b->z)))
1393      return FALSE;
1394    mpz_t  bb;
1395    mpz_init_set(bb,b->n);
1396    mpz_mul_si(bb,bb,(long)SR_TO_INT(a));
1397    bo=(mpz_cmp(bb,b->z)==0);
1398    mpz_clear(bb);
1399    return bo;
1400  }
1401// long - long
1402  if (((a->s==1) && (b->s==3))
1403  ||  ((b->s==1) && (a->s==3)))
1404    return FALSE;
1405  if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
1406    return FALSE;
1407  if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
1408    return FALSE;
1409  mpz_t  aa;
1410  mpz_t  bb;
1411  mpz_init_set(aa,a->z);
1412  mpz_init_set(bb,b->z);
1413  if (a->s<2) mpz_mul(bb,bb,a->n);
1414  if (b->s<2) mpz_mul(aa,aa,b->n);
1415  bo=(mpz_cmp(aa,bb)==0);
1416  mpz_clear(aa);
1417  mpz_clear(bb);
1418  return bo;
1419}
1420
1421// copy not immediate number a
1422number _nlCopy_NoImm(number a)
1423{
1424  assume(!((SR_HDL(a) & SR_INT)||(a==NULL)));
1425  //nlTest(a, r);
1426  number b=ALLOC_RNUMBER();
1427#if defined(LDEBUG)
1428  b->debug=123456;
1429#endif
1430  switch (a->s)
1431  {
1432    case 0:
1433    case 1:
1434            mpz_init_set(b->n,a->n);
1435    case 3:
1436            mpz_init_set(b->z,a->z);
1437            break;
1438  }
1439  b->s = a->s;
1440  return b;
1441}
1442
1443void _nlDelete_NoImm(number *a)
1444{
1445  {
1446    switch ((*a)->s)
1447    {
1448      case 0:
1449      case 1:
1450        mpz_clear((*a)->n);
1451      case 3:
1452        mpz_clear((*a)->z);
1453#ifdef LDEBUG
1454        (*a)->s=2;
1455#endif
1456    }
1457    FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
1458  }
1459}
1460
1461number _nlNeg_NoImm(number a)
1462{
1463  {
1464    mpz_neg(a->z,a->z);
1465    if (a->s==3)
1466    {
1467      a=nlShort3(a);
1468    }
1469  }
1470  return a;
1471}
1472
1473number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1474{
1475  number u=ALLOC_RNUMBER();
1476#if defined(LDEBUG)
1477  u->debug=123456;
1478#endif
1479  mpz_init(u->z);
1480  if (SR_HDL(b) & SR_INT)
1481  {
1482    number x=a;
1483    a=b;
1484    b=x;
1485  }
1486  if (SR_HDL(a) & SR_INT)
1487  {
1488    switch (b->s)
1489    {
1490      case 0:
1491      case 1:/* a:short, b:1 */
1492      {
1493        mpz_t x;
1494        mpz_init(x);
1495        mpz_mul_si(x,b->n,SR_TO_INT(a));
1496        mpz_add(u->z,b->z,x);
1497        mpz_clear(x);
1498        if (mpz_cmp_ui(u->z,(long)0)==0)
1499        {
1500          mpz_clear(u->z);
1501          FREE_RNUMBER(u);
1502          return INT_TO_SR(0);
1503        }
1504        if (mpz_cmp(u->z,b->n)==0)
1505        {
1506          mpz_clear(u->z);
1507          FREE_RNUMBER(u);
1508          return INT_TO_SR(1);
1509        }
1510        mpz_init_set(u->n,b->n);
1511        u->s = 0;
1512        break;
1513      }
1514      case 3:
1515      {
1516        if ((long)a>0L)
1517          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1518        else
1519          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1520        if (mpz_cmp_ui(u->z,(long)0)==0)
1521        {
1522          mpz_clear(u->z);
1523          FREE_RNUMBER(u);
1524          return INT_TO_SR(0);
1525        }
1526        u->s = 3;
1527        u=nlShort3(u);
1528        break;
1529      }
1530    }
1531  }
1532  else
1533  {
1534    switch (a->s)
1535    {
1536      case 0:
1537      case 1:
1538      {
1539        switch(b->s)
1540        {
1541          case 0:
1542          case 1:
1543          {
1544            mpz_t x;
1545            mpz_init(x);
1546
1547            mpz_mul(x,b->z,a->n);
1548            mpz_mul(u->z,a->z,b->n);
1549            mpz_add(u->z,u->z,x);
1550            mpz_clear(x);
1551
1552            if (mpz_cmp_ui(u->z,(long)0)==0)
1553            {
1554              mpz_clear(u->z);
1555              FREE_RNUMBER(u);
1556              return INT_TO_SR(0);
1557            }
1558            mpz_init(u->n);
1559            mpz_mul(u->n,a->n,b->n);
1560            if (mpz_cmp(u->z,u->n)==0)
1561            {
1562               mpz_clear(u->z);
1563               mpz_clear(u->n);
1564               FREE_RNUMBER(u);
1565               return INT_TO_SR(1);
1566            }
1567            u->s = 0;
1568            break;
1569          }
1570          case 3: /* a:1 b:3 */
1571          {
1572            mpz_mul(u->z,b->z,a->n);
1573            mpz_add(u->z,u->z,a->z);
1574            if (mpz_cmp_ui(u->z,(long)0)==0)
1575            {
1576              mpz_clear(u->z);
1577              FREE_RNUMBER(u);
1578              return INT_TO_SR(0);
1579            }
1580            if (mpz_cmp(u->z,a->n)==0)
1581            {
1582              mpz_clear(u->z);
1583              FREE_RNUMBER(u);
1584              return INT_TO_SR(1);
1585            }
1586            mpz_init_set(u->n,a->n);
1587            u->s = 0;
1588            break;
1589          }
1590        } /*switch (b->s) */
1591        break;
1592      }
1593      case 3:
1594      {
1595        switch(b->s)
1596        {
1597          case 0:
1598          case 1:/* a:3, b:1 */
1599          {
1600            mpz_mul(u->z,a->z,b->n);
1601            mpz_add(u->z,u->z,b->z);
1602            if (mpz_cmp_ui(u->z,(long)0)==0)
1603            {
1604              mpz_clear(u->z);
1605              FREE_RNUMBER(u);
1606              return INT_TO_SR(0);
1607            }
1608            if (mpz_cmp(u->z,b->n)==0)
1609            {
1610              mpz_clear(u->z);
1611              FREE_RNUMBER(u);
1612              return INT_TO_SR(1);
1613            }
1614            mpz_init_set(u->n,b->n);
1615            u->s = 0;
1616            break;
1617          }
1618          case 3:
1619          {
1620            mpz_add(u->z,a->z,b->z);
1621            if (mpz_cmp_ui(u->z,(long)0)==0)
1622            {
1623              mpz_clear(u->z);
1624              FREE_RNUMBER(u);
1625              return INT_TO_SR(0);
1626            }
1627            u->s = 3;
1628            u=nlShort3(u);
1629            break;
1630          }
1631        }
1632        break;
1633      }
1634    }
1635  }
1636  return u;
1637}
1638
1639void _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b)
1640{
1641  if (SR_HDL(b) & SR_INT)
1642  {
1643    switch (a->s)
1644    {
1645      case 0:
1646      case 1:/* b:short, a:1 */
1647      {
1648        mpz_t x;
1649        mpz_init(x);
1650        mpz_mul_si(x,a->n,SR_TO_INT(b));
1651        mpz_add(a->z,a->z,x);
1652        mpz_clear(x);
1653        a->s = 0;
1654        a=nlShort1(a);
1655        break;
1656      }
1657      case 3:
1658      {
1659        if ((long)b>0L)
1660          mpz_add_ui(a->z,a->z,SR_TO_INT(b));
1661        else
1662          mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
1663        a->s = 3;
1664        a=nlShort3_noinline(a);
1665        break;
1666      }
1667    }
1668    return;
1669  }
1670  else if (SR_HDL(a) & SR_INT)
1671  {
1672    number u=ALLOC_RNUMBER();
1673    #if defined(LDEBUG)
1674    u->debug=123456;
1675    #endif
1676    mpz_init(u->z);
1677    switch (b->s)
1678    {
1679      case 0:
1680      case 1:/* a:short, b:1 */
1681      {
1682        mpz_t x;
1683        mpz_init(x);
1684
1685        mpz_mul_si(x,b->n,SR_TO_INT(a));
1686        mpz_add(u->z,b->z,x);
1687        mpz_clear(x);
1688        // result cannot be 0, if coeffs are normalized
1689        mpz_init_set(u->n,b->n);
1690        u->s = 0;
1691        u=nlShort1(u);
1692        break;
1693      }
1694      case 3:
1695      {
1696        if ((long)a>0L)
1697          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1698        else
1699          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1700        // result cannot be 0, if coeffs are normalized
1701        u->s = 3;
1702        u=nlShort3_noinline(u);
1703        break;
1704      }
1705    }
1706    a=u;
1707  }
1708  else
1709  {
1710    switch (a->s)
1711    {
1712      case 0:
1713      case 1:
1714      {
1715        switch(b->s)
1716        {
1717          case 0:
1718          case 1: /* a:1 b:1 */
1719          {
1720            mpz_t x;
1721            mpz_t y;
1722            mpz_init(x);
1723            mpz_init(y);
1724            mpz_mul(x,b->z,a->n);
1725            mpz_mul(y,a->z,b->n);
1726            mpz_add(a->z,x,y);
1727            mpz_clear(x);
1728            mpz_clear(y);
1729            mpz_mul(a->n,a->n,b->n);
1730            a->s = 0;
1731            break;
1732          }
1733          case 3: /* a:1 b:3 */
1734          {
1735            mpz_t x;
1736            mpz_init(x);
1737            mpz_mul(x,b->z,a->n);
1738            mpz_add(a->z,a->z,x);
1739            mpz_clear(x);
1740            a->s = 0;
1741            break;
1742          }
1743        } /*switch (b->s) */
1744        a=nlShort1(a);
1745        break;
1746      }
1747      case 3:
1748      {
1749        switch(b->s)
1750        {
1751          case 0:
1752          case 1:/* a:3, b:1 */
1753          {
1754            mpz_t x;
1755            mpz_init(x);
1756            mpz_mul(x,a->z,b->n);
1757            mpz_add(a->z,b->z,x);
1758            mpz_clear(x);
1759            mpz_init_set(a->n,b->n);
1760            a->s = 0;
1761            a=nlShort1(a);
1762            break;
1763          }
1764          case 3:
1765          {
1766            mpz_add(a->z,a->z,b->z);
1767            a->s = 3;
1768            a=nlShort3_noinline(a);
1769            break;
1770          }
1771        }
1772        break;
1773      }
1774    }
1775  }
1776}
1777
1778number _nlSub_aNoImm_OR_bNoImm(number a, number b)
1779{
1780  number u=ALLOC_RNUMBER();
1781#if defined(LDEBUG)
1782  u->debug=123456;
1783#endif
1784  mpz_init(u->z);
1785  if (SR_HDL(a) & SR_INT)
1786  {
1787    switch (b->s)
1788    {
1789      case 0:
1790      case 1:/* a:short, b:1 */
1791      {
1792        mpz_t x;
1793        mpz_init(x);
1794        mpz_mul_si(x,b->n,SR_TO_INT(a));
1795        mpz_sub(u->z,x,b->z);
1796        mpz_clear(x);
1797        if (mpz_cmp_ui(u->z,(long)0)==0)
1798        {
1799          mpz_clear(u->z);
1800          FREE_RNUMBER(u);
1801          return INT_TO_SR(0);
1802        }
1803        if (mpz_cmp(u->z,b->n)==0)
1804        {
1805          mpz_clear(u->z);
1806          FREE_RNUMBER(u);
1807          return INT_TO_SR(1);
1808        }
1809        mpz_init_set(u->n,b->n);
1810        u->s = 0;
1811        break;
1812      }
1813      case 3:
1814      {
1815        if ((long)a>0L)
1816        {
1817          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
1818          mpz_neg(u->z,u->z);
1819        }
1820        else
1821        {
1822          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
1823          mpz_neg(u->z,u->z);
1824        }
1825        if (mpz_cmp_ui(u->z,(long)0)==0)
1826        {
1827          mpz_clear(u->z);
1828          FREE_RNUMBER(u);
1829          return INT_TO_SR(0);
1830        }
1831        u->s = 3;
1832        u=nlShort3(u);
1833        break;
1834      }
1835    }
1836  }
1837  else if (SR_HDL(b) & SR_INT)
1838  {
1839    switch (a->s)
1840    {
1841      case 0:
1842      case 1:/* b:short, a:1 */
1843      {
1844        mpz_t x;
1845        mpz_init(x);
1846        mpz_mul_si(x,a->n,SR_TO_INT(b));
1847        mpz_sub(u->z,a->z,x);
1848        mpz_clear(x);
1849        if (mpz_cmp_ui(u->z,(long)0)==0)
1850        {
1851          mpz_clear(u->z);
1852          FREE_RNUMBER(u);
1853          return INT_TO_SR(0);
1854        }
1855        if (mpz_cmp(u->z,a->n)==0)
1856        {
1857          mpz_clear(u->z);
1858          FREE_RNUMBER(u);
1859          return INT_TO_SR(1);
1860        }
1861        mpz_init_set(u->n,a->n);
1862        u->s = 0;
1863        break;
1864      }
1865      case 3:
1866      {
1867        if ((long)b>0L)
1868        {
1869          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
1870        }
1871        else
1872        {
1873          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
1874        }
1875        if (mpz_cmp_ui(u->z,(long)0)==0)
1876        {
1877          mpz_clear(u->z);
1878          FREE_RNUMBER(u);
1879          return INT_TO_SR(0);
1880        }
1881        u->s = 3;
1882        u=nlShort3(u);
1883        break;
1884      }
1885    }
1886  }
1887  else
1888  {
1889    switch (a->s)
1890    {
1891      case 0:
1892      case 1:
1893      {
1894        switch(b->s)
1895        {
1896          case 0:
1897          case 1:
1898          {
1899            mpz_t x;
1900            mpz_t y;
1901            mpz_init(x);
1902            mpz_init(y);
1903            mpz_mul(x,b->z,a->n);
1904            mpz_mul(y,a->z,b->n);
1905            mpz_sub(u->z,y,x);
1906            mpz_clear(x);
1907            mpz_clear(y);
1908            if (mpz_cmp_ui(u->z,(long)0)==0)
1909            {
1910              mpz_clear(u->z);
1911              FREE_RNUMBER(u);
1912              return INT_TO_SR(0);
1913            }
1914            mpz_init(u->n);
1915            mpz_mul(u->n,a->n,b->n);
1916            if (mpz_cmp(u->z,u->n)==0)
1917            {
1918              mpz_clear(u->z);
1919              mpz_clear(u->n);
1920              FREE_RNUMBER(u);
1921              return INT_TO_SR(1);
1922            }
1923            u->s = 0;
1924            break;
1925          }
1926          case 3: /* a:1, b:3 */
1927          {
1928            mpz_t x;
1929            mpz_init(x);
1930            mpz_mul(x,b->z,a->n);
1931            mpz_sub(u->z,a->z,x);
1932            mpz_clear(x);
1933            if (mpz_cmp_ui(u->z,(long)0)==0)
1934            {
1935              mpz_clear(u->z);
1936              FREE_RNUMBER(u);
1937              return INT_TO_SR(0);
1938            }
1939            if (mpz_cmp(u->z,a->n)==0)
1940            {
1941              mpz_clear(u->z);
1942              FREE_RNUMBER(u);
1943              return INT_TO_SR(1);
1944            }
1945            mpz_init_set(u->n,a->n);
1946            u->s = 0;
1947            break;
1948          }
1949        }
1950        break;
1951      }
1952      case 3:
1953      {
1954        switch(b->s)
1955        {
1956          case 0:
1957          case 1: /* a:3, b:1 */
1958          {
1959            mpz_t x;
1960            mpz_init(x);
1961            mpz_mul(x,a->z,b->n);
1962            mpz_sub(u->z,x,b->z);
1963            mpz_clear(x);
1964            if (mpz_cmp_ui(u->z,(long)0)==0)
1965            {
1966              mpz_clear(u->z);
1967              FREE_RNUMBER(u);
1968              return INT_TO_SR(0);
1969            }
1970            if (mpz_cmp(u->z,b->n)==0)
1971            {
1972              mpz_clear(u->z);
1973              FREE_RNUMBER(u);
1974              return INT_TO_SR(1);
1975            }
1976            mpz_init_set(u->n,b->n);
1977            u->s = 0;
1978            break;
1979          }
1980          case 3: /* a:3 , b:3 */
1981          {
1982            mpz_sub(u->z,a->z,b->z);
1983            if (mpz_cmp_ui(u->z,(long)0)==0)
1984            {
1985              mpz_clear(u->z);
1986              FREE_RNUMBER(u);
1987              return INT_TO_SR(0);
1988            }
1989            u->s = 3;
1990            u=nlShort3(u);
1991            break;
1992          }
1993        }
1994        break;
1995      }
1996    }
1997  }
1998  return u;
1999}
2000
2001// a and b are intermediate, but a*b not
2002number _nlMult_aImm_bImm_rNoImm(number a, number b)
2003{
2004  number u=ALLOC_RNUMBER();
2005#if defined(LDEBUG)
2006  u->debug=123456;
2007#endif
2008  u->s=3;
2009  mpz_init_set_si(u->z,SR_TO_INT(a));
2010  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
2011  return u;
2012}
2013
2014// a or b are not immediate
2015number _nlMult_aNoImm_OR_bNoImm(number a, number b)
2016{
2017  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2018  number u=ALLOC_RNUMBER();
2019#if defined(LDEBUG)
2020  u->debug=123456;
2021#endif
2022  mpz_init(u->z);
2023  if (SR_HDL(b) & SR_INT)
2024  {
2025    number x=a;
2026    a=b;
2027    b=x;
2028  }
2029  if (SR_HDL(a) & SR_INT)
2030  {
2031    u->s=b->s;
2032    if (u->s==1) u->s=0;
2033    if ((long)a>0L)
2034    {
2035      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2036    }
2037    else
2038    {
2039      if (a==INT_TO_SR(-1))
2040      {
2041        mpz_set(u->z,b->z);
2042        mpz_neg(u->z,u->z);
2043        u->s=b->s;
2044      }
2045      else
2046      {
2047        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2048        mpz_neg(u->z,u->z);
2049      }
2050    }
2051    if (u->s<2)
2052    {
2053      if (mpz_cmp(u->z,b->n)==0)
2054      {
2055        mpz_clear(u->z);
2056        FREE_RNUMBER(u);
2057        return INT_TO_SR(1);
2058      }
2059      mpz_init_set(u->n,b->n);
2060    }
2061    else //u->s==3
2062    {
2063      u=nlShort3(u);
2064    }
2065  }
2066  else
2067  {
2068    mpz_mul(u->z,a->z,b->z);
2069    u->s = 0;
2070    if(a->s==3)
2071    {
2072      if(b->s==3)
2073      {
2074        u->s = 3;
2075      }
2076      else
2077      {
2078        if (mpz_cmp(u->z,b->n)==0)
2079        {
2080          mpz_clear(u->z);
2081          FREE_RNUMBER(u);
2082          return INT_TO_SR(1);
2083        }
2084        mpz_init_set(u->n,b->n);
2085      }
2086    }
2087    else
2088    {
2089      if(b->s==3)
2090      {
2091        if (mpz_cmp(u->z,a->n)==0)
2092        {
2093          mpz_clear(u->z);
2094          FREE_RNUMBER(u);
2095          return INT_TO_SR(1);
2096        }
2097        mpz_init_set(u->n,a->n);
2098      }
2099      else
2100      {
2101        mpz_init(u->n);
2102        mpz_mul(u->n,a->n,b->n);
2103        if (mpz_cmp(u->z,u->n)==0)
2104        {
2105          mpz_clear(u->z);
2106          mpz_clear(u->n);
2107          FREE_RNUMBER(u);
2108          return INT_TO_SR(1);
2109        }
2110      }
2111    }
2112  }
2113  return u;
2114}
2115
2116/*2
2117* copy a to b for mapping
2118*/
2119number nlCopyMap(number a, const coeffs src, const coeffs dst)
2120{
2121  assume( getCoeffType(dst) == ID );
2122  assume( getCoeffType(src) == ID );
2123
2124  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2125  {
2126    return a;
2127  }
2128  return _nlCopy_NoImm(a);
2129}
2130
2131nMapFunc nlSetMap(const coeffs src, const coeffs dst)
2132{
2133  assume( getCoeffType(dst) == ID );
2134//  assume( getCoeffType(src) == ID );
2135
2136  if (nCoeff_is_Q(src))
2137  {
2138    return ndCopyMap;
2139  }
2140  if (nCoeff_is_Zp(src))
2141  {
2142    return nlMapP;
2143  }
2144  if (nCoeff_is_R(src))
2145  {
2146    return nlMapR;
2147  }
2148  if (nCoeff_is_long_R(src))
2149  {
2150    return nlMapLongR; /* long R -> Q */
2151  }
2152#ifdef HAVE_RINGS
2153  if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
2154  {
2155    return nlMapGMP;
2156  }
2157  if (nCoeff_is_Ring_2toM(src))
2158  {
2159    return nlMapMachineInt;
2160  }
2161#endif
2162  return NULL;
2163}
2164/*2
2165* z := i
2166*/
2167number nlRInit (long i)
2168{
2169  number z=ALLOC_RNUMBER();
2170#if defined(LDEBUG)
2171  z->debug=123456;
2172#endif
2173  mpz_init_set_si(z->z,i);
2174  z->s = 3;
2175  return z;
2176}
2177
2178/*2
2179* z := i/j
2180*/
2181number nlInit2 (int i, int j, const coeffs r)
2182{
2183  number z=ALLOC_RNUMBER();
2184#if defined(LDEBUG)
2185  z->debug=123456;
2186#endif
2187  mpz_init_set_si(z->z,(long)i);
2188  mpz_init_set_si(z->n,(long)j);
2189  z->s = 0;
2190  nlNormalize(z,r);
2191  return z;
2192}
2193
2194number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
2195{
2196  number z=ALLOC_RNUMBER();
2197#if defined(LDEBUG)
2198  z->debug=123456;
2199#endif
2200  mpz_init_set(z->z,i);
2201  mpz_init_set(z->n,j);
2202  z->s = 0;
2203  nlNormalize(z,r);
2204  return z;
2205}
2206
2207#else // DO_LINLINE
2208
2209// declare immedate routines
2210number nlRInit (long i);
2211BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2212number  _nlCopy_NoImm(number a);
2213number  _nlNeg_NoImm(number a);
2214number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2215void    _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b);
2216number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2217number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2218number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2219
2220#endif
2221
2222/***************************************************************
2223 *
2224 * Routines which might be inlined by p_Numbers.h
2225 *
2226 *******************************************************************/
2227#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2228
2229// routines which are always inlined/static
2230
2231/*2
2232* a = b ?
2233*/
2234LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
2235{
2236  nlTest(a, r);
2237  nlTest(b, r);
2238// short - short
2239  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2240  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2241}
2242
2243LINLINE number nlInit (long i, const coeffs r)
2244{
2245  number n;
2246  LONG ii=(LONG)i;
2247  if ( ((ii << 3) >> 3) == ii ) n=INT_TO_SR(ii);
2248  else                          n=nlRInit(ii);
2249  nlTest(n, r);
2250  return n;
2251}
2252
2253/*2
2254* a == 1 ?
2255*/
2256LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
2257{
2258#ifdef LDEBUG
2259  if (a==NULL) return FALSE;
2260  nlTest(a, r);
2261#endif
2262  return (a==INT_TO_SR(1));
2263}
2264
2265LINLINE BOOLEAN nlIsZero (number a, const coeffs)
2266{
2267  #if 0
2268  if (a==INT_TO_SR(0)) return TRUE;
2269  if ((SR_HDL(a) & SR_INT)||(a==NULL)) return FALSE;
2270  if (mpz_cmp_si(a->z,(long)0)==0)
2271  {
2272    printf("gmp-0 in nlIsZero\n");
2273    dErrorBreak();
2274    return TRUE;
2275  }
2276  return FALSE;
2277  #else
2278  return (a==INT_TO_SR(0));
2279  #endif
2280}
2281
2282/*2
2283* copy a to b
2284*/
2285LINLINE number nlCopy(number a, const coeffs)
2286{
2287  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2288  {
2289    return a;
2290  }
2291  return _nlCopy_NoImm(a);
2292}
2293
2294
2295/*2
2296* delete a
2297*/
2298LINLINE void nlDelete (number * a, const coeffs r)
2299{
2300  if (*a!=NULL)
2301  {
2302    nlTest(*a, r);
2303    if ((SR_HDL(*a) & SR_INT)==0)
2304    {
2305      _nlDelete_NoImm(a);
2306    }
2307    *a=NULL;
2308  }
2309}
2310
2311/*2
2312* za:= - za
2313*/
2314LINLINE number nlNeg (number a, const coeffs R)
2315{
2316  nlTest(a, R);
2317  if(SR_HDL(a) &SR_INT)
2318  {
2319    LONG r=SR_TO_INT(a);
2320    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2321    else               a=INT_TO_SR(-r);
2322    return a;
2323  }
2324  a = _nlNeg_NoImm(a);
2325  nlTest(a, R);
2326  return a;
2327
2328}
2329
2330/*2
2331* u:= a + b
2332*/
2333LINLINE number nlAdd (number a, number b, const coeffs R)
2334{
2335  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2336  {
2337    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2338    if ( ((r << 1) >> 1) == r )
2339      return (number)(long)r;
2340    else
2341      return nlRInit(SR_TO_INT(r));
2342  }
2343  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2344  nlTest(u, R);
2345  return u;
2346}
2347
2348number nlShort1(number a);
2349number nlShort3_noinline(number x);
2350
2351LINLINE void nlInpAdd(number &a, number b, const coeffs r)
2352{
2353  // a=a+b
2354  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2355  {
2356    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2357    if ( ((r << 1) >> 1) == r )
2358      a=(number)(long)r;
2359    else
2360      a=nlRInit(SR_TO_INT(r));
2361  }
2362  else
2363  {
2364    _nlInpAdd_aNoImm_OR_bNoImm(a,b);
2365    nlTest(a,r);
2366  }
2367}
2368
2369LINLINE number nlMult (number a, number b, const coeffs R)
2370{
2371  nlTest(a, R);
2372  nlTest(b, R);
2373  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2374  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2375  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2376  {
2377    LONG r=(LONG)((unsigned LONG)(SR_HDL(a)-1L))*((unsigned LONG)(SR_HDL(b)>>1));
2378    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2379    {
2380      number u=((number) ((r>>1)+SR_INT));
2381      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2382      return nlRInit(SR_HDL(u)>>2);
2383    }
2384    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2385    nlTest(u, R);
2386    return u;
2387
2388  }
2389  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2390  nlTest(u, R);
2391  return u;
2392
2393}
2394
2395
2396/*2
2397* u:= a - b
2398*/
2399LINLINE number nlSub (number a, number b, const coeffs r)
2400{
2401  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2402  {
2403    LONG r=SR_HDL(a)-SR_HDL(b)+1;
2404    if ( ((r << 1) >> 1) == r )
2405    {
2406      return (number)(long)r;
2407    }
2408    else
2409      return nlRInit(SR_TO_INT(r));
2410  }
2411  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2412  nlTest(u, r);
2413  return u;
2414
2415}
2416
2417LINLINE void nlInpMult(number &a, number b, const coeffs r)
2418{
2419  number aa=a;
2420  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2421  {
2422    number n=nlMult(aa,b,r);
2423    nlDelete(&a,r);
2424    a=n;
2425  }
2426  else
2427  {
2428    mpz_mul(aa->z,a->z,b->z);
2429    if (aa->s==3)
2430    {
2431      if(b->s!=3)
2432      {
2433        mpz_init_set(a->n,b->n);
2434        a->s=0;
2435      }
2436    }
2437    else
2438    {
2439      if(b->s!=3)
2440      {
2441        mpz_mul(a->n,a->n,b->n);
2442      }
2443      a->s=0;
2444    }
2445  }
2446}
2447#endif // DO_LINLINE
2448
2449#ifndef P_NUMBERS_H
2450
2451static void nlMPZ(mpz_t m, number &n, const coeffs r)
2452{
2453  assume( getCoeffType(r) == ID );
2454
2455  nlTest(n, r);
2456  nlNormalize(n, r);
2457  if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n));     /* n fits in an int */
2458  else             mpz_init_set(m, (mpz_ptr)n->z);
2459}
2460
2461
2462static number nlInitMPZ(mpz_t m, const coeffs)
2463{
2464  number z = ALLOC_RNUMBER();
2465  mpz_init_set(z->z, m);
2466  mpz_init_set_ui(z->n, 1);
2467  z->s = 3;
2468  return z;
2469}
2470
2471
2472void nlInpGcd(number &a, number b, const coeffs r)
2473{
2474  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2475  {
2476    number n=nlGcd(a,b,r);
2477    nlDelete(&a,r);
2478    a=n;
2479  }
2480  else
2481  {
2482    mpz_gcd(a->z,a->z,b->z);
2483    a=nlShort3_noinline(a);
2484  }
2485}
2486
2487void nlInpIntDiv(number &a, number b, const coeffs r)
2488{
2489  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2490  {
2491    number n=nlIntDiv(a,b, r);
2492    nlDelete(&a,r);
2493    a=n;
2494  }
2495  else
2496  {
2497    if (mpz_isNeg(a->z))
2498    {
2499      if (mpz_isNeg(b->z))
2500      {
2501        mpz_add(a->z,a->z,b->z);
2502      }
2503      else
2504      {
2505        mpz_sub(a->z,a->z,b->z);
2506      }
2507      mpz_add_ui(a->z,a->z,1);
2508    }
2509    else
2510    {
2511      if (mpz_isNeg(b->z))
2512      {
2513        mpz_sub(a->z,a->z,b->z);
2514      }
2515      else
2516      {
2517        mpz_add(a->z,a->z,b->z);
2518      }
2519      mpz_sub_ui(a->z,a->z,1);
2520    }
2521    MPZ_DIV(a->z,a->z,b->z);
2522    a=nlShort3_noinline(a);
2523  }
2524}
2525
2526number nlFarey(number nN, number nP, const coeffs r)
2527{
2528  mpz_t tmp; mpz_init(tmp);
2529  mpz_t A,B,C,D,E,N,P;
2530  if (SR_HDL(nN) & SR_INT) mpz_init_set_si(N,SR_TO_INT(nN));
2531  else                     mpz_init_set(N,nN->z);
2532  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2533  else                     mpz_init_set(P,nP->z);
2534  assume(!mpz_isNeg(P));
2535  if (mpz_isNeg(N))  mpz_add(N,N,P);
2536  mpz_init_set_si(A,(long)0);
2537  mpz_init_set_ui(B,(unsigned long)1);
2538  mpz_init_set_si(C,(long)0);
2539  mpz_init(D);
2540  mpz_init_set(E,P);
2541  number z=INT_TO_SR(0);
2542  while(mpz_cmp_si(N,(long)0)!=0)
2543  {
2544    mpz_mul(tmp,N,N);
2545    mpz_add(tmp,tmp,tmp);
2546    if (mpz_cmp(tmp,P)<0)
2547    {
2548       if (mpz_isNeg(B))
2549       {
2550         mpz_neg(B,B);
2551         mpz_neg(N,N);
2552       }
2553       // check for gcd(N,B)==1
2554       mpz_gcd(tmp,N,B);
2555       if (mpz_cmp_ui(tmp,1)==0)
2556       {
2557         // return N/B
2558         z=ALLOC_RNUMBER();
2559         #ifdef LDEBUG
2560         z->debug=123456;
2561         #endif
2562         mpz_init_set(z->z,N);
2563         mpz_init_set(z->n,B);
2564         z->s = 0;
2565         nlNormalize(z,r);
2566       }
2567       else
2568       {
2569         // return nN (the input) instead of "fail"
2570         z=nlCopy(nN,r);
2571       }
2572       break;
2573    }
2574    //mpz_mod(D,E,N);
2575    //mpz_div(tmp,E,N);
2576    mpz_divmod(tmp,D,E,N);
2577    mpz_mul(tmp,tmp,B);
2578    mpz_sub(C,A,tmp);
2579    mpz_set(E,N);
2580    mpz_set(N,D);
2581    mpz_set(A,B);
2582    mpz_set(B,C);
2583  }
2584  mpz_clear(tmp);
2585  mpz_clear(A);
2586  mpz_clear(B);
2587  mpz_clear(C);
2588  mpz_clear(D);
2589  mpz_clear(E);
2590  mpz_clear(N);
2591  mpz_clear(P);
2592  return z;
2593}
2594
2595number nlExtGcd(number a, number b, number *s, number *t, const coeffs)
2596{
2597  mpz_t aa,bb;
2598  *s=ALLOC_RNUMBER();
2599  mpz_init((*s)->z); (*s)->s=3;
2600  (*t)=ALLOC_RNUMBER();
2601  mpz_init((*t)->z); (*t)->s=3;
2602  number g=ALLOC_RNUMBER();
2603  mpz_init(g->z); g->s=3;
2604  if (SR_HDL(a) & SR_INT)
2605  {
2606    mpz_init_set_si(aa,SR_TO_INT(a));
2607  }
2608  else
2609  {
2610    mpz_init_set(aa,a->z);
2611  }
2612  if (SR_HDL(b) & SR_INT)
2613  {
2614    mpz_init_set_si(bb,SR_TO_INT(b));
2615  }
2616  else
2617  {
2618    mpz_init_set(bb,b->z);
2619  }
2620  mpz_gcdext(g->z,(*s)->z,(*t)->z,aa,bb);
2621  mpz_clear(aa);
2622  mpz_clear(bb);
2623  (*s)=nlShort3((*s));
2624  (*t)=nlShort3((*t));
2625  g=nlShort3(g);
2626  return g;
2627}
2628
2629void    nlCoeffWrite  (const coeffs, BOOLEAN /*details*/)
2630{
2631  PrintS("//   characteristic : 0\n");
2632}
2633
2634number   nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
2635// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
2636{
2637#ifdef HAVE_FACTORY
2638  setCharacteristic( 0 ); // only in char 0
2639  CFArray X(rl), Q(rl);
2640  int i;
2641  for(i=rl-1;i>=0;i--)
2642  {
2643    X[i]=C->convSingNFactoryN(x[i],FALSE,C); // may be larger MAX_INT
2644    Q[i]=C->convSingNFactoryN(q[i],FALSE,C); // may be larger MAX_INT
2645  }
2646  CanonicalForm xnew,qnew;
2647  chineseRemainder(X,Q,xnew,qnew);
2648  number n=C->convFactoryNSingN(xnew,C);
2649  number p=C->convFactoryNSingN(qnew,C);
2650  number p2=nlIntDiv(p,nlInit(2, C),C);
2651  if (nlGreater(n,p2,C))
2652  {
2653     number n2=nlSub(n,p,C);
2654     nlDelete(&n,C);
2655     n=n2;
2656  }
2657  nlDelete(&p,C);
2658  nlDelete(&p2,C);
2659  return n;
2660#else
2661  WerrorS("not implemented");
2662  return nlInit(0, C);
2663#endif
2664}
2665
2666static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2667{
2668  assume(cf != NULL);
2669  assume(getCoeffType(cf) == ID);
2670
2671  numberCollectionEnumerator.Reset();
2672
2673  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2674  {
2675    c = n_Init(1, cf);
2676    return;
2677  }
2678 
2679  // all coeffs are given by integers!!!
2680
2681  // part 1, find a small candidate for gcd
2682  number cand1,cand;
2683  int s1,s;
2684  s=2147483647; // max. int
2685
2686  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
2687
2688  int normalcount = 0;
2689  do
2690  {
2691    number& n = numberCollectionEnumerator.Current();
2692    nlNormalize(n, cf); ++normalcount;
2693    cand1 = n;
2694   
2695    if (SR_HDL(cand1)&SR_INT) { cand=cand1; break; }
2696    assume(cand1->s==3); // all coeffs should be integers // ==0?!! after printing
2697    s1=mpz_size1(cand1->z);
2698    if (s>s1)
2699    {
2700      cand=cand1;
2701      s=s1;
2702    }
2703  } while (numberCollectionEnumerator.MoveNext() );
2704
2705//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
2706
2707  cand=nlCopy(cand,cf);
2708  // part 2: compute gcd(cand,all coeffs)
2709
2710  numberCollectionEnumerator.Reset();
2711 
2712  while (numberCollectionEnumerator.MoveNext() )
2713  {
2714    number& n = numberCollectionEnumerator.Current();
2715
2716    if( (--normalcount) <= 0)
2717      nlNormalize(n, cf);
2718
2719    nlInpGcd(cand, n, cf);
2720
2721    assume( nlGreaterZero(cand,cf) );
2722   
2723    if(nlIsOne(cand,cf))
2724    {
2725      c = cand;
2726
2727      if(!lc_is_pos)
2728      {
2729        // make the leading coeff positive
2730        c = nlNeg(c, cf);
2731        numberCollectionEnumerator.Reset();
2732       
2733        while (numberCollectionEnumerator.MoveNext() )
2734        {
2735          number& nn = numberCollectionEnumerator.Current();
2736          nn = nlNeg(nn, cf);
2737        }
2738      }
2739      return;
2740    }
2741  } 
2742
2743  // part3: all coeffs = all coeffs / cand
2744  if (!lc_is_pos)
2745    cand = nlNeg(cand,cf);
2746 
2747  c = cand;
2748  numberCollectionEnumerator.Reset();
2749
2750  while (numberCollectionEnumerator.MoveNext() )
2751  {
2752    number& n = numberCollectionEnumerator.Current();
2753    number t=nlIntDiv(n, cand, cf); // simple integer exact division, no ratios to remain
2754    nlDelete(&n, cf);
2755    n = t;
2756  }
2757}
2758
2759static void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2760{
2761  assume(cf != NULL);
2762  assume(getCoeffType(cf) == ID);
2763
2764  numberCollectionEnumerator.Reset();
2765 
2766  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2767  {
2768    c = n_Init(1, cf);
2769//    assume( n_GreaterZero(c, cf) );
2770    return;
2771  }
2772
2773  // all coeffs are given by integers after returning from this routine
2774
2775  // part 1, collect product of all denominators /gcds
2776  number cand;
2777  cand=ALLOC_RNUMBER();
2778#if defined(LDEBUG)
2779  cand->debug=123456;
2780#endif
2781  cand->s=3;
2782
2783  int s=0;
2784 
2785  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
2786
2787  do
2788  {
2789    number& cand1 = numberCollectionEnumerator.Current();
2790
2791    if (!(SR_HDL(cand1)&SR_INT))
2792    {
2793      nlNormalize(cand1, cf);
2794      if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
2795      && (cand1->s==1))             // and is a normalised rational
2796      {
2797        if (s==0) // first denom, we meet
2798        {
2799          mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
2800          s=1;
2801        }
2802        else // we have already something
2803        {
2804          mpz_lcm(cand->z, cand->z, cand1->n);
2805        }
2806      }
2807    }
2808  }
2809  while (numberCollectionEnumerator.MoveNext() );
2810 
2811
2812  if (s==0) // nothing to do, all coeffs are already integers
2813  {
2814//    mpz_clear(tmp);
2815    FREE_RNUMBER(cand);
2816    if (lc_is_pos)
2817      c=nlInit(1,cf);
2818    else
2819    {
2820      // make the leading coeff positive
2821      c=nlInit(-1,cf);
2822
2823      // TODO: incorporate the following into the loop below?
2824      numberCollectionEnumerator.Reset();
2825      while (numberCollectionEnumerator.MoveNext() )
2826      {
2827        number& n = numberCollectionEnumerator.Current();
2828        n = nlNeg(n, cf);
2829      } 
2830    }
2831//    assume( n_GreaterZero(c, cf) );
2832    return;
2833  }
2834
2835  cand = nlShort3(cand);
2836
2837  // part2: all coeffs = all coeffs * cand
2838  // make the lead coeff positive
2839  numberCollectionEnumerator.Reset();
2840 
2841  if (!lc_is_pos)
2842    cand = nlNeg(cand, cf);
2843 
2844  c = cand;
2845 
2846  while (numberCollectionEnumerator.MoveNext() )
2847  {
2848    number &n = numberCollectionEnumerator.Current();
2849    n_InpMult(n, cand, cf);
2850  }
2851
2852}
2853
2854BOOLEAN nlInitChar(coeffs r, void*)
2855{
2856  assume( getCoeffType(r) == ID );
2857  //const int ch = (int)(long)(p);
2858
2859  r->cfKillChar=NULL;
2860  r->nCoeffIsEqual=ndCoeffIsEqual;
2861  r->cfKillChar = ndKillChar; /* dummy */
2862
2863  r->cfInitMPZ = nlInitMPZ;
2864  r->cfMPZ  = nlMPZ;
2865
2866  r->cfMult  = nlMult;
2867  r->cfSub   = nlSub;
2868  r->cfAdd   = nlAdd;
2869  r->cfDiv   = nlDiv;
2870  r->cfIntDiv= nlIntDiv;
2871  r->cfIntMod= nlIntMod;
2872  r->cfExactDiv= nlExactDiv;
2873  r->cfInit = nlInit;
2874  r->cfSize  = nlSize;
2875  r->cfInt  = nlInt;
2876
2877  r->cfChineseRemainder=nlChineseRemainder;
2878  r->cfFarey=nlFarey;
2879  #ifdef HAVE_RINGS
2880  //r->cfDivComp = NULL; // only for ring stuff
2881  //r->cfIsUnit = NULL; // only for ring stuff
2882  //r->cfGetUnit = NULL; // only for ring stuff
2883  //r->cfDivBy = NULL; // only for ring stuff
2884  #endif
2885  r->cfNeg   = nlNeg;
2886  r->cfInvers= nlInvers;
2887  r->cfCopy  = nlCopy;
2888  r->cfRePart = nlCopy;
2889  //r->cfImPart = ndReturn0;
2890  r->cfWriteLong = nlWrite;
2891  r->cfRead = nlRead;
2892  r->cfNormalize=nlNormalize;
2893  r->cfGreater = nlGreater;
2894  r->cfEqual = nlEqual;
2895  r->cfIsZero = nlIsZero;
2896  r->cfIsOne = nlIsOne;
2897  r->cfIsMOne = nlIsMOne;
2898  r->cfGreaterZero = nlGreaterZero;
2899  r->cfPower = nlPower;
2900  r->cfGetDenom = nlGetDenom;
2901  r->cfGetNumerator = nlGetNumerator;
2902  r->cfGcd  = nlGcd;
2903  r->cfExtGcd = nlExtGcd; // only for ring stuff and Z
2904  r->cfLcm  = nlLcm;
2905  r->cfDelete= nlDelete;
2906  r->cfSetMap = nlSetMap;
2907  //r->cfName = ndName;
2908  r->cfInpMult=nlInpMult;
2909  r->cfInit_bigint=nlCopyMap;
2910  r->cfCoeffWrite=nlCoeffWrite;
2911
2912  r->cfClearContent = nlClearContent;
2913  r->cfClearDenominators = nlClearDenominators;
2914
2915#ifdef LDEBUG
2916  // debug stuff
2917  r->cfDBTest=nlDBTest;
2918#endif
2919#ifdef HAVE_FACTORY
2920  r->convSingNFactoryN=nlConvSingNFactoryN;
2921  r->convFactoryNSingN=nlConvFactoryNSingN;
2922#endif
2923
2924  // the variables: general stuff (required)
2925  r->nNULL = INT_TO_SR(0);
2926  //r->type = n_Q;
2927  r->ch = 0;
2928  r->has_simple_Alloc=FALSE;
2929  r->has_simple_Inverse=FALSE;
2930
2931  // variables for this type of coeffs:
2932  // (none)
2933  return FALSE;
2934}
2935#if 0
2936number nlMod(number a, number b)
2937{
2938  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
2939  {
2940    int bi=SR_TO_INT(b);
2941    int ai=SR_TO_INT(a);
2942    int bb=ABS(bi);
2943    int c=ai%bb;
2944    if (c<0)  c+=bb;
2945    return (INT_TO_SR(c));
2946  }
2947  number al;
2948  number bl;
2949  if (SR_HDL(a))&SR_INT)
2950    al=nlRInit(SR_TO_INT(a));
2951  else
2952    al=nlCopy(a);
2953  if (SR_HDL(b))&SR_INT)
2954    bl=nlRInit(SR_TO_INT(b));
2955  else
2956    bl=nlCopy(b);
2957  number r=nlRInit(0);
2958  mpz_mod(r->z,al->z,bl->z);
2959  nlDelete(&al);
2960  nlDelete(&bl);
2961  if (mpz_size1(&r->z)<=MP_SMALL)
2962  {
2963    LONG ui=(int)mpz_get_si(&r->z);
2964    if ((((ui<<3)>>3)==ui)
2965    && (mpz_cmp_si(x->z,(long)ui)==0))
2966    {
2967      mpz_clear(&r->z);
2968      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
2969      r=INT_TO_SR(ui);
2970    }
2971  }
2972  return r;
2973}
2974#endif
2975#endif // not P_NUMBERS_H
2976#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.