source: git/Singular/longrat.cc @ 01cd520

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