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