source: git/kernel/longrat.cc @ 1a97da

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