source: git/kernel/longrat.cc @ b1c0a9

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