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

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