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

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