source: git/coeffs/longrat.cc @ 7d90aa

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