My Project
longrat.cc
Go to the documentation of this file.
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 
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 
159 number nlShort3_noinline(number x) // assume x->s==3
160 {
161  return nlShort3(x);
162 }
163 
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)
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 
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 */
206 static inline number nlMapGMP(number from, const coeffs /*src*/, const coeffs dst)
207 {
208  return nlInitMPZ((mpz_ptr)from,dst);
209 }
210 
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 */
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
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 
330 static CanonicalForm nlConvSingNFactoryN( number n, const BOOLEAN setChar, const coeffs /*r*/ )
331 {
332  if (setChar) setCharacteristic( 0 );
333 
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 
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 
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 
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 
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 
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 */
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 */
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 */
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 */
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 */
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 */
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 
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 
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 
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 
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 
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 */
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 */
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 */
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 */
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 */
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 */
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)
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)
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)
1359  else
1360  result=INT_TO_SR(t);
1361  }
1362  else
1363  {
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);
1372  }
1373  nlTest(result, r);
1374  return result;
1375 }
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 
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 */
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  {
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;
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 */
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  }
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);
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)
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 */
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 */
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 */
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  *******************************************************************/
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
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 
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 
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 
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;
1780  }
1781  }
1782  mpz_clear(gcd);
1783 }
1784 
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 
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 
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
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
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 */
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 
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 
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 */
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 */
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 
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 */
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 
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 */
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 
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 */
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 */
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 */
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 */
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 
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  {
2693  nlTest(a,r);
2694  }
2695 }
2696 
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 */
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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);
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 
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 
3282 char * nlCoeffName(const coeffs r)
3283 {
3284  if (r->cfDiv==nlDiv) return (char*)"QQ";
3285  else return (char*)"ZZ";
3286 }
3287 
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 
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 
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 
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 
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 
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
All the auxiliary stuff.
#define SSI_BASE
Definition: auxiliary.h:135
static int ABS(int v)
Definition: auxiliary.h:112
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void On(int sw)
switches
void Off(int sw)
switches
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
Variable x
Definition: cfModGcd.cc:4084
int p
Definition: cfModGcd.cc:4080
g
Definition: cfModGcd.cc:4092
CanonicalForm cf
Definition: cfModGcd.cc:4085
CanonicalForm b
Definition: cfModGcd.cc:4105
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition: cf_chinese.cc:57
void FACTORY_PUBLIC chineseRemainderCached(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew, CFArray &inv)
Definition: cf_chinese.cc:308
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:30
FILE * f
Definition: checklibs.c:9
factory's main class
Definition: canonicalform.h:86
virtual reference Current()=0
Gets the current element in the collection (read and write).
virtual void Reset()=0
Sets the enumerator to its initial position: -1, which is before the first element in the collection.
virtual bool MoveNext()=0
Advances the enumerator to the next element of the collection. returns true if the enumerator was suc...
Templated enumerator interface for simple iteration over a generic collection of T's.
Definition: Enumerator.h:125
gmp_complex numbers based on
Definition: mpr_complex.h:179
mpf_t * _mpfp()
Definition: mpr_complex.h:134
Definition: int_poly.h:33
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE BOOLEAN nCoeff_is_long_R(const coeffs r)
Definition: coeffs.h:915
n_coeffType
Definition: coeffs.h:28
@ n_R
single prescision (6,6) real numbers
Definition: coeffs.h:32
@ n_Q
rational (GMP) numbers
Definition: coeffs.h:31
@ n_Zn
only used if HAVE_RINGS is defined
Definition: coeffs.h:45
@ n_long_R
real floating point (GMP) numbers
Definition: coeffs.h:34
@ n_Zp
\F{p < 2^31}
Definition: coeffs.h:30
@ n_long_C
complex floating point (GMP) numbers
Definition: coeffs.h:42
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:616
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:358
#define ALLOC_RNUMBER()
Definition: coeffs.h:88
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:422
static FORCE_INLINE int n_GetChar(const coeffs r)
Return the characteristic of the coeff. domain.
Definition: coeffs.h:445
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:456
static FORCE_INLINE BOOLEAN nCoeff_is_Zp(const coeffs r)
Definition: coeffs.h:824
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:539
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_2toM(const coeffs r)
Definition: coeffs.h:748
#define FREE_RNUMBER(x)
Definition: coeffs.h:87
@ n_rep_gap_rat
(number), see longrat.h
Definition: coeffs.h:112
@ n_rep_gap_gmp
(), see rinteger.h, new impl.
Definition: coeffs.h:113
@ n_rep_float
(float), see shortfl.h
Definition: coeffs.h:117
@ n_rep_int
(int), see modulop.h
Definition: coeffs.h:111
@ n_rep_gmp_float
(gmp_float), see
Definition: coeffs.h:118
@ n_rep_gmp
(mpz_ptr), see rmodulon,h
Definition: coeffs.h:116
#define ALLOC0_RNUMBER()
Definition: coeffs.h:89
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:74
static FORCE_INLINE BOOLEAN nCoeff_is_R(const coeffs r)
Definition: coeffs.h:860
static FORCE_INLINE BOOLEAN nCoeff_is_long_C(const coeffs r)
Definition: coeffs.h:918
#define Print
Definition: emacs.cc:80
#define WarnS
Definition: emacs.cc:78
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
REvaluation E(1, terms.length(), IntRandom(25))
b *CanonicalForm B
Definition: facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
const ExtensionInfo & info
< [in] sqrfree poly
int j
Definition: facHensel.cc:110
bool isZero(const CFArray &A)
checks if entries of A are zero
CanonicalForm FACTORY_PUBLIC make_cf(const mpz_ptr n)
Definition: singext.cc:66
void FACTORY_PUBLIC gmp_numerator(const CanonicalForm &f, mpz_ptr result)
Definition: singext.cc:20
void FACTORY_PUBLIC gmp_denominator(const CanonicalForm &f, mpz_ptr result)
Definition: singext.cc:40
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define D(A)
Definition: gentable.cc:131
#define VAR
Definition: globaldefs.h:5
STATIC_VAR jList * Q
Definition: janet.cc:30
static number nlMapP(number from, const coeffs src, const coeffs dst)
Definition: longrat.cc:189
#define nlTest(a, r)
Definition: longrat.cc:87
void nlWriteFd(number n, const ssiInfo *d, const coeffs)
Definition: longrat.cc:3288
LINLINE void nlInpMult(number &a, number b, const coeffs r)
Definition: longrat.cc:2745
LINLINE BOOLEAN nlEqual(number a, number b, const coeffs r)
Definition: longrat.cc:2557
LINLINE number nlAdd(number la, number li, const coeffs r)
Definition: longrat.cc:2661
long nlInt(number &n, const coeffs r)
Definition: longrat.cc:701
static number nlLcm(number a, number b, const coeffs r)
Definition: longrat.cc:3409
number nlInit2(int i, int j, const coeffs r)
create a rational i/j (implicitly) over Q NOTE: make sure to use correct Q in debug mode
Definition: longrat.cc:2504
#define POW_2_28
Definition: longrat.cc:103
LINLINE number nl_Copy(number a, const coeffs r)
number nlInit2gmp(mpz_t i, mpz_t j, const coeffs r)
create a rational i/j (implicitly) over Q NOTE: make sure to use correct Q in debug mode
Definition: longrat.cc:2517
const char * nlRead(const char *s, number *a, const coeffs r)
Definition: longrat0.cc:31
void _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b)
Definition: longrat.cc:1943
LINLINE number nlSub(number la, number li, const coeffs r)
Definition: longrat.cc:2727
number nlIntMod(number a, number b, const coeffs r)
Definition: longrat.cc:979
number _nlCopy_NoImm(number a)
Definition: longrat.cc:1708
number _nlSub_aNoImm_OR_bNoImm(number a, number b)
Definition: longrat.cc:2086
LINLINE number nlCopy(number a, const coeffs r)
Definition: longrat.cc:2613
LINLINE number nlNeg(number za, const coeffs r)
Definition: longrat.cc:2642
number nlXExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
Definition: longrat.cc:2788
void nlPower(number x, int exp, number *lu, const coeffs r)
Definition: longrat.cc:1215
number nlQuotRem(number a, number b, number *r, const coeffs R)
Definition: longrat.cc:2840
number nlFarey(number nN, number nP, const coeffs CF)
Definition: longrat.cc:2926
LINLINE BOOLEAN nlIsOne(number a, const coeffs r)
Definition: longrat.cc:2584
#define mpz_isNeg(A)
Definition: longrat.cc:146
static number nlMapC(number from, const coeffs src, const coeffs dst)
Definition: longrat.cc:506
number nlNormalizeHelper(number a, number b, const coeffs r)
Definition: longrat.cc:1491
LINLINE void nlDelete(number *a, const coeffs r)
Definition: longrat.cc:2626
number nlMapZ(number from, const coeffs src, const coeffs dst)
Definition: longrat.cc:211
BOOLEAN nlGreaterZero(number za, const coeffs r)
Definition: longrat.cc:1268
number _nlNeg_NoImm(number a)
Definition: longrat.cc:1750
number nlModP(number q, const coeffs, const coeffs Zp)
Definition: longrat.cc:1538
LINLINE void nlInpAdd(number &a, number b, const coeffs r)
Definition: longrat.cc:2679
number nlExactDiv(number a, number b, const coeffs r)
Definition: longrat.cc:831
void mpz_mul_si(mpz_ptr r, mpz_srcptr s, long int si)
Definition: longrat.cc:177
VAR int n_SwitchChinRem
Definition: longrat.cc:3052
void nlMPZ(mpz_t m, number &n, const coeffs r)
Definition: longrat.cc:2779
number nlInvers(number a, const coeffs r)
Definition: longrat.cc:751
BOOLEAN nlIsUnit(number a, const coeffs)
Definition: longrat.cc:1096
void nlInpIntDiv(number &a, number b, const coeffs r)
Definition: longrat.cc:2908
static void nlNormalize_Gcd(number &x)
Definition: longrat.cc:1765
static number nlConvFactoryNSingN(const CanonicalForm f, const coeffs r)
Definition: longrat.cc:368
number nlChineseRemainderSym(number *x, number *q, int rl, BOOLEAN sym, CFArray &inv_cache, const coeffs CF)
Definition: longrat.cc:3053
int nlDivComp(number a, number b, const coeffs r)
Definition: longrat.cc:1054
void _nlDelete_NoImm(number *a)
Definition: longrat.cc:1729
#define LINLINE
Definition: longrat.cc:31
#define POW_2_28_32
Definition: longrat.cc:104
BOOLEAN nlInitChar(coeffs r, void *p)
Definition: longrat.cc:3433
number nlCopyMap(number a, const coeffs, const coeffs)
Definition: longrat.cc:2418
number nlExtGcd(number a, number b, number *s, number *t, const coeffs)
Definition: longrat.cc:2997
static number nlMapGMP(number from, const coeffs, const coeffs dst)
Definition: longrat.cc:206
LINLINE number nlMult(number a, number b, const coeffs r)
Definition: longrat.cc:2697
static number nlInitMPZ(mpz_t m, const coeffs)
Definition: longrat.cc:164
number nlIntDiv(number a, number b, const coeffs r)
Definition: longrat.cc:898
static void nlClearDenominators(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs cf)
Definition: longrat.cc:3188
static number nlMapLongR(number from, const coeffs src, const coeffs dst)
Definition: longrat.cc:425
LINLINE BOOLEAN nlIsZero(number za, const coeffs r)
Definition: longrat.cc:2593
number nlGetDenom(number &n, const coeffs r)
Definition: longrat.cc:1601
number nlGcd(number a, number b, const coeffs r)
Definition: longrat.cc:1305
number _nlMult_aImm_bImm_rNoImm(number a, number b)
Definition: longrat.cc:2297
number nlReadFd(const ssiInfo *d, const coeffs)
Definition: longrat.cc:3334
int nlSize(number a, const coeffs)
Definition: longrat.cc:672
number nlMapMachineInt(number from, const coeffs, const coeffs)
Definition: longrat.cc:223
nMapFunc nlSetMap(const coeffs src, const coeffs dst)
Definition: longrat.cc:2446
number nlBigInt(number &n)
char * nlCoeffName(const coeffs r)
Definition: longrat.cc:3282
static number nlShort3(number x)
Definition: longrat.cc:109
#define GCD_NORM_COND(OLD, NEW)
Definition: longrat.cc:1763
BOOLEAN nlDBTest(number a, const char *f, const int l)
number nlDiv(number a, number b, const coeffs r)
Definition: longrat.cc:1105
number nlRInit(long i)
Definition: longrat.cc:2490
BOOLEAN nlIsMOne(number a, const coeffs r)
Definition: longrat.cc:1293
static void nlClearContent(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs cf)
Definition: longrat.cc:3097
number _nlMult_aNoImm_OR_bNoImm(number a, number b)
Definition: longrat.cc:2310
LINLINE number nlInit(long i, const coeffs r)
Definition: longrat.cc:2566
number nlShort3_noinline(number x)
Definition: longrat.cc:159
number nlGetNumerator(number &n, const coeffs r)
Definition: longrat.cc:1630
number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
Definition: longrat.cc:1785
#define LONG
Definition: longrat.cc:105
BOOLEAN nlCoeffIsEqual(const coeffs r, n_coeffType n, void *p)
Definition: longrat.cc:3397
static CanonicalForm nlConvSingNFactoryN(number n, const BOOLEAN setChar, const coeffs)
Definition: longrat.cc:330
static number nlMapR(number from, const coeffs src, const coeffs dst)
Definition: longrat.cc:395
number nlGetUnit(number n, const coeffs cf)
Definition: longrat.cc:1065
coeffs nlQuot1(number c, const coeffs r)
Definition: longrat.cc:1071
BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
Definition: longrat.cc:1661
number nlShort1(number x)
Definition: longrat.cc:1426
#define MP_SMALL
Definition: longrat.cc:144
BOOLEAN nlGreater(number a, number b, const coeffs r)
Definition: longrat.cc:1278
void nlGMP(number &i, mpz_t n, const coeffs r)
Definition: longrat.cc:1580
void nlNormalize(number &x, const coeffs r)
Definition: longrat.cc:1447
BOOLEAN nlDivBy(number a, number b, const coeffs)
Definition: longrat.cc:1040
static int int_extgcd(int a, int b, int *u, int *x, int *v, int *y)
Definition: longrat.cc:1376
void nlWrite(number a, const coeffs r)
Definition: longrat0.cc:90
void nlInpGcd(number &a, number b, const coeffs r)
Definition: longrat.cc:2893
static number nlRandom(siRandProc p, number v2, number, const coeffs cf)
Definition: longrat.cc:3419
number nlMapQtoZ(number a, const coeffs src, const coeffs dst)
Definition: longrat.cc:2427
#define SR_INT
Definition: longrat.h:67
#define INT_TO_SR(INT)
Definition: longrat.h:68
#define SR_TO_INT(SR)
Definition: longrat.h:69
#define assume(x)
Definition: mod2.h:387
void dErrorBreak()
Definition: dError.cc:139
long npInt(number &n, const coeffs r)
Definition: modulop.cc:85
char * floatToStr(const gmp_float &r, const unsigned int oprec)
Definition: mpr_complex.cc:578
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
bool negative(N n)
Definition: ValueTraits.h:119
The main handler for Singular numbers which are suitable for Singular polynomials.
char * nEatLong(char *s, mpz_ptr i)
extracts a long integer from s, returns the rest
Definition: numbers.cc:667
const char *const nDivBy0
Definition: numbers.h:87
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omCheckIf(cond, test)
Definition: omAllocDecl.h:323
#define omCheckAddrSize(addr, size)
Definition: omAllocDecl.h:327
#define omFree(addr)
Definition: omAllocDecl.h:261
#define NULL
Definition: omList.c:12
int IsPrime(int p)
Definition: prime.cc:61
void Werror(const char *fmt,...)
Definition: reporter.cc:189
void s_readmpz(s_buff F, mpz_t a)
Definition: s_buff.cc:184
void s_readmpz_base(s_buff F, mpz_ptr a, int base)
Definition: s_buff.cc:209
int s_readint(s_buff F)
Definition: s_buff.cc:112
long s_readlong(s_buff F)
Definition: s_buff.cc:140
s_buff f_read
Definition: s_buff.h:22
FILE * f_write
Definition: s_buff.h:23
Definition: s_buff.h:21
SI_FLOAT nrFloat(number n)
Converts a n_R number into a float. Needed by Maps.
Definition: shortfl.cc:56
#define mpz_size1(A)
Definition: si_gmp.h:12
#define mpz_sgn1(A)
Definition: si_gmp.h:13
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
int(* siRandProc)()
Definition: sirandom.h:9
#define SR_HDL(A)
Definition: tgb.cc:35
int gcd(int a, int b)
Definition: walkSupport.cc:836