source: git/kernel/longrat.cc @ 815d74

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