source: git/kernel/longrat.cc @ bc669c

fieker-DuValspielwiese
Last change on this file since bc669c was ea66a97, checked in by Hans Schönemann <hannes@…>, 14 years ago
memory leak in nlGCD fixed git-svn-id: file:///usr/local/Singular/svn/trunk@12585 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 56.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
7*/
8
9#ifndef LONGRAT_CC
10#define LONGRAT_CC
11
12#include <string.h>
13#include <float.h>
14#include "mod2.h"
15#include "structs.h"
16#include "omalloc.h"
17#include "febase.h"
18#include "numbers.h"
19#include "modulop.h"
20#include "ring.h"
21#include "shortfl.h"
22#include "mpr_complex.h"
23#include "longrat.h"
24
25
26#ifndef BYTES_PER_MP_LIMB
27#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
28#endif
29
30/*-----------------------------------------------------------------*/
31/*3
32* parameter s in number:
33* 0 (or FALSE): not normalised rational
34* 1 (or TRUE):  normalised rational
35* 3          :  integer with n==NULL
36*/
37/*3
38**  'SR_INT' is the type of those integers small enough to fit into  29  bits.
39**  Therefor the value range of this small integers is: $-2^{28}...2^{28}-1$.
40**
41**  Small integers are represented by an immediate integer handle, containing
42**  the value instead of pointing  to  it,  which  has  the  following  form:
43**
44**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
45**      | guard | sign  | bit   | bit   |       | bit   | tag   | tag   |
46**      | bit   | bit   | 27    | 26    |       | 0     | 0     | 1     |
47**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
48**
49**  Immediate integers handles carry the tag 'SR_INT', i.e. the last bit is 1.
50**  This distuingishes immediate integers from other handles which  point  to
51**  structures aligned on 4 byte boundaries and therefor have last bit  zero.
52**  (The second bit is reserved as tag to allow extensions of  this  scheme.)
53**  Using immediates as pointers and dereferencing them gives address errors.
54**
55**  To aid overflow check the most significant two bits must always be equal,
56**  that is to say that the sign bit of immediate integers has a  guard  bit.
57**
58**  The macros 'INT_TO_SR' and 'SR_TO_INT' should be used to convert  between
59**  a small integer value and its representation as immediate integer handle.
60**
61**  Large integers and rationals are represented by z and n
62**  where n may be undefined (if s==3)
63**  NULL represents only deleted values
64*/
65#define SR_HDL(A) ((long)(A))
66/*#define SR_INT    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,(MP_INT*) 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._mp_d,a->z._mp_alloc*BYTES_PER_MP_LIMB);
238  if (a->z._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._mp_d,a->n._mp_alloc*BYTES_PER_MP_LIMB), return FALSE);
244    if (a->z._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  lint *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._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._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  lint 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}
1100
1101/*2
1102* result =gcd(a,b)
1103*/
1104number nlGcd(number a, number b, const ring r)
1105{
1106  number result;
1107#ifdef LDEBUG
1108  nlTest(a);
1109  nlTest(b);
1110#endif
1111  //nlNormalize(a);
1112  //nlNormalize(b);
1113  if ((a==INT_TO_SR(1))||(a==INT_TO_SR(-1))
1114  ||  (b==INT_TO_SR(1))||(b==INT_TO_SR(-1)))
1115    return INT_TO_SR(1);
1116  if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
1117    return nlCopy(b);
1118  if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
1119    return nlCopy(a);
1120  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1121  {
1122    int i=SR_TO_INT(a);
1123    int j=SR_TO_INT(b);
1124    if((i==0)||(j==0))
1125      return INT_TO_SR(1);
1126    int l;
1127    i=ABS(i);
1128    j=ABS(j);
1129    do
1130    {
1131      l=i%j;
1132      i=j;
1133      j=l;
1134    } while (l!=0);
1135    return INT_TO_SR(i);
1136  }
1137  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1138  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1139  if (SR_HDL(a) & SR_INT)
1140  {
1141    unsigned long t=mpz_gcd_ui(NULL,&b->z,ABS(SR_TO_INT(a)));
1142    return INT_TO_SR((int)t);
1143  }
1144  else
1145  if (SR_HDL(b) & SR_INT)
1146  {
1147    unsigned long t=mpz_gcd_ui(NULL,&a->z,ABS(SR_TO_INT(b)));
1148    return INT_TO_SR((int)t);
1149  }
1150  else
1151  {
1152    result=(number)omAllocBin(rnumber_bin);
1153    mpz_init(&result->z);
1154    mpz_gcd(&result->z,&a->z,&b->z);
1155    if (mpz_size1(&result->z)<=MP_SMALL)
1156    {
1157      int ui=(int)mpz_get_si(&result->z);
1158      if ((((ui<<3)>>3)==ui)
1159      && (mpz_cmp_si(&result->z,(long)ui)==0))
1160      {
1161        mpz_clear(&result->z);
1162        omFreeBin((ADDRESS)result, rnumber_bin);
1163        return INT_TO_SR(ui);
1164      }
1165    }
1166    result->s = 3;
1167  #ifdef LDEBUG
1168    result->debug=123456;
1169    nlTest(result);
1170  #endif
1171  }
1172  return result;
1173}
1174
1175number nlShort1(number x) // assume x->s==0/1
1176{
1177  assume(x->s<2);
1178  if (mpz_cmp_ui(&x->z,(long)0)==0)
1179  {
1180    _nlDelete_NoImm(&x);
1181    return INT_TO_SR(0);
1182  }
1183  if (x->s<2)
1184  {
1185    if (mpz_cmp(&x->z,&x->n)==0)
1186    {
1187      _nlDelete_NoImm(&x);
1188      return INT_TO_SR(1);
1189    }
1190  }
1191  return x;
1192}
1193number nlShort3(number x) // assume x->s==3
1194{
1195  assume(x->s==3);
1196  if ((mpz_cmp_ui(&x->z,(long)0)==0)
1197  || (mpz_size1(&x->z)<=MP_SMALL))
1198  {
1199    int ui=(int)mpz_get_si(&x->z);
1200    if ((((ui<<3)>>3)==ui)
1201    && (mpz_cmp_si(&x->z,(long)ui)==0))
1202    {
1203      mpz_clear(&x->z);
1204      omFreeBin((ADDRESS)x, rnumber_bin);
1205      return INT_TO_SR(ui);
1206    }
1207  }
1208  return x;
1209}
1210/*2
1211* simplify x
1212*/
1213void nlNormalize (number &x)
1214{
1215  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1216    return;
1217#ifdef LDEBUG
1218  if (!nlTest(x)) { if (x->s!=3) x->s=1; return; }
1219#endif
1220  if (x->s==3)
1221  {
1222    x=nlShort3(x);
1223    return;
1224  }
1225  else if (x->s==0)
1226  {
1227    if (mpz_cmp_si(&x->n,(long)1)==0)
1228    {
1229      mpz_clear(&x->n);
1230      if (mpz_size1(&x->z)<=MP_SMALL)
1231      {
1232        int ui=(int)mpz_get_si(&x->z);
1233        if ((((ui<<3)>>3)==ui)
1234        && (mpz_cmp_si(&x->z,(long)ui)==0))
1235        {
1236          mpz_clear(&x->z);
1237#if defined(LDEBUG)
1238          x->debug=654324;
1239#endif
1240          omFreeBin((ADDRESS)x, rnumber_bin);
1241          x=INT_TO_SR(ui);
1242          return;
1243        }
1244      }
1245      x->s=3;
1246    }
1247    else
1248    {
1249      MP_INT gcd;
1250      mpz_init(&gcd);
1251      mpz_gcd(&gcd,&x->z,&x->n);
1252      x->s=1;
1253      if (mpz_cmp_si(&gcd,(long)1)!=0)
1254      {
1255        MP_INT r;
1256        mpz_init(&r);
1257        MPZ_EXACTDIV(&r,&x->z,&gcd);
1258        mpz_set(&x->z,&r);
1259        MPZ_EXACTDIV(&r,&x->n,&gcd);
1260        mpz_set(&x->n,&r);
1261        mpz_clear(&r);
1262        if (mpz_cmp_si(&x->n,(long)1)==0)
1263        {
1264          mpz_clear(&x->n);
1265          if (mpz_size1(&x->z)<=MP_SMALL)
1266          {
1267            int ui=(int)mpz_get_si(&x->z);
1268            if ((((ui<<3)>>3)==ui)
1269            && (mpz_cmp_si(&x->z,(long)ui)==0))
1270            {
1271              mpz_clear(&x->z);
1272              mpz_clear(&gcd);
1273#if defined(LDEBUG)
1274              x->debug=654324;
1275#endif
1276              omFreeBin((ADDRESS)x, rnumber_bin);
1277              x=INT_TO_SR(ui);
1278              return;
1279            }
1280          }
1281          x->s=3;
1282        }
1283      }
1284      mpz_clear(&gcd);
1285    }
1286  }
1287#ifdef LDEBUG
1288  nlTest(x);
1289#endif
1290}
1291
1292/*2
1293* returns in result->z the lcm(a->z,b->n)
1294*/
1295number nlLcm(number a, number b, const ring r)
1296{
1297  number result;
1298#ifdef LDEBUG
1299  nlTest(a);
1300  nlTest(b);
1301#endif
1302  if ((SR_HDL(b) & SR_INT)
1303  || (b->s==3))
1304  {
1305    // b is 1/(b->n) => b->n is 1 => result is a
1306    return nlCopy(a);
1307  }
1308  result=(number)omAllocBin(rnumber_bin);
1309#if defined(LDEBUG)
1310  result->debug=123456;
1311#endif
1312  result->s=3;
1313  MP_INT gcd;
1314  mpz_init(&gcd);
1315  mpz_init(&result->z);
1316  if (SR_HDL(a) & SR_INT)
1317    mpz_gcd_ui(&gcd,&b->n,ABS(SR_TO_INT(a)));
1318  else
1319    mpz_gcd(&gcd,&a->z,&b->n);
1320  if (mpz_cmp_si(&gcd,(long)1)!=0)
1321  {
1322    MP_INT bt;
1323    mpz_init_set(&bt,&b->n);
1324    MPZ_EXACTDIV(&bt,&bt,&gcd);
1325    if (SR_HDL(a) & SR_INT)
1326      mpz_mul_si(&result->z,&bt,SR_TO_INT(a));
1327    else
1328      mpz_mul(&result->z,&bt,&a->z);
1329    mpz_clear(&bt);
1330  }
1331  else
1332    if (SR_HDL(a) & SR_INT)
1333      mpz_mul_si(&result->z,&b->n,SR_TO_INT(a));
1334    else
1335      mpz_mul(&result->z,&b->n,&a->z);
1336  mpz_clear(&gcd);
1337  if (mpz_size1(&result->z)<=MP_SMALL)
1338  {
1339    int ui=(int)mpz_get_si(&result->z);
1340    if ((((ui<<3)>>3)==ui)
1341    && (mpz_cmp_si(&result->z,(long)ui)==0))
1342    {
1343      mpz_clear(&result->z);
1344      omFreeBin((ADDRESS)result, rnumber_bin);
1345      return INT_TO_SR(ui);
1346    }
1347  }
1348#ifdef LDEBUG
1349  nlTest(result);
1350#endif
1351  return result;
1352}
1353
1354int nlModP(number n, int p)
1355{
1356  if (SR_HDL(n) & SR_INT)
1357  {
1358    int i=SR_TO_INT(n);
1359    if (i<0) return (p-((-i)%p));
1360    return i%p;
1361  }
1362  int iz=(int)mpz_mmod_ui(NULL,&n->z,(unsigned long)p);
1363  if (n->s!=3)
1364  {
1365    int in=mpz_mmod_ui(NULL,&n->n,(unsigned long)p);
1366    #ifdef NV_OPS
1367    if (npPrimeM>NV_MAX_PRIME)
1368    return (int)((long)nvDiv((number)iz,(number)in));
1369    #endif
1370    return (int)((long)npDiv((number)iz,(number)in));
1371  }
1372  return iz;
1373}
1374
1375#ifdef HAVE_RINGS
1376/*2
1377* convert number to GMP and warn if denom != 1
1378*/
1379void nlGMP(number &i, number n)
1380{
1381  // Hier brauche ich einfach die GMP Zahl
1382#ifdef LDEBUG
1383  nlTest(i);
1384#endif
1385  nlNormalize(i);
1386  if (SR_HDL(i) & SR_INT)
1387  {
1388    mpz_set_si((MP_INT*) n, (long) SR_TO_INT(i));
1389    return;
1390  }
1391  if (i->s!=3)
1392  {
1393     WarnS("Omitted denominator during coefficient mapping !");
1394  }
1395  mpz_set((MP_INT*) n, &i->z);
1396}
1397#endif
1398
1399/*2
1400* acces to denominator, other 1 for integers
1401*/
1402number   nlGetDenom(number &n, const ring r)
1403{
1404  if (!(SR_HDL(n) & SR_INT))
1405  {
1406    if (n->s==0)
1407    {
1408      nlNormalize(n);
1409    }
1410    if (!(SR_HDL(n) & SR_INT))
1411    {
1412      if (n->s!=3)
1413      {
1414        number u=(number)omAllocBin(rnumber_bin);
1415        u->s=3;
1416#if defined(LDEBUG)
1417        u->debug=123456;
1418#endif
1419        mpz_init_set(&u->z,&n->n);
1420        {
1421          int ui=(int)mpz_get_si(&u->z);
1422          if ((((ui<<3)>>3)==ui)
1423          && (mpz_cmp_si(&u->z,(long)ui)==0))
1424          {
1425            mpz_clear(&u->z);
1426            omFreeBin((ADDRESS)u, rnumber_bin);
1427            return INT_TO_SR(ui);
1428          }
1429        }
1430        return u;
1431      }
1432    }
1433  }
1434  return INT_TO_SR(1);
1435}
1436
1437/*2
1438* acces to Nominator, nlCopy(n) for integers
1439*/
1440number   nlGetNumerator(number &n, const ring r)
1441{
1442  if (!(SR_HDL(n) & SR_INT))
1443  {
1444    if (n->s==0)
1445    {
1446      nlNormalize(n);
1447    }
1448    if (!(SR_HDL(n) & SR_INT))
1449    {
1450      number u=(number)omAllocBin(rnumber_bin);
1451#if defined(LDEBUG)
1452      u->debug=123456;
1453#endif
1454      u->s=3;
1455      mpz_init_set(&u->z,&n->z);
1456      if (n->s!=3)
1457      {
1458        int ui=(int)mpz_get_si(&u->z);
1459        if ((((ui<<3)>>3)==ui)
1460        && (mpz_cmp_si(&u->z,(long)ui)==0))
1461        {
1462          mpz_clear(&u->z);
1463          omFreeBin((ADDRESS)u, rnumber_bin);
1464          return INT_TO_SR(ui);
1465        }
1466      }
1467      return u;
1468    }
1469  }
1470  return n; // imm. int
1471}
1472
1473/***************************************************************
1474 *
1475 * routines which are needed by Inline(d) routines
1476 *
1477 *******************************************************************/
1478BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
1479{
1480  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1481//  long - short
1482  BOOLEAN bo;
1483  if (SR_HDL(b) & SR_INT)
1484  {
1485    if (a->s!=0) return FALSE;
1486    number n=b; b=a; a=n;
1487  }
1488//  short - long
1489  if (SR_HDL(a) & SR_INT)
1490  {
1491    if (b->s!=0)
1492      return FALSE;
1493    if (((long)a > 0L) && (mpz_isNeg(&b->z)))
1494      return FALSE;
1495    if (((long)a < 0L) && (!mpz_isNeg(&b->z)))
1496      return FALSE;
1497    MP_INT  bb;
1498    mpz_init_set(&bb,&b->n);
1499    mpz_mul_si(&bb,&bb,(long)SR_TO_INT(a));
1500    bo=(mpz_cmp(&bb,&b->z)==0);
1501    mpz_clear(&bb);
1502    return bo;
1503  }
1504// long - long
1505  if (((a->s==1) && (b->s==3))
1506  ||  ((b->s==1) && (a->s==3)))
1507    return FALSE;
1508  if (mpz_isNeg(&a->z)&&(!mpz_isNeg(&b->z)))
1509    return FALSE;
1510  if (mpz_isNeg(&b->z)&&(!mpz_isNeg(&a->z)))
1511    return FALSE;
1512  MP_INT  aa;
1513  MP_INT  bb;
1514  mpz_init_set(&aa,&a->z);
1515  mpz_init_set(&bb,&b->z);
1516  if (a->s<2) mpz_mul(&bb,&bb,&a->n);
1517  if (b->s<2) mpz_mul(&aa,&aa,&b->n);
1518  bo=(mpz_cmp(&aa,&bb)==0);
1519  mpz_clear(&aa);
1520  mpz_clear(&bb);
1521  return bo;
1522}
1523
1524// copy not immediate number a
1525number _nlCopy_NoImm(number a)
1526{
1527  assume(!((SR_HDL(a) & SR_INT)||(a==NULL)));
1528#ifdef LDEBUG
1529  nlTest(a);
1530#endif
1531  number b=(number)omAllocBin(rnumber_bin);
1532#if defined(LDEBUG)
1533  b->debug=123456;
1534#endif
1535  switch (a->s)
1536  {
1537    case 0:
1538    case 1:
1539            mpz_init_set(&b->n,&a->n);
1540    case 3:
1541            mpz_init_set(&b->z,&a->z);
1542            break;
1543  }
1544  b->s = a->s;
1545#ifdef LDEBUG
1546  nlTest(b);
1547#endif
1548  return b;
1549}
1550
1551void _nlDelete_NoImm(number *a)
1552{
1553  {
1554    switch ((*a)->s)
1555    {
1556      case 0:
1557      case 1:
1558        mpz_clear(&(*a)->n);
1559      case 3:
1560        mpz_clear(&(*a)->z);
1561#ifdef LDEBUG
1562        (*a)->s=2;
1563#endif
1564    }
1565    omFreeBin((ADDRESS) *a, rnumber_bin);
1566  }
1567}
1568
1569number _nlNeg_NoImm(number a)
1570{
1571  {
1572    mpz_neg(&a->z,&a->z);
1573    if ((a->s==3) && (mpz_size1(&a->z)<=MP_SMALL))
1574    {
1575      int ui=(int)mpz_get_si(&a->z);
1576      if ((((ui<<3)>>3)==ui)
1577        && (mpz_cmp_si(&a->z,(long)ui)==0))
1578      {
1579        mpz_clear(&a->z);
1580        omFreeBin((ADDRESS)a, rnumber_bin);
1581        a=INT_TO_SR(ui);
1582      }
1583    }
1584  }
1585#ifdef LDEBUG
1586  nlTest(a);
1587#endif
1588  return a;
1589}
1590
1591number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1592{
1593  number u=(number)omAllocBin(rnumber_bin);
1594#if defined(LDEBUG)
1595  u->debug=123456;
1596#endif
1597  mpz_init(&u->z);
1598  if (SR_HDL(b) & SR_INT)
1599  {
1600    number x=a;
1601    a=b;
1602    b=x;
1603  }
1604  if (SR_HDL(a) & SR_INT)
1605  {
1606    switch (b->s)
1607    {
1608      case 0:
1609      case 1:/* a:short, b:1 */
1610      {
1611        MP_INT x;
1612        mpz_init(&x);
1613        mpz_mul_si(&x,&b->n,SR_TO_INT(a));
1614        mpz_add(&u->z,&b->z,&x);
1615        mpz_clear(&x);
1616        if (mpz_cmp_ui(&u->z,(long)0)==0)
1617        {
1618          mpz_clear(&u->z);
1619          omFreeBin((ADDRESS)u, rnumber_bin);
1620          return INT_TO_SR(0);
1621        }
1622        if (mpz_cmp(&u->z,&b->n)==0)
1623        {
1624          mpz_clear(&u->z);
1625          omFreeBin((ADDRESS)u, rnumber_bin);
1626          return INT_TO_SR(1);
1627        }
1628        mpz_init_set(&u->n,&b->n);
1629        u->s = 0;
1630        break;
1631      }
1632      case 3:
1633      {
1634        if ((long)a>0L)
1635          mpz_add_ui(&u->z,&b->z,SR_TO_INT(a));
1636        else
1637          mpz_sub_ui(&u->z,&b->z,-SR_TO_INT(a));
1638        if (mpz_cmp_ui(&u->z,(long)0)==0)
1639        {
1640          mpz_clear(&u->z);
1641          omFreeBin((ADDRESS)u, rnumber_bin);
1642          return INT_TO_SR(0);
1643        }
1644        //u->s = 3;
1645        if (mpz_size1(&u->z)<=MP_SMALL)
1646        {
1647          int ui=(int)mpz_get_si(&u->z);
1648          if ((((ui<<3)>>3)==ui)
1649          && (mpz_cmp_si(&u->z,(long)ui)==0))
1650          {
1651            mpz_clear(&u->z);
1652            omFreeBin((ADDRESS)u, rnumber_bin);
1653            return INT_TO_SR(ui);
1654          }
1655        }
1656        u->s = 3;
1657        break;
1658      }
1659    }
1660  }
1661  else
1662  {
1663    switch (a->s)
1664    {
1665      case 0:
1666      case 1:
1667      {
1668        switch(b->s)
1669        {
1670          case 0:
1671          case 1:
1672          {
1673            MP_INT x;
1674            mpz_init(&x);
1675
1676            mpz_mul(&x,&b->z,&a->n);
1677            mpz_mul(&u->z,&a->z,&b->n);
1678            mpz_add(&u->z,&u->z,&x);
1679            mpz_clear(&x);
1680
1681            if (mpz_cmp_ui(&u->z,(long)0)==0)
1682            {
1683              mpz_clear(&u->z);
1684              omFreeBin((ADDRESS)u, rnumber_bin);
1685              return INT_TO_SR(0);
1686            }
1687            mpz_init(&u->n);
1688            mpz_mul(&u->n,&a->n,&b->n);
1689            if (mpz_cmp(&u->z,&u->n)==0)
1690            {
1691               mpz_clear(&u->z);
1692               mpz_clear(&u->n);
1693               omFreeBin((ADDRESS)u, rnumber_bin);
1694               return INT_TO_SR(1);
1695            }
1696            u->s = 0;
1697            break;
1698          }
1699          case 3: /* a:1 b:3 */
1700          {
1701            mpz_mul(&u->z,&b->z,&a->n);
1702            mpz_add(&u->z,&u->z,&a->z);
1703            if (mpz_cmp_ui(&u->z,(long)0)==0)
1704            {
1705              mpz_clear(&u->z);
1706              omFreeBin((ADDRESS)u, rnumber_bin);
1707              return INT_TO_SR(0);
1708            }
1709            if (mpz_cmp(&u->z,&a->n)==0)
1710            {
1711              mpz_clear(&u->z);
1712              omFreeBin((ADDRESS)u, rnumber_bin);
1713              return INT_TO_SR(1);
1714            }
1715            mpz_init_set(&u->n,&a->n);
1716            u->s = 0;
1717            break;
1718          }
1719        } /*switch (b->s) */
1720        break;
1721      }
1722      case 3:
1723      {
1724        switch(b->s)
1725        {
1726          case 0:
1727          case 1:/* a:3, b:1 */
1728          {
1729            mpz_mul(&u->z,&a->z,&b->n);
1730            mpz_add(&u->z,&u->z,&b->z);
1731            if (mpz_cmp_ui(&u->z,(long)0)==0)
1732            {
1733              mpz_clear(&u->z);
1734              omFreeBin((ADDRESS)u, rnumber_bin);
1735              return INT_TO_SR(0);
1736            }
1737            if (mpz_cmp(&u->z,&b->n)==0)
1738            {
1739              mpz_clear(&u->z);
1740              omFreeBin((ADDRESS)u, rnumber_bin);
1741              return INT_TO_SR(1);
1742            }
1743            mpz_init_set(&u->n,&b->n);
1744            u->s = 0;
1745            break;
1746          }
1747          case 3:
1748          {
1749            mpz_add(&u->z,&a->z,&b->z);
1750            if (mpz_cmp_ui(&u->z,(long)0)==0)
1751            {
1752              mpz_clear(&u->z);
1753              omFreeBin((ADDRESS)u, rnumber_bin);
1754              return INT_TO_SR(0);
1755            }
1756            if (mpz_size1(&u->z)<=MP_SMALL)
1757            {
1758              int ui=(int)mpz_get_si(&u->z);
1759              if ((((ui<<3)>>3)==ui)
1760              && (mpz_cmp_si(&u->z,(long)ui)==0))
1761              {
1762                mpz_clear(&u->z);
1763                omFreeBin((ADDRESS)u, rnumber_bin);
1764                return INT_TO_SR(ui);
1765              }
1766            }
1767            u->s = 3;
1768            break;
1769          }
1770        }
1771        break;
1772      }
1773    }
1774  }
1775#ifdef LDEBUG
1776  nlTest(u);
1777#endif
1778  return u;
1779}
1780
1781number _nlSub_aNoImm_OR_bNoImm(number a, number b)
1782{
1783  number u=(number)omAllocBin(rnumber_bin);
1784#if defined(LDEBUG)
1785  u->debug=123456;
1786#endif
1787  mpz_init(&u->z);
1788  if (SR_HDL(a) & SR_INT)
1789  {
1790    switch (b->s)
1791    {
1792      case 0:
1793      case 1:/* a:short, b:1 */
1794      {
1795        MP_INT x;
1796        mpz_init(&x);
1797        mpz_mul_si(&x,&b->n,SR_TO_INT(a));
1798        mpz_sub(&u->z,&x,&b->z);
1799        mpz_clear(&x);
1800        if (mpz_cmp_ui(&u->z,(long)0)==0)
1801        {
1802          mpz_clear(&u->z);
1803          omFreeBin((ADDRESS)u, rnumber_bin);
1804          return INT_TO_SR(0);
1805        }
1806        if (mpz_cmp(&u->z,&b->n)==0)
1807        {
1808          mpz_clear(&u->z);
1809          omFreeBin((ADDRESS)u, rnumber_bin);
1810          return INT_TO_SR(1);
1811        }
1812        mpz_init_set(&u->n,&b->n);
1813        u->s = 0;
1814        break;
1815      }
1816      case 3:
1817      {
1818        if ((long)a>0L)
1819        {
1820          mpz_sub_ui(&u->z,&b->z,SR_TO_INT(a));
1821          mpz_neg(&u->z,&u->z);
1822        }
1823        else
1824        {
1825          mpz_add_ui(&u->z,&b->z,-SR_TO_INT(a));
1826          mpz_neg(&u->z,&u->z);
1827        }
1828        if (mpz_cmp_ui(&u->z,(long)0)==0)
1829        {
1830          mpz_clear(&u->z);
1831          omFreeBin((ADDRESS)u, rnumber_bin);
1832          return INT_TO_SR(0);
1833        }
1834        if (mpz_size1(&u->z)<=MP_SMALL)
1835        {
1836          int ui=(int)mpz_get_si(&u->z);
1837          if ((((ui<<3)>>3)==ui)
1838          && (mpz_cmp_si(&u->z,(long)ui)==0))
1839          {
1840            mpz_clear(&u->z);
1841            omFreeBin((ADDRESS)u, rnumber_bin);
1842            return INT_TO_SR(ui);
1843          }
1844        }
1845        u->s = 3;
1846        break;
1847      }
1848    }
1849  }
1850  else if (SR_HDL(b) & SR_INT)
1851  {
1852    switch (a->s)
1853    {
1854      case 0:
1855      case 1:/* b:short, a:1 */
1856      {
1857        MP_INT x;
1858        mpz_init(&x);
1859        mpz_mul_si(&x,&a->n,SR_TO_INT(b));
1860        mpz_sub(&u->z,&a->z,&x);
1861        mpz_clear(&x);
1862        if (mpz_cmp_ui(&u->z,(long)0)==0)
1863        {
1864          mpz_clear(&u->z);
1865          omFreeBin((ADDRESS)u, rnumber_bin);
1866          return INT_TO_SR(0);
1867        }
1868        if (mpz_cmp(&u->z,&a->n)==0)
1869        {
1870          mpz_clear(&u->z);
1871          omFreeBin((ADDRESS)u, rnumber_bin);
1872          return INT_TO_SR(1);
1873        }
1874        mpz_init_set(&u->n,&a->n);
1875        u->s = 0;
1876        break;
1877      }
1878      case 3:
1879      {
1880        if ((long)b>0L)
1881        {
1882          mpz_sub_ui(&u->z,&a->z,SR_TO_INT(b));
1883        }
1884        else
1885        {
1886          mpz_add_ui(&u->z,&a->z,-SR_TO_INT(b));
1887        }
1888        if (mpz_cmp_ui(&u->z,(long)0)==0)
1889        {
1890          mpz_clear(&u->z);
1891          omFreeBin((ADDRESS)u, rnumber_bin);
1892          return INT_TO_SR(0);
1893        }
1894        //u->s = 3;
1895        if (mpz_size1(&u->z)<=MP_SMALL)
1896        {
1897          int ui=(int)mpz_get_si(&u->z);
1898          if ((((ui<<3)>>3)==ui)
1899          && (mpz_cmp_si(&u->z,(long)ui)==0))
1900          {
1901            mpz_clear(&u->z);
1902            omFreeBin((ADDRESS)u, rnumber_bin);
1903            return INT_TO_SR(ui);
1904          }
1905        }
1906        u->s = 3;
1907        break;
1908      }
1909    }
1910  }
1911  else
1912  {
1913    switch (a->s)
1914    {
1915      case 0:
1916      case 1:
1917      {
1918        switch(b->s)
1919        {
1920          case 0:
1921          case 1:
1922          {
1923            MP_INT x;
1924            MP_INT y;
1925            mpz_init(&x);
1926            mpz_init(&y);
1927            mpz_mul(&x,&b->z,&a->n);
1928            mpz_mul(&y,&a->z,&b->n);
1929            mpz_sub(&u->z,&y,&x);
1930            mpz_clear(&x);
1931            mpz_clear(&y);
1932            if (mpz_cmp_ui(&u->z,(long)0)==0)
1933            {
1934              mpz_clear(&u->z);
1935              omFreeBin((ADDRESS)u, rnumber_bin);
1936              return INT_TO_SR(0);
1937            }
1938            mpz_init(&u->n);
1939            mpz_mul(&u->n,&a->n,&b->n);
1940            if (mpz_cmp(&u->z,&u->n)==0)
1941            {
1942              mpz_clear(&u->z);
1943              mpz_clear(&u->n);
1944              omFreeBin((ADDRESS)u, rnumber_bin);
1945              return INT_TO_SR(1);
1946            }
1947            u->s = 0;
1948            break;
1949          }
1950          case 3: /* a:1, b:3 */
1951          {
1952            MP_INT x;
1953            mpz_init(&x);
1954            mpz_mul(&x,&b->z,&a->n);
1955            mpz_sub(&u->z,&a->z,&x);
1956            mpz_clear(&x);
1957            if (mpz_cmp_ui(&u->z,(long)0)==0)
1958            {
1959              mpz_clear(&u->z);
1960              omFreeBin((ADDRESS)u, rnumber_bin);
1961              return INT_TO_SR(0);
1962            }
1963            if (mpz_cmp(&u->z,&a->n)==0)
1964            {
1965              mpz_clear(&u->z);
1966              omFreeBin((ADDRESS)u, rnumber_bin);
1967              return INT_TO_SR(1);
1968            }
1969            mpz_init_set(&u->n,&a->n);
1970            u->s = 0;
1971            break;
1972          }
1973        }
1974        break;
1975      }
1976      case 3:
1977      {
1978        switch(b->s)
1979        {
1980          case 0:
1981          case 1: /* a:3, b:1 */
1982          {
1983            MP_INT x;
1984            mpz_init(&x);
1985            mpz_mul(&x,&a->z,&b->n);
1986            mpz_sub(&u->z,&x,&b->z);
1987            mpz_clear(&x);
1988            if (mpz_cmp_ui(&u->z,(long)0)==0)
1989            {
1990              mpz_clear(&u->z);
1991              omFreeBin((ADDRESS)u, rnumber_bin);
1992              return INT_TO_SR(0);
1993            }
1994            if (mpz_cmp(&u->z,&b->n)==0)
1995            {
1996              mpz_clear(&u->z);
1997              omFreeBin((ADDRESS)u, rnumber_bin);
1998              return INT_TO_SR(1);
1999            }
2000            mpz_init_set(&u->n,&b->n);
2001            u->s = 0;
2002            break;
2003          }
2004          case 3: /* a:3 , b:3 */
2005          {
2006            mpz_sub(&u->z,&a->z,&b->z);
2007            if (mpz_cmp_ui(&u->z,(long)0)==0)
2008            {
2009              mpz_clear(&u->z);
2010              omFreeBin((ADDRESS)u, rnumber_bin);
2011              return INT_TO_SR(0);
2012            }
2013            //u->s = 3;
2014            if (mpz_size1(&u->z)<=MP_SMALL)
2015            {
2016              int ui=(int)mpz_get_si(&u->z);
2017              if ((((ui<<3)>>3)==ui)
2018              && (mpz_cmp_si(&u->z,(long)ui)==0))
2019              {
2020                mpz_clear(&u->z);
2021                omFreeBin((ADDRESS)u, rnumber_bin);
2022                return INT_TO_SR(ui);
2023              }
2024            }
2025            u->s = 3;
2026            break;
2027          }
2028        }
2029        break;
2030      }
2031    }
2032  }
2033#ifdef LDEBUG
2034  nlTest(u);
2035#endif
2036  return u;
2037}
2038
2039// a and b are intermediate, but a*b not
2040number _nlMult_aImm_bImm_rNoImm(number a, number b)
2041{
2042  number u=(number)omAllocBin(rnumber_bin);
2043#if defined(LDEBUG)
2044  u->debug=123456;
2045#endif
2046  u->s=3;
2047  mpz_init_set_si(&u->z,(long)SR_TO_INT(a));
2048  mpz_mul_si(&u->z,&u->z,(long)SR_TO_INT(b));
2049#ifdef LDEBUG
2050  nlTest(u);
2051#endif
2052  return u;
2053}
2054
2055// a or b are not immediate
2056number _nlMult_aNoImm_OR_bNoImm(number a, number b)
2057{
2058  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2059  number u=(number)omAllocBin(rnumber_bin);
2060#if defined(LDEBUG)
2061  u->debug=123456;
2062#endif
2063  mpz_init(&u->z);
2064  if (SR_HDL(b) & SR_INT)
2065  {
2066    number x=a;
2067    a=b;
2068    b=x;
2069  }
2070  if (SR_HDL(a) & SR_INT)
2071  {
2072    u->s=b->s;
2073    if (u->s==1) u->s=0;
2074    if ((long)a>0L)
2075    {
2076      mpz_mul_ui(&u->z,&b->z,(unsigned long)SR_TO_INT(a));
2077    }
2078    else
2079    {
2080      if (a==INT_TO_SR(-1))
2081      {
2082        mpz_set(&u->z,&b->z);
2083        mpz_neg(&u->z,&u->z);
2084        u->s=b->s;
2085      }
2086      else
2087      {
2088        mpz_mul_ui(&u->z,&b->z,(unsigned long)-SR_TO_INT(a));
2089        mpz_neg(&u->z,&u->z);
2090      }
2091    }
2092    if (u->s<2)
2093    {
2094      if (mpz_cmp(&u->z,&b->n)==0)
2095      {
2096        mpz_clear(&u->z);
2097        omFreeBin((ADDRESS)u, rnumber_bin);
2098        return INT_TO_SR(1);
2099      }
2100      mpz_init_set(&u->n,&b->n);
2101    }
2102    else //u->s==3
2103    {
2104      if (mpz_size1(&u->z)<=MP_SMALL)
2105      {
2106        int ui=(int)mpz_get_si(&u->z);
2107        if ((((ui<<3)>>3)==ui)
2108            && (mpz_cmp_si(&u->z,(long)ui)==0))
2109        {
2110          mpz_clear(&u->z);
2111          omFreeBin((ADDRESS)u, rnumber_bin);
2112          return INT_TO_SR(ui);
2113        }
2114      }
2115    }
2116  }
2117  else
2118  {
2119    mpz_mul(&u->z,&a->z,&b->z);
2120    u->s = 0;
2121    if(a->s==3)
2122    {
2123      if(b->s==3)
2124      {
2125        u->s = 3;
2126      }
2127      else
2128      {
2129        if (mpz_cmp(&u->z,&b->n)==0)
2130        {
2131          mpz_clear(&u->z);
2132          omFreeBin((ADDRESS)u, rnumber_bin);
2133          return INT_TO_SR(1);
2134        }
2135        mpz_init_set(&u->n,&b->n);
2136      }
2137    }
2138    else
2139    {
2140      if(b->s==3)
2141      {
2142        if (mpz_cmp(&u->z,&a->n)==0)
2143        {
2144          mpz_clear(&u->z);
2145          omFreeBin((ADDRESS)u, rnumber_bin);
2146          return INT_TO_SR(1);
2147        }
2148        mpz_init_set(&u->n,&a->n);
2149      }
2150      else
2151      {
2152        mpz_init(&u->n);
2153        mpz_mul(&u->n,&a->n,&b->n);
2154        if (mpz_cmp(&u->z,&u->n)==0)
2155        {
2156          mpz_clear(&u->z);
2157          mpz_clear(&u->n);
2158          omFreeBin((ADDRESS)u, rnumber_bin);
2159          return INT_TO_SR(1);
2160        }
2161      }
2162    }
2163  }
2164#ifdef LDEBUG
2165  nlTest(u);
2166#endif
2167  return u;
2168}
2169
2170/*2
2171* z := i
2172*/
2173number nlRInit (int i)
2174{
2175  number z=(number)omAllocBin(rnumber_bin);
2176#if defined(LDEBUG)
2177  z->debug=123456;
2178#endif
2179  mpz_init_set_si(&z->z,(long)i);
2180  z->s = 3;
2181  return z;
2182}
2183
2184/*2
2185* z := i/j
2186*/
2187number nlInit2 (int i, int j)
2188{
2189  number z=(number)omAllocBin(rnumber_bin);
2190#if defined(LDEBUG)
2191  z->debug=123456;
2192#endif
2193  mpz_init_set_si(&z->z,(long)i);
2194  mpz_init_set_si(&z->n,(long)j);
2195  z->s = 0;
2196  nlNormalize(z);
2197  return z;
2198}
2199
2200number nlInit2gmp (mpz_t i, mpz_t j)
2201{
2202  number z=(number)omAllocBin(rnumber_bin);
2203#if defined(LDEBUG)
2204  z->debug=123456;
2205#endif
2206  mpz_init_set(&z->z,i);
2207  mpz_init_set(&z->n,j);
2208  z->s = 0;
2209  nlNormalize(z);
2210  return z;
2211}
2212
2213#else // DO_LINLINE
2214
2215// declare immedate routines
2216number nlRInit (int i);
2217BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2218number  _nlCopy_NoImm(number a);
2219number  _nlNeg_NoImm(number a);
2220number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2221number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2222number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2223number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2224
2225#endif
2226
2227/***************************************************************
2228 *
2229 * Routines which might be inlined by p_Numbers.h
2230 *
2231 *******************************************************************/
2232#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2233
2234// routines which are always inlined/static
2235
2236/*2
2237* a = b ?
2238*/
2239LINLINE BOOLEAN nlEqual (number a, number b)
2240{
2241#ifdef LDEBUG
2242  nlTest(a);
2243  nlTest(b);
2244#endif
2245// short - short
2246  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2247  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2248}
2249
2250
2251LINLINE number nlInit (int i, const ring r)
2252{
2253  number n;
2254  if ( ((i << 3) >> 3) == i ) n=INT_TO_SR(i);
2255  else                        n=nlRInit(i);
2256#ifdef LDEBUG
2257  nlTest(n);
2258#endif
2259  return n;
2260}
2261
2262
2263/*2
2264* a == 1 ?
2265*/
2266LINLINE BOOLEAN nlIsOne (number a)
2267{
2268#ifdef LDEBUG
2269  if (a==NULL) return FALSE;
2270  nlTest(a);
2271#endif
2272  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(1));
2273  return FALSE;
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        MP_INT 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        MP_INT 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            MP_INT x;
2446            MP_INT 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            MP_INT 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            MP_INT 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  MP_INT tmp; mpz_init(&tmp);
2658  MP_INT 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.