source: git/kernel/longrat.cc @ 967a9d

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