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

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