source: git/libpolys/coeffs/longrat.cc @ 863333

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