source: git/Singular/longrat.cc @ 512a2b

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