source: git/kernel/longrat.cc @ 1919112

spielwiese
Last change on this file since 1919112 was 1919112, checked in by Hans Schoenemann <hannes@…>, 14 years ago
simplified comments etc git-svn-id: file:///usr/local/Singular/svn/trunk@13041 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 55.4 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.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  nlNormalize(tmp);
534  return tmp;
535}
536
537/*
538* 1/a
539*/
540number nlInvers(number a)
541{
542#ifdef LDEBUG
543  nlTest(a);
544#endif
545  number n;
546  if (SR_HDL(a) & SR_INT)
547  {
548    if ((a==INT_TO_SR(1)) || (a==INT_TO_SR(-1)))
549    {
550      return a;
551    }
552    if (nlIsZero(a))
553    {
554      WerrorS(nDivBy0);
555      return INT_TO_SR(0);
556    }
557    n=(number)omAllocBin(rnumber_bin);
558#if defined(LDEBUG)
559    n->debug=123456;
560#endif
561    n->s=1;
562    if ((long)a>0L)
563    {
564      mpz_init_set_si(n->z,(long)1);
565      mpz_init_set_si(n->n,(long)SR_TO_INT(a));
566    }
567    else
568    {
569      mpz_init_set_si(n->z,(long)-1);
570      mpz_init_set_si(n->n,(long)-SR_TO_INT(a));
571    }
572#ifdef LDEBUG
573    nlTest(n);
574#endif
575    return n;
576  }
577  n=(number)omAllocBin(rnumber_bin);
578#if defined(LDEBUG)
579  n->debug=123456;
580#endif
581  {
582    n->s=a->s;
583    mpz_init_set(n->n,a->z);
584    switch (a->s)
585    {
586      case 0:
587      case 1:
588              mpz_init_set(n->z,a->n);
589              if (mpz_isNeg(n->n)) /* && n->s<2*/
590              {
591                mpz_neg(n->z,n->z);
592                mpz_neg(n->n,n->n);
593              }
594              if (mpz_cmp_si(n->n,(long)1)==0)
595              {
596                mpz_clear(n->n);
597                n->s=3;
598                if (mpz_size1(n->z)<=MP_SMALL)
599                {
600                  int ui=(int)mpz_get_si(n->z);
601                  if ((((ui<<3)>>3)==ui)
602                  && (mpz_cmp_si(n->z,(long)ui)==0))
603                  {
604                    mpz_clear(n->z);
605#if defined(LDEBUG)
606                    n->debug=654324;
607#endif
608                    omFreeBin((ADDRESS)n, rnumber_bin);
609                    n=INT_TO_SR(ui);
610                  }
611                }
612              }
613              break;
614      case 3:
615              n->s=1;
616              if (mpz_isNeg(n->n)) /* && n->s<2*/
617              {
618                mpz_neg(n->n,n->n);
619                mpz_init_set_si(n->z,(long)-1);
620              }
621              else
622              {
623                mpz_init_set_si(n->z,(long)1);
624              }
625              break;
626    }
627  }
628#ifdef LDEBUG
629  nlTest(n);
630#endif
631  return n;
632}
633
634
635/*2
636* u := a / b in Z, if b | a (else undefined)
637*/
638number   nlExactDiv(number a, number b)
639{
640  if (b==INT_TO_SR(0))
641  {
642    WerrorS(nDivBy0);
643    return INT_TO_SR(0);
644  }
645  if (a==INT_TO_SR(0))
646    return INT_TO_SR(0);
647  number u;
648  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
649  {
650    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
651    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
652    {
653      return nlRInit(1<<28);
654    }
655    int aa=SR_TO_INT(a);
656    int bb=SR_TO_INT(b);
657    return INT_TO_SR(aa/bb);
658  }
659  number bb=NULL;
660  if (SR_HDL(b) & SR_INT)
661  {
662    bb=nlRInit(SR_TO_INT(b));
663    b=bb;
664  }
665  u=(number)omAllocBin(rnumber_bin);
666#if defined(LDEBUG)
667  u->debug=123456;
668#endif
669  mpz_init(u->z);
670  /* u=a/b */
671  u->s = 3;
672  MPZ_EXACTDIV(u->z,a->z,b->z);
673  if (bb!=NULL)
674  {
675    mpz_clear(bb->z);
676#if defined(LDEBUG)
677    bb->debug=654324;
678#endif
679    omFreeBin((ADDRESS)bb, rnumber_bin);
680  }
681  if (mpz_size1(u->z)<=MP_SMALL)
682  {
683    int ui=(int)mpz_get_si(u->z);
684    if ((((ui<<3)>>3)==ui)
685    && (mpz_cmp_si(u->z,(long)ui)==0))
686    {
687      mpz_clear(u->z);
688#if defined(LDEBUG)
689      u->debug=654324;
690#endif
691      omFreeBin((ADDRESS)u, rnumber_bin);
692      u=INT_TO_SR(ui);
693    }
694  }
695#ifdef LDEBUG
696  nlTest(u);
697#endif
698  return u;
699}
700
701/*2
702* u := a / b in Z
703*/
704number nlIntDiv (number a, number b)
705{
706  if (b==INT_TO_SR(0))
707  {
708    WerrorS(nDivBy0);
709    return INT_TO_SR(0);
710  }
711  if (a==INT_TO_SR(0))
712    return INT_TO_SR(0);
713  number u;
714  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
715  {
716    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
717    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
718    {
719      return nlRInit(1<<28);
720    }
721    int aa=SR_TO_INT(a);
722    int bb=SR_TO_INT(b);
723    return INT_TO_SR(aa/bb);
724  }
725  if (SR_HDL(a) & SR_INT)
726  {
727    /* the small int -(1<<28) divided by 2^28 is 1   */
728    if (a==INT_TO_SR(-(1<<28)))
729    {
730      if(mpz_cmp_si(b->z,(long)(1<<28))==0)
731      {
732        return INT_TO_SR(-1);
733      }
734    }
735    /* a is a small and b is a large int: -> 0 */
736    return INT_TO_SR(0);
737  }
738  number bb=NULL;
739  if (SR_HDL(b) & SR_INT)
740  {
741    bb=nlRInit(SR_TO_INT(b));
742    b=bb;
743  }
744  u=(number)omAllocBin(rnumber_bin);
745#if defined(LDEBUG)
746  u->debug=123456;
747#endif
748  assume(a->s==3);
749  assume(b->s==3);
750  mpz_init_set(u->z,a->z);
751  /* u=u/b */
752  u->s = 3;
753  MPZ_DIV(u->z,u->z,b->z);
754  if (bb!=NULL)
755  {
756    mpz_clear(bb->z);
757#if defined(LDEBUG)
758    bb->debug=654324;
759#endif
760    omFreeBin((ADDRESS)bb, rnumber_bin);
761  }
762  if (mpz_size1(u->z)<=MP_SMALL)
763  {
764    int ui=(int)mpz_get_si(u->z);
765    if ((((ui<<3)>>3)==ui)
766    && (mpz_cmp_si(u->z,(long)ui)==0))
767    {
768      mpz_clear(u->z);
769#if defined(LDEBUG)
770      u->debug=654324;
771#endif
772      omFreeBin((ADDRESS)u, rnumber_bin);
773      u=INT_TO_SR(ui);
774    }
775  }
776#ifdef LDEBUG
777  nlTest(u);
778#endif
779  return u;
780}
781
782/*2
783* u := a mod b in Z, u>=0
784*/
785number nlIntMod (number a, number b)
786{
787  if (b==INT_TO_SR(0))
788  {
789    WerrorS(nDivBy0);
790    return INT_TO_SR(0);
791  }
792  if (a==INT_TO_SR(0))
793    return INT_TO_SR(0);
794  number u;
795  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
796  {
797    if ((long)a>0L)
798    {
799      if ((long)b>0L)
800        return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
801      else
802        return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
803    }
804    else
805    {
806      if ((long)b>0L)
807      {
808        int i=(-SR_TO_INT(a))%SR_TO_INT(b);
809        if ( i != 0 ) i = (SR_TO_INT(b))-i;
810        return INT_TO_SR(i);
811      }
812      else
813      {
814        int i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
815        if ( i != 0 ) i = (-SR_TO_INT(b))-i;
816        return INT_TO_SR(i);
817      }
818    }
819  }
820  if (SR_HDL(a) & SR_INT)
821  {
822    /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
823    if ((long)a<0L)
824    {
825      if (mpz_isNeg(b->z))
826        return nlSub(a,b);
827      /*else*/
828        return nlAdd(a,b);
829    }
830    /*else*/
831      return a;
832  }
833  number bb=NULL;
834  if (SR_HDL(b) & SR_INT)
835  {
836    bb=nlRInit(SR_TO_INT(b));
837    b=bb;
838  }
839  u=(number)omAllocBin(rnumber_bin);
840#if defined(LDEBUG)
841  u->debug=123456;
842#endif
843  mpz_init(u->z);
844  u->s = 3;
845  mpz_mod(u->z,a->z,b->z);
846  if (bb!=NULL)
847  {
848    mpz_clear(bb->z);
849#if defined(LDEBUG)
850    bb->debug=654324;
851#endif
852    omFreeBin((ADDRESS)bb, rnumber_bin);
853  }
854  if (mpz_isNeg(u->z))
855  {
856    if (mpz_isNeg(b->z))
857      mpz_sub(u->z,u->z,b->z);
858    else
859      mpz_add(u->z,u->z,b->z);
860  }
861  if (mpz_size1(u->z)<=MP_SMALL)
862  {
863    int ui=(int)mpz_get_si(u->z);
864    if ((((ui<<3)>>3)==ui)
865    && (mpz_cmp_si(u->z,(long)ui)==0))
866    {
867      mpz_clear(u->z);
868#if defined(LDEBUG)
869      u->debug=654324;
870#endif
871      omFreeBin((ADDRESS)u, rnumber_bin);
872      u=INT_TO_SR(ui);
873    }
874  }
875#ifdef LDEBUG
876  nlTest(u);
877#endif
878  return u;
879}
880
881/*2
882* u := a / b
883*/
884number nlDiv (number a, number b)
885{
886  number u;
887  if (nlIsZero(b))
888  {
889    WerrorS(nDivBy0);
890    return INT_TO_SR(0);
891  }
892  u=(number)omAllocBin(rnumber_bin);
893  u->s=0;
894#if defined(LDEBUG)
895  u->debug=123456;
896#endif
897// ---------- short / short ------------------------------------
898  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
899  {
900    int i=SR_TO_INT(a);
901    int j=SR_TO_INT(b);
902    int r=i%j;
903    if (r==0)
904    {
905      omFreeBin((ADDRESS)u, rnumber_bin);
906      return INT_TO_SR(i/j);
907    }
908    mpz_init_set_si(u->z,(long)i);
909    mpz_init_set_si(u->n,(long)j);
910  }
911  else
912  {
913    mpz_init(u->z);
914// ---------- short / long ------------------------------------
915    if (SR_HDL(a) & SR_INT)
916    {
917      // short a / (z/n) -> (a*n)/z
918      if (b->s<2)
919      {
920        mpz_mul_si(u->z,b->n,SR_TO_INT(a));
921      }
922      else
923      // short a / long z -> a/z
924      {
925        mpz_set_si(u->z,SR_TO_INT(a));
926      }
927      if (mpz_cmp(u->z,b->z)==0)
928      {
929        mpz_clear(u->z);
930        omFreeBin((ADDRESS)u, rnumber_bin);
931        return INT_TO_SR(1);
932      }
933      mpz_init_set(u->n,b->z);
934    }
935// ---------- long / short ------------------------------------
936    else if (SR_HDL(b) & SR_INT)
937    {
938      mpz_set(u->z,a->z);
939      // (z/n) / b -> z/(n*b)
940      if (a->s<2)
941      {
942        mpz_init_set(u->n,a->n);
943        if ((long)b>0L)
944          mpz_mul_ui(u->n,u->n,SR_TO_INT(b));
945        else
946        {
947          mpz_mul_ui(u->n,u->n,-SR_TO_INT(b));
948          mpz_neg(u->z,u->z);
949        }
950      }
951      else
952      // long z / short b -> z/b
953      {
954        //mpz_set(u->z,a->z);
955        mpz_init_set_si(u->n,SR_TO_INT(b));
956      }
957    }
958// ---------- long / long ------------------------------------
959    else
960    {
961      mpz_set(u->z,a->z);
962      mpz_init_set(u->n,b->z);
963      if (a->s<2) mpz_mul(u->n,u->n,a->n);
964      if (b->s<2) mpz_mul(u->z,u->z,b->n);
965    }
966  }
967  if (mpz_isNeg(u->n))
968  {
969    mpz_neg(u->z,u->z);
970    mpz_neg(u->n,u->n);
971  }
972  if (mpz_cmp_si(u->n,(long)1)==0)
973  {
974    mpz_clear(u->n);
975    if (mpz_size1(u->z)<=MP_SMALL)
976    {
977      int ui=(int)mpz_get_si(u->z);
978      if ((((ui<<3)>>3)==ui)
979      && (mpz_cmp_si(u->z,(long)ui)==0))
980      {
981        mpz_clear(u->z);
982        omFreeBin((ADDRESS)u, rnumber_bin);
983        return INT_TO_SR(ui);
984      }
985    }
986    u->s=3;
987  }
988#ifdef LDEBUG
989  nlTest(u);
990#endif
991  return u;
992}
993
994/*2
995* u:= x ^ exp
996*/
997void nlPower (number x,int exp,number * u)
998{
999  *u = INT_TO_SR(0); // 0^e, e!=0
1000  if (!nlIsZero(x))
1001  {
1002#ifdef LDEBUG
1003    nlTest(x);
1004#endif
1005    number aa=NULL;
1006    if (SR_HDL(x) & SR_INT)
1007    {
1008      aa=nlRInit(SR_TO_INT(x));
1009      x=aa;
1010    }
1011    else if (x->s==0)
1012      nlNormalize(x);
1013    *u=(number)omAllocBin(rnumber_bin);
1014#if defined(LDEBUG)
1015    (*u)->debug=123456;
1016#endif
1017    mpz_init((*u)->z);
1018    mpz_pow_ui((*u)->z,x->z,(unsigned long)exp);
1019    if (x->s<2)
1020    {
1021      if (mpz_cmp_si(x->n,(long)1)==0)
1022      {
1023        x->s=3;
1024        mpz_clear(x->n);
1025      }
1026      else
1027      {
1028        mpz_init((*u)->n);
1029        mpz_pow_ui((*u)->n,x->n,(unsigned long)exp);
1030      }
1031    }
1032    (*u)->s = x->s;
1033    if (((*u)->s==3) && (mpz_size1((*u)->z)<=MP_SMALL))
1034    {
1035      int ui=(int)mpz_get_si((*u)->z);
1036      if ((((ui<<3)>>3)==ui)
1037      && (mpz_cmp_si((*u)->z,(long)ui)==0))
1038      {
1039        mpz_clear((*u)->z);
1040        omFreeBin((ADDRESS)*u, rnumber_bin);
1041        *u=INT_TO_SR(ui);
1042      }
1043    }
1044    if (aa!=NULL)
1045    {
1046      mpz_clear(aa->z);
1047      omFreeBin((ADDRESS)aa, rnumber_bin);
1048    }
1049  }
1050  else if (exp==0)
1051    *u = INT_TO_SR(1); // 0^0
1052#ifdef LDEBUG
1053  if (exp<0) Print("nlPower: neg. exp. %d\n",exp);
1054  nlTest(*u);
1055#endif
1056}
1057
1058
1059/*2
1060* za >= 0 ?
1061*/
1062BOOLEAN nlGreaterZero (number a)
1063{
1064#ifdef LDEBUG
1065  nlTest(a);
1066#endif
1067  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>1 /* represents number(0) */;
1068  return (!mpz_isNeg(a->z));
1069}
1070
1071/*2
1072* a > b ?
1073*/
1074BOOLEAN nlGreater (number a, number b)
1075{
1076#ifdef LDEBUG
1077  nlTest(a);
1078  nlTest(b);
1079#endif
1080  number r;
1081  BOOLEAN rr;
1082  r=nlSub(a,b);
1083  rr=(!nlIsZero(r)) && (nlGreaterZero(r));
1084  nlDelete(&r,currRing);
1085  return rr;
1086}
1087
1088/*2
1089* a == -1 ?
1090*/
1091BOOLEAN nlIsMOne (number a)
1092{
1093#ifdef LDEBUG
1094  if (a==NULL) return FALSE;
1095  nlTest(a);
1096#endif
1097  //if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1));
1098  //return FALSE;
1099  return  (a==INT_TO_SR(-1));
1100}
1101
1102/*2
1103* result =gcd(a,b)
1104*/
1105number nlGcd(number a, number b, const ring r)
1106{
1107  number result;
1108#ifdef LDEBUG
1109  nlTest(a);
1110  nlTest(b);
1111#endif
1112  //nlNormalize(a);
1113  //nlNormalize(b);
1114  if ((a==INT_TO_SR(1))||(a==INT_TO_SR(-1))
1115  ||  (b==INT_TO_SR(1))||(b==INT_TO_SR(-1)))
1116    return INT_TO_SR(1);
1117  if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
1118    return nlCopy(b);
1119  if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
1120    return nlCopy(a);
1121  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1122  {
1123    int i=SR_TO_INT(a);
1124    int j=SR_TO_INT(b);
1125    if((i==0)||(j==0))
1126      return INT_TO_SR(1);
1127    int l;
1128    i=ABS(i);
1129    j=ABS(j);
1130    do
1131    {
1132      l=i%j;
1133      i=j;
1134      j=l;
1135    } while (l!=0);
1136    return INT_TO_SR(i);
1137  }
1138  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1139  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1140  if (SR_HDL(a) & SR_INT)
1141  {
1142    unsigned long t=mpz_gcd_ui(NULL,b->z,ABS(SR_TO_INT(a)));
1143    return INT_TO_SR((int)t);
1144  }
1145  else
1146  if (SR_HDL(b) & SR_INT)
1147  {
1148    unsigned long t=mpz_gcd_ui(NULL,a->z,ABS(SR_TO_INT(b)));
1149    return INT_TO_SR((int)t);
1150  }
1151  else
1152  {
1153    result=(number)omAllocBin(rnumber_bin);
1154    mpz_init(result->z);
1155    mpz_gcd(result->z,a->z,b->z);
1156    if (mpz_size1(result->z)<=MP_SMALL)
1157    {
1158      int ui=(int)mpz_get_si(result->z);
1159      if ((((ui<<3)>>3)==ui)
1160      && (mpz_cmp_si(result->z,(long)ui)==0))
1161      {
1162        mpz_clear(result->z);
1163        omFreeBin((ADDRESS)result, rnumber_bin);
1164        return INT_TO_SR(ui);
1165      }
1166    }
1167    result->s = 3;
1168  #ifdef LDEBUG
1169    result->debug=123456;
1170    nlTest(result);
1171  #endif
1172  }
1173  return result;
1174}
1175
1176number nlShort1(number x) // assume x->s==0/1
1177{
1178  assume(x->s<2);
1179  if (mpz_cmp_ui(x->z,(long)0)==0)
1180  {
1181    _nlDelete_NoImm(&x);
1182    return INT_TO_SR(0);
1183  }
1184  if (x->s<2)
1185  {
1186    if (mpz_cmp(x->z,x->n)==0)
1187    {
1188      _nlDelete_NoImm(&x);
1189      return INT_TO_SR(1);
1190    }
1191  }
1192  return x;
1193}
1194number nlShort3(number x) // assume x->s==3
1195{
1196  assume(x->s==3);
1197  if ((mpz_cmp_ui(x->z,(long)0)==0)
1198  || (mpz_size1(x->z)<=MP_SMALL))
1199  {
1200    int ui=(int)mpz_get_si(x->z);
1201    if ((((ui<<3)>>3)==ui)
1202    && (mpz_cmp_si(x->z,(long)ui)==0))
1203    {
1204      mpz_clear(x->z);
1205      omFreeBin((ADDRESS)x, rnumber_bin);
1206      return INT_TO_SR(ui);
1207    }
1208  }
1209  return x;
1210}
1211/*2
1212* simplify x
1213*/
1214void nlNormalize (number &x)
1215{
1216  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1217    return;
1218#ifdef LDEBUG
1219  if (!nlTest(x)) { if (x->s!=3) x->s=1; return; }
1220#endif
1221  if (x->s==3)
1222  {
1223    x=nlShort3(x);
1224    return;
1225  }
1226  else if (x->s==0)
1227  {
1228    if (mpz_cmp_si(x->n,(long)1)==0)
1229    {
1230      mpz_clear(x->n);
1231      if (mpz_size1(x->z)<=MP_SMALL)
1232      {
1233        int ui=(int)mpz_get_si(x->z);
1234        if ((((ui<<3)>>3)==ui)
1235        && (mpz_cmp_si(x->z,(long)ui)==0))
1236        {
1237          mpz_clear(x->z);
1238#if defined(LDEBUG)
1239          x->debug=654324;
1240#endif
1241          omFreeBin((ADDRESS)x, rnumber_bin);
1242          x=INT_TO_SR(ui);
1243          return;
1244        }
1245      }
1246      x->s=3;
1247    }
1248    else
1249    {
1250      mpz_t gcd;
1251      mpz_init(gcd);
1252      mpz_gcd(gcd,x->z,x->n);
1253      x->s=1;
1254      if (mpz_cmp_si(gcd,(long)1)!=0)
1255      {
1256        mpz_t r;
1257        mpz_init(r);
1258        MPZ_EXACTDIV(r,x->z,gcd);
1259        mpz_set(x->z,r);
1260        MPZ_EXACTDIV(r,x->n,gcd);
1261        mpz_set(x->n,r);
1262        mpz_clear(r);
1263        if (mpz_cmp_si(x->n,(long)1)==0)
1264        {
1265          mpz_clear(x->n);
1266          if (mpz_size1(x->z)<=MP_SMALL)
1267          {
1268            int ui=(int)mpz_get_si(x->z);
1269            if ((((ui<<3)>>3)==ui)
1270            && (mpz_cmp_si(x->z,(long)ui)==0))
1271            {
1272              mpz_clear(x->z);
1273              mpz_clear(gcd);
1274#if defined(LDEBUG)
1275              x->debug=654324;
1276#endif
1277              omFreeBin((ADDRESS)x, rnumber_bin);
1278              x=INT_TO_SR(ui);
1279              return;
1280            }
1281          }
1282          x->s=3;
1283        }
1284      }
1285      mpz_clear(gcd);
1286    }
1287  }
1288#ifdef LDEBUG
1289  nlTest(x);
1290#endif
1291}
1292
1293/*2
1294* returns in result->z the lcm(a->z,b->n)
1295*/
1296number nlLcm(number a, number b, const ring r)
1297{
1298  number result;
1299#ifdef LDEBUG
1300  nlTest(a);
1301  nlTest(b);
1302#endif
1303  if ((SR_HDL(b) & SR_INT)
1304  || (b->s==3))
1305  {
1306    // b is 1/(b->n) => b->n is 1 => result is a
1307    return nlCopy(a);
1308  }
1309  result=(number)omAllocBin(rnumber_bin);
1310#if defined(LDEBUG)
1311  result->debug=123456;
1312#endif
1313  result->s=3;
1314  mpz_t gcd;
1315  mpz_init(gcd);
1316  mpz_init(result->z);
1317  if (SR_HDL(a) & SR_INT)
1318    mpz_gcd_ui(gcd,b->n,ABS(SR_TO_INT(a)));
1319  else
1320    mpz_gcd(gcd,a->z,b->n);
1321  if (mpz_cmp_si(gcd,(long)1)!=0)
1322  {
1323    mpz_t bt;
1324    mpz_init_set(bt,b->n);
1325    MPZ_EXACTDIV(bt,bt,gcd);
1326    if (SR_HDL(a) & SR_INT)
1327      mpz_mul_si(result->z,bt,SR_TO_INT(a));
1328    else
1329      mpz_mul(result->z,bt,a->z);
1330    mpz_clear(bt);
1331  }
1332  else
1333    if (SR_HDL(a) & SR_INT)
1334      mpz_mul_si(result->z,b->n,SR_TO_INT(a));
1335    else
1336      mpz_mul(result->z,b->n,a->z);
1337  mpz_clear(gcd);
1338  if (mpz_size1(result->z)<=MP_SMALL)
1339  {
1340    int ui=(int)mpz_get_si(result->z);
1341    if ((((ui<<3)>>3)==ui)
1342    && (mpz_cmp_si(result->z,(long)ui)==0))
1343    {
1344      mpz_clear(result->z);
1345      omFreeBin((ADDRESS)result, rnumber_bin);
1346      return INT_TO_SR(ui);
1347    }
1348  }
1349#ifdef LDEBUG
1350  nlTest(result);
1351#endif
1352  return result;
1353}
1354
1355int nlModP(number n, int p)
1356{
1357  if (SR_HDL(n) & SR_INT)
1358  {
1359    int i=SR_TO_INT(n);
1360    if (i<0) return (p-((-i)%p));
1361    return i%p;
1362  }
1363  int iz=(int)mpz_fdiv_ui(n->z,(unsigned long)p);
1364  if (n->s!=3)
1365  {
1366    int in=mpz_fdiv_ui(n->n,(unsigned long)p);
1367    #ifdef NV_OPS
1368    if (npPrimeM>NV_MAX_PRIME)
1369    return (int)((long)nvDiv((number)iz,(number)in));
1370    #endif
1371    return (int)((long)npDiv((number)iz,(number)in));
1372  }
1373  return iz;
1374}
1375
1376#ifdef HAVE_RINGS
1377/*2
1378* convert number i (from Q) to GMP and warn if denom != 1
1379*/
1380void nlGMP(number &i, number n)
1381{
1382  // Hier brauche ich einfach die GMP Zahl
1383#ifdef LDEBUG
1384  nlTest(i);
1385#endif
1386  nlNormalize(i);
1387  if (SR_HDL(i) & SR_INT)
1388  {
1389    mpz_set_si((mpz_ptr) n, (long) SR_TO_INT(i));
1390    return;
1391  }
1392  if (i->s!=3)
1393  {
1394     WarnS("Omitted denominator during coefficient mapping !");
1395  }
1396  mpz_set((mpz_ptr) n, i->z);
1397}
1398#endif
1399
1400/*2
1401* acces to denominator, other 1 for integers
1402*/
1403number   nlGetDenom(number &n, const ring r)
1404{
1405  if (!(SR_HDL(n) & SR_INT))
1406  {
1407    if (n->s==0)
1408    {
1409      nlNormalize(n);
1410    }
1411    if (!(SR_HDL(n) & SR_INT))
1412    {
1413      if (n->s!=3)
1414      {
1415        number u=(number)omAllocBin(rnumber_bin);
1416        u->s=3;
1417#if defined(LDEBUG)
1418        u->debug=123456;
1419#endif
1420        mpz_init_set(u->z,n->n);
1421        {
1422          int ui=(int)mpz_get_si(u->z);
1423          if ((((ui<<3)>>3)==ui)
1424          && (mpz_cmp_si(u->z,(long)ui)==0))
1425          {
1426            mpz_clear(u->z);
1427            omFreeBin((ADDRESS)u, rnumber_bin);
1428            return INT_TO_SR(ui);
1429          }
1430        }
1431        return u;
1432      }
1433    }
1434  }
1435  return INT_TO_SR(1);
1436}
1437
1438/*2
1439* acces to Nominator, nlCopy(n) for integers
1440*/
1441number   nlGetNumerator(number &n, const ring r)
1442{
1443  if (!(SR_HDL(n) & SR_INT))
1444  {
1445    if (n->s==0)
1446    {
1447      nlNormalize(n);
1448    }
1449    if (!(SR_HDL(n) & SR_INT))
1450    {
1451      number u=(number)omAllocBin(rnumber_bin);
1452#if defined(LDEBUG)
1453      u->debug=123456;
1454#endif
1455      u->s=3;
1456      mpz_init_set(u->z,n->z);
1457      if (n->s!=3)
1458      {
1459        int ui=(int)mpz_get_si(u->z);
1460        if ((((ui<<3)>>3)==ui)
1461        && (mpz_cmp_si(u->z,(long)ui)==0))
1462        {
1463          mpz_clear(u->z);
1464          omFreeBin((ADDRESS)u, rnumber_bin);
1465          return INT_TO_SR(ui);
1466        }
1467      }
1468      return u;
1469    }
1470  }
1471  return n; // imm. int
1472}
1473
1474/***************************************************************
1475 *
1476 * routines which are needed by Inline(d) routines
1477 *
1478 *******************************************************************/
1479BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
1480{
1481  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1482//  long - short
1483  BOOLEAN bo;
1484  if (SR_HDL(b) & SR_INT)
1485  {
1486    if (a->s!=0) return FALSE;
1487    number n=b; b=a; a=n;
1488  }
1489//  short - long
1490  if (SR_HDL(a) & SR_INT)
1491  {
1492    if (b->s!=0)
1493      return FALSE;
1494    if (((long)a > 0L) && (mpz_isNeg(b->z)))
1495      return FALSE;
1496    if (((long)a < 0L) && (!mpz_isNeg(b->z)))
1497      return FALSE;
1498    mpz_t  bb;
1499    mpz_init_set(bb,b->n);
1500    mpz_mul_si(bb,bb,(long)SR_TO_INT(a));
1501    bo=(mpz_cmp(bb,b->z)==0);
1502    mpz_clear(bb);
1503    return bo;
1504  }
1505// long - long
1506  if (((a->s==1) && (b->s==3))
1507  ||  ((b->s==1) && (a->s==3)))
1508    return FALSE;
1509  if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
1510    return FALSE;
1511  if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
1512    return FALSE;
1513  mpz_t  aa;
1514  mpz_t  bb;
1515  mpz_init_set(aa,a->z);
1516  mpz_init_set(bb,b->z);
1517  if (a->s<2) mpz_mul(bb,bb,a->n);
1518  if (b->s<2) mpz_mul(aa,aa,b->n);
1519  bo=(mpz_cmp(aa,bb)==0);
1520  mpz_clear(aa);
1521  mpz_clear(bb);
1522  return bo;
1523}
1524
1525// copy not immediate number a
1526number _nlCopy_NoImm(number a)
1527{
1528  assume(!((SR_HDL(a) & SR_INT)||(a==NULL)));
1529#ifdef LDEBUG
1530  nlTest(a);
1531#endif
1532  number b=(number)omAllocBin(rnumber_bin);
1533#if defined(LDEBUG)
1534  b->debug=123456;
1535#endif
1536  switch (a->s)
1537  {
1538    case 0:
1539    case 1:
1540            mpz_init_set(b->n,a->n);
1541    case 3:
1542            mpz_init_set(b->z,a->z);
1543            break;
1544  }
1545  b->s = a->s;
1546#ifdef LDEBUG
1547  nlTest(b);
1548#endif
1549  return b;
1550}
1551
1552void _nlDelete_NoImm(number *a)
1553{
1554  {
1555    switch ((*a)->s)
1556    {
1557      case 0:
1558      case 1:
1559        mpz_clear((*a)->n);
1560      case 3:
1561        mpz_clear((*a)->z);
1562#ifdef LDEBUG
1563        (*a)->s=2;
1564#endif
1565    }
1566    omFreeBin((ADDRESS) *a, rnumber_bin);
1567  }
1568}
1569
1570number _nlNeg_NoImm(number a)
1571{
1572  {
1573    mpz_neg(a->z,a->z);
1574    if ((a->s==3) && (mpz_size1(a->z)<=MP_SMALL))
1575    {
1576      int ui=(int)mpz_get_si(a->z);
1577      if ((((ui<<3)>>3)==ui)
1578        && (mpz_cmp_si(a->z,(long)ui)==0))
1579      {
1580        mpz_clear(a->z);
1581        omFreeBin((ADDRESS)a, rnumber_bin);
1582        a=INT_TO_SR(ui);
1583      }
1584    }
1585  }
1586#ifdef LDEBUG
1587  nlTest(a);
1588#endif
1589  return a;
1590}
1591
1592number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1593{
1594  number u=(number)omAllocBin(rnumber_bin);
1595#if defined(LDEBUG)
1596  u->debug=123456;
1597#endif
1598  mpz_init(u->z);
1599  if (SR_HDL(b) & SR_INT)
1600  {
1601    number x=a;
1602    a=b;
1603    b=x;
1604  }
1605  if (SR_HDL(a) & SR_INT)
1606  {
1607    switch (b->s)
1608    {
1609      case 0:
1610      case 1:/* a:short, b:1 */
1611      {
1612        mpz_t x;
1613        mpz_init(x);
1614        mpz_mul_si(x,b->n,SR_TO_INT(a));
1615        mpz_add(u->z,b->z,x);
1616        mpz_clear(x);
1617        if (mpz_cmp_ui(u->z,(long)0)==0)
1618        {
1619          mpz_clear(u->z);
1620          omFreeBin((ADDRESS)u, rnumber_bin);
1621          return INT_TO_SR(0);
1622        }
1623        if (mpz_cmp(u->z,b->n)==0)
1624        {
1625          mpz_clear(u->z);
1626          omFreeBin((ADDRESS)u, rnumber_bin);
1627          return INT_TO_SR(1);
1628        }
1629        mpz_init_set(u->n,b->n);
1630        u->s = 0;
1631        break;
1632      }
1633      case 3:
1634      {
1635        if ((long)a>0L)
1636          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1637        else
1638          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1639        if (mpz_cmp_ui(u->z,(long)0)==0)
1640        {
1641          mpz_clear(u->z);
1642          omFreeBin((ADDRESS)u, rnumber_bin);
1643          return INT_TO_SR(0);
1644        }
1645        //u->s = 3;
1646        if (mpz_size1(u->z)<=MP_SMALL)
1647        {
1648          int ui=(int)mpz_get_si(u->z);
1649          if ((((ui<<3)>>3)==ui)
1650          && (mpz_cmp_si(u->z,(long)ui)==0))
1651          {
1652            mpz_clear(u->z);
1653            omFreeBin((ADDRESS)u, rnumber_bin);
1654            return INT_TO_SR(ui);
1655          }
1656        }
1657        u->s = 3;
1658        break;
1659      }
1660    }
1661  }
1662  else
1663  {
1664    switch (a->s)
1665    {
1666      case 0:
1667      case 1:
1668      {
1669        switch(b->s)
1670        {
1671          case 0:
1672          case 1:
1673          {
1674            mpz_t x;
1675            mpz_init(x);
1676
1677            mpz_mul(x,b->z,a->n);
1678            mpz_mul(u->z,a->z,b->n);
1679            mpz_add(u->z,u->z,x);
1680            mpz_clear(x);
1681
1682            if (mpz_cmp_ui(u->z,(long)0)==0)
1683            {
1684              mpz_clear(u->z);
1685              omFreeBin((ADDRESS)u, rnumber_bin);
1686              return INT_TO_SR(0);
1687            }
1688            mpz_init(u->n);
1689            mpz_mul(u->n,a->n,b->n);
1690            if (mpz_cmp(u->z,u->n)==0)
1691            {
1692               mpz_clear(u->z);
1693               mpz_clear(u->n);
1694               omFreeBin((ADDRESS)u, rnumber_bin);
1695               return INT_TO_SR(1);
1696            }
1697            u->s = 0;
1698            break;
1699          }
1700          case 3: /* a:1 b:3 */
1701          {
1702            mpz_mul(u->z,b->z,a->n);
1703            mpz_add(u->z,u->z,a->z);
1704            if (mpz_cmp_ui(u->z,(long)0)==0)
1705            {
1706              mpz_clear(u->z);
1707              omFreeBin((ADDRESS)u, rnumber_bin);
1708              return INT_TO_SR(0);
1709            }
1710            if (mpz_cmp(u->z,a->n)==0)
1711            {
1712              mpz_clear(u->z);
1713              omFreeBin((ADDRESS)u, rnumber_bin);
1714              return INT_TO_SR(1);
1715            }
1716            mpz_init_set(u->n,a->n);
1717            u->s = 0;
1718            break;
1719          }
1720        } /*switch (b->s) */
1721        break;
1722      }
1723      case 3:
1724      {
1725        switch(b->s)
1726        {
1727          case 0:
1728          case 1:/* a:3, b:1 */
1729          {
1730            mpz_mul(u->z,a->z,b->n);
1731            mpz_add(u->z,u->z,b->z);
1732            if (mpz_cmp_ui(u->z,(long)0)==0)
1733            {
1734              mpz_clear(u->z);
1735              omFreeBin((ADDRESS)u, rnumber_bin);
1736              return INT_TO_SR(0);
1737            }
1738            if (mpz_cmp(u->z,b->n)==0)
1739            {
1740              mpz_clear(u->z);
1741              omFreeBin((ADDRESS)u, rnumber_bin);
1742              return INT_TO_SR(1);
1743            }
1744            mpz_init_set(u->n,b->n);
1745            u->s = 0;
1746            break;
1747          }
1748          case 3:
1749          {
1750            mpz_add(u->z,a->z,b->z);
1751            if (mpz_cmp_ui(u->z,(long)0)==0)
1752            {
1753              mpz_clear(u->z);
1754              omFreeBin((ADDRESS)u, rnumber_bin);
1755              return INT_TO_SR(0);
1756            }
1757            if (mpz_size1(u->z)<=MP_SMALL)
1758            {
1759              int ui=(int)mpz_get_si(u->z);
1760              if ((((ui<<3)>>3)==ui)
1761              && (mpz_cmp_si(u->z,(long)ui)==0))
1762              {
1763                mpz_clear(u->z);
1764                omFreeBin((ADDRESS)u, rnumber_bin);
1765                return INT_TO_SR(ui);
1766              }
1767            }
1768            u->s = 3;
1769            break;
1770          }
1771        }
1772        break;
1773      }
1774    }
1775  }
1776#ifdef LDEBUG
1777  nlTest(u);
1778#endif
1779  return u;
1780}
1781
1782number _nlSub_aNoImm_OR_bNoImm(number a, number b)
1783{
1784  number u=(number)omAllocBin(rnumber_bin);
1785#if defined(LDEBUG)
1786  u->debug=123456;
1787#endif
1788  mpz_init(u->z);
1789  if (SR_HDL(a) & SR_INT)
1790  {
1791    switch (b->s)
1792    {
1793      case 0:
1794      case 1:/* a:short, b:1 */
1795      {
1796        mpz_t x;
1797        mpz_init(x);
1798        mpz_mul_si(x,b->n,SR_TO_INT(a));
1799        mpz_sub(u->z,x,b->z);
1800        mpz_clear(x);
1801        if (mpz_cmp_ui(u->z,(long)0)==0)
1802        {
1803          mpz_clear(u->z);
1804          omFreeBin((ADDRESS)u, rnumber_bin);
1805          return INT_TO_SR(0);
1806        }
1807        if (mpz_cmp(u->z,b->n)==0)
1808        {
1809          mpz_clear(u->z);
1810          omFreeBin((ADDRESS)u, rnumber_bin);
1811          return INT_TO_SR(1);
1812        }
1813        mpz_init_set(u->n,b->n);
1814        u->s = 0;
1815        break;
1816      }
1817      case 3:
1818      {
1819        if ((long)a>0L)
1820        {
1821          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
1822          mpz_neg(u->z,u->z);
1823        }
1824        else
1825        {
1826          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
1827          mpz_neg(u->z,u->z);
1828        }
1829        if (mpz_cmp_ui(u->z,(long)0)==0)
1830        {
1831          mpz_clear(u->z);
1832          omFreeBin((ADDRESS)u, rnumber_bin);
1833          return INT_TO_SR(0);
1834        }
1835        if (mpz_size1(u->z)<=MP_SMALL)
1836        {
1837          int ui=(int)mpz_get_si(u->z);
1838          if ((((ui<<3)>>3)==ui)
1839          && (mpz_cmp_si(u->z,(long)ui)==0))
1840          {
1841            mpz_clear(u->z);
1842            omFreeBin((ADDRESS)u, rnumber_bin);
1843            return INT_TO_SR(ui);
1844          }
1845        }
1846        u->s = 3;
1847        break;
1848      }
1849    }
1850  }
1851  else if (SR_HDL(b) & SR_INT)
1852  {
1853    switch (a->s)
1854    {
1855      case 0:
1856      case 1:/* b:short, a:1 */
1857      {
1858        mpz_t x;
1859        mpz_init(x);
1860        mpz_mul_si(x,a->n,SR_TO_INT(b));
1861        mpz_sub(u->z,a->z,x);
1862        mpz_clear(x);
1863        if (mpz_cmp_ui(u->z,(long)0)==0)
1864        {
1865          mpz_clear(u->z);
1866          omFreeBin((ADDRESS)u, rnumber_bin);
1867          return INT_TO_SR(0);
1868        }
1869        if (mpz_cmp(u->z,a->n)==0)
1870        {
1871          mpz_clear(u->z);
1872          omFreeBin((ADDRESS)u, rnumber_bin);
1873          return INT_TO_SR(1);
1874        }
1875        mpz_init_set(u->n,a->n);
1876        u->s = 0;
1877        break;
1878      }
1879      case 3:
1880      {
1881        if ((long)b>0L)
1882        {
1883          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
1884        }
1885        else
1886        {
1887          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
1888        }
1889        if (mpz_cmp_ui(u->z,(long)0)==0)
1890        {
1891          mpz_clear(u->z);
1892          omFreeBin((ADDRESS)u, rnumber_bin);
1893          return INT_TO_SR(0);
1894        }
1895        //u->s = 3;
1896        if (mpz_size1(u->z)<=MP_SMALL)
1897        {
1898          int ui=(int)mpz_get_si(u->z);
1899          if ((((ui<<3)>>3)==ui)
1900          && (mpz_cmp_si(u->z,(long)ui)==0))
1901          {
1902            mpz_clear(u->z);
1903            omFreeBin((ADDRESS)u, rnumber_bin);
1904            return INT_TO_SR(ui);
1905          }
1906        }
1907        u->s = 3;
1908        break;
1909      }
1910    }
1911  }
1912  else
1913  {
1914    switch (a->s)
1915    {
1916      case 0:
1917      case 1:
1918      {
1919        switch(b->s)
1920        {
1921          case 0:
1922          case 1:
1923          {
1924            mpz_t x;
1925            mpz_t y;
1926            mpz_init(x);
1927            mpz_init(y);
1928            mpz_mul(x,b->z,a->n);
1929            mpz_mul(y,a->z,b->n);
1930            mpz_sub(u->z,y,x);
1931            mpz_clear(x);
1932            mpz_clear(y);
1933            if (mpz_cmp_ui(u->z,(long)0)==0)
1934            {
1935              mpz_clear(u->z);
1936              omFreeBin((ADDRESS)u, rnumber_bin);
1937              return INT_TO_SR(0);
1938            }
1939            mpz_init(u->n);
1940            mpz_mul(u->n,a->n,b->n);
1941            if (mpz_cmp(u->z,u->n)==0)
1942            {
1943              mpz_clear(u->z);
1944              mpz_clear(u->n);
1945              omFreeBin((ADDRESS)u, rnumber_bin);
1946              return INT_TO_SR(1);
1947            }
1948            u->s = 0;
1949            break;
1950          }
1951          case 3: /* a:1, b:3 */
1952          {
1953            mpz_t x;
1954            mpz_init(x);
1955            mpz_mul(x,b->z,a->n);
1956            mpz_sub(u->z,a->z,x);
1957            mpz_clear(x);
1958            if (mpz_cmp_ui(u->z,(long)0)==0)
1959            {
1960              mpz_clear(u->z);
1961              omFreeBin((ADDRESS)u, rnumber_bin);
1962              return INT_TO_SR(0);
1963            }
1964            if (mpz_cmp(u->z,a->n)==0)
1965            {
1966              mpz_clear(u->z);
1967              omFreeBin((ADDRESS)u, rnumber_bin);
1968              return INT_TO_SR(1);
1969            }
1970            mpz_init_set(u->n,a->n);
1971            u->s = 0;
1972            break;
1973          }
1974        }
1975        break;
1976      }
1977      case 3:
1978      {
1979        switch(b->s)
1980        {
1981          case 0:
1982          case 1: /* a:3, b:1 */
1983          {
1984            mpz_t x;
1985            mpz_init(x);
1986            mpz_mul(x,a->z,b->n);
1987            mpz_sub(u->z,x,b->z);
1988            mpz_clear(x);
1989            if (mpz_cmp_ui(u->z,(long)0)==0)
1990            {
1991              mpz_clear(u->z);
1992              omFreeBin((ADDRESS)u, rnumber_bin);
1993              return INT_TO_SR(0);
1994            }
1995            if (mpz_cmp(u->z,b->n)==0)
1996            {
1997              mpz_clear(u->z);
1998              omFreeBin((ADDRESS)u, rnumber_bin);
1999              return INT_TO_SR(1);
2000            }
2001            mpz_init_set(u->n,b->n);
2002            u->s = 0;
2003            break;
2004          }
2005          case 3: /* a:3 , b:3 */
2006          {
2007            mpz_sub(u->z,a->z,b->z);
2008            if (mpz_cmp_ui(u->z,(long)0)==0)
2009            {
2010              mpz_clear(u->z);
2011              omFreeBin((ADDRESS)u, rnumber_bin);
2012              return INT_TO_SR(0);
2013            }
2014            //u->s = 3;
2015            if (mpz_size1(u->z)<=MP_SMALL)
2016            {
2017              int ui=(int)mpz_get_si(u->z);
2018              if ((((ui<<3)>>3)==ui)
2019              && (mpz_cmp_si(u->z,(long)ui)==0))
2020              {
2021                mpz_clear(u->z);
2022                omFreeBin((ADDRESS)u, rnumber_bin);
2023                return INT_TO_SR(ui);
2024              }
2025            }
2026            u->s = 3;
2027            break;
2028          }
2029        }
2030        break;
2031      }
2032    }
2033  }
2034#ifdef LDEBUG
2035  nlTest(u);
2036#endif
2037  return u;
2038}
2039
2040// a and b are intermediate, but a*b not
2041number _nlMult_aImm_bImm_rNoImm(number a, number b)
2042{
2043  number u=(number)omAllocBin(rnumber_bin);
2044#if defined(LDEBUG)
2045  u->debug=123456;
2046#endif
2047  u->s=3;
2048  mpz_init_set_si(u->z,(long)SR_TO_INT(a));
2049  mpz_mul_si(u->z,u->z,(long)SR_TO_INT(b));
2050#ifdef LDEBUG
2051  nlTest(u);
2052#endif
2053  return u;
2054}
2055
2056// a or b are not immediate
2057number _nlMult_aNoImm_OR_bNoImm(number a, number b)
2058{
2059  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2060  number u=(number)omAllocBin(rnumber_bin);
2061#if defined(LDEBUG)
2062  u->debug=123456;
2063#endif
2064  mpz_init(u->z);
2065  if (SR_HDL(b) & SR_INT)
2066  {
2067    number x=a;
2068    a=b;
2069    b=x;
2070  }
2071  if (SR_HDL(a) & SR_INT)
2072  {
2073    u->s=b->s;
2074    if (u->s==1) u->s=0;
2075    if ((long)a>0L)
2076    {
2077      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2078    }
2079    else
2080    {
2081      if (a==INT_TO_SR(-1))
2082      {
2083        mpz_set(u->z,b->z);
2084        mpz_neg(u->z,u->z);
2085        u->s=b->s;
2086      }
2087      else
2088      {
2089        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2090        mpz_neg(u->z,u->z);
2091      }
2092    }
2093    if (u->s<2)
2094    {
2095      if (mpz_cmp(u->z,b->n)==0)
2096      {
2097        mpz_clear(u->z);
2098        omFreeBin((ADDRESS)u, rnumber_bin);
2099        return INT_TO_SR(1);
2100      }
2101      mpz_init_set(u->n,b->n);
2102    }
2103    else //u->s==3
2104    {
2105      if (mpz_size1(u->z)<=MP_SMALL)
2106      {
2107        int ui=(int)mpz_get_si(u->z);
2108        if ((((ui<<3)>>3)==ui)
2109            && (mpz_cmp_si(u->z,(long)ui)==0))
2110        {
2111          mpz_clear(u->z);
2112          omFreeBin((ADDRESS)u, rnumber_bin);
2113          return INT_TO_SR(ui);
2114        }
2115      }
2116    }
2117  }
2118  else
2119  {
2120    mpz_mul(u->z,a->z,b->z);
2121    u->s = 0;
2122    if(a->s==3)
2123    {
2124      if(b->s==3)
2125      {
2126        u->s = 3;
2127      }
2128      else
2129      {
2130        if (mpz_cmp(u->z,b->n)==0)
2131        {
2132          mpz_clear(u->z);
2133          omFreeBin((ADDRESS)u, rnumber_bin);
2134          return INT_TO_SR(1);
2135        }
2136        mpz_init_set(u->n,b->n);
2137      }
2138    }
2139    else
2140    {
2141      if(b->s==3)
2142      {
2143        if (mpz_cmp(u->z,a->n)==0)
2144        {
2145          mpz_clear(u->z);
2146          omFreeBin((ADDRESS)u, rnumber_bin);
2147          return INT_TO_SR(1);
2148        }
2149        mpz_init_set(u->n,a->n);
2150      }
2151      else
2152      {
2153        mpz_init(u->n);
2154        mpz_mul(u->n,a->n,b->n);
2155        if (mpz_cmp(u->z,u->n)==0)
2156        {
2157          mpz_clear(u->z);
2158          mpz_clear(u->n);
2159          omFreeBin((ADDRESS)u, rnumber_bin);
2160          return INT_TO_SR(1);
2161        }
2162      }
2163    }
2164  }
2165#ifdef LDEBUG
2166  nlTest(u);
2167#endif
2168  return u;
2169}
2170
2171/*2
2172* z := i
2173*/
2174number nlRInit (int i)
2175{
2176  number z=(number)omAllocBin(rnumber_bin);
2177#if defined(LDEBUG)
2178  z->debug=123456;
2179#endif
2180  mpz_init_set_si(z->z,(long)i);
2181  z->s = 3;
2182  return z;
2183}
2184
2185/*2
2186* z := i/j
2187*/
2188number nlInit2 (int i, int j)
2189{
2190  number z=(number)omAllocBin(rnumber_bin);
2191#if defined(LDEBUG)
2192  z->debug=123456;
2193#endif
2194  mpz_init_set_si(z->z,(long)i);
2195  mpz_init_set_si(z->n,(long)j);
2196  z->s = 0;
2197  nlNormalize(z);
2198  return z;
2199}
2200
2201number nlInit2gmp (mpz_t i, mpz_t j)
2202{
2203  number z=(number)omAllocBin(rnumber_bin);
2204#if defined(LDEBUG)
2205  z->debug=123456;
2206#endif
2207  mpz_init_set(z->z,i);
2208  mpz_init_set(z->n,j);
2209  z->s = 0;
2210  nlNormalize(z);
2211  return z;
2212}
2213
2214#else // DO_LINLINE
2215
2216// declare immedate routines
2217number nlRInit (int i);
2218BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2219number  _nlCopy_NoImm(number a);
2220number  _nlNeg_NoImm(number a);
2221number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2222number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2223number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2224number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2225
2226#endif
2227
2228/***************************************************************
2229 *
2230 * Routines which might be inlined by p_Numbers.h
2231 *
2232 *******************************************************************/
2233#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2234
2235// routines which are always inlined/static
2236
2237/*2
2238* a = b ?
2239*/
2240LINLINE BOOLEAN nlEqual (number a, number b)
2241{
2242#ifdef LDEBUG
2243  nlTest(a);
2244  nlTest(b);
2245#endif
2246// short - short
2247  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2248  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2249}
2250
2251
2252LINLINE number nlInit (int i, const ring r)
2253{
2254  number n;
2255  if ( ((i << 3) >> 3) == i ) n=INT_TO_SR(i);
2256  else                        n=nlRInit(i);
2257#ifdef LDEBUG
2258  nlTest(n);
2259#endif
2260  return n;
2261}
2262
2263
2264/*2
2265* a == 1 ?
2266*/
2267LINLINE BOOLEAN nlIsOne (number a)
2268{
2269#ifdef LDEBUG
2270  if (a==NULL) return FALSE;
2271  nlTest(a);
2272#endif
2273  return (a==INT_TO_SR(1));
2274}
2275
2276LINLINE BOOLEAN nlIsZero (number a)
2277{
2278  return (a==INT_TO_SR(0));
2279  //return (mpz_cmp_si(a->z,(long)0)==0);
2280}
2281
2282/*2
2283* copy a to b
2284*/
2285LINLINE number nlCopy(number a)
2286{
2287  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2288  {
2289    return a;
2290  }
2291  return _nlCopy_NoImm(a);
2292}
2293
2294
2295/*2
2296* delete a
2297*/
2298LINLINE void nlDelete (number * a, const ring r)
2299{
2300  if (*a!=NULL)
2301  {
2302#ifdef LDEBUG
2303    nlTest(*a);
2304#endif
2305    if ((SR_HDL(*a) & SR_INT)==0)
2306    {
2307      _nlDelete_NoImm(a);
2308    }
2309    *a=NULL;
2310  }
2311}
2312
2313/*2
2314* za:= - za
2315*/
2316LINLINE number nlNeg (number a)
2317{
2318#ifdef LDEBUG
2319  nlTest(a);
2320#endif
2321  if(SR_HDL(a) &SR_INT)
2322  {
2323    int r=SR_TO_INT(a);
2324    if (r==(-(1<<28))) a=nlRInit(1<<28);
2325    else               a=INT_TO_SR(-r);
2326    return a;
2327  }
2328  return _nlNeg_NoImm(a);
2329}
2330
2331/*2
2332* u:= a + b
2333*/
2334LINLINE number nlAdd (number a, number b)
2335{
2336  number u;
2337  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2338  {
2339    int r=SR_HDL(a)+SR_HDL(b)-1;
2340    if ( ((r << 1) >> 1) == r )
2341      return (number)(long)r;
2342    else
2343      return nlRInit(SR_TO_INT(r));
2344  }
2345  return _nlAdd_aNoImm_OR_bNoImm(a, b);
2346}
2347
2348number nlShort1(number a);
2349number nlShort3(number a);
2350
2351LINLINE number nlInpAdd (number a, number b, const ring r)
2352{
2353  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2354  {
2355    int r=SR_HDL(a)+SR_HDL(b)-1;
2356    if ( ((r << 1) >> 1) == r )
2357      return (number)(long)r;
2358    else
2359      return nlRInit(SR_TO_INT(r));
2360  }
2361  // a=a+b
2362  if (SR_HDL(b) & SR_INT)
2363  {
2364    switch (a->s)
2365    {
2366      case 0:
2367      case 1:/* b:short, a:1 */
2368      {
2369        mpz_t x;
2370        mpz_init(x);
2371        mpz_mul_si(x,a->n,SR_TO_INT(b));
2372        mpz_add(a->z,a->z,x);
2373        mpz_clear(x);
2374        a->s = 0;
2375        a=nlShort1(a);
2376        break;
2377      }
2378      case 3:
2379      {
2380        if ((long)b>0L)
2381          mpz_add_ui(a->z,a->z,SR_TO_INT(b));
2382        else
2383          mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
2384        a->s = 3;
2385        a=nlShort3(a);
2386        //nlNormalize(a);
2387        break;
2388      }
2389    }
2390    return a;
2391  }
2392  if (SR_HDL(a) & SR_INT)
2393  {
2394    number u=(number)omAllocBin(rnumber_bin);
2395    #if defined(LDEBUG)
2396    u->debug=123456;
2397    #endif
2398    mpz_init(u->z);
2399    switch (b->s)
2400    {
2401      case 0:
2402      case 1:/* a:short, b:1 */
2403      {
2404        mpz_t x;
2405        mpz_init(x);
2406
2407        mpz_mul_si(x,b->n,SR_TO_INT(a));
2408        mpz_add(u->z,b->z,x);
2409        mpz_clear(x);
2410        // result cannot be 0, if coeffs are normalized
2411        mpz_init_set(u->n,b->n);
2412        u->s = 0;
2413        u=nlShort1(u);
2414        break;
2415      }
2416      case 3:
2417      {
2418        if ((long)a>0L)
2419          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
2420        else
2421          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
2422        // result cannot be 0, if coeffs are normalized
2423        u->s = 3;
2424        u=nlShort3(u);
2425        break;
2426      }
2427    }
2428    #ifdef LDEBUG
2429    nlTest(u);
2430    #endif
2431    return u;
2432  }
2433  else
2434  {
2435    switch (a->s)
2436    {
2437      case 0:
2438      case 1:
2439      {
2440        switch(b->s)
2441        {
2442          case 0:
2443          case 1: /* a:1 b:1 */
2444          {
2445            mpz_t x;
2446            mpz_t y;
2447            mpz_init(x);
2448            mpz_init(y);
2449            mpz_mul(x,b->z,a->n);
2450            mpz_mul(y,a->z,b->n);
2451            mpz_add(a->z,x,y);
2452            mpz_clear(x);
2453            mpz_clear(y);
2454            mpz_mul(a->n,a->n,b->n);
2455            a->s = 0;
2456            break;
2457          }
2458          case 3: /* a:1 b:3 */
2459          {
2460            mpz_t x;
2461            mpz_init(x);
2462            mpz_mul(x,b->z,a->n);
2463            mpz_add(a->z,a->z,x);
2464            mpz_clear(x);
2465            a->s = 0;
2466            break;
2467          }
2468        } /*switch (b->s) */
2469        a=nlShort1(a);
2470        break;
2471      }
2472      case 3:
2473      {
2474        switch(b->s)
2475        {
2476          case 0:
2477          case 1:/* a:3, b:1 */
2478          {
2479            mpz_t x;
2480            mpz_init(x);
2481            mpz_mul(x,a->z,b->n);
2482            mpz_add(a->z,b->z,x);
2483            mpz_clear(x);
2484            mpz_init_set(a->n,b->n);
2485            a->s = 0;
2486            a=nlShort1(a);
2487            break;
2488          }
2489          case 3:
2490          {
2491            mpz_add(a->z,a->z,b->z);
2492            a->s = 3;
2493            a=nlShort3(a);
2494            break;
2495          }
2496        }
2497        break;
2498      }
2499    }
2500    #ifdef LDEBUG
2501    nlTest(a);
2502    #endif
2503    return a;
2504  }
2505}
2506
2507LINLINE number nlMult (number a, number b)
2508{
2509#ifdef LDEBUG
2510  nlTest(a);
2511  nlTest(b);
2512#endif
2513  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2514  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2515  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2516  {
2517    int r=(SR_HDL(a)-1)*(SR_HDL(b)>>1);
2518    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1))
2519    {
2520      number u=((number) ((r>>1)+SR_INT));
2521      if (((SR_HDL(u)<<1)>>1)==SR_HDL(u)) return (u);
2522      return nlRInit(SR_HDL(u)>>2);
2523    }
2524    return _nlMult_aImm_bImm_rNoImm(a, b);
2525  }
2526  return _nlMult_aNoImm_OR_bNoImm(a, b);
2527}
2528
2529
2530/*2
2531* u:= a - b
2532*/
2533LINLINE number nlSub (number a, number b)
2534{
2535  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2536  {
2537    int r=SR_HDL(a)-SR_HDL(b)+1;
2538    if ( ((r << 1) >> 1) == r )
2539    {
2540      return (number)r;
2541    }
2542    else
2543      return nlRInit(SR_TO_INT(r));
2544  }
2545  return _nlSub_aNoImm_OR_bNoImm(a, b);
2546}
2547
2548#endif // DO_LINLINE
2549
2550#ifndef P_NUMBERS_H
2551
2552void nlInpGcd(number &a, number b, const ring r)
2553{
2554  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2555  {
2556    number n=nlGcd(a,b,r);
2557    nlDelete(&a,r);
2558    a=n;
2559  }
2560  else
2561  {
2562    mpz_gcd(a->z,a->z,b->z);
2563    if (mpz_size1(a->z)<=MP_SMALL)
2564    {
2565      int ui=(int)mpz_get_si(a->z);
2566      if ((((ui<<3)>>3)==ui)
2567      && (mpz_cmp_si(a->z,(long)ui)==0))
2568      {
2569        mpz_clear(a->z);
2570        omFreeBin((ADDRESS)a, rnumber_bin);
2571        a=INT_TO_SR(ui);
2572      }
2573    }
2574  }
2575}
2576void nlInpIntDiv(number &a, number b, const ring r)
2577{
2578  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2579  {
2580    number n=nlIntDiv(a,b);
2581    nlDelete(&a,r);
2582    a=n;
2583  }
2584  else
2585  {
2586    if (mpz_isNeg(a->z))
2587    {
2588      if (mpz_isNeg(b->z))
2589      {
2590        mpz_add(a->z,a->z,b->z);
2591      }
2592      else
2593      {
2594        mpz_sub(a->z,a->z,b->z);
2595      }
2596      mpz_add_ui(a->z,a->z,1);
2597    }
2598    else
2599    {
2600      if (mpz_isNeg(b->z))
2601      {
2602        mpz_sub(a->z,a->z,b->z);
2603      }
2604      else
2605      {
2606        mpz_add(a->z,a->z,b->z);
2607      }
2608      mpz_sub_ui(a->z,a->z,1);
2609    }
2610    MPZ_DIV(a->z,a->z,b->z);
2611    if (mpz_size1(a->z)<=MP_SMALL)
2612    {
2613      int ui=(int)mpz_get_si(a->z);
2614      if ((((ui<<3)>>3)==ui)
2615      && (mpz_cmp_si(a->z,(long)ui)==0))
2616      {
2617        mpz_clear(a->z);
2618        omFreeBin((ADDRESS)a, rnumber_bin);
2619        a=INT_TO_SR(ui);
2620      }
2621    }
2622  }
2623}
2624void nlInpMult(number &a, number b, const ring r)
2625{
2626  if (((SR_HDL(b)|SR_HDL(a))&SR_INT)
2627  )
2628  {
2629    number n=nlMult(a,b);
2630    nlDelete(&a,r);
2631    a=n;
2632  }
2633  else
2634  {
2635    mpz_mul(a->z,a->z,b->z);
2636    if (a->s==3)
2637    {
2638      if(b->s!=3)
2639      {
2640        mpz_init_set(a->n,b->n);
2641        a->s=0;
2642      }
2643    }
2644    else
2645    {
2646      if(b->s!=3)
2647      {
2648        mpz_mul(a->n,a->n,b->n);
2649      }
2650      a->s=0;
2651    }
2652  }
2653}
2654
2655number nlFarey(number nN, number nP)
2656{
2657  mpz_t tmp; mpz_init(tmp);
2658  mpz_t A,B,C,D,E,N,P;
2659  if (SR_HDL(nN) & SR_INT) mpz_init_set_si(N,SR_TO_INT(nN));
2660  else                     mpz_init_set(N,nN->z);
2661  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2662  else                     mpz_init_set(P,nP->z);
2663  assume(!mpz_isNeg(P));
2664  if (mpz_isNeg(N))  mpz_add(N,N,P);
2665  mpz_init_set_si(A,(long)0);
2666  mpz_init_set_ui(B,(unsigned long)1);
2667  mpz_init_set_si(C,(long)0);
2668  mpz_init(D);
2669  mpz_init_set(E,P);
2670  number z=INT_TO_SR(0);
2671  while(mpz_cmp_si(N,(long)0)!=0)
2672  {
2673    mpz_mul(tmp,N,N);
2674    mpz_add(tmp,tmp,tmp);
2675    if (mpz_cmp(tmp,P)<0)
2676    {
2677       // return N/B
2678       z=(number)omAllocBin(rnumber_bin);
2679       if (mpz_isNeg(B))
2680       {
2681         mpz_neg(B,B);
2682         mpz_neg(N,N);
2683       }
2684       mpz_init_set(z->z,N);
2685       mpz_init_set(z->n,B);
2686       z->s = 0;
2687       nlNormalize(z);
2688       break;
2689    }
2690    mpz_mod(D,E,N);
2691    mpz_div(tmp,E,N);
2692    mpz_mul(tmp,tmp,B);
2693    mpz_sub(C,A,tmp);
2694    mpz_set(E,N);
2695    mpz_set(N,D);
2696    mpz_set(A,B);
2697    mpz_set(B,C);
2698  }
2699  mpz_clear(tmp);
2700  mpz_clear(A);
2701  mpz_clear(B);
2702  mpz_clear(C);
2703  mpz_clear(D);
2704  mpz_clear(E);
2705  mpz_clear(N);
2706  mpz_clear(P);
2707  return z;
2708}
2709#if 0
2710number nlMod(number a, number b)
2711{
2712  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
2713  {
2714    int bi=SR_TO_INT(b);
2715    int ai=SR_TO_INT(a);
2716    int bb=ABS(bi);
2717    int c=ai%bb;
2718    if (c<0)  c+=bb;
2719    return (INT_TO_SR(c));
2720  }
2721  number al;
2722  number bl;
2723  if (SR_HDL(a))&SR_INT)
2724    al=nlRInit(SR_TO_INT(a));
2725  else
2726    al=nlCopy(a);
2727  if (SR_HDL(b))&SR_INT)
2728    bl=nlRInit(SR_TO_INT(b));
2729  else
2730    bl=nlCopy(b);
2731  number r=nlRInit(0);
2732  mpz_mod(r->z,al->z,bl->z);
2733  nlDelete(&al);
2734  nlDelete(&bl);
2735  if (mpz_size1(&r->z)<=MP_SMALL)
2736  {
2737    int ui=(int)mpz_get_si(&r->z);
2738    if ((((ui<<3)>>3)==ui)
2739    && (mpz_cmp_si(x->z,(long)ui)==0))
2740    {
2741      mpz_clear(&r->z);
2742      omFreeBin((ADDRESS)r, rnumber_bin);
2743      r=INT_TO_SR(ui);
2744    }
2745  }
2746  return r;
2747}
2748#endif
2749#endif // not P_NUMBERS_H
2750#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.