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