source: git/libpolys/coeffs/longrat.cc @ 3d752a

spielwiese
Last change on this file since 3d752a was 3d752a, checked in by Hans Schoenemann <hannes@…>, 5 years ago
debug helper: kTst_T, kTest_L, _nlDelete_NoImm
  • Property mode set to 100644
File size: 70.8 KB
Line 
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
36LINLINE BOOLEAN  nlEqual(number a, number b, const coeffs r);
37LINLINE number   nlInit(long i, const coeffs r);
38LINLINE BOOLEAN  nlIsOne(number a, const coeffs r);
39LINLINE BOOLEAN  nlIsZero(number za, const coeffs r);
40LINLINE number   nlCopy(number a, const coeffs r);
41LINLINE number   nl_Copy(number a, const coeffs r);
42LINLINE void     nlDelete(number *a, const coeffs r);
43LINLINE number   nlNeg(number za, const coeffs r);
44LINLINE number   nlAdd(number la, number li, const coeffs r);
45LINLINE number   nlSub(number la, number li, const coeffs r);
46LINLINE number   nlMult(number a, number b, const coeffs r);
47LINLINE void nlInpAdd(number &a, number b, const coeffs r);
48LINLINE void nlInpMult(number &a, number b, const coeffs r);
49
50number 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
56void     nlNormalize(number &x, const coeffs r);
57
58number   nlGcd(number a, number b, const coeffs r);
59number nlExtGcd(number a, number b, number *s, number *t, const coeffs);
60number   nlNormalizeHelper(number a, number b, const coeffs r);   /*special routine !*/
61BOOLEAN  nlGreater(number a, number b, const coeffs r);
62BOOLEAN  nlIsMOne(number a, const coeffs r);
63long     nlInt(number &n, const coeffs r);
64number   nlBigInt(number &n);
65
66#ifdef HAVE_RINGS
67number nlMapGMP(number from, const coeffs src, const coeffs dst);
68#endif
69
70BOOLEAN  nlGreaterZero(number za, const coeffs r);
71number   nlInvers(number a, const coeffs r);
72number   nlDiv(number a, number b, const coeffs r);
73number   nlExactDiv(number a, number b, const coeffs r);
74number   nlIntDiv(number a, number b, const coeffs r);
75number   nlIntMod(number a, number b, const coeffs r);
76void     nlPower(number x, int exp, number *lu, const coeffs r);
77const char *   nlRead (const char *s, number *a, const coeffs r);
78void     nlWrite(number a, const coeffs r);
79
80void     nlCoeffWrite(const coeffs r, BOOLEAN details);
81number   nlFarey(number nN, number nP, const coeffs CF);
82
83#ifdef LDEBUG
84BOOLEAN  nlDBTest(number a, const char *f, const int l);
85#endif
86
87nMapFunc nlSetMap(const coeffs src, const coeffs dst);
88
89// in-place operations
90void nlInpIntDiv(number &a, number b, const coeffs r);
91
92#ifdef LDEBUG
93#define nlTest(a, r) nlDBTest(a,__FILE__,__LINE__, r)
94BOOLEAN 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
115static 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
156void    _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
165number 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)
172void 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
184static 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
193static number nlMapLongR(number from, const coeffs src, const coeffs dst);
194static 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*/
201number 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
213number 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*/
225number 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
240BOOLEAN 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
332static CanonicalForm nlConvSingNFactoryN( number n, const BOOLEAN setChar, const coeffs /*r*/ )
333{
334  if (setChar) setCharacteristic( 0 );
335
336  CanonicalForm term;
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
368number nlRInit (long i);
369
370static 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
397static 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
427static 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
570int 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*/
599long 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*/
631number 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*/
649number 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*/
729number   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*/
796number 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*/
877number 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
938BOOLEAN 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
952int 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
963number  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
969coeffs 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
994BOOLEAN 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*/
1003number 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*/
1113void 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*/
1166BOOLEAN 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*/
1176BOOLEAN 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*/
1191BOOLEAN 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*/
1203number 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}
1274static 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
1324number 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*/
1345void 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*/
1389number 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)
1436number 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*/
1478void 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*/
1499number   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*/
1528number   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 *******************************************************************/
1559BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
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
1606number _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
1627void _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    #ifdef LDEBUG
1642    memset(*a,0,sizeof(**a));
1643    #endif
1644    FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
1645  }
1646}
1647
1648number _nlNeg_NoImm(number a)
1649{
1650  {
1651    mpz_neg(a->z,a->z);
1652    if (a->s==3)
1653    {
1654      a=nlShort3(a);
1655    }
1656  }
1657  return a;
1658}
1659
1660// conditio to use nlNormalize_Gcd in intermediate computations:
1661#define GCD_NORM_COND(OLD,NEW) (mpz_size1(NEW->z)>mpz_size1(OLD->z))
1662
1663static void nlNormalize_Gcd(number &x)
1664{
1665  mpz_t gcd;
1666  mpz_init(gcd);
1667  mpz_gcd(gcd,x->z,x->n);
1668  x->s=1;
1669  if (mpz_cmp_si(gcd,1L)!=0)
1670  {
1671    mpz_divexact(x->z,x->z,gcd);
1672    mpz_divexact(x->n,x->n,gcd);
1673    if (mpz_cmp_si(x->n,1L)==0)
1674    {
1675      mpz_clear(x->n);
1676      x->s=3;
1677      x=nlShort3_noinline(x);
1678    }
1679  }
1680  mpz_clear(gcd);
1681}
1682
1683number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1684{
1685  number u=ALLOC_RNUMBER();
1686#if defined(LDEBUG)
1687  u->debug=123456;
1688#endif
1689  mpz_init(u->z);
1690  if (SR_HDL(b) & SR_INT)
1691  {
1692    number x=a;
1693    a=b;
1694    b=x;
1695  }
1696  if (SR_HDL(a) & SR_INT)
1697  {
1698    switch (b->s)
1699    {
1700      case 0:
1701      case 1:/* a:short, b:1 */
1702      {
1703        mpz_t x;
1704        mpz_init(x);
1705        mpz_mul_si(x,b->n,SR_TO_INT(a));
1706        mpz_add(u->z,b->z,x);
1707        mpz_clear(x);
1708        if (mpz_sgn1(u->z)==0)
1709        {
1710          mpz_clear(u->z);
1711          FREE_RNUMBER(u);
1712          return INT_TO_SR(0);
1713        }
1714        if (mpz_cmp(u->z,b->n)==0)
1715        {
1716          mpz_clear(u->z);
1717          FREE_RNUMBER(u);
1718          return INT_TO_SR(1);
1719        }
1720        mpz_init_set(u->n,b->n);
1721        u->s = 0;
1722        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1723        break;
1724      }
1725      case 3:
1726      {
1727        if (((long)a)>0L)
1728          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1729        else
1730          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1731        u->s = 3;
1732        u=nlShort3(u);
1733        break;
1734      }
1735    }
1736  }
1737  else
1738  {
1739    switch (a->s)
1740    {
1741      case 0:
1742      case 1:
1743      {
1744        switch(b->s)
1745        {
1746          case 0:
1747          case 1:
1748          {
1749            mpz_t x;
1750            mpz_init(x);
1751
1752            mpz_mul(x,b->z,a->n);
1753            mpz_mul(u->z,a->z,b->n);
1754            mpz_add(u->z,u->z,x);
1755            mpz_clear(x);
1756
1757            if (mpz_sgn1(u->z)==0)
1758            {
1759              mpz_clear(u->z);
1760              FREE_RNUMBER(u);
1761              return INT_TO_SR(0);
1762            }
1763            mpz_init(u->n);
1764            mpz_mul(u->n,a->n,b->n);
1765            if (mpz_cmp(u->z,u->n)==0)
1766            {
1767               mpz_clear(u->z);
1768               mpz_clear(u->n);
1769               FREE_RNUMBER(u);
1770               return INT_TO_SR(1);
1771            }
1772            u->s = 0;
1773            if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1774            break;
1775          }
1776          case 3: /* a:1 b:3 */
1777          {
1778            mpz_mul(u->z,b->z,a->n);
1779            mpz_add(u->z,u->z,a->z);
1780            if (mpz_sgn1(u->z)==0)
1781            {
1782              mpz_clear(u->z);
1783              FREE_RNUMBER(u);
1784              return INT_TO_SR(0);
1785            }
1786            if (mpz_cmp(u->z,a->n)==0)
1787            {
1788              mpz_clear(u->z);
1789              FREE_RNUMBER(u);
1790              return INT_TO_SR(1);
1791            }
1792            mpz_init_set(u->n,a->n);
1793            u->s = 0;
1794            if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
1795            break;
1796          }
1797        } /*switch (b->s) */
1798        break;
1799      }
1800      case 3:
1801      {
1802        switch(b->s)
1803        {
1804          case 0:
1805          case 1:/* a:3, b:1 */
1806          {
1807            mpz_mul(u->z,a->z,b->n);
1808            mpz_add(u->z,u->z,b->z);
1809            if (mpz_sgn1(u->z)==0)
1810            {
1811              mpz_clear(u->z);
1812              FREE_RNUMBER(u);
1813              return INT_TO_SR(0);
1814            }
1815            if (mpz_cmp(u->z,b->n)==0)
1816            {
1817              mpz_clear(u->z);
1818              FREE_RNUMBER(u);
1819              return INT_TO_SR(1);
1820            }
1821            mpz_init_set(u->n,b->n);
1822            u->s = 0;
1823            if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1824            break;
1825          }
1826          case 3:
1827          {
1828            mpz_add(u->z,a->z,b->z);
1829            u->s = 3;
1830            u=nlShort3(u);
1831            break;
1832          }
1833        }
1834        break;
1835      }
1836    }
1837  }
1838  return u;
1839}
1840
1841void _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b)
1842{
1843  if (SR_HDL(b) & SR_INT)
1844  {
1845    switch (a->s)
1846    {
1847      case 0:
1848      case 1:/* b:short, a:1 */
1849      {
1850        mpz_t x;
1851        mpz_init(x);
1852        mpz_mul_si(x,a->n,SR_TO_INT(b));
1853        mpz_add(a->z,a->z,x);
1854        mpz_clear(x);
1855        nlNormalize_Gcd(a);
1856        break;
1857      }
1858      case 3:
1859      {
1860        if (((long)b)>0L)
1861          mpz_add_ui(a->z,a->z,SR_TO_INT(b));
1862        else
1863          mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
1864        a->s = 3;
1865        a=nlShort3_noinline(a);
1866        break;
1867      }
1868    }
1869    return;
1870  }
1871  else if (SR_HDL(a) & SR_INT)
1872  {
1873    number u=ALLOC_RNUMBER();
1874    #if defined(LDEBUG)
1875    u->debug=123456;
1876    #endif
1877    mpz_init(u->z);
1878    switch (b->s)
1879    {
1880      case 0:
1881      case 1:/* a:short, b:1 */
1882      {
1883        mpz_t x;
1884        mpz_init(x);
1885
1886        mpz_mul_si(x,b->n,SR_TO_INT(a));
1887        mpz_add(u->z,b->z,x);
1888        mpz_clear(x);
1889        // result cannot be 0, if coeffs are normalized
1890        mpz_init_set(u->n,b->n);
1891        u->s=0;
1892        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1893        else { u=nlShort1(u); }
1894        break;
1895      }
1896      case 3:
1897      {
1898        if (((long)a)>0L)
1899          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1900        else
1901          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1902        // result cannot be 0, if coeffs are normalized
1903        u->s = 3;
1904        u=nlShort3_noinline(u);
1905        break;
1906      }
1907    }
1908    a=u;
1909  }
1910  else
1911  {
1912    switch (a->s)
1913    {
1914      case 0:
1915      case 1:
1916      {
1917        switch(b->s)
1918        {
1919          case 0:
1920          case 1: /* a:1 b:1 */
1921          {
1922            mpz_t x;
1923            mpz_t y;
1924            mpz_init(x);
1925            mpz_init(y);
1926            mpz_mul(x,b->z,a->n);
1927            mpz_mul(y,a->z,b->n);
1928            mpz_add(a->z,x,y);
1929            mpz_clear(x);
1930            mpz_clear(y);
1931            mpz_mul(a->n,a->n,b->n);
1932            a->s=0;
1933            if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
1934            else { a=nlShort1(a);}
1935            break;
1936          }
1937          case 3: /* a:1 b:3 */
1938          {
1939            mpz_t x;
1940            mpz_init(x);
1941            mpz_mul(x,b->z,a->n);
1942            mpz_add(a->z,a->z,x);
1943            mpz_clear(x);
1944            a->s=0;
1945            if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
1946            else { a=nlShort1(a);}
1947            break;
1948          }
1949        } /*switch (b->s) */
1950        break;
1951      }
1952      case 3:
1953      {
1954        switch(b->s)
1955        {
1956          case 0:
1957          case 1:/* a:3, b:1 */
1958          {
1959            mpz_t x;
1960            mpz_init(x);
1961            mpz_mul(x,a->z,b->n);
1962            mpz_add(a->z,b->z,x);
1963            mpz_clear(x);
1964            mpz_init_set(a->n,b->n);
1965            a->s=0;
1966            if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
1967            else { a=nlShort1(a);}
1968            break;
1969          }
1970          case 3:
1971          {
1972            mpz_add(a->z,a->z,b->z);
1973            a->s = 3;
1974            a=nlShort3_noinline(a);
1975            break;
1976          }
1977        }
1978        break;
1979      }
1980    }
1981  }
1982}
1983
1984number _nlSub_aNoImm_OR_bNoImm(number a, number b)
1985{
1986  number u=ALLOC_RNUMBER();
1987#if defined(LDEBUG)
1988  u->debug=123456;
1989#endif
1990  mpz_init(u->z);
1991  if (SR_HDL(a) & SR_INT)
1992  {
1993    switch (b->s)
1994    {
1995      case 0:
1996      case 1:/* a:short, b:1 */
1997      {
1998        mpz_t x;
1999        mpz_init(x);
2000        mpz_mul_si(x,b->n,SR_TO_INT(a));
2001        mpz_sub(u->z,x,b->z);
2002        mpz_clear(x);
2003        if (mpz_sgn1(u->z)==0)
2004        {
2005          mpz_clear(u->z);
2006          FREE_RNUMBER(u);
2007          return INT_TO_SR(0);
2008        }
2009        if (mpz_cmp(u->z,b->n)==0)
2010        {
2011          mpz_clear(u->z);
2012          FREE_RNUMBER(u);
2013          return INT_TO_SR(1);
2014        }
2015        mpz_init_set(u->n,b->n);
2016        u->s=0;
2017        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2018        break;
2019      }
2020      case 3:
2021      {
2022        if (((long)a)>0L)
2023        {
2024          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
2025          mpz_neg(u->z,u->z);
2026        }
2027        else
2028        {
2029          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
2030          mpz_neg(u->z,u->z);
2031        }
2032        u->s = 3;
2033        u=nlShort3(u);
2034        break;
2035      }
2036    }
2037  }
2038  else if (SR_HDL(b) & SR_INT)
2039  {
2040    switch (a->s)
2041    {
2042      case 0:
2043      case 1:/* b:short, a:1 */
2044      {
2045        mpz_t x;
2046        mpz_init(x);
2047        mpz_mul_si(x,a->n,SR_TO_INT(b));
2048        mpz_sub(u->z,a->z,x);
2049        mpz_clear(x);
2050        if (mpz_sgn1(u->z)==0)
2051        {
2052          mpz_clear(u->z);
2053          FREE_RNUMBER(u);
2054          return INT_TO_SR(0);
2055        }
2056        if (mpz_cmp(u->z,a->n)==0)
2057        {
2058          mpz_clear(u->z);
2059          FREE_RNUMBER(u);
2060          return INT_TO_SR(1);
2061        }
2062        mpz_init_set(u->n,a->n);
2063        u->s=0;
2064        if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2065        break;
2066      }
2067      case 3:
2068      {
2069        if (((long)b)>0L)
2070        {
2071          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
2072        }
2073        else
2074        {
2075          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
2076        }
2077        u->s = 3;
2078        u=nlShort3(u);
2079        break;
2080      }
2081    }
2082  }
2083  else
2084  {
2085    switch (a->s)
2086    {
2087      case 0:
2088      case 1:
2089      {
2090        switch(b->s)
2091        {
2092          case 0:
2093          case 1:
2094          {
2095            mpz_t x;
2096            mpz_t y;
2097            mpz_init(x);
2098            mpz_init(y);
2099            mpz_mul(x,b->z,a->n);
2100            mpz_mul(y,a->z,b->n);
2101            mpz_sub(u->z,y,x);
2102            mpz_clear(x);
2103            mpz_clear(y);
2104            if (mpz_sgn1(u->z)==0)
2105            {
2106              mpz_clear(u->z);
2107              FREE_RNUMBER(u);
2108              return INT_TO_SR(0);
2109            }
2110            mpz_init(u->n);
2111            mpz_mul(u->n,a->n,b->n);
2112            if (mpz_cmp(u->z,u->n)==0)
2113            {
2114              mpz_clear(u->z);
2115              mpz_clear(u->n);
2116              FREE_RNUMBER(u);
2117              return INT_TO_SR(1);
2118            }
2119            u->s=0;
2120            if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2121            break;
2122          }
2123          case 3: /* a:1, b:3 */
2124          {
2125            mpz_t x;
2126            mpz_init(x);
2127            mpz_mul(x,b->z,a->n);
2128            mpz_sub(u->z,a->z,x);
2129            mpz_clear(x);
2130            if (mpz_sgn1(u->z)==0)
2131            {
2132              mpz_clear(u->z);
2133              FREE_RNUMBER(u);
2134              return INT_TO_SR(0);
2135            }
2136            if (mpz_cmp(u->z,a->n)==0)
2137            {
2138              mpz_clear(u->z);
2139              FREE_RNUMBER(u);
2140              return INT_TO_SR(1);
2141            }
2142            mpz_init_set(u->n,a->n);
2143            u->s=0;
2144            if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2145            break;
2146          }
2147        }
2148        break;
2149      }
2150      case 3:
2151      {
2152        switch(b->s)
2153        {
2154          case 0:
2155          case 1: /* a:3, b:1 */
2156          {
2157            mpz_t x;
2158            mpz_init(x);
2159            mpz_mul(x,a->z,b->n);
2160            mpz_sub(u->z,x,b->z);
2161            mpz_clear(x);
2162            if (mpz_sgn1(u->z)==0)
2163            {
2164              mpz_clear(u->z);
2165              FREE_RNUMBER(u);
2166              return INT_TO_SR(0);
2167            }
2168            if (mpz_cmp(u->z,b->n)==0)
2169            {
2170              mpz_clear(u->z);
2171              FREE_RNUMBER(u);
2172              return INT_TO_SR(1);
2173            }
2174            mpz_init_set(u->n,b->n);
2175            u->s=0;
2176            if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2177            break;
2178          }
2179          case 3: /* a:3 , b:3 */
2180          {
2181            mpz_sub(u->z,a->z,b->z);
2182            u->s = 3;
2183            u=nlShort3(u);
2184            break;
2185          }
2186        }
2187        break;
2188      }
2189    }
2190  }
2191  return u;
2192}
2193
2194// a and b are intermediate, but a*b not
2195number _nlMult_aImm_bImm_rNoImm(number a, number b)
2196{
2197  number u=ALLOC_RNUMBER();
2198#if defined(LDEBUG)
2199  u->debug=123456;
2200#endif
2201  u->s=3;
2202  mpz_init_set_si(u->z,SR_TO_INT(a));
2203  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
2204  return u;
2205}
2206
2207// a or b are not immediate
2208number _nlMult_aNoImm_OR_bNoImm(number a, number b)
2209{
2210  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2211  number u=ALLOC_RNUMBER();
2212#if defined(LDEBUG)
2213  u->debug=123456;
2214#endif
2215  mpz_init(u->z);
2216  if (SR_HDL(b) & SR_INT)
2217  {
2218    number x=a;
2219    a=b;
2220    b=x;
2221  }
2222  if (SR_HDL(a) & SR_INT)
2223  {
2224    u->s=b->s;
2225    if (u->s==1) u->s=0;
2226    if (((long)a)>0L)
2227    {
2228      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2229    }
2230    else
2231    {
2232      if (a==INT_TO_SR(-1))
2233      {
2234        mpz_set(u->z,b->z);
2235        mpz_neg(u->z,u->z);
2236        u->s=b->s;
2237      }
2238      else
2239      {
2240        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2241        mpz_neg(u->z,u->z);
2242      }
2243    }
2244    if (u->s<2)
2245    {
2246      if (mpz_cmp(u->z,b->n)==0)
2247      {
2248        mpz_clear(u->z);
2249        FREE_RNUMBER(u);
2250        return INT_TO_SR(1);
2251      }
2252      mpz_init_set(u->n,b->n);
2253      if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2254    }
2255    else //u->s==3
2256    {
2257      u=nlShort3(u);
2258    }
2259  }
2260  else
2261  {
2262    mpz_mul(u->z,a->z,b->z);
2263    u->s = 0;
2264    if(a->s==3)
2265    {
2266      if(b->s==3)
2267      {
2268        u->s = 3;
2269      }
2270      else
2271      {
2272        if (mpz_cmp(u->z,b->n)==0)
2273        {
2274          mpz_clear(u->z);
2275          FREE_RNUMBER(u);
2276          return INT_TO_SR(1);
2277        }
2278        mpz_init_set(u->n,b->n);
2279        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2280      }
2281    }
2282    else
2283    {
2284      if(b->s==3)
2285      {
2286        if (mpz_cmp(u->z,a->n)==0)
2287        {
2288          mpz_clear(u->z);
2289          FREE_RNUMBER(u);
2290          return INT_TO_SR(1);
2291        }
2292        mpz_init_set(u->n,a->n);
2293        if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2294      }
2295      else
2296      {
2297        mpz_init(u->n);
2298        mpz_mul(u->n,a->n,b->n);
2299        if (mpz_cmp(u->z,u->n)==0)
2300        {
2301          mpz_clear(u->z);
2302          mpz_clear(u->n);
2303          FREE_RNUMBER(u);
2304          return INT_TO_SR(1);
2305        }
2306        if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2307      }
2308    }
2309  }
2310  return u;
2311}
2312
2313/*2
2314* copy a to b for mapping
2315*/
2316number nlCopyMap(number a, const coeffs /*src*/, const coeffs /*dst*/)
2317{
2318  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2319  {
2320    return a;
2321  }
2322  return _nlCopy_NoImm(a);
2323}
2324
2325nMapFunc nlSetMap(const coeffs src, const coeffs /*dst*/)
2326{
2327  if (src->rep==n_rep_gap_rat)  /*Q, coeffs_BIGINT */
2328  {
2329    return ndCopyMap;
2330  }
2331  if ((src->rep==n_rep_int) && nCoeff_is_Zp(src))
2332  {
2333    return nlMapP;
2334  }
2335  if ((src->rep==n_rep_float) && nCoeff_is_R(src))
2336  {
2337    return nlMapR;
2338  }
2339  if ((src->rep==n_rep_gmp_float) && nCoeff_is_long_R(src))
2340  {
2341    return nlMapLongR; /* long R -> Q */
2342  }
2343#ifdef HAVE_RINGS
2344  if (src->rep==n_rep_gmp) // nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
2345  {
2346    return nlMapGMP;
2347  }
2348  if (src->rep==n_rep_gap_gmp)
2349  {
2350    return nlMapZ;
2351  }
2352  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src))
2353  {
2354    return nlMapMachineInt;
2355  }
2356#endif
2357  return NULL;
2358}
2359/*2
2360* z := i
2361*/
2362number nlRInit (long i)
2363{
2364  number z=ALLOC_RNUMBER();
2365#if defined(LDEBUG)
2366  z->debug=123456;
2367#endif
2368  mpz_init_set_si(z->z,i);
2369  z->s = 3;
2370  return z;
2371}
2372
2373/*2
2374* z := i/j
2375*/
2376number nlInit2 (int i, int j, const coeffs r)
2377{
2378  number z=ALLOC_RNUMBER();
2379#if defined(LDEBUG)
2380  z->debug=123456;
2381#endif
2382  mpz_init_set_si(z->z,(long)i);
2383  mpz_init_set_si(z->n,(long)j);
2384  z->s = 0;
2385  nlNormalize(z,r);
2386  return z;
2387}
2388
2389number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
2390{
2391  number z=ALLOC_RNUMBER();
2392#if defined(LDEBUG)
2393  z->debug=123456;
2394#endif
2395  mpz_init_set(z->z,i);
2396  mpz_init_set(z->n,j);
2397  z->s = 0;
2398  nlNormalize(z,r);
2399  return z;
2400}
2401
2402#else // DO_LINLINE
2403
2404// declare immedate routines
2405number nlRInit (long i);
2406BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2407number  _nlCopy_NoImm(number a);
2408number  _nlNeg_NoImm(number a);
2409number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2410void    _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b);
2411number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2412number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2413number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2414
2415#endif
2416
2417/***************************************************************
2418 *
2419 * Routines which might be inlined by p_Numbers.h
2420 *
2421 *******************************************************************/
2422#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2423
2424// routines which are always inlined/static
2425
2426/*2
2427* a = b ?
2428*/
2429LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
2430{
2431  nlTest(a, r);
2432  nlTest(b, r);
2433// short - short
2434  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2435  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2436}
2437
2438LINLINE number nlInit (long i, const coeffs r)
2439{
2440  number n;
2441  #if MAX_NUM_SIZE == 60
2442  if (((i << 3) >> 3) == i) n=INT_TO_SR(i);
2443  else                      n=nlRInit(i);
2444  #else
2445  LONG ii=(LONG)i;
2446  if ( ((((long)ii)==i) && ((ii << 3) >> 3) == ii )) n=INT_TO_SR(ii);
2447  else                                             n=nlRInit(i);
2448  #endif
2449  nlTest(n, r);
2450  return n;
2451}
2452
2453/*2
2454* a == 1 ?
2455*/
2456LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
2457{
2458#ifdef LDEBUG
2459  if (a==NULL) return FALSE;
2460  nlTest(a, r);
2461#endif
2462  return (a==INT_TO_SR(1));
2463}
2464
2465LINLINE BOOLEAN nlIsZero (number a, const coeffs)
2466{
2467  #if 0
2468  if (a==INT_TO_SR(0)) return TRUE;
2469  if ((SR_HDL(a) & SR_INT)||(a==NULL)) return FALSE;
2470  if (mpz_cmp_si(a->z,0L)==0)
2471  {
2472    printf("gmp-0 in nlIsZero\n");
2473    dErrorBreak();
2474    return TRUE;
2475  }
2476  return FALSE;
2477  #else
2478  return (a==NULL)|| (a==INT_TO_SR(0));
2479  #endif
2480}
2481
2482/*2
2483* copy a to b
2484*/
2485LINLINE number nlCopy(number a, const coeffs)
2486{
2487  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2488  {
2489    return a;
2490  }
2491  return _nlCopy_NoImm(a);
2492}
2493
2494
2495/*2
2496* delete a
2497*/
2498LINLINE void nlDelete (number * a, const coeffs r)
2499{
2500  if (*a!=NULL)
2501  {
2502    nlTest(*a, r);
2503    if ((SR_HDL(*a) & SR_INT)==0)
2504    {
2505      _nlDelete_NoImm(a);
2506    }
2507    *a=NULL;
2508  }
2509}
2510
2511/*2
2512* za:= - za
2513*/
2514LINLINE number nlNeg (number a, const coeffs R)
2515{
2516  nlTest(a, R);
2517  if(SR_HDL(a) &SR_INT)
2518  {
2519    LONG r=SR_TO_INT(a);
2520    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2521    else               a=INT_TO_SR(-r);
2522    return a;
2523  }
2524  a = _nlNeg_NoImm(a);
2525  nlTest(a, R);
2526  return a;
2527
2528}
2529
2530/*2
2531* u:= a + b
2532*/
2533LINLINE number nlAdd (number a, number b, const coeffs R)
2534{
2535  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2536  {
2537    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2538    if ( ((r << 1) >> 1) == r )
2539      return (number)(long)r;
2540    else
2541      return nlRInit(SR_TO_INT(r));
2542  }
2543  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2544  nlTest(u, R);
2545  return u;
2546}
2547
2548number nlShort1(number a);
2549number nlShort3_noinline(number x);
2550
2551LINLINE void nlInpAdd(number &a, number b, const coeffs r)
2552{
2553  // a=a+b
2554  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2555  {
2556    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2557    if ( ((r << 1) >> 1) == r )
2558      a=(number)(long)r;
2559    else
2560      a=nlRInit(SR_TO_INT(r));
2561  }
2562  else
2563  {
2564    _nlInpAdd_aNoImm_OR_bNoImm(a,b);
2565    nlTest(a,r);
2566  }
2567}
2568
2569LINLINE number nlMult (number a, number b, const coeffs R)
2570{
2571  nlTest(a, R);
2572  nlTest(b, R);
2573  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2574  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2575  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2576  {
2577    LONG r=(LONG)((unsigned LONG)(SR_HDL(a)-1L))*((unsigned LONG)(SR_HDL(b)>>1));
2578    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2579    {
2580      number u=((number) ((r>>1)+SR_INT));
2581      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2582      return nlRInit(SR_HDL(u)>>2);
2583    }
2584    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2585    nlTest(u, R);
2586    return u;
2587
2588  }
2589  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2590  nlTest(u, R);
2591  return u;
2592
2593}
2594
2595
2596/*2
2597* u:= a - b
2598*/
2599LINLINE number nlSub (number a, number b, const coeffs r)
2600{
2601  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2602  {
2603    LONG r=SR_HDL(a)-SR_HDL(b)+1;
2604    if ( ((r << 1) >> 1) == r )
2605    {
2606      return (number)(long)r;
2607    }
2608    else
2609      return nlRInit(SR_TO_INT(r));
2610  }
2611  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2612  nlTest(u, r);
2613  return u;
2614
2615}
2616
2617LINLINE void nlInpMult(number &a, number b, const coeffs r)
2618{
2619  number aa=a;
2620  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2621  {
2622    number n=nlMult(aa,b,r);
2623    nlDelete(&a,r);
2624    a=n;
2625  }
2626  else
2627  {
2628    mpz_mul(aa->z,a->z,b->z);
2629    if (aa->s==3)
2630    {
2631      if(b->s!=3)
2632      {
2633        mpz_init_set(a->n,b->n);
2634        a->s=0;
2635      }
2636    }
2637    else
2638    {
2639      if(b->s!=3)
2640      {
2641        mpz_mul(a->n,a->n,b->n);
2642      }
2643      a->s=0;
2644    }
2645  }
2646}
2647#endif // DO_LINLINE
2648
2649#ifndef P_NUMBERS_H
2650
2651static void nlMPZ(mpz_t m, number &n, const coeffs r)
2652{
2653  nlTest(n, r);
2654  nlNormalize(n, r);
2655  if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n));     /* n fits in an int */
2656  else             mpz_init_set(m, (mpz_ptr)n->z);
2657}
2658
2659
2660static number nlInitMPZ(mpz_t m, const coeffs)
2661{
2662  number z = ALLOC_RNUMBER();
2663  z->s = 3;
2664  #ifdef LDEBUG
2665  z->debug=123456;
2666  #endif
2667  mpz_init_set(z->z, m);
2668  z=nlShort3(z);
2669  return z;
2670}
2671
2672number  nlXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
2673{
2674  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2675  {
2676    int uu, vv, x, y;
2677    int g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &uu, &vv, &x, &y);
2678    *s = INT_TO_SR(uu);
2679    *t = INT_TO_SR(vv);
2680    *u = INT_TO_SR(x);
2681    *v = INT_TO_SR(y);
2682    return INT_TO_SR(g);
2683  }
2684  else
2685  {
2686    mpz_t aa, bb;
2687    if (SR_HDL(a) & SR_INT)
2688    {
2689      mpz_init_set_si(aa, SR_TO_INT(a));
2690    }
2691    else
2692    {
2693      mpz_init_set(aa, a->z);
2694    }
2695    if (SR_HDL(b) & SR_INT)
2696    {
2697      mpz_init_set_si(bb, SR_TO_INT(b));
2698    }
2699    else
2700    {
2701      mpz_init_set(bb, b->z);
2702    }
2703    mpz_t erg; mpz_t bs; mpz_t bt;
2704    mpz_init(erg);
2705    mpz_init(bs);
2706    mpz_init(bt);
2707
2708    mpz_gcdext(erg, bs, bt, aa, bb);
2709
2710    mpz_div(aa, aa, erg);
2711    *u=nlInitMPZ(bb,r);
2712    *u=nlNeg(*u,r);
2713    *v=nlInitMPZ(aa,r);
2714
2715    mpz_clear(aa);
2716    mpz_clear(bb);
2717
2718    *s = nlInitMPZ(bs,r);
2719    *t = nlInitMPZ(bt,r);
2720    return nlInitMPZ(erg,r);
2721  }
2722}
2723
2724number nlQuotRem (number a, number b, number * r, const coeffs R)
2725{
2726  assume(SR_TO_INT(b)!=0);
2727  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2728  {
2729    if (r!=NULL)
2730      *r = INT_TO_SR(SR_TO_INT(a) % SR_TO_INT(b));
2731    return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
2732  }
2733  else if (SR_HDL(a) & SR_INT)
2734  {
2735    // -2^xx / 2^xx
2736    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
2737    {
2738      if (r!=NULL) *r=INT_TO_SR(0);
2739      return nlRInit(POW_2_28);
2740    }
2741    //a is small, b is not, so q=0, r=a
2742    if (r!=NULL)
2743      *r = a;
2744    return INT_TO_SR(0);
2745  }
2746  else if (SR_HDL(b) & SR_INT)
2747  {
2748    unsigned long rr;
2749    mpz_t qq;
2750    mpz_init(qq);
2751    mpz_t rrr;
2752    mpz_init(rrr);
2753    rr = mpz_divmod_ui(qq, rrr, a->z, (unsigned long)ABS(SR_TO_INT(b)));
2754    mpz_clear(rrr);
2755
2756    if (r!=NULL)
2757      *r = INT_TO_SR(rr);
2758    if (SR_TO_INT(b)<0)
2759    {
2760      mpz_neg(qq, qq);
2761    }
2762    return nlInitMPZ(qq,R);
2763  }
2764  mpz_t qq,rr;
2765  mpz_init(qq);
2766  mpz_init(rr);
2767  mpz_divmod(qq, rr, a->z, b->z);
2768  if (r!=NULL)
2769    *r = nlInitMPZ(rr,R);
2770  else
2771  {
2772    mpz_clear(rr);
2773  }
2774  return nlInitMPZ(qq,R);
2775}
2776
2777void nlInpGcd(number &a, number b, const coeffs r)
2778{
2779  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2780  {
2781    number n=nlGcd(a,b,r);
2782    nlDelete(&a,r);
2783    a=n;
2784  }
2785  else
2786  {
2787    mpz_gcd(a->z,a->z,b->z);
2788    a=nlShort3_noinline(a);
2789  }
2790}
2791
2792void nlInpIntDiv(number &a, number b, const coeffs r)
2793{
2794  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2795  {
2796    number n=nlIntDiv(a,b, r);
2797    nlDelete(&a,r);
2798    a=n;
2799  }
2800  else
2801  {
2802    number rr=nlIntMod(a,b,r);
2803    if (SR_HDL(rr) &  SR_INT) mpz_sub_ui(a->z,a->z,SR_TO_INT(rr));
2804    else                      mpz_sub(a->z,a->z,rr->z);
2805    mpz_divexact(a->z,a->z,b->z);
2806    a=nlShort3_noinline(a);
2807  }
2808}
2809
2810number nlFarey(number nN, number nP, const coeffs r)
2811{
2812  mpz_t A,B,C,D,E,N,P,tmp;
2813  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2814  else                     mpz_init_set(P,nP->z);
2815  const mp_bitcnt_t bits=2*(mpz_size1(P)+1)*GMP_LIMB_BITS;
2816  mpz_init2(N,bits);
2817  if (SR_HDL(nN) & SR_INT) mpz_set_si(N,SR_TO_INT(nN));
2818  else                     mpz_set(N,nN->z);
2819  assume(!mpz_isNeg(P));
2820  if (mpz_isNeg(N))  mpz_add(N,N,P);
2821  mpz_init2(A,bits); mpz_set_ui(A,0L);
2822  mpz_init2(B,bits); mpz_set_ui(B,1L);
2823  mpz_init2(C,bits); mpz_set_ui(C,0L);
2824  mpz_init2(D,bits);
2825  mpz_init2(E,bits); mpz_set(E,P);
2826  mpz_init2(tmp,bits);
2827  number z=INT_TO_SR(0);
2828  while(mpz_sgn1(N)!=0)
2829  {
2830    mpz_mul(tmp,N,N);
2831    mpz_add(tmp,tmp,tmp);
2832    if (mpz_cmp(tmp,P)<0)
2833    {
2834       if (mpz_isNeg(B))
2835       {
2836         mpz_neg(B,B);
2837         mpz_neg(N,N);
2838       }
2839       // check for gcd(N,B)==1
2840       mpz_gcd(tmp,N,B);
2841       if (mpz_cmp_ui(tmp,1)==0)
2842       {
2843         // return N/B
2844         z=ALLOC_RNUMBER();
2845         #ifdef LDEBUG
2846         z->debug=123456;
2847         #endif
2848         mpz_init_set(z->z,N);
2849         mpz_init_set(z->n,B);
2850         z->s = 0;
2851         nlNormalize(z,r);
2852       }
2853       else
2854       {
2855         // return nN (the input) instead of "fail"
2856         z=nlCopy(nN,r);
2857       }
2858       break;
2859    }
2860    //mpz_mod(D,E,N);
2861    //mpz_div(tmp,E,N);
2862    mpz_divmod(tmp,D,E,N);
2863    mpz_mul(tmp,tmp,B);
2864    mpz_sub(C,A,tmp);
2865    mpz_set(E,N);
2866    mpz_set(N,D);
2867    mpz_set(A,B);
2868    mpz_set(B,C);
2869  }
2870  mpz_clear(tmp);
2871  mpz_clear(A);
2872  mpz_clear(B);
2873  mpz_clear(C);
2874  mpz_clear(D);
2875  mpz_clear(E);
2876  mpz_clear(N);
2877  mpz_clear(P);
2878  return z;
2879}
2880
2881number nlExtGcd(number a, number b, number *s, number *t, const coeffs)
2882{
2883  mpz_ptr aa,bb;
2884  *s=ALLOC_RNUMBER();
2885  mpz_init((*s)->z); (*s)->s=3;
2886  (*t)=ALLOC_RNUMBER();
2887  mpz_init((*t)->z); (*t)->s=3;
2888  number g=ALLOC_RNUMBER();
2889  mpz_init(g->z); g->s=3;
2890  #ifdef LDEBUG
2891  g->debug=123456;
2892  (*s)->debug=123456;
2893  (*t)->debug=123456;
2894  #endif
2895  if (SR_HDL(a) & SR_INT)
2896  {
2897    aa=(mpz_ptr)omAlloc(sizeof(mpz_t));
2898    mpz_init_set_si(aa,SR_TO_INT(a));
2899  }
2900  else
2901  {
2902    aa=a->z;
2903  }
2904  if (SR_HDL(b) & SR_INT)
2905  {
2906    bb=(mpz_ptr)omAlloc(sizeof(mpz_t));
2907    mpz_init_set_si(bb,SR_TO_INT(b));
2908  }
2909  else
2910  {
2911    bb=b->z;
2912  }
2913  mpz_gcdext(g->z,(*s)->z,(*t)->z,aa,bb);
2914  g=nlShort3(g);
2915  (*s)=nlShort3((*s));
2916  (*t)=nlShort3((*t));
2917  if (SR_HDL(a) & SR_INT)
2918  {
2919    mpz_clear(aa);
2920    omFreeSize(aa, sizeof(mpz_t));
2921  }
2922  if (SR_HDL(b) & SR_INT)
2923  {
2924    mpz_clear(bb);
2925    omFreeSize(bb, sizeof(mpz_t));
2926  }
2927  return g;
2928}
2929
2930void    nlCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
2931{
2932  if (r->is_field)
2933  PrintS("QQ");
2934  else
2935  PrintS("ZZ");
2936}
2937
2938int n_SwitchChinRem=0;
2939number   nlChineseRemainderSym(number *x, number *q,int rl, BOOLEAN sym, CFArray &inv_cache,const coeffs CF)
2940// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
2941{
2942  setCharacteristic( 0 ); // only in char 0
2943  Off(SW_RATIONAL);
2944  CFArray X(rl), Q(rl);
2945  int i;
2946  for(i=rl-1;i>=0;i--)
2947  {
2948    X[i]=CF->convSingNFactoryN(x[i],FALSE,CF); // may be larger MAX_INT
2949    Q[i]=CF->convSingNFactoryN(q[i],FALSE,CF); // may be larger MAX_INT
2950  }
2951  CanonicalForm xnew,qnew;
2952  if (n_SwitchChinRem)
2953    chineseRemainder(X,Q,xnew,qnew);
2954  else
2955    chineseRemainderCached(X,Q,xnew,qnew,inv_cache);
2956  number n=CF->convFactoryNSingN(xnew,CF);
2957  if (sym)
2958  {
2959    number p=CF->convFactoryNSingN(qnew,CF);
2960    number p2;
2961    if (getCoeffType(CF) == n_Q) p2=nlIntDiv(p,nlInit(2, CF),CF);
2962    else                         p2=CF->cfDiv(p,CF->cfInit(2, CF),CF);
2963    if (CF->cfGreater(n,p2,CF))
2964    {
2965       number n2=CF->cfSub(n,p,CF);
2966       CF->cfDelete(&n,CF);
2967       n=n2;
2968    }
2969    CF->cfDelete(&p2,CF);
2970    CF->cfDelete(&p,CF);
2971  }
2972  CF->cfNormalize(n,CF);
2973  return n;
2974}
2975#if 0
2976number nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
2977{
2978  CFArray inv(rl);
2979  return nlChineseRemainderSym(x,q,rl,TRUE,inv,C);
2980}
2981#endif
2982
2983static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2984{
2985  assume(cf != NULL);
2986
2987  numberCollectionEnumerator.Reset();
2988
2989  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2990  {
2991    c = nlInit(1, cf);
2992    return;
2993  }
2994
2995  // all coeffs are given by integers!!!
2996
2997  // part 1, find a small candidate for gcd
2998  number cand1,cand;
2999  int s1,s;
3000  s=2147483647; // max. int
3001
3002  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
3003
3004  int normalcount = 0;
3005  do
3006  {
3007    number& n = numberCollectionEnumerator.Current();
3008    nlNormalize(n, cf); ++normalcount;
3009    cand1 = n;
3010
3011    if (SR_HDL(cand1)&SR_INT) { cand=cand1; break; }
3012    assume(cand1->s==3); // all coeffs should be integers // ==0?!! after printing
3013    s1=mpz_size1(cand1->z);
3014    if (s>s1)
3015    {
3016      cand=cand1;
3017      s=s1;
3018    }
3019  } while (numberCollectionEnumerator.MoveNext() );
3020
3021//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
3022
3023  cand=nlCopy(cand,cf);
3024  // part 2: compute gcd(cand,all coeffs)
3025
3026  numberCollectionEnumerator.Reset();
3027
3028  while (numberCollectionEnumerator.MoveNext() )
3029  {
3030    number& n = numberCollectionEnumerator.Current();
3031
3032    if( (--normalcount) <= 0)
3033      nlNormalize(n, cf);
3034
3035    nlInpGcd(cand, n, cf);
3036    assume( nlGreaterZero(cand,cf) );
3037
3038    if(nlIsOne(cand,cf))
3039    {
3040      c = cand;
3041
3042      if(!lc_is_pos)
3043      {
3044        // make the leading coeff positive
3045        c = nlNeg(c, cf);
3046        numberCollectionEnumerator.Reset();
3047
3048        while (numberCollectionEnumerator.MoveNext() )
3049        {
3050          number& nn = numberCollectionEnumerator.Current();
3051          nn = nlNeg(nn, cf);
3052        }
3053      }
3054      return;
3055    }
3056  }
3057
3058  // part3: all coeffs = all coeffs / cand
3059  if (!lc_is_pos)
3060    cand = nlNeg(cand,cf);
3061
3062  c = cand;
3063  numberCollectionEnumerator.Reset();
3064
3065  while (numberCollectionEnumerator.MoveNext() )
3066  {
3067    number& n = numberCollectionEnumerator.Current();
3068    number t=nlIntDiv(n, cand, cf); // simple integer exact division, no ratios to remain
3069    nlDelete(&n, cf);
3070    n = t;
3071  }
3072}
3073
3074static void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
3075{
3076  assume(cf != NULL);
3077
3078  numberCollectionEnumerator.Reset();
3079
3080  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
3081  {
3082    c = nlInit(1, cf);
3083//    assume( n_GreaterZero(c, cf) );
3084    return;
3085  }
3086
3087  // all coeffs are given by integers after returning from this routine
3088
3089  // part 1, collect product of all denominators /gcds
3090  number cand;
3091  cand=ALLOC_RNUMBER();
3092#if defined(LDEBUG)
3093  cand->debug=123456;
3094#endif
3095  cand->s=3;
3096
3097  int s=0;
3098
3099  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
3100
3101  do
3102  {
3103    number& cand1 = numberCollectionEnumerator.Current();
3104
3105    if (!(SR_HDL(cand1)&SR_INT))
3106    {
3107      nlNormalize(cand1, cf);
3108      if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
3109      && (cand1->s==1))             // and is a normalised rational
3110      {
3111        if (s==0) // first denom, we meet
3112        {
3113          mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
3114          s=1;
3115        }
3116        else // we have already something
3117        {
3118          mpz_lcm(cand->z, cand->z, cand1->n);
3119        }
3120      }
3121    }
3122  }
3123  while (numberCollectionEnumerator.MoveNext() );
3124
3125
3126  if (s==0) // nothing to do, all coeffs are already integers
3127  {
3128//    mpz_clear(tmp);
3129    FREE_RNUMBER(cand);
3130    if (lc_is_pos)
3131      c=nlInit(1,cf);
3132    else
3133    {
3134      // make the leading coeff positive
3135      c=nlInit(-1,cf);
3136
3137      // TODO: incorporate the following into the loop below?
3138      numberCollectionEnumerator.Reset();
3139      while (numberCollectionEnumerator.MoveNext() )
3140      {
3141        number& n = numberCollectionEnumerator.Current();
3142        n = nlNeg(n, cf);
3143      }
3144    }
3145//    assume( n_GreaterZero(c, cf) );
3146    return;
3147  }
3148
3149  cand = nlShort3(cand);
3150
3151  // part2: all coeffs = all coeffs * cand
3152  // make the lead coeff positive
3153  numberCollectionEnumerator.Reset();
3154
3155  if (!lc_is_pos)
3156    cand = nlNeg(cand, cf);
3157
3158  c = cand;
3159
3160  while (numberCollectionEnumerator.MoveNext() )
3161  {
3162    number &n = numberCollectionEnumerator.Current();
3163    nlInpMult(n, cand, cf);
3164  }
3165
3166}
3167
3168char * nlCoeffName(const coeffs r)
3169{
3170  if (r->cfDiv==nlDiv) return (char*)"QQ";
3171  else                 return (char*)"ZZ";
3172}
3173
3174static char* nlCoeffString(const coeffs r)
3175{
3176  //return omStrDup(nlCoeffName(r));
3177  if (r->cfDiv==nlDiv) return omStrDup("QQ");
3178  else                 return omStrDup("ZZ");
3179}
3180
3181void nlWriteFd(number n, const ssiInfo* d, const coeffs)
3182{
3183  if(SR_HDL(n) & SR_INT)
3184  {
3185    #if SIZEOF_LONG == 4
3186    fprintf(d->f_write,"4 %ld ",SR_TO_INT(n));
3187    #else
3188    long nn=SR_TO_INT(n);
3189    if ((nn<POW_2_28_32)&&(nn>= -POW_2_28_32))
3190    {
3191      int nnn=(int)nn;
3192      fprintf(d->f_write,"4 %d ",nnn);
3193    }
3194    else
3195    {
3196      mpz_t tmp;
3197      mpz_init_set_si(tmp,nn);
3198      fputs("8 ",d->f_write);
3199      mpz_out_str (d->f_write,SSI_BASE, tmp);
3200      fputc(' ',d->f_write);
3201      mpz_clear(tmp);
3202    }
3203    #endif
3204  }
3205  else if (n->s<2)
3206  {
3207    //gmp_fprintf(f,"%d %Zd %Zd ",n->s,n->z,n->n);
3208    fprintf(d->f_write,"%d ",n->s+5);
3209    mpz_out_str (d->f_write,SSI_BASE, n->z);
3210    fputc(' ',d->f_write);
3211    mpz_out_str (d->f_write,SSI_BASE, n->n);
3212    fputc(' ',d->f_write);
3213
3214    //if (d->f_debug!=NULL) gmp_fprintf(d->f_debug,"number: s=%d gmp/gmp \"%Zd %Zd\" ",n->s,n->z,n->n);
3215  }
3216  else /*n->s==3*/
3217  {
3218    //gmp_fprintf(d->f_write,"3 %Zd ",n->z);
3219    fputs("8 ",d->f_write);
3220    mpz_out_str (d->f_write,SSI_BASE, n->z);
3221    fputc(' ',d->f_write);
3222
3223    //if (d->f_debug!=NULL) gmp_fprintf(d->f_debug,"number: gmp \"%Zd\" ",n->z);
3224  }
3225}
3226
3227number nlReadFd(const ssiInfo *d, const coeffs)
3228{
3229  int sub_type=-1;
3230  sub_type=s_readint(d->f_read);
3231  switch(sub_type)
3232  {
3233     case 0:
3234     case 1:
3235       {// read mpz_t, mpz_t
3236         number n=nlRInit(0);
3237         mpz_init(n->n);
3238         s_readmpz(d->f_read,n->z);
3239         s_readmpz(d->f_read,n->n);
3240         n->s=sub_type;
3241         return n;
3242       }
3243
3244     case 3:
3245       {// read mpz_t
3246         number n=nlRInit(0);
3247         s_readmpz(d->f_read,n->z);
3248         n->s=3; /*sub_type*/
3249         #if SIZEOF_LONG == 8
3250         n=nlShort3(n);
3251         #endif
3252         return n;
3253       }
3254     case 4:
3255       {
3256         LONG dd=s_readlong(d->f_read);
3257         //#if SIZEOF_LONG == 8
3258         return INT_TO_SR(dd);
3259         //#else
3260         //return nlInit(dd,NULL);
3261         //#endif
3262       }
3263     case 5:
3264     case 6:
3265       {// read raw mpz_t, mpz_t
3266         number n=nlRInit(0);
3267         mpz_init(n->n);
3268         s_readmpz_base (d->f_read,n->z, SSI_BASE);
3269         s_readmpz_base (d->f_read,n->n, SSI_BASE);
3270         n->s=sub_type-5;
3271         return n;
3272       }
3273     case 8:
3274       {// read raw mpz_t
3275         number n=nlRInit(0);
3276         s_readmpz_base (d->f_read,n->z, SSI_BASE);
3277         n->s=sub_type=3; /*subtype-5*/
3278         #if SIZEOF_LONG == 8
3279         n=nlShort3(n);
3280         #endif
3281         return n;
3282       }
3283
3284     default: Werror("error in reading number: invalid subtype %d",sub_type);
3285              return NULL;
3286  }
3287  return NULL;
3288}
3289
3290BOOLEAN nlCoeffIsEqual(const coeffs r, n_coeffType n, void *p)
3291{
3292  /* test, if r is an instance of nInitCoeffs(n,parameter) */
3293  /* if parameter is not needed */
3294  if (n==r->type)
3295  {
3296    if ((p==NULL)&&(r->cfDiv==nlDiv)) return TRUE;
3297    if ((p!=NULL)&&(r->cfDiv!=nlDiv)) return TRUE;
3298  }
3299  return FALSE;
3300}
3301
3302static number nlLcm(number a,number b,const coeffs r)
3303{
3304  number g=nlGcd(a,b,r);
3305  number n1=nlMult(a,b,r);
3306  number n2=nlIntDiv(n1,g,r);
3307  nlDelete(&g,r);
3308  nlDelete(&n1,r);
3309  return n2;
3310}
3311
3312static number nlRandom(siRandProc p, number v2, number, const coeffs cf)
3313{
3314  number a=nlInit(p(),cf);
3315  if (v2!=NULL)
3316  {
3317    number b=nlInit(p(),cf);
3318    number c=nlDiv(a,b,cf);
3319    nlDelete(&b,cf);
3320    nlDelete(&a,cf);
3321    a=c;
3322  }
3323  return a;
3324}
3325
3326BOOLEAN nlInitChar(coeffs r, void*p)
3327{
3328  r->is_domain=TRUE;
3329  r->rep=n_rep_gap_rat;
3330
3331  r->nCoeffIsEqual=nlCoeffIsEqual;
3332  //r->cfKillChar = ndKillChar; /* dummy */
3333  r->cfCoeffString=nlCoeffString;
3334  r->cfCoeffName=nlCoeffName;
3335
3336  r->cfInitMPZ = nlInitMPZ;
3337  r->cfMPZ  = nlMPZ;
3338
3339  r->cfMult  = nlMult;
3340  r->cfSub   = nlSub;
3341  r->cfAdd   = nlAdd;
3342  r->cfExactDiv= nlExactDiv;
3343  if (p==NULL) /* Q */
3344  {
3345    r->is_field=TRUE;
3346    r->cfDiv   = nlDiv;
3347    //r->cfGcd  = ndGcd_dummy;
3348    r->cfSubringGcd  = nlGcd;
3349  }
3350  else /* Z: coeffs_BIGINT */
3351  {
3352    r->is_field=FALSE;
3353    r->cfDiv   = nlIntDiv;
3354    r->cfIntMod= nlIntMod;
3355    r->cfGcd  = nlGcd;
3356    r->cfDivBy=nlDivBy;
3357    r->cfDivComp = nlDivComp;
3358    r->cfIsUnit = nlIsUnit;
3359    r->cfGetUnit = nlGetUnit;
3360    r->cfQuot1 = nlQuot1;
3361    r->cfLcm = nlLcm;
3362    r->cfXExtGcd=nlXExtGcd;
3363    r->cfQuotRem=nlQuotRem;
3364  }
3365  r->cfInit = nlInit;
3366  r->cfSize  = nlSize;
3367  r->cfInt  = nlInt;
3368
3369  r->cfChineseRemainder=nlChineseRemainderSym;
3370  r->cfFarey=nlFarey;
3371  r->cfInpNeg   = nlNeg;
3372  r->cfInvers= nlInvers;
3373  r->cfCopy  = nlCopy;
3374  r->cfRePart = nlCopy;
3375  //r->cfImPart = ndReturn0;
3376  r->cfWriteLong = nlWrite;
3377  r->cfRead = nlRead;
3378  r->cfNormalize=nlNormalize;
3379  r->cfGreater = nlGreater;
3380  r->cfEqual = nlEqual;
3381  r->cfIsZero = nlIsZero;
3382  r->cfIsOne = nlIsOne;
3383  r->cfIsMOne = nlIsMOne;
3384  r->cfGreaterZero = nlGreaterZero;
3385  r->cfPower = nlPower;
3386  r->cfGetDenom = nlGetDenom;
3387  r->cfGetNumerator = nlGetNumerator;
3388  r->cfExtGcd = nlExtGcd; // only for ring stuff and Z
3389  r->cfNormalizeHelper  = nlNormalizeHelper;
3390  r->cfDelete= nlDelete;
3391  r->cfSetMap = nlSetMap;
3392  //r->cfName = ndName;
3393  r->cfInpMult=nlInpMult;
3394  r->cfInpAdd=nlInpAdd;
3395  r->cfCoeffWrite=nlCoeffWrite;
3396
3397  r->cfClearContent = nlClearContent;
3398  r->cfClearDenominators = nlClearDenominators;
3399
3400#ifdef LDEBUG
3401  // debug stuff
3402  r->cfDBTest=nlDBTest;
3403#endif
3404  r->convSingNFactoryN=nlConvSingNFactoryN;
3405  r->convFactoryNSingN=nlConvFactoryNSingN;
3406
3407  r->cfRandom=nlRandom;
3408
3409  // io via ssi
3410  r->cfWriteFd=nlWriteFd;
3411  r->cfReadFd=nlReadFd;
3412
3413  //r->type = n_Q;
3414  r->ch = 0;
3415  r->has_simple_Alloc=FALSE;
3416  r->has_simple_Inverse=FALSE;
3417
3418  // variables for this type of coeffs:
3419  // (none)
3420  return FALSE;
3421}
3422#if 0
3423number nlMod(number a, number b)
3424{
3425  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
3426  {
3427    int bi=SR_TO_INT(b);
3428    int ai=SR_TO_INT(a);
3429    int bb=ABS(bi);
3430    int c=ai%bb;
3431    if (c<0)  c+=bb;
3432    return (INT_TO_SR(c));
3433  }
3434  number al;
3435  number bl;
3436  if (SR_HDL(a))&SR_INT)
3437    al=nlRInit(SR_TO_INT(a));
3438  else
3439    al=nlCopy(a);
3440  if (SR_HDL(b))&SR_INT)
3441    bl=nlRInit(SR_TO_INT(b));
3442  else
3443    bl=nlCopy(b);
3444  number r=nlRInit(0);
3445  mpz_mod(r->z,al->z,bl->z);
3446  nlDelete(&al);
3447  nlDelete(&bl);
3448  if (mpz_size1(&r->z)<=MP_SMALL)
3449  {
3450    LONG ui=(int)mpz_get_si(&r->z);
3451    if ((((ui<<3)>>3)==ui)
3452    && (mpz_cmp_si(x->z,(long)ui)==0))
3453    {
3454      mpz_clear(&r->z);
3455      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
3456      r=INT_TO_SR(ui);
3457    }
3458  }
3459  return r;
3460}
3461#endif
3462#endif // not P_NUMBERS_H
3463#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.