source: git/libpolys/coeffs/longrat.cc @ 8167af

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