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