source: git/libpolys/coeffs/longrat.cc @ 4f9652

spielwiese
Last change on this file since 4f9652 was 4f9652, checked in by Hans Schoenemann <hannes@…>, 12 years ago
add: nlClearContent
  • Property mode set to 100644
File size: 57.1 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/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    if ((long)a>0L)
806    {
807      if ((long)b>0L)
808        return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
809      else
810        return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
811    }
812    else
813    {
814      if ((long)b>0L)
815      {
816        long i=(-SR_TO_INT(a))%SR_TO_INT(b);
817        if ( i != 0L ) i = (SR_TO_INT(b))-i;
818        return INT_TO_SR(i);
819      }
820      else
821      {
822        long i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
823        if ( i != 0L ) i = (-SR_TO_INT(b))-i;
824        return INT_TO_SR(i);
825      }
826    }
827  }
828  if (SR_HDL(a) & SR_INT)
829  {
830    /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
831    if ((long)a<0L)
832    {
833      if (mpz_isNeg(b->z))
834        return nlSub(a,b,r);
835      /*else*/
836        return nlAdd(a,b,r);
837    }
838    /*else*/
839      return a;
840  }
841  number bb=NULL;
842  if (SR_HDL(b) & SR_INT)
843  {
844    bb=nlRInit(SR_TO_INT(b));
845    b=bb;
846  }
847  u=ALLOC_RNUMBER();
848#if defined(LDEBUG)
849  u->debug=123456;
850#endif
851  mpz_init(u->z);
852  u->s = 3;
853  mpz_mod(u->z,a->z,b->z);
854  if (bb!=NULL)
855  {
856    mpz_clear(bb->z);
857#if defined(LDEBUG)
858    bb->debug=654324;
859#endif
860    FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
861  }
862  if (mpz_isNeg(u->z))
863  {
864    if (mpz_isNeg(b->z))
865      mpz_sub(u->z,u->z,b->z);
866    else
867      mpz_add(u->z,u->z,b->z);
868  }
869  u=nlShort3(u);
870  nlTest(u, r);
871  return u;
872}
873
874/*2
875* u := a / b
876*/
877number nlDiv (number a, number b, const coeffs r)
878{
879  number u;
880  if (nlIsZero(b,r))
881  {
882    WerrorS(nDivBy0);
883    return INT_TO_SR(0);
884  }
885  u=ALLOC_RNUMBER();
886  u->s=0;
887#if defined(LDEBUG)
888  u->debug=123456;
889#endif
890// ---------- short / short ------------------------------------
891  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
892  {
893    long i=SR_TO_INT(a);
894    long j=SR_TO_INT(b);
895    if ((i==-POW_2_28) && (j== -1L))
896    { 
897      FREE_RNUMBER(u);
898      return nlRInit(POW_2_28);
899    }
900    long r=i%j;
901    if (r==0)
902    {
903      FREE_RNUMBER(u); // omFreeBin((void *)u, rnumber_bin);
904      return INT_TO_SR(i/j);
905    }
906    mpz_init_set_si(u->z,(long)i);
907    mpz_init_set_si(u->n,(long)j);
908  }
909  else
910  {
911    mpz_init(u->z);
912// ---------- short / long ------------------------------------
913    if (SR_HDL(a) & SR_INT)
914    {
915      // short a / (z/n) -> (a*n)/z
916      if (b->s<2)
917      {
918        mpz_mul_si(u->z,b->n,SR_TO_INT(a));
919      }
920      else
921      // short a / long z -> a/z
922      {
923        mpz_set_si(u->z,SR_TO_INT(a));
924      }
925      if (mpz_cmp(u->z,b->z)==0)
926      {
927        mpz_clear(u->z);
928        FREE_RNUMBER(u);
929        return INT_TO_SR(1);
930      }
931      mpz_init_set(u->n,b->z);
932    }
933// ---------- long / short ------------------------------------
934    else if (SR_HDL(b) & SR_INT)
935    {
936      mpz_set(u->z,a->z);
937      // (z/n) / b -> z/(n*b)
938      if (a->s<2)
939      {
940        mpz_init_set(u->n,a->n);
941        if ((long)b>0L)
942          mpz_mul_ui(u->n,u->n,SR_TO_INT(b));
943        else
944        {
945          mpz_mul_ui(u->n,u->n,-SR_TO_INT(b));
946          mpz_neg(u->z,u->z);
947        }
948      }
949      else
950      // long z / short b -> z/b
951      {
952        //mpz_set(u->z,a->z);
953        mpz_init_set_si(u->n,SR_TO_INT(b));
954      }
955    }
956// ---------- long / long ------------------------------------
957    else
958    {
959      mpz_set(u->z,a->z);
960      mpz_init_set(u->n,b->z);
961      if (a->s<2) mpz_mul(u->n,u->n,a->n);
962      if (b->s<2) mpz_mul(u->z,u->z,b->n);
963    }
964  }
965  if (mpz_isNeg(u->n))
966  {
967    mpz_neg(u->z,u->z);
968    mpz_neg(u->n,u->n);
969  }
970  if (mpz_cmp_si(u->n,(long)1)==0)
971  {
972    mpz_clear(u->n);
973    u->s=3;
974    u=nlShort3(u);
975  }
976  nlTest(u, r);
977  return u;
978}
979
980/*2
981* u:= x ^ exp
982*/
983void nlPower (number x,int exp,number * u, const coeffs r)
984{
985  *u = INT_TO_SR(0); // 0^e, e!=0
986  if (!nlIsZero(x,r))
987  {
988    nlTest(x, r);
989    number aa=NULL;
990    if (SR_HDL(x) & SR_INT)
991    {
992      aa=nlRInit(SR_TO_INT(x));
993      x=aa;
994    }
995    else if (x->s==0)
996      nlNormalize(x,r);
997    *u=ALLOC_RNUMBER();
998#if defined(LDEBUG)
999    (*u)->debug=123456;
1000#endif
1001    mpz_init((*u)->z);
1002    mpz_pow_ui((*u)->z,x->z,(unsigned long)exp);
1003    if (x->s<2)
1004    {
1005      if (mpz_cmp_si(x->n,(long)1)==0)
1006      {
1007        x->s=3;
1008        mpz_clear(x->n);
1009      }
1010      else
1011      {
1012        mpz_init((*u)->n);
1013        mpz_pow_ui((*u)->n,x->n,(unsigned long)exp);
1014      }
1015    }
1016    (*u)->s = x->s;
1017    if ((*u)->s==3) *u=nlShort3(*u);
1018    if (aa!=NULL)
1019    {
1020      mpz_clear(aa->z);
1021      FREE_RNUMBER(aa);
1022    }
1023  }
1024  else if (exp==0)
1025    *u = INT_TO_SR(1); // 0^0
1026#ifdef LDEBUG
1027  if (exp<0) Print("nlPower: neg. exp. %d\n",exp);
1028  nlTest(*u, r);
1029#endif
1030}
1031
1032
1033/*2
1034* za >= 0 ?
1035*/
1036BOOLEAN nlGreaterZero (number a, const coeffs r)
1037{
1038  nlTest(a, r);
1039  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>1L /* represents number(0) */;
1040  return (!mpz_isNeg(a->z));
1041}
1042
1043/*2
1044* a > b ?
1045*/
1046BOOLEAN nlGreater (number a, number b, const coeffs r)
1047{
1048  nlTest(a, r);
1049  nlTest(b, r);
1050  number re;
1051  BOOLEAN rr;
1052  re=nlSub(a,b,r);
1053  rr=(!nlIsZero(re,r)) && (nlGreaterZero(re,r));
1054  nlDelete(&re,r);
1055  return rr;
1056}
1057
1058/*2
1059* a == -1 ?
1060*/
1061BOOLEAN nlIsMOne (number a, const coeffs r)
1062{
1063#ifdef LDEBUG
1064  if (a==NULL) return FALSE;
1065  nlTest(a, r);
1066#endif
1067  //if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1L));
1068  //return FALSE;
1069  return  (a==INT_TO_SR(-1L));
1070}
1071
1072/*2
1073* result =gcd(a,b)
1074*/
1075number nlGcd(number a, number b, const coeffs r)
1076{
1077  number result;
1078  nlTest(a, r);
1079  nlTest(b, r);
1080  //nlNormalize(a);
1081  //nlNormalize(b);
1082  if ((a==INT_TO_SR(1L))||(a==INT_TO_SR(-1L))
1083  ||  (b==INT_TO_SR(1L))||(b==INT_TO_SR(-1L)))
1084    return INT_TO_SR(1L);
1085  if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
1086    return nlCopy(b,r);
1087  if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
1088    return nlCopy(a,r);
1089  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1090  {
1091    long i=SR_TO_INT(a);
1092    long j=SR_TO_INT(b);
1093    if((i==0L)||(j==0L))
1094      return INT_TO_SR(1);
1095    long l;
1096    i=ABS(i);
1097    j=ABS(j);
1098    do
1099    {
1100      l=i%j;
1101      i=j;
1102      j=l;
1103    } while (l!=0L);
1104    if (i==POW_2_28)
1105      result=nlRInit(POW_2_28);
1106    else
1107     result=INT_TO_SR(i);
1108    nlTest(result,r);
1109    return result;
1110  }
1111  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1112  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1113  if (SR_HDL(a) & SR_INT)
1114  {
1115    LONG aa=ABS(SR_TO_INT(a));
1116    unsigned long t=mpz_gcd_ui(NULL,b->z,(long)aa);
1117    if (t==POW_2_28)
1118      result=nlRInit(POW_2_28);
1119    else
1120     result=INT_TO_SR(t);
1121  }
1122  else
1123  if (SR_HDL(b) & SR_INT)
1124  {
1125    LONG bb=ABS(SR_TO_INT(b));
1126    unsigned long t=mpz_gcd_ui(NULL,a->z,(long)bb);
1127    if (t==POW_2_28)
1128      result=nlRInit(POW_2_28);
1129    else
1130     result=INT_TO_SR(t);
1131  }
1132  else
1133  {
1134    result=ALLOC_RNUMBER();
1135    mpz_init(result->z);
1136    mpz_gcd(result->z,a->z,b->z);
1137    result->s = 3;
1138  #ifdef LDEBUG
1139    result->debug=123456;
1140  #endif
1141    result=nlShort3(result);
1142  }
1143  nlTest(result, r);
1144  return result;
1145}
1146
1147number nlShort1(number x) // assume x->s==0/1
1148{
1149  assume(x->s<2);
1150  if (mpz_cmp_ui(x->z,(long)0)==0)
1151  {
1152    _nlDelete_NoImm(&x);
1153    return INT_TO_SR(0);
1154  }
1155  if (x->s<2)
1156  {
1157    if (mpz_cmp(x->z,x->n)==0)
1158    {
1159      _nlDelete_NoImm(&x);
1160      return INT_TO_SR(1);
1161    }
1162  }
1163  return x;
1164}
1165/*2
1166* simplify x
1167*/
1168void nlNormalize (number &x, const coeffs r)
1169{
1170  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1171    return;
1172  if (x->s==3)
1173  {
1174    x=nlShort3_noinline(x);
1175    nlTest(x,r);
1176    return;
1177  }
1178  else if (x->s==0)
1179  {
1180    if (mpz_cmp_si(x->n,(long)1)==0)
1181    {
1182      mpz_clear(x->n);
1183      x->s=3;
1184      x=nlShort3(x);
1185    }
1186    else
1187    {
1188      mpz_t gcd;
1189      mpz_init(gcd);
1190      mpz_gcd(gcd,x->z,x->n);
1191      x->s=1;
1192      if (mpz_cmp_si(gcd,(long)1)!=0)
1193      {
1194        mpz_t r;
1195        mpz_init(r);
1196        MPZ_EXACTDIV(r,x->z,gcd);
1197        mpz_set(x->z,r);
1198        MPZ_EXACTDIV(r,x->n,gcd);
1199        mpz_set(x->n,r);
1200        mpz_clear(r);
1201        if (mpz_cmp_si(x->n,(long)1)==0)
1202        {
1203          mpz_clear(x->n);
1204          x->s=3;
1205          x=nlShort3_noinline(x);
1206        }
1207      }
1208      mpz_clear(gcd);
1209    }
1210  }
1211  nlTest(x, r);
1212}
1213
1214/*2
1215* returns in result->z the lcm(a->z,b->n)
1216*/
1217number nlLcm(number a, number b, const coeffs r)
1218{
1219  number result;
1220  nlTest(a, r);
1221  nlTest(b, r);
1222  if ((SR_HDL(b) & SR_INT)
1223  || (b->s==3))
1224  {
1225    // b is 1/(b->n) => b->n is 1 => result is a
1226    return nlCopy(a,r);
1227  }
1228  result=ALLOC_RNUMBER();
1229#if defined(LDEBUG)
1230  result->debug=123456;
1231#endif
1232  result->s=3;
1233  mpz_t gcd;
1234  mpz_init(gcd);
1235  mpz_init(result->z);
1236  if (SR_HDL(a) & SR_INT)
1237    mpz_gcd_ui(gcd,b->n,ABS(SR_TO_INT(a)));
1238  else
1239    mpz_gcd(gcd,a->z,b->n);
1240  if (mpz_cmp_si(gcd,(long)1)!=0)
1241  {
1242    mpz_t bt;
1243    mpz_init_set(bt,b->n);
1244    MPZ_EXACTDIV(bt,bt,gcd);
1245    if (SR_HDL(a) & SR_INT)
1246      mpz_mul_si(result->z,bt,SR_TO_INT(a));
1247    else
1248      mpz_mul(result->z,bt,a->z);
1249    mpz_clear(bt);
1250  }
1251  else
1252    if (SR_HDL(a) & SR_INT)
1253      mpz_mul_si(result->z,b->n,SR_TO_INT(a));
1254    else
1255      mpz_mul(result->z,b->n,a->z);
1256  mpz_clear(gcd);
1257  result=nlShort3(result);
1258  nlTest(result, r);
1259  return result;
1260}
1261
1262// Map q \in QQ \to Zp
1263// src = Q, dst = Zp (or an extension of Zp?)
1264number nlModP(number q, const coeffs Q, const coeffs Zp)
1265{
1266  assume( getCoeffType(Q) == ID );
1267
1268  const int p = n_GetChar(Zp);
1269  assume( p > 0 );
1270
1271  const long P = p;
1272  assume( P > 0 );
1273
1274  // embedded long within q => only long numerator has to be converted
1275  // to int (modulo char.)
1276  if (SR_HDL(q) & SR_INT)
1277  {
1278    long i = SR_TO_INT(q);
1279    if (i<0L)
1280      return n_Init( static_cast<int>( P - ((-i)%P) ), Zp);
1281    else
1282      return n_Init( static_cast<int>( i % P ), Zp );
1283  }
1284
1285  const unsigned long PP = p;
1286
1287  // numerator modulo char. should fit into int
1288  number z = n_Init( static_cast<int>(mpz_fdiv_ui(q->z, PP)), Zp ); 
1289
1290  // denominator != 1?
1291  if (q->s!=3)
1292  {
1293    // denominator modulo char. should fit into int
1294    number n = n_Init( static_cast<int>(mpz_fdiv_ui(q->n, PP)), Zp );
1295
1296    number res = n_Div( z, n, Zp );
1297
1298    n_Delete(&z, Zp);
1299    n_Delete(&n, Zp);
1300
1301    return res;
1302  }
1303
1304  return z;
1305}
1306
1307#ifdef HAVE_RINGS
1308/*2
1309* convert number i (from Q) to GMP and warn if denom != 1
1310*/
1311void nlGMP(number &i, number n, const coeffs r)
1312{
1313  // Hier brauche ich einfach die GMP Zahl
1314  nlTest(i, r);
1315  nlNormalize(i, r);
1316  if (SR_HDL(i) & SR_INT)
1317  {
1318    mpz_set_si((mpz_ptr) n, SR_TO_INT(i));
1319    return;
1320  }
1321  if (i->s!=3)
1322  {
1323     WarnS("Omitted denominator during coefficient mapping !");
1324  }
1325  mpz_set((mpz_ptr) n, i->z);
1326}
1327#endif
1328
1329/*2
1330* acces to denominator, other 1 for integers
1331*/
1332number   nlGetDenom(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      if (n->s!=3)
1343      {
1344        number u=ALLOC_RNUMBER();
1345        u->s=3;
1346#if defined(LDEBUG)
1347        u->debug=123456;
1348#endif
1349        mpz_init_set(u->z,n->n);
1350        u=nlShort3_noinline(u);
1351        return u;
1352      }
1353    }
1354  }
1355  return INT_TO_SR(1);
1356}
1357
1358/*2
1359* acces to Nominator, nlCopy(n) for integers
1360*/
1361number   nlGetNumerator(number &n, const coeffs r)
1362{
1363  if (!(SR_HDL(n) & SR_INT))
1364  {
1365    if (n->s==0)
1366    {
1367      nlNormalize(n,r);
1368    }
1369    if (!(SR_HDL(n) & SR_INT))
1370    {
1371      number u=ALLOC_RNUMBER();
1372#if defined(LDEBUG)
1373      u->debug=123456;
1374#endif
1375      u->s=3;
1376      mpz_init_set(u->z,n->z);
1377      if (n->s!=3)
1378      {
1379        u=nlShort3_noinline(u);
1380      }
1381      return u;
1382    }
1383  }
1384  return n; // imm. int
1385}
1386
1387/***************************************************************
1388 *
1389 * routines which are needed by Inline(d) routines
1390 *
1391 *******************************************************************/
1392BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
1393{
1394  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1395//  long - short
1396  BOOLEAN bo;
1397  if (SR_HDL(b) & SR_INT)
1398  {
1399    if (a->s!=0) return FALSE;
1400    number n=b; b=a; a=n;
1401  }
1402//  short - long
1403  if (SR_HDL(a) & SR_INT)
1404  {
1405    if (b->s!=0)
1406      return FALSE;
1407    if (((long)a > 0L) && (mpz_isNeg(b->z)))
1408      return FALSE;
1409    if (((long)a < 0L) && (!mpz_isNeg(b->z)))
1410      return FALSE;
1411    mpz_t  bb;
1412    mpz_init_set(bb,b->n);
1413    mpz_mul_si(bb,bb,(long)SR_TO_INT(a));
1414    bo=(mpz_cmp(bb,b->z)==0);
1415    mpz_clear(bb);
1416    return bo;
1417  }
1418// long - long
1419  if (((a->s==1) && (b->s==3))
1420  ||  ((b->s==1) && (a->s==3)))
1421    return FALSE;
1422  if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
1423    return FALSE;
1424  if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
1425    return FALSE;
1426  mpz_t  aa;
1427  mpz_t  bb;
1428  mpz_init_set(aa,a->z);
1429  mpz_init_set(bb,b->z);
1430  if (a->s<2) mpz_mul(bb,bb,a->n);
1431  if (b->s<2) mpz_mul(aa,aa,b->n);
1432  bo=(mpz_cmp(aa,bb)==0);
1433  mpz_clear(aa);
1434  mpz_clear(bb);
1435  return bo;
1436}
1437
1438// copy not immediate number a
1439number _nlCopy_NoImm(number a)
1440{
1441  assume(!((SR_HDL(a) & SR_INT)||(a==NULL)));
1442  //nlTest(a, r);
1443  number b=ALLOC_RNUMBER();
1444#if defined(LDEBUG)
1445  b->debug=123456;
1446#endif
1447  switch (a->s)
1448  {
1449    case 0:
1450    case 1:
1451            mpz_init_set(b->n,a->n);
1452    case 3:
1453            mpz_init_set(b->z,a->z);
1454            break;
1455  }
1456  b->s = a->s;
1457  return b;
1458}
1459
1460void _nlDelete_NoImm(number *a)
1461{
1462  {
1463    switch ((*a)->s)
1464    {
1465      case 0:
1466      case 1:
1467        mpz_clear((*a)->n);
1468      case 3:
1469        mpz_clear((*a)->z);
1470#ifdef LDEBUG
1471        (*a)->s=2;
1472#endif
1473    }
1474    FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
1475  }
1476}
1477
1478number _nlNeg_NoImm(number a)
1479{
1480  {
1481    mpz_neg(a->z,a->z);
1482    if (a->s==3)
1483    {
1484      a=nlShort3(a);
1485    }
1486  }
1487  return a;
1488}
1489
1490number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1491{
1492  number u=ALLOC_RNUMBER();
1493#if defined(LDEBUG)
1494  u->debug=123456;
1495#endif
1496  mpz_init(u->z);
1497  if (SR_HDL(b) & SR_INT)
1498  {
1499    number x=a;
1500    a=b;
1501    b=x;
1502  }
1503  if (SR_HDL(a) & SR_INT)
1504  {
1505    switch (b->s)
1506    {
1507      case 0:
1508      case 1:/* a:short, b:1 */
1509      {
1510        mpz_t x;
1511        mpz_init(x);
1512        mpz_mul_si(x,b->n,SR_TO_INT(a));
1513        mpz_add(u->z,b->z,x);
1514        mpz_clear(x);
1515        if (mpz_cmp_ui(u->z,(long)0)==0)
1516        {
1517          mpz_clear(u->z);
1518          FREE_RNUMBER(u);
1519          return INT_TO_SR(0);
1520        }
1521        if (mpz_cmp(u->z,b->n)==0)
1522        {
1523          mpz_clear(u->z);
1524          FREE_RNUMBER(u);
1525          return INT_TO_SR(1);
1526        }
1527        mpz_init_set(u->n,b->n);
1528        u->s = 0;
1529        break;
1530      }
1531      case 3:
1532      {
1533        if ((long)a>0L)
1534          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1535        else
1536          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1537        if (mpz_cmp_ui(u->z,(long)0)==0)
1538        {
1539          mpz_clear(u->z);
1540          FREE_RNUMBER(u);
1541          return INT_TO_SR(0);
1542        }
1543        u->s = 3;
1544        u=nlShort3(u);
1545        break;
1546      }
1547    }
1548  }
1549  else
1550  {
1551    switch (a->s)
1552    {
1553      case 0:
1554      case 1:
1555      {
1556        switch(b->s)
1557        {
1558          case 0:
1559          case 1:
1560          {
1561            mpz_t x;
1562            mpz_init(x);
1563
1564            mpz_mul(x,b->z,a->n);
1565            mpz_mul(u->z,a->z,b->n);
1566            mpz_add(u->z,u->z,x);
1567            mpz_clear(x);
1568
1569            if (mpz_cmp_ui(u->z,(long)0)==0)
1570            {
1571              mpz_clear(u->z);
1572              FREE_RNUMBER(u);
1573              return INT_TO_SR(0);
1574            }
1575            mpz_init(u->n);
1576            mpz_mul(u->n,a->n,b->n);
1577            if (mpz_cmp(u->z,u->n)==0)
1578            {
1579               mpz_clear(u->z);
1580               mpz_clear(u->n);
1581               FREE_RNUMBER(u);
1582               return INT_TO_SR(1);
1583            }
1584            u->s = 0;
1585            break;
1586          }
1587          case 3: /* a:1 b:3 */
1588          {
1589            mpz_mul(u->z,b->z,a->n);
1590            mpz_add(u->z,u->z,a->z);
1591            if (mpz_cmp_ui(u->z,(long)0)==0)
1592            {
1593              mpz_clear(u->z);
1594              FREE_RNUMBER(u);
1595              return INT_TO_SR(0);
1596            }
1597            if (mpz_cmp(u->z,a->n)==0)
1598            {
1599              mpz_clear(u->z);
1600              FREE_RNUMBER(u);
1601              return INT_TO_SR(1);
1602            }
1603            mpz_init_set(u->n,a->n);
1604            u->s = 0;
1605            break;
1606          }
1607        } /*switch (b->s) */
1608        break;
1609      }
1610      case 3:
1611      {
1612        switch(b->s)
1613        {
1614          case 0:
1615          case 1:/* a:3, b:1 */
1616          {
1617            mpz_mul(u->z,a->z,b->n);
1618            mpz_add(u->z,u->z,b->z);
1619            if (mpz_cmp_ui(u->z,(long)0)==0)
1620            {
1621              mpz_clear(u->z);
1622              FREE_RNUMBER(u);
1623              return INT_TO_SR(0);
1624            }
1625            if (mpz_cmp(u->z,b->n)==0)
1626            {
1627              mpz_clear(u->z);
1628              FREE_RNUMBER(u);
1629              return INT_TO_SR(1);
1630            }
1631            mpz_init_set(u->n,b->n);
1632            u->s = 0;
1633            break;
1634          }
1635          case 3:
1636          {
1637            mpz_add(u->z,a->z,b->z);
1638            if (mpz_cmp_ui(u->z,(long)0)==0)
1639            {
1640              mpz_clear(u->z);
1641              FREE_RNUMBER(u);
1642              return INT_TO_SR(0);
1643            }
1644            u->s = 3;
1645            u=nlShort3(u);
1646            break;
1647          }
1648        }
1649        break;
1650      }
1651    }
1652  }
1653  return u;
1654}
1655
1656void _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b)
1657{
1658  if (SR_HDL(b) & SR_INT)
1659  {
1660    switch (a->s)
1661    {
1662      case 0:
1663      case 1:/* b:short, a:1 */
1664      {
1665        mpz_t x;
1666        mpz_init(x);
1667        mpz_mul_si(x,a->n,SR_TO_INT(b));
1668        mpz_add(a->z,a->z,x);
1669        mpz_clear(x);
1670        a->s = 0;
1671        a=nlShort1(a);
1672        break;
1673      }
1674      case 3:
1675      {
1676        if ((long)b>0L)
1677          mpz_add_ui(a->z,a->z,SR_TO_INT(b));
1678        else
1679          mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
1680        a->s = 3;
1681        a=nlShort3_noinline(a);
1682        break;
1683      }
1684    }
1685    return;
1686  }
1687  else if (SR_HDL(a) & SR_INT)
1688  {
1689    number u=ALLOC_RNUMBER();
1690    #if defined(LDEBUG)
1691    u->debug=123456;
1692    #endif
1693    mpz_init(u->z);
1694    switch (b->s)
1695    {
1696      case 0:
1697      case 1:/* a:short, b:1 */
1698      {
1699        mpz_t x;
1700        mpz_init(x);
1701
1702        mpz_mul_si(x,b->n,SR_TO_INT(a));
1703        mpz_add(u->z,b->z,x);
1704        mpz_clear(x);
1705        // result cannot be 0, if coeffs are normalized
1706        mpz_init_set(u->n,b->n);
1707        u->s = 0;
1708        u=nlShort1(u);
1709        break;
1710      }
1711      case 3:
1712      {
1713        if ((long)a>0L)
1714          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1715        else
1716          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1717        // result cannot be 0, if coeffs are normalized
1718        u->s = 3;
1719        u=nlShort3_noinline(u);
1720        break;
1721      }
1722    }
1723    a=u;
1724  }
1725  else
1726  {
1727    switch (a->s)
1728    {
1729      case 0:
1730      case 1:
1731      {
1732        switch(b->s)
1733        {
1734          case 0:
1735          case 1: /* a:1 b:1 */
1736          {
1737            mpz_t x;
1738            mpz_t y;
1739            mpz_init(x);
1740            mpz_init(y);
1741            mpz_mul(x,b->z,a->n);
1742            mpz_mul(y,a->z,b->n);
1743            mpz_add(a->z,x,y);
1744            mpz_clear(x);
1745            mpz_clear(y);
1746            mpz_mul(a->n,a->n,b->n);
1747            a->s = 0;
1748            break;
1749          }
1750          case 3: /* a:1 b:3 */
1751          {
1752            mpz_t x;
1753            mpz_init(x);
1754            mpz_mul(x,b->z,a->n);
1755            mpz_add(a->z,a->z,x);
1756            mpz_clear(x);
1757            a->s = 0;
1758            break;
1759          }
1760        } /*switch (b->s) */
1761        a=nlShort1(a);
1762        break;
1763      }
1764      case 3:
1765      {
1766        switch(b->s)
1767        {
1768          case 0:
1769          case 1:/* a:3, b:1 */
1770          {
1771            mpz_t x;
1772            mpz_init(x);
1773            mpz_mul(x,a->z,b->n);
1774            mpz_add(a->z,b->z,x);
1775            mpz_clear(x);
1776            mpz_init_set(a->n,b->n);
1777            a->s = 0;
1778            a=nlShort1(a);
1779            break;
1780          }
1781          case 3:
1782          {
1783            mpz_add(a->z,a->z,b->z);
1784            a->s = 3;
1785            a=nlShort3_noinline(a);
1786            break;
1787          }
1788        }
1789        break;
1790      }
1791    }
1792  }
1793}
1794
1795number _nlSub_aNoImm_OR_bNoImm(number a, number b)
1796{
1797  number u=ALLOC_RNUMBER();
1798#if defined(LDEBUG)
1799  u->debug=123456;
1800#endif
1801  mpz_init(u->z);
1802  if (SR_HDL(a) & SR_INT)
1803  {
1804    switch (b->s)
1805    {
1806      case 0:
1807      case 1:/* a:short, b:1 */
1808      {
1809        mpz_t x;
1810        mpz_init(x);
1811        mpz_mul_si(x,b->n,SR_TO_INT(a));
1812        mpz_sub(u->z,x,b->z);
1813        mpz_clear(x);
1814        if (mpz_cmp_ui(u->z,(long)0)==0)
1815        {
1816          mpz_clear(u->z);
1817          FREE_RNUMBER(u);
1818          return INT_TO_SR(0);
1819        }
1820        if (mpz_cmp(u->z,b->n)==0)
1821        {
1822          mpz_clear(u->z);
1823          FREE_RNUMBER(u);
1824          return INT_TO_SR(1);
1825        }
1826        mpz_init_set(u->n,b->n);
1827        u->s = 0;
1828        break;
1829      }
1830      case 3:
1831      {
1832        if ((long)a>0L)
1833        {
1834          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
1835          mpz_neg(u->z,u->z);
1836        }
1837        else
1838        {
1839          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
1840          mpz_neg(u->z,u->z);
1841        }
1842        if (mpz_cmp_ui(u->z,(long)0)==0)
1843        {
1844          mpz_clear(u->z);
1845          FREE_RNUMBER(u);
1846          return INT_TO_SR(0);
1847        }
1848        u->s = 3;
1849        u=nlShort3(u);
1850        break;
1851      }
1852    }
1853  }
1854  else if (SR_HDL(b) & SR_INT)
1855  {
1856    switch (a->s)
1857    {
1858      case 0:
1859      case 1:/* b:short, a:1 */
1860      {
1861        mpz_t x;
1862        mpz_init(x);
1863        mpz_mul_si(x,a->n,SR_TO_INT(b));
1864        mpz_sub(u->z,a->z,x);
1865        mpz_clear(x);
1866        if (mpz_cmp_ui(u->z,(long)0)==0)
1867        {
1868          mpz_clear(u->z);
1869          FREE_RNUMBER(u);
1870          return INT_TO_SR(0);
1871        }
1872        if (mpz_cmp(u->z,a->n)==0)
1873        {
1874          mpz_clear(u->z);
1875          FREE_RNUMBER(u);
1876          return INT_TO_SR(1);
1877        }
1878        mpz_init_set(u->n,a->n);
1879        u->s = 0;
1880        break;
1881      }
1882      case 3:
1883      {
1884        if ((long)b>0L)
1885        {
1886          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
1887        }
1888        else
1889        {
1890          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
1891        }
1892        if (mpz_cmp_ui(u->z,(long)0)==0)
1893        {
1894          mpz_clear(u->z);
1895          FREE_RNUMBER(u);
1896          return INT_TO_SR(0);
1897        }
1898        u->s = 3;
1899        u=nlShort3(u);
1900        break;
1901      }
1902    }
1903  }
1904  else
1905  {
1906    switch (a->s)
1907    {
1908      case 0:
1909      case 1:
1910      {
1911        switch(b->s)
1912        {
1913          case 0:
1914          case 1:
1915          {
1916            mpz_t x;
1917            mpz_t y;
1918            mpz_init(x);
1919            mpz_init(y);
1920            mpz_mul(x,b->z,a->n);
1921            mpz_mul(y,a->z,b->n);
1922            mpz_sub(u->z,y,x);
1923            mpz_clear(x);
1924            mpz_clear(y);
1925            if (mpz_cmp_ui(u->z,(long)0)==0)
1926            {
1927              mpz_clear(u->z);
1928              FREE_RNUMBER(u);
1929              return INT_TO_SR(0);
1930            }
1931            mpz_init(u->n);
1932            mpz_mul(u->n,a->n,b->n);
1933            if (mpz_cmp(u->z,u->n)==0)
1934            {
1935              mpz_clear(u->z);
1936              mpz_clear(u->n);
1937              FREE_RNUMBER(u);
1938              return INT_TO_SR(1);
1939            }
1940            u->s = 0;
1941            break;
1942          }
1943          case 3: /* a:1, b:3 */
1944          {
1945            mpz_t x;
1946            mpz_init(x);
1947            mpz_mul(x,b->z,a->n);
1948            mpz_sub(u->z,a->z,x);
1949            mpz_clear(x);
1950            if (mpz_cmp_ui(u->z,(long)0)==0)
1951            {
1952              mpz_clear(u->z);
1953              FREE_RNUMBER(u);
1954              return INT_TO_SR(0);
1955            }
1956            if (mpz_cmp(u->z,a->n)==0)
1957            {
1958              mpz_clear(u->z);
1959              FREE_RNUMBER(u);
1960              return INT_TO_SR(1);
1961            }
1962            mpz_init_set(u->n,a->n);
1963            u->s = 0;
1964            break;
1965          }
1966        }
1967        break;
1968      }
1969      case 3:
1970      {
1971        switch(b->s)
1972        {
1973          case 0:
1974          case 1: /* a:3, b:1 */
1975          {
1976            mpz_t x;
1977            mpz_init(x);
1978            mpz_mul(x,a->z,b->n);
1979            mpz_sub(u->z,x,b->z);
1980            mpz_clear(x);
1981            if (mpz_cmp_ui(u->z,(long)0)==0)
1982            {
1983              mpz_clear(u->z);
1984              FREE_RNUMBER(u);
1985              return INT_TO_SR(0);
1986            }
1987            if (mpz_cmp(u->z,b->n)==0)
1988            {
1989              mpz_clear(u->z);
1990              FREE_RNUMBER(u);
1991              return INT_TO_SR(1);
1992            }
1993            mpz_init_set(u->n,b->n);
1994            u->s = 0;
1995            break;
1996          }
1997          case 3: /* a:3 , b:3 */
1998          {
1999            mpz_sub(u->z,a->z,b->z);
2000            if (mpz_cmp_ui(u->z,(long)0)==0)
2001            {
2002              mpz_clear(u->z);
2003              FREE_RNUMBER(u);
2004              return INT_TO_SR(0);
2005            }
2006            u->s = 3;
2007            u=nlShort3(u);
2008            break;
2009          }
2010        }
2011        break;
2012      }
2013    }
2014  }
2015  return u;
2016}
2017
2018// a and b are intermediate, but a*b not
2019number _nlMult_aImm_bImm_rNoImm(number a, number b)
2020{
2021  number u=ALLOC_RNUMBER();
2022#if defined(LDEBUG)
2023  u->debug=123456;
2024#endif
2025  u->s=3;
2026  mpz_init_set_si(u->z,SR_TO_INT(a));
2027  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
2028  return u;
2029}
2030
2031// a or b are not immediate
2032number _nlMult_aNoImm_OR_bNoImm(number a, number b)
2033{
2034  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2035  number u=ALLOC_RNUMBER();
2036#if defined(LDEBUG)
2037  u->debug=123456;
2038#endif
2039  mpz_init(u->z);
2040  if (SR_HDL(b) & SR_INT)
2041  {
2042    number x=a;
2043    a=b;
2044    b=x;
2045  }
2046  if (SR_HDL(a) & SR_INT)
2047  {
2048    u->s=b->s;
2049    if (u->s==1) u->s=0;
2050    if ((long)a>0L)
2051    {
2052      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2053    }
2054    else
2055    {
2056      if (a==INT_TO_SR(-1))
2057      {
2058        mpz_set(u->z,b->z);
2059        mpz_neg(u->z,u->z);
2060        u->s=b->s;
2061      }
2062      else
2063      {
2064        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2065        mpz_neg(u->z,u->z);
2066      }
2067    }
2068    if (u->s<2)
2069    {
2070      if (mpz_cmp(u->z,b->n)==0)
2071      {
2072        mpz_clear(u->z);
2073        FREE_RNUMBER(u);
2074        return INT_TO_SR(1);
2075      }
2076      mpz_init_set(u->n,b->n);
2077    }
2078    else //u->s==3
2079    {
2080      u=nlShort3(u);
2081    }
2082  }
2083  else
2084  {
2085    mpz_mul(u->z,a->z,b->z);
2086    u->s = 0;
2087    if(a->s==3)
2088    {
2089      if(b->s==3)
2090      {
2091        u->s = 3;
2092      }
2093      else
2094      {
2095        if (mpz_cmp(u->z,b->n)==0)
2096        {
2097          mpz_clear(u->z);
2098          FREE_RNUMBER(u);
2099          return INT_TO_SR(1);
2100        }
2101        mpz_init_set(u->n,b->n);
2102      }
2103    }
2104    else
2105    {
2106      if(b->s==3)
2107      {
2108        if (mpz_cmp(u->z,a->n)==0)
2109        {
2110          mpz_clear(u->z);
2111          FREE_RNUMBER(u);
2112          return INT_TO_SR(1);
2113        }
2114        mpz_init_set(u->n,a->n);
2115      }
2116      else
2117      {
2118        mpz_init(u->n);
2119        mpz_mul(u->n,a->n,b->n);
2120        if (mpz_cmp(u->z,u->n)==0)
2121        {
2122          mpz_clear(u->z);
2123          mpz_clear(u->n);
2124          FREE_RNUMBER(u);
2125          return INT_TO_SR(1);
2126        }
2127      }
2128    }
2129  }
2130  return u;
2131}
2132
2133/*2
2134* copy a to b for mapping
2135*/
2136number nlCopyMap(number a, const coeffs src, const coeffs dst)
2137{
2138  assume( getCoeffType(dst) == ID );
2139  assume( getCoeffType(src) == ID );
2140
2141  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2142  {
2143    return a;
2144  }
2145  return _nlCopy_NoImm(a);
2146}
2147
2148nMapFunc nlSetMap(const coeffs src, const coeffs dst)
2149{
2150  assume( getCoeffType(dst) == ID );
2151//  assume( getCoeffType(src) == ID );
2152
2153  if (nCoeff_is_Q(src))
2154  {
2155    return ndCopyMap;
2156  }
2157  if (nCoeff_is_Zp(src))
2158  {
2159    return nlMapP;
2160  }
2161  if (nCoeff_is_R(src))
2162  {
2163    return nlMapR;
2164  }
2165  if (nCoeff_is_long_R(src))
2166  {
2167    return nlMapLongR; /* long R -> Q */
2168  }
2169#ifdef HAVE_RINGS
2170  if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
2171  {
2172    return nlMapGMP;
2173  }
2174  if (nCoeff_is_Ring_2toM(src))
2175  {
2176    return nlMapMachineInt;
2177  }
2178#endif
2179  return NULL;
2180}
2181/*2
2182* z := i
2183*/
2184number nlRInit (long i)
2185{
2186  number z=ALLOC_RNUMBER();
2187#if defined(LDEBUG)
2188  z->debug=123456;
2189#endif
2190  mpz_init_set_si(z->z,i);
2191  z->s = 3;
2192  return z;
2193}
2194
2195/*2
2196* z := i/j
2197*/
2198number nlInit2 (int i, int j, const coeffs r)
2199{
2200  number z=ALLOC_RNUMBER();
2201#if defined(LDEBUG)
2202  z->debug=123456;
2203#endif
2204  mpz_init_set_si(z->z,(long)i);
2205  mpz_init_set_si(z->n,(long)j);
2206  z->s = 0;
2207  nlNormalize(z,r);
2208  return z;
2209}
2210
2211number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
2212{
2213  number z=ALLOC_RNUMBER();
2214#if defined(LDEBUG)
2215  z->debug=123456;
2216#endif
2217  mpz_init_set(z->z,i);
2218  mpz_init_set(z->n,j);
2219  z->s = 0;
2220  nlNormalize(z,r);
2221  return z;
2222}
2223
2224#else // DO_LINLINE
2225
2226// declare immedate routines
2227number nlRInit (long i);
2228BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2229number  _nlCopy_NoImm(number a);
2230number  _nlNeg_NoImm(number a);
2231number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2232void    _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b);
2233number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2234number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2235number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2236
2237#endif
2238
2239/***************************************************************
2240 *
2241 * Routines which might be inlined by p_Numbers.h
2242 *
2243 *******************************************************************/
2244#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2245
2246// routines which are always inlined/static
2247
2248/*2
2249* a = b ?
2250*/
2251LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
2252{
2253  nlTest(a, r);
2254  nlTest(b, r);
2255// short - short
2256  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2257  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2258}
2259
2260LINLINE number nlInit (long i, const coeffs r)
2261{
2262  number n;
2263  LONG ii=(LONG)i;
2264  if ( ((ii << 3) >> 3) == ii ) n=INT_TO_SR(ii);
2265  else                          n=nlRInit(ii);
2266  nlTest(n, r);
2267  return n;
2268}
2269
2270/*2
2271* a == 1 ?
2272*/
2273LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
2274{
2275#ifdef LDEBUG
2276  if (a==NULL) return FALSE;
2277  nlTest(a, r);
2278#endif
2279  return (a==INT_TO_SR(1));
2280}
2281
2282LINLINE BOOLEAN nlIsZero (number a, const coeffs)
2283{
2284  #if 0
2285  if (a==INT_TO_SR(0)) return TRUE;
2286  if ((SR_HDL(a) & SR_INT)||(a==NULL)) return FALSE;
2287  if (mpz_cmp_si(a->z,(long)0)==0)
2288  {
2289    printf("gmp-0 in nlIsZero\n");
2290    dErrorBreak();
2291    return TRUE;
2292  }
2293  return FALSE;
2294  #else
2295  return (a==INT_TO_SR(0));
2296  #endif
2297}
2298
2299/*2
2300* copy a to b
2301*/
2302LINLINE number nlCopy(number a, const coeffs)
2303{
2304  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2305  {
2306    return a;
2307  }
2308  return _nlCopy_NoImm(a);
2309}
2310
2311
2312/*2
2313* delete a
2314*/
2315LINLINE void nlDelete (number * a, const coeffs r)
2316{
2317  if (*a!=NULL)
2318  {
2319    nlTest(*a, r);
2320    if ((SR_HDL(*a) & SR_INT)==0)
2321    {
2322      _nlDelete_NoImm(a);
2323    }
2324    *a=NULL;
2325  }
2326}
2327
2328/*2
2329* za:= - za
2330*/
2331LINLINE number nlNeg (number a, const coeffs R)
2332{
2333  nlTest(a, R);
2334  if(SR_HDL(a) &SR_INT)
2335  {
2336    int r=SR_TO_INT(a);
2337    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2338    else               a=INT_TO_SR(-r);
2339    return a;
2340  }
2341  a = _nlNeg_NoImm(a);
2342  nlTest(a, R);
2343  return a;
2344
2345}
2346
2347/*2
2348* u:= a + b
2349*/
2350LINLINE number nlAdd (number a, number b, const coeffs R)
2351{
2352  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2353  {
2354    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2355    if ( ((r << 1) >> 1) == r )
2356      return (number)(long)r;
2357    else
2358      return nlRInit(SR_TO_INT(r));
2359  }
2360  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2361  nlTest(u, R);
2362  return u;
2363}
2364
2365number nlShort1(number a);
2366number nlShort3_noinline(number x);
2367
2368LINLINE void nlInpAdd(number &a, number b, const coeffs r)
2369{
2370  // a=a+b
2371  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2372  {
2373    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2374    if ( ((r << 1) >> 1) == r )
2375      a=(number)(long)r;
2376    else
2377      a=nlRInit(SR_TO_INT(r));
2378  }
2379  else
2380  {
2381    _nlInpAdd_aNoImm_OR_bNoImm(a,b);
2382    nlTest(a,r);
2383  }
2384}
2385
2386LINLINE number nlMult (number a, number b, const coeffs R)
2387{
2388  nlTest(a, R);
2389  nlTest(b, R);
2390  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2391  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2392  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2393  {
2394    LONG r=(SR_HDL(a)-1L)*(SR_HDL(b)>>1);
2395    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2396    {
2397      number u=((number) ((r>>1)+SR_INT));
2398      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2399      return nlRInit(SR_HDL(u)>>2);
2400    }
2401    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2402    nlTest(u, R);
2403    return u;
2404
2405  }
2406  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2407  nlTest(u, R);
2408  return u;
2409
2410}
2411
2412
2413/*2
2414* u:= a - b
2415*/
2416LINLINE number nlSub (number a, number b, const coeffs r)
2417{
2418  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2419  {
2420    LONG r=SR_HDL(a)-SR_HDL(b)+1;
2421    if ( ((r << 1) >> 1) == r )
2422    {
2423      return (number)(long)r;
2424    }
2425    else
2426      return nlRInit(SR_TO_INT(r));
2427  }
2428  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2429  nlTest(u, r);
2430  return u;
2431
2432}
2433
2434LINLINE void nlInpMult(number &a, number b, const coeffs r)
2435{
2436  number aa=a;
2437  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2438  {
2439    number n=nlMult(aa,b,r);
2440    nlDelete(&a,r);
2441    a=n;
2442  }
2443  else
2444  {
2445    mpz_mul(aa->z,a->z,b->z);
2446    if (aa->s==3)
2447    {
2448      if(b->s!=3)
2449      {
2450        mpz_init_set(a->n,b->n);
2451        a->s=0;
2452      }
2453    }
2454    else
2455    {
2456      if(b->s!=3)
2457      {
2458        mpz_mul(a->n,a->n,b->n);
2459      }
2460      a->s=0;
2461    }
2462  }
2463}
2464#endif // DO_LINLINE
2465
2466#ifndef P_NUMBERS_H
2467
2468static void nlMPZ(mpz_t m, number &n, const coeffs r)
2469{
2470  assume( getCoeffType(r) == ID );
2471   
2472  nlTest(n, r);
2473  nlNormalize(n, r);
2474  if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n));     /* n fits in an int */
2475  else             mpz_init_set(m, (mpz_ptr)n->z);
2476}
2477
2478
2479static number nlInitMPZ(mpz_t m, const coeffs)
2480{
2481  number z = ALLOC_RNUMBER();
2482  mpz_init_set(z->z, m);
2483  mpz_init_set_ui(z->n, 1);
2484  z->s = 3;
2485  return z;   
2486}
2487
2488
2489void nlInpGcd(number &a, number b, const coeffs r)
2490{
2491  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2492  {
2493    number n=nlGcd(a,b,r);
2494    nlDelete(&a,r);
2495    a=n;
2496  }
2497  else
2498  {
2499    mpz_gcd(a->z,a->z,b->z);
2500    a=nlShort3_noinline(a);
2501  }
2502}
2503
2504void nlInpIntDiv(number &a, number b, const coeffs r)
2505{
2506  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2507  {
2508    number n=nlIntDiv(a,b, r);
2509    nlDelete(&a,r);
2510    a=n;
2511  }
2512  else
2513  {
2514    if (mpz_isNeg(a->z))
2515    {
2516      if (mpz_isNeg(b->z))
2517      {
2518        mpz_add(a->z,a->z,b->z);
2519      }
2520      else
2521      {
2522        mpz_sub(a->z,a->z,b->z);
2523      }
2524      mpz_add_ui(a->z,a->z,1);
2525    }
2526    else
2527    {
2528      if (mpz_isNeg(b->z))
2529      {
2530        mpz_sub(a->z,a->z,b->z);
2531      }
2532      else
2533      {
2534        mpz_add(a->z,a->z,b->z);
2535      }
2536      mpz_sub_ui(a->z,a->z,1);
2537    }
2538    MPZ_DIV(a->z,a->z,b->z);
2539    a=nlShort3_noinline(a);
2540  }
2541}
2542
2543number nlFarey(number nN, number nP, const coeffs r)
2544{
2545  mpz_t tmp; mpz_init(tmp);
2546  mpz_t A,B,C,D,E,N,P;
2547  if (SR_HDL(nN) & SR_INT) mpz_init_set_si(N,SR_TO_INT(nN));
2548  else                     mpz_init_set(N,nN->z);
2549  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2550  else                     mpz_init_set(P,nP->z);
2551  assume(!mpz_isNeg(P));
2552  if (mpz_isNeg(N))  mpz_add(N,N,P);
2553  mpz_init_set_si(A,(long)0);
2554  mpz_init_set_ui(B,(unsigned long)1);
2555  mpz_init_set_si(C,(long)0);
2556  mpz_init(D);
2557  mpz_init_set(E,P);
2558  number z=INT_TO_SR(0);
2559  while(mpz_cmp_si(N,(long)0)!=0)
2560  {
2561    mpz_mul(tmp,N,N);
2562    mpz_add(tmp,tmp,tmp);
2563    if (mpz_cmp(tmp,P)<0)
2564    {
2565       // return N/B
2566       z=ALLOC_RNUMBER();
2567       #ifdef LDEBUG
2568       z->debug=123456;
2569       #endif
2570       if (mpz_isNeg(B))
2571       {
2572         mpz_neg(B,B);
2573         mpz_neg(N,N);
2574       }
2575       mpz_init_set(z->z,N);
2576       mpz_init_set(z->n,B);
2577       z->s = 0;
2578       nlNormalize(z,r);
2579       break;
2580    }
2581    //mpz_mod(D,E,N);
2582    //mpz_div(tmp,E,N);
2583    mpz_divmod(tmp,D,E,N);
2584    mpz_mul(tmp,tmp,B);
2585    mpz_sub(C,A,tmp);
2586    mpz_set(E,N);
2587    mpz_set(N,D);
2588    mpz_set(A,B);
2589    mpz_set(B,C);
2590  }
2591  mpz_clear(tmp);
2592  mpz_clear(A);
2593  mpz_clear(B);
2594  mpz_clear(C);
2595  mpz_clear(D);
2596  mpz_clear(E);
2597  mpz_clear(N);
2598  mpz_clear(P);
2599  return z;
2600}
2601
2602void    nlCoeffWrite  (const coeffs, BOOLEAN /*details*/)
2603{
2604  PrintS("//   characteristic : 0\n");
2605}
2606
2607number   nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
2608// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
2609{
2610#ifdef HAVE_FACTORY
2611  setCharacteristic( 0 ); // only in char 0
2612  CFArray X(rl), Q(rl);
2613  int i;
2614  for(i=rl-1;i>=0;i--)
2615  {
2616    X[i]=C->convSingNFactoryN(x[i],FALSE,C); // may be larger MAX_INT
2617    Q[i]=C->convSingNFactoryN(q[i],FALSE,C); // may be larger MAX_INT
2618  }
2619  CanonicalForm xnew,qnew;
2620  chineseRemainder(X,Q,xnew,qnew);
2621  number n=C->convFactoryNSingN(xnew,C);
2622  number p=C->convFactoryNSingN(qnew,C);
2623  number p2=nlIntDiv(p,nlInit(2, C),C);
2624  if (nlGreater(n,p2,C))
2625  {
2626     number n2=nlSub(n,p,C);
2627     nlDelete(&n,C);
2628     n=n2;
2629  }
2630  nlDelete(&p,C);
2631  nlDelete(&p2,C);
2632  return n;
2633#else
2634  WerrorS("not implemented");
2635  return nlInit(0, C);
2636#endif
2637}
2638
2639static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2640{
2641  assume(cf != NULL);
2642  assume(getCoeffType(cf) == ID);
2643  // all coeffs are given by integers!!!
2644
2645  // part 1, find a small candidate for gcd
2646  number cand1,cand;
2647  int s1,s;
2648  s=2147483647; // max. int
2649  numberCollectionEnumerator.Reset();
2650  do
2651  {
2652    cand1= numberCollectionEnumerator.Current();
2653    if (SR_HDL(cand1)&SR_INT) { cand=cand1;break;}
2654    s1=mpz_size1(cand1->z);
2655    if (s>s1)
2656    {
2657      cand=cand1;
2658      s=s1;
2659    }
2660  } while (numberCollectionEnumerator.MoveNext() ); 
2661 
2662  // part 2: compute gcd(cand,all coeffs)
2663  numberCollectionEnumerator.Reset();
2664  do
2665  {
2666    nlNormalize(numberCollectionEnumerator.Current(),cf);
2667    cand1= nlGcd(cand,numberCollectionEnumerator.Current(),cf);
2668    nlDelete(&cand,cf);
2669    cand=cand1;
2670    if(nlIsOne(cand,cf)) { c=cand; return; }
2671  } while (numberCollectionEnumerator.MoveNext() ); 
2672 
2673  // part3: all coeffs = all coeffs / cand
2674  c=cand;
2675  numberCollectionEnumerator.Reset();
2676  do
2677  {
2678    number t=nlIntDiv(numberCollectionEnumerator.Current(),cand,cf);
2679    nlDelete(&numberCollectionEnumerator.Current(),cf);
2680    numberCollectionEnumerator.Current()=t;
2681  } while (numberCollectionEnumerator.MoveNext() ); 
2682 
2683}
2684
2685static void nlClearDenominators(ICoeffsEnumerator& /*numberCollectionEnumerator*/, number& c, const coeffs cf)
2686{
2687  assume(cf != NULL);
2688  assume(getCoeffType(cf) == ID);
2689  // all coeffs are given by integers!!!
2690
2691  c = n_Init(1, cf);
2692  assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
2693}
2694
2695BOOLEAN nlInitChar(coeffs r, void*)
2696{
2697  assume( getCoeffType(r) == ID );
2698  //const int ch = (int)(long)(p);
2699
2700  r->cfKillChar=NULL;
2701  r->nCoeffIsEqual=ndCoeffIsEqual;
2702  r->cfKillChar = ndKillChar; /* dummy */
2703
2704  r->cfInitMPZ = nlInitMPZ;
2705  r->cfMPZ  = nlMPZ;
2706   
2707  r->cfMult  = nlMult;
2708  r->cfSub   = nlSub;
2709  r->cfAdd   = nlAdd;
2710  r->cfDiv   = nlDiv;
2711  r->cfIntDiv= nlIntDiv;
2712  r->cfIntMod= nlIntMod;
2713  r->cfExactDiv= nlExactDiv;
2714  r->cfInit = nlInit;
2715  r->cfSize  = nlSize;
2716  r->cfInt  = nlInt;
2717   
2718  r->cfChineseRemainder=nlChineseRemainder;
2719  r->cfFarey=nlFarey;
2720  #ifdef HAVE_RINGS
2721  //r->cfDivComp = NULL; // only for ring stuff
2722  //r->cfIsUnit = NULL; // only for ring stuff
2723  //r->cfGetUnit = NULL; // only for ring stuff
2724  //r->cfExtGcd = NULL; // only for ring stuff
2725  //r->cfDivBy = NULL; // only for ring stuff
2726  #endif
2727  r->cfNeg   = nlNeg;
2728  r->cfInvers= nlInvers;
2729  r->cfCopy  = nlCopy;
2730  r->cfRePart = nlCopy;
2731  //r->cfImPart = ndReturn0;
2732  r->cfWriteLong = nlWrite;
2733  r->cfRead = nlRead;
2734  r->cfNormalize=nlNormalize;
2735  r->cfGreater = nlGreater;
2736  r->cfEqual = nlEqual;
2737  r->cfIsZero = nlIsZero;
2738  r->cfIsOne = nlIsOne;
2739  r->cfIsMOne = nlIsMOne;
2740  r->cfGreaterZero = nlGreaterZero;
2741  r->cfPower = nlPower;
2742  r->cfGetDenom = nlGetDenom;
2743  r->cfGetNumerator = nlGetNumerator;
2744  r->cfGcd  = nlGcd;
2745  r->cfLcm  = nlLcm;
2746  r->cfDelete= nlDelete;
2747  r->cfSetMap = nlSetMap;
2748  //r->cfName = ndName;
2749  r->cfInpMult=nlInpMult;
2750  r->cfInit_bigint=nlCopyMap;
2751  r->cfCoeffWrite=nlCoeffWrite;
2752
2753  r->cfClearContent = nlClearContent;
2754  r->cfClearDenominators = nlClearDenominators;
2755 
2756#ifdef LDEBUG
2757  // debug stuff
2758  r->cfDBTest=nlDBTest;
2759#endif
2760#ifdef HAVE_FACTORY
2761  r->convSingNFactoryN=nlConvSingNFactoryN;
2762  r->convFactoryNSingN=nlConvFactoryNSingN;
2763#endif
2764
2765  // the variables: general stuff (required)
2766  r->nNULL = INT_TO_SR(0);
2767  //r->type = n_Q;
2768  r->ch = 0;
2769  r->has_simple_Alloc=FALSE;
2770  r->has_simple_Inverse=FALSE;
2771
2772  // variables for this type of coeffs:
2773  // (none)
2774  return FALSE;
2775}
2776#if 0
2777number nlMod(number a, number b)
2778{
2779  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
2780  {
2781    int bi=SR_TO_INT(b);
2782    int ai=SR_TO_INT(a);
2783    int bb=ABS(bi);
2784    int c=ai%bb;
2785    if (c<0)  c+=bb;
2786    return (INT_TO_SR(c));
2787  }
2788  number al;
2789  number bl;
2790  if (SR_HDL(a))&SR_INT)
2791    al=nlRInit(SR_TO_INT(a));
2792  else
2793    al=nlCopy(a);
2794  if (SR_HDL(b))&SR_INT)
2795    bl=nlRInit(SR_TO_INT(b));
2796  else
2797    bl=nlCopy(b);
2798  number r=nlRInit(0);
2799  mpz_mod(r->z,al->z,bl->z);
2800  nlDelete(&al);
2801  nlDelete(&bl);
2802  if (mpz_size1(&r->z)<=MP_SMALL)
2803  {
2804    LONG ui=(int)mpz_get_si(&r->z);
2805    if ((((ui<<3)>>3)==ui)
2806    && (mpz_cmp_si(x->z,(long)ui)==0))
2807    {
2808      mpz_clear(&r->z);
2809      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
2810      r=INT_TO_SR(ui);
2811    }
2812  }
2813  return r;
2814}
2815#endif
2816#endif // not P_NUMBERS_H
2817#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.