# source:git/kernel/longrat.cc@37a003

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