source: git/Singular/longrat.cc @ a86d4c

spielwiese
Last change on this file since a86d4c was a86d4c, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: minor optimization:nlMult git-svn-id: file:///usr/local/Singular/svn/trunk@4287 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 45.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: longrat.cc,v 1.28 2000-05-02 14:55: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 (a==INT_TO_SR(0)) return INT_TO_SR(0);
1192  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
1193  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1194  {
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    u=(number)AllocSizeOf(rnumber);
1230#if defined(LDEBUG) && ! defined(HAVE_ASO)
1231    u->debug=123456;
1232#endif
1233    mpz_init(&u->z);
1234    if (SR_HDL(b) & SR_INT)
1235    {
1236      number x=a;
1237      a=b;
1238      b=x;
1239    }
1240    if (SR_HDL(a) & SR_INT)
1241    {
1242      u->s=b->s;
1243      if (u->s==1) u->s=0;
1244      if ((int)a>0)
1245      {
1246        mpz_mul_ui(&u->z,&b->z,(unsigned long)SR_TO_INT(a));
1247      }
1248      else
1249      {
1250        if (a==INT_TO_SR(-1))
1251        {
1252          mpz_set(&u->z,&b->z);
1253          mpz_neg(&u->z,&u->z);
1254          u->s=b->s;
1255        }
1256        else
1257        {
1258          mpz_mul_ui(&u->z,&b->z,(unsigned long)-SR_TO_INT(a));
1259          mpz_neg(&u->z,&u->z);
1260        }
1261      }
1262      nlGmpSimple(&u->z);
1263      if (u->s<2)
1264      {
1265        if (mpz_cmp(&u->z,&b->n)==0)
1266        {
1267          mpz_clear(&u->z);
1268          FreeSizeOf((ADDRESS)u,rnumber);
1269          return INT_TO_SR(1);
1270        }
1271        mpz_init_set(&u->n,&b->n);
1272      }
1273      else //u->s==3
1274      {
1275        if (mpz_size1(&u->z)<=MP_SMALL)
1276        {
1277          int ui=(int)mpz_get_si(&u->z);
1278          if ((((ui<<3)>>3)==ui)
1279          && (mpz_cmp_si(&u->z,(long)ui)==0))
1280          {
1281            mpz_clear(&u->z);
1282            FreeSizeOf((ADDRESS)u,rnumber);
1283            return INT_TO_SR(ui);
1284          }
1285        }
1286      }
1287    }
1288    else
1289    {
1290      mpz_mul(&u->z,&a->z,&b->z);
1291      u->s = 0;
1292      if(a->s==3)
1293      {
1294        if(b->s==3)
1295        {
1296          u->s = 3;
1297        }
1298        else
1299        {
1300          if (mpz_cmp(&u->z,&b->n)==0)
1301          {
1302            mpz_clear(&u->z);
1303            FreeSizeOf((ADDRESS)u,rnumber);
1304            return INT_TO_SR(1);
1305          }
1306          mpz_init_set(&u->n,&b->n);
1307        }
1308      }
1309      else
1310      {
1311        if(b->s==3)
1312        {
1313          if (mpz_cmp(&u->z,&a->n)==0)
1314          {
1315            mpz_clear(&u->z);
1316            FreeSizeOf((ADDRESS)u,rnumber);
1317            return INT_TO_SR(1);
1318          }
1319          mpz_init_set(&u->n,&a->n);
1320        }
1321        else
1322        {
1323          mpz_init(&u->n);
1324          mpz_mul(&u->n,&a->n,&b->n);
1325          if (mpz_cmp(&u->z,&u->n)==0)
1326          {
1327            mpz_clear(&u->z);
1328            mpz_clear(&u->n);
1329            FreeSizeOf((ADDRESS)u,rnumber);
1330            return INT_TO_SR(1);
1331          }
1332        }
1333      }
1334    }
1335  }
1336#ifdef LDEBUG
1337  nlTest(u);
1338#endif
1339  return u;
1340}
1341
1342/*2
1343* u := a / b in Z, if b | a (else undefined)
1344*/
1345number   nlExactDiv(number a, number b)
1346{
1347  if (b==INT_TO_SR(0))
1348  {
1349    WerrorS("div. by 0");
1350    return INT_TO_SR(0);
1351  }
1352  if (a==INT_TO_SR(0))
1353    return INT_TO_SR(0);
1354  number u;
1355  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1356  {
1357    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
1358    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
1359    {
1360      return nlRInit(1<<28);
1361    }
1362    int aa=SR_TO_INT(a);
1363    int bb=SR_TO_INT(b);
1364    return INT_TO_SR(aa/bb);
1365  }
1366  number bb=NULL;
1367  if (SR_HDL(b) & SR_INT)
1368  {
1369    bb=nlRInit(SR_TO_INT(b));
1370    b=bb;
1371  }
1372  u=(number)AllocSizeOf(rnumber);
1373#if defined(LDEBUG) && ! defined(HAVE_ASO)
1374  u->debug=123456;
1375#endif
1376  mpz_init(&u->z);
1377  /* u=a/b */
1378  u->s = 3;
1379  MPZ_EXACTDIV(&u->z,&a->z,&b->z);
1380  if (bb!=NULL)
1381  {
1382    mpz_clear(&bb->z);
1383    FreeSizeOf((ADDRESS)bb,rnumber);
1384  }
1385  if (mpz_size1(&u->z)<=MP_SMALL)
1386  {
1387    int ui=(int)mpz_get_si(&u->z);
1388    if ((((ui<<3)>>3)==ui)
1389    && (mpz_cmp_si(&u->z,(long)ui)==0))
1390    {
1391      mpz_clear(&u->z);
1392      FreeSizeOf((ADDRESS)u,rnumber);
1393      u=INT_TO_SR(ui);
1394    }
1395  }
1396#ifdef LDEBUG
1397  nlTest(u);
1398#endif
1399  return u;
1400}
1401
1402/*2
1403* u := a / b in Z
1404*/
1405number nlIntDiv (number a, number b)
1406{
1407  if (b==INT_TO_SR(0))
1408  {
1409    WerrorS("div. by 0");
1410    return INT_TO_SR(0);
1411  }
1412  if (a==INT_TO_SR(0))
1413    return INT_TO_SR(0);
1414  number u;
1415  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1416  {
1417    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
1418    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
1419    {
1420      return nlRInit(1<<28);
1421    }
1422    /* consider the signs of a and b:
1423    *  + + -> u=a+b-1
1424    *  + - -> u=a-b-1
1425    *  - + -> u=a-b+1
1426    *  - - -> u=a+b+1
1427    */
1428    int aa=SR_TO_INT(a);
1429    int bb=SR_TO_INT(b);
1430    if (aa>0)
1431    {
1432      if (bb>0)
1433        return INT_TO_SR((aa+bb-1)/bb);
1434      else
1435        return INT_TO_SR((aa-bb-1)/bb);
1436    }
1437    else
1438    {
1439      if (bb>0)
1440      {
1441        return INT_TO_SR((aa-bb+1)/bb);
1442      }
1443      else
1444      {
1445        return INT_TO_SR((aa+bb+1)/bb);
1446      }
1447    }
1448  }
1449  if (SR_HDL(a) & SR_INT)
1450  {
1451    /* the small int -(1<<28) divided by 2^28 is 1   */
1452    if (a==INT_TO_SR(-(1<<28)))
1453    {
1454      if(mpz_cmp_si(&b->z,(long)(1<<28))==0)
1455      {
1456        return INT_TO_SR(-1);
1457      }
1458    }
1459    /* a is a small and b is a large int: -> 0 */
1460    return INT_TO_SR(0);
1461  }
1462  number bb=NULL;
1463  if (SR_HDL(b) & SR_INT)
1464  {
1465    bb=nlRInit(SR_TO_INT(b));
1466    b=bb;
1467  }
1468  u=(number)AllocSizeOf(rnumber);
1469#if defined(LDEBUG) && ! defined(HAVE_ASO)
1470  u->debug=123456;
1471#endif
1472  mpz_init_set(&u->z,&a->z);
1473  /* consider the signs of a and b:
1474  *  + + -> u=a+b-1
1475  *  + - -> u=a-b-1
1476  *  - + -> u=a-b+1
1477  *  - - -> u=a+b+1
1478  */
1479  if (mpz_isNeg(&a->z))
1480  {
1481    if (mpz_isNeg(&b->z))
1482    {
1483      mpz_add(&u->z,&u->z,&b->z);
1484    }
1485    else
1486    {
1487      mpz_sub(&u->z,&u->z,&b->z);
1488    }
1489    mpz_add_ui(&u->z,&u->z,1);
1490  }
1491  else
1492  {
1493    if (mpz_isNeg(&b->z))
1494    {
1495      mpz_sub(&u->z,&u->z,&b->z);
1496    }
1497    else
1498    {
1499      mpz_add(&u->z,&u->z,&b->z);
1500    }
1501    mpz_sub_ui(&u->z,&u->z,1);
1502  }
1503  /* u=u/b */
1504  u->s = 3;
1505  MPZ_DIV(&u->z,&u->z,&b->z);
1506  nlGmpSimple(&u->z);
1507  if (bb!=NULL)
1508  {
1509    mpz_clear(&bb->z);
1510    FreeSizeOf((ADDRESS)bb,rnumber);
1511  }
1512  if (mpz_size1(&u->z)<=MP_SMALL)
1513  {
1514    int ui=(int)mpz_get_si(&u->z);
1515    if ((((ui<<3)>>3)==ui)
1516    && (mpz_cmp_si(&u->z,(long)ui)==0))
1517    {
1518      mpz_clear(&u->z);
1519      FreeSizeOf((ADDRESS)u,rnumber);
1520      u=INT_TO_SR(ui);
1521    }
1522  }
1523#ifdef LDEBUG
1524  nlTest(u);
1525#endif
1526  return u;
1527}
1528
1529/*2
1530* u := a mod b in Z, u>=0
1531*/
1532number nlIntMod (number a, number b)
1533{
1534  if (b==INT_TO_SR(0))
1535  {
1536    WerrorS("div. by 0");
1537    return INT_TO_SR(0);
1538  }
1539  if (a==INT_TO_SR(0))
1540    return INT_TO_SR(0);
1541  number u;
1542  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1543  {
1544    if ((int)a>0)
1545    {
1546      if ((int)b>0)
1547        return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
1548      else
1549        return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
1550    }
1551    else
1552    {
1553      if ((int)b>0)
1554      {
1555        int i=(-SR_TO_INT(a))%SR_TO_INT(b);
1556        if ( i != 0 ) i = (SR_TO_INT(b))-i;
1557        return INT_TO_SR(i);
1558      }
1559      else
1560      {
1561        int i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
1562        if ( i != 0 ) i = (-SR_TO_INT(b))-i;
1563        return INT_TO_SR(i);
1564      }
1565    }
1566  }
1567  if (SR_HDL(a) & SR_INT)
1568  {
1569    /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
1570    if ((int)a<0)
1571    {
1572      if (mpz_isNeg(&b->z))
1573        return nlSub(a,b);
1574      /*else*/
1575        return nlAdd(a,b);
1576    }
1577    /*else*/
1578      return a;
1579  }
1580  number bb=NULL;
1581  if (SR_HDL(b) & SR_INT)
1582  {
1583    bb=nlRInit(SR_TO_INT(b));
1584    b=bb;
1585  }
1586  u=(number)AllocSizeOf(rnumber);
1587#if defined(LDEBUG) && ! defined(HAVE_ASO)
1588  u->debug=123456;
1589#endif
1590  mpz_init(&u->z);
1591  u->s = 3;
1592  mpz_mod(&u->z,&a->z,&b->z);
1593  if (bb!=NULL)
1594  {
1595    mpz_clear(&bb->z);
1596    FreeSizeOf((ADDRESS)bb,rnumber);
1597  }
1598  if (mpz_isNeg(&u->z))
1599  {
1600    if (mpz_isNeg(&b->z))
1601      mpz_sub(&u->z,&u->z,&b->z);
1602    else
1603      mpz_add(&u->z,&u->z,&b->z);
1604  }
1605  nlGmpSimple(&u->z);
1606  if (mpz_size1(&u->z)<=MP_SMALL)
1607  {
1608    int ui=(int)mpz_get_si(&u->z);
1609    if ((((ui<<3)>>3)==ui)
1610    && (mpz_cmp_si(&u->z,(long)ui)==0))
1611    {
1612      mpz_clear(&u->z);
1613      FreeSizeOf((ADDRESS)u,rnumber);
1614      u=INT_TO_SR(ui);
1615    }
1616  }
1617#ifdef LDEBUG
1618  nlTest(u);
1619#endif
1620  return u;
1621}
1622
1623/*2
1624* u := a / b
1625*/
1626number nlDiv (number a, number b)
1627{
1628  number u;
1629  if (nlIsZero(b))
1630  {
1631    WerrorS("div. by 0");
1632    return INT_TO_SR(0);
1633  }
1634  u=(number)AllocSizeOf(rnumber);
1635  u->s=0;
1636#if defined(LDEBUG) && ! defined(HAVE_ASO)
1637  u->debug=123456;
1638#endif
1639// ---------- short / short ------------------------------------
1640  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1641  {
1642    int i=SR_TO_INT(a);
1643    int j=SR_TO_INT(b);
1644    int r=i%j;
1645    if (r==0)
1646    {
1647      FreeSizeOf((ADDRESS)u,rnumber);
1648      return INT_TO_SR(i/j);
1649    }
1650    mpz_init_set_si(&u->z,(long)i);
1651    mpz_init_set_si(&u->n,(long)j);
1652  }
1653  else
1654  {
1655    mpz_init(&u->z);
1656// ---------- short / long ------------------------------------
1657    if (SR_HDL(a) & SR_INT)
1658    {
1659      // short a / (z/n) -> (a*n)/z
1660      if (b->s<2)
1661      {
1662        if ((int)a>=0)
1663          mpz_mul_ui(&u->z,&b->n,SR_TO_INT(a));
1664        else
1665        {
1666          mpz_mul_ui(&u->z,&b->n,-SR_TO_INT(a));
1667          mpz_neg(&u->z,&u->z);
1668        }
1669      }
1670      else
1671      // short a / long z -> a/z
1672      {
1673        mpz_set_si(&u->z,SR_TO_INT(a));
1674      }
1675      if (mpz_cmp(&u->z,&b->z)==0)
1676      {
1677        mpz_clear(&u->z);
1678        FreeSizeOf((ADDRESS)u,rnumber);
1679        return INT_TO_SR(1);
1680      }
1681      mpz_init_set(&u->n,&b->z);
1682    }
1683// ---------- long / short ------------------------------------
1684    else if (SR_HDL(b) & SR_INT)
1685    {
1686      mpz_set(&u->z,&a->z);
1687      // (z/n) / b -> z/(n*b)
1688      if (a->s<2)
1689      {
1690        mpz_init_set(&u->n,&a->n);
1691        if ((int)b>0)
1692          mpz_mul_ui(&u->n,&u->n,SR_TO_INT(b));
1693        else
1694        {
1695          mpz_mul_ui(&u->n,&u->n,-SR_TO_INT(b));
1696          mpz_neg(&u->z,&u->z);
1697        }
1698      }
1699      else
1700      // long z / short b -> z/b
1701      {
1702        //mpz_set(&u->z,&a->z);
1703        mpz_init_set_si(&u->n,SR_TO_INT(b));
1704      }
1705    }
1706// ---------- long / long ------------------------------------
1707    else
1708    {
1709      //u->s=0;
1710      mpz_set(&u->z,&a->z);
1711      mpz_init_set(&u->n,&b->z);
1712      if (a->s<2) mpz_mul(&u->n,&u->n,&a->n);
1713      if (b->s<2) mpz_mul(&u->z,&u->z,&b->n);
1714    }
1715  }
1716  if (mpz_isNeg(&u->n))
1717  {
1718    mpz_neg(&u->z,&u->z);
1719    mpz_neg(&u->n,&u->n);
1720  }
1721  if (mpz_cmp_si(&u->n,(long)1)==0)
1722  {
1723    mpz_clear(&u->n);
1724    if (mpz_size1(&u->z)<=MP_SMALL)
1725    {
1726      int ui=(int)mpz_get_si(&u->z);
1727      if ((((ui<<3)>>3)==ui)
1728      && (mpz_cmp_si(&u->z,(long)ui)==0))
1729      {
1730        mpz_clear(&u->z);
1731        FreeSizeOf((ADDRESS)u,rnumber);
1732        return INT_TO_SR(ui);
1733      }
1734    }
1735    u->s=3;
1736  }
1737#ifdef LDEBUG
1738  nlTest(u);
1739#endif
1740  return u;
1741}
1742
1743#if (defined(i386)) && (defined(HAVE_LIBGMP1))
1744/*
1745* compute x^exp for x in Z
1746* there seems to be an bug in mpz_pow_ui for i386
1747*/
1748static inline void nlPow (MP_INT * res,MP_INT * x,int exp)
1749{
1750  if (exp==0)
1751  {
1752    mpz_set_si(res,(long)1);
1753  }
1754  else
1755  {
1756    MP_INT xh;
1757    mpz_init(&xh);
1758    mpz_set(res,x);
1759    exp--;
1760    while (exp!=0)
1761    {
1762      mpz_mul(&xh,res,x);
1763      mpz_set(res,&xh);
1764      exp--;
1765    }
1766    mpz_clear(&xh);
1767  }
1768}
1769#endif
1770
1771/*2
1772* u:= x ^ exp
1773*/
1774void nlPower (number x,int exp,number * u)
1775{
1776  if (!nlIsZero(x))
1777  {
1778#ifdef LDEBUG
1779    nlTest(x);
1780#endif
1781    number aa=NULL;
1782    if (SR_HDL(x) & SR_INT)
1783    {
1784      aa=nlRInit(SR_TO_INT(x));
1785      x=aa;
1786    }
1787    *u=(number)AllocSizeOf(rnumber);
1788#if defined(LDEBUG) && ! defined(HAVE_ASO)
1789    (*u)->debug=123456;
1790#endif
1791    mpz_init(&(*u)->z);
1792    (*u)->s = x->s;
1793#if (!defined(i386)) || (defined(HAVE_LIBGMP2))
1794    mpz_pow_ui(&(*u)->z,&x->z,(unsigned long)exp);
1795#else
1796    nlPow(&(*u)->z,&x->z,exp);
1797#endif
1798    if (x->s<2)
1799    {
1800      mpz_init(&(*u)->n);
1801#if (!defined(i386)) || (defined(HAVE_LIBGMP2))
1802      mpz_pow_ui(&(*u)->n,&x->n,(unsigned long)exp);
1803#else
1804      nlPow(&(*u)->n,&x->n,exp);
1805#endif
1806    }
1807    if (aa!=NULL)
1808    {
1809      mpz_clear(&aa->z);
1810      FreeSizeOf((ADDRESS)aa,rnumber);
1811    }
1812    if (((*u)->s<=2) && (mpz_cmp_si(&(*u)->n,(long)1)==0))
1813    {
1814      mpz_clear(&(*u)->n);
1815      (*u)->s=3;
1816    }
1817    if (((*u)->s==3) && (mpz_size1(&(*u)->z)<=MP_SMALL))
1818    {
1819      int ui=(int)mpz_get_si(&(*u)->z);
1820      if ((((ui<<3)>>3)==ui)
1821      && (mpz_cmp_si(&(*u)->z,(long)ui)==0))
1822      {
1823        mpz_clear(&(*u)->z);
1824        FreeSizeOf((ADDRESS)*u,rnumber);
1825        *u=INT_TO_SR(ui);
1826      }
1827    }
1828  }
1829  else
1830    *u = INT_TO_SR(0);
1831#ifdef LDEBUG
1832  nlTest(*u);
1833#endif
1834}
1835
1836BOOLEAN nlIsZero (number a)
1837{
1838// #ifdef LDEBUG
1839//   nlTest(a);
1840// #endif
1841  //if (a==INT_TO_SR(0)) return TRUE;
1842  //if (SR_HDL(a) & SR_INT) return FALSE;
1843  //number aa=nlCopy(a);
1844  //nlNormalize(aa);
1845  return (a==INT_TO_SR(0));
1846}
1847
1848/*2
1849* za >= 0 ?
1850*/
1851BOOLEAN nlGreaterZero (number a)
1852{
1853#ifdef LDEBUG
1854  nlTest(a);
1855#endif
1856  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>=0;
1857  return (!mpz_isNeg(&a->z));
1858}
1859
1860/*2
1861* a > b ?
1862*/
1863BOOLEAN nlGreater (number a, number b)
1864{
1865#ifdef LDEBUG
1866  nlTest(a);
1867  nlTest(b);
1868#endif
1869  number r;
1870  BOOLEAN rr;
1871  r=nlSub(a,b);
1872  rr=(!nlIsZero(r)) && (nlGreaterZero(r));
1873  nlDelete(&r);
1874  return rr;
1875}
1876
1877/*2
1878* a = b ?
1879*/
1880BOOLEAN nlEqual (number a, number b)
1881{
1882#ifdef LDEBUG
1883  nlTest(a);
1884  nlTest(b);
1885#endif
1886// short - short
1887  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
1888//  long - short
1889  BOOLEAN bo;
1890  if (SR_HDL(b) & SR_INT)
1891  {
1892    if (a->s!=0) return FALSE;
1893    number n=b; b=a; a=n;
1894  }
1895//  short - long
1896  if (SR_HDL(a) & SR_INT)
1897  {
1898    if (b->s!=0)
1899      return FALSE;
1900    if (((int)a > 0) && (mpz_isNeg(&b->z)))
1901      return FALSE;
1902    if (((int)a < 0) && (!mpz_isNeg(&b->z)))
1903      return FALSE;
1904    MP_INT  bb;
1905    mpz_init_set(&bb,&b->n);
1906    if ((int)a<0)
1907    {
1908      mpz_neg(&bb,&bb);
1909      mpz_mul_ui(&bb,&bb,(long)-SR_TO_INT(a));
1910    }
1911    else
1912    {
1913      mpz_mul_ui(&bb,&bb,(long)SR_TO_INT(a));
1914    }
1915    bo=(mpz_cmp(&bb,&b->z)==0);
1916    mpz_clear(&bb);
1917    return bo;
1918  }
1919// long - long
1920  if (((a->s==1) && (b->s==3))
1921  ||  ((b->s==1) && (a->s==3)))
1922    return FALSE;
1923  if (mpz_isNeg(&a->z)&&(!mpz_isNeg(&b->z)))
1924    return FALSE;
1925  if (mpz_isNeg(&b->z)&&(!mpz_isNeg(&a->z)))
1926    return FALSE;
1927  nlGmpSimple(&a->z);
1928  nlGmpSimple(&b->z);
1929  if (a->s<2)
1930    nlGmpSimple(&a->n);
1931  if (b->s<2)
1932    nlGmpSimple(&b->n);
1933  MP_INT  aa;
1934  MP_INT  bb;
1935  mpz_init_set(&aa,&a->z);
1936  mpz_init_set(&bb,&b->z);
1937  if (a->s<2) mpz_mul(&bb,&bb,&a->n);
1938  if (b->s<2) mpz_mul(&aa,&aa,&b->n);
1939  bo=(mpz_cmp(&aa,&bb)==0);
1940  mpz_clear(&aa);
1941  mpz_clear(&bb);
1942  return bo;
1943}
1944
1945/*2
1946* a == 1 ?
1947*/
1948BOOLEAN nlIsOne (number a)
1949{
1950#ifdef LDEBUG
1951  if (a==NULL) return FALSE;
1952  nlTest(a);
1953#endif
1954  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(1));
1955  return FALSE;
1956}
1957
1958/*2
1959* a == -1 ?
1960*/
1961BOOLEAN nlIsMOne (number a)
1962{
1963#ifdef LDEBUG
1964  if (a==NULL) return FALSE;
1965  nlTest(a);
1966#endif
1967  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1));
1968  return FALSE;
1969}
1970
1971/*2
1972* result =gcd(a,b)
1973*/
1974number nlGcd(number a, number b)
1975{
1976  number result;
1977#ifdef LDEBUG
1978  nlTest(a);
1979  nlTest(b);
1980#endif
1981  //nlNormalize(a);
1982  //nlNormalize(b);
1983  if ((SR_HDL(a)==5)||(a==INT_TO_SR(-1))
1984  ||  (SR_HDL(b)==5)||(b==INT_TO_SR(-1))) return INT_TO_SR(1);
1985  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1986  {
1987    int i=SR_TO_INT(a);
1988    int j=SR_TO_INT(b);
1989    if((i==0)||(j==0))
1990      return INT_TO_SR(1);
1991    int l;
1992    i=ABS(i);
1993    j=ABS(j);
1994    do
1995    {
1996      l=i%j;
1997      i=j;
1998      j=l;
1999    } while (l!=0);
2000    return INT_TO_SR(i);
2001  }
2002  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
2003  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
2004  number aa=NULL;
2005  if (SR_HDL(a) & SR_INT)
2006  {
2007    aa=nlRInit(SR_TO_INT(a));
2008    a=aa;
2009  }
2010  else
2011  if (SR_HDL(b) & SR_INT)
2012  {
2013    aa=nlRInit(SR_TO_INT(b));
2014    b=aa;
2015  }
2016  result=(number)AllocSizeOf(rnumber);
2017#if defined(LDEBUG) && ! defined(HAVE_ASO)
2018  result->debug=123456;
2019#endif
2020  mpz_init(&result->z);
2021  mpz_gcd(&result->z,&a->z,&b->z);
2022  nlGmpSimple(&result->z);
2023  result->s = 3;
2024  if (mpz_size1(&result->z)<=MP_SMALL)
2025  {
2026    int ui=(int)mpz_get_si(&result->z);
2027    if ((((ui<<3)>>3)==ui)
2028    && (mpz_cmp_si(&result->z,(long)ui)==0))
2029    {
2030      mpz_clear(&result->z);
2031      FreeSizeOf((ADDRESS)result,rnumber);
2032      result=INT_TO_SR(ui);
2033    }
2034  }
2035  if (aa!=NULL)
2036  {
2037    mpz_clear(&aa->z);
2038    FreeSizeOf((ADDRESS)aa,rnumber);
2039  }
2040#ifdef LDEBUG
2041  nlTest(result);
2042#endif
2043  return result;
2044}
2045
2046/*2
2047* simplify x
2048*/
2049void nlNormalize (number &x)
2050{
2051  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
2052    return;
2053#ifdef LDEBUG
2054  if (!nlTest(x)) return;
2055#endif
2056  if (x->s==3)
2057  {
2058    if (mpz_cmp_ui(&x->z,(long)0)==0)
2059    {
2060      nlDelete(&x);
2061      x=nlInit(0);
2062      return;
2063    }
2064    if (mpz_size1(&x->z)<=MP_SMALL)
2065    {
2066      int ui=(int)mpz_get_si(&x->z);
2067      if ((((ui<<3)>>3)==ui)
2068      && (mpz_cmp_si(&x->z,(long)ui)==0))
2069      {
2070        mpz_clear(&x->z);
2071        FreeSizeOf((ADDRESS)x,rnumber);
2072        x=INT_TO_SR(ui);
2073        return;
2074      }
2075    }
2076  }
2077  if (x->s!=1) nlGmpSimple(&x->z);
2078  if (x->s==0)
2079  {
2080    BOOLEAN divided=FALSE;
2081    MP_INT gcd;
2082    mpz_init(&gcd);
2083    mpz_gcd(&gcd,&x->z,&x->n);
2084    x->s=1;
2085    if (mpz_cmp_si(&gcd,(long)1)!=0)
2086    {
2087      MP_INT r;
2088      mpz_init(&r);
2089      MPZ_EXACTDIV(&r,&x->z,&gcd);
2090      mpz_set(&x->z,&r);
2091      MPZ_EXACTDIV(&r,&x->n,&gcd);
2092      mpz_set(&x->n,&r);
2093      mpz_clear(&r);
2094      nlGmpSimple(&x->z);
2095      divided=TRUE;
2096    }
2097    mpz_clear(&gcd);
2098    nlGmpSimple(&x->n);
2099    if (mpz_cmp_si(&x->n,(long)1)==0)
2100    {
2101      mpz_clear(&x->n);
2102      if (mpz_size1(&x->z)<=MP_SMALL)
2103      {
2104        int ui=(int)mpz_get_si(&x->z);
2105        if ((((ui<<3)>>3)==ui)
2106        && (mpz_cmp_si(&x->z,(long)ui)==0))
2107        {
2108          mpz_clear(&x->z);
2109#if defined(LDEBUG) && ! defined(HAVE_ASO)
2110          x->debug=654324;
2111#endif
2112          FreeSizeOf((ADDRESS)x,rnumber);
2113          x=INT_TO_SR(ui);
2114          return;
2115        }
2116      }
2117      x->s=3;
2118    }
2119    else if (divided)
2120    {
2121      _mpz_realloc(&x->n,mpz_size1(&x->n));
2122    }
2123    if (divided) _mpz_realloc(&x->z,mpz_size1(&x->z));
2124  }
2125#ifdef LDEBUG
2126  nlTest(x);
2127#endif
2128}
2129
2130/*2
2131* returns in result->z the lcm(a->z,b->n)
2132*/
2133number nlLcm(number a, number b)
2134{
2135  number result;
2136#ifdef LDEBUG
2137  nlTest(a);
2138  nlTest(b);
2139#endif
2140  if ((SR_HDL(b) & SR_INT)
2141  || (b->s==3))
2142  {
2143    // b is 1/(b->n) => b->n is 1 => result is a
2144    return nlCopy(a);
2145  }
2146  number aa=NULL;
2147  if (SR_HDL(a) & SR_INT)
2148  {
2149    aa=nlRInit(SR_TO_INT(a));
2150    a=aa;
2151  }
2152  result=(number)AllocSizeOf(rnumber);
2153#if defined(LDEBUG) && ! defined(HAVE_ASO)
2154  result->debug=123456;
2155#endif
2156  result->s=3;
2157  MP_INT gcd;
2158  mpz_init(&gcd);
2159  mpz_init(&result->z);
2160  mpz_gcd(&gcd,&a->z,&b->n);
2161  if (mpz_cmp_si(&gcd,(long)1)!=0)
2162  {
2163    MP_INT bt;
2164    mpz_init_set(&bt,&b->n);
2165    MPZ_EXACTDIV(&bt,&bt,&gcd);
2166    mpz_mul(&result->z,&a->z,&bt);
2167    mpz_clear(&bt);
2168  }
2169  else
2170    mpz_mul(&result->z,&a->z,&b->n);
2171  mpz_clear(&gcd);
2172  if (aa!=NULL)
2173  {
2174    mpz_clear(&aa->z);
2175    FreeSizeOf((ADDRESS)aa,rnumber);
2176  }
2177  nlGmpSimple(&result->z);
2178  if (mpz_size1(&result->z)<=MP_SMALL)
2179  {
2180    int ui=(int)mpz_get_si(&result->z);
2181    if ((((ui<<3)>>3)==ui)
2182    && (mpz_cmp_si(&result->z,(long)ui)==0))
2183    {
2184      mpz_clear(&result->z);
2185      FreeSizeOf((ADDRESS)result,rnumber);
2186      return INT_TO_SR(ui);
2187    }
2188  }
2189#ifdef LDEBUG
2190  nlTest(result);
2191#endif
2192  return result;
2193}
2194
2195int nlModP(number n, int p)
2196{
2197  if (SR_HDL(n) & SR_INT)
2198  {
2199    int i=SR_TO_INT(n);
2200    if (i<0) return (p-((-i)%p));
2201    return i%p;
2202  }
2203  int iz=mpz_mmod_ui(NULL,&n->z,(unsigned long)p);
2204  if (n->s!=3)
2205  {
2206    int in=mpz_mmod_ui(NULL,&n->n,(unsigned long)p);
2207    return (int)npDiv((number)iz,(number)in);
2208  }
2209  return iz;
2210}
2211
2212/*2
2213* acces to denominator, other 1 for integers
2214*/
2215number   nlGetDenom(number &n)
2216{
2217  if (!(SR_HDL(n) & SR_INT))
2218  {
2219    if (n->s==0)
2220    {
2221      nlNormalize(n);
2222    }
2223    if (!(SR_HDL(n) & SR_INT))
2224    {
2225      if (n->s!=3)
2226      {
2227        number u=(number)AllocSizeOf(rnumber);
2228        u->s=3;
2229#if defined(LDEBUG) && ! defined(HAVE_ASO)
2230        u->debug=123456;
2231#endif
2232
2233        mpz_init_set(&u->z,&n->n);
2234        {
2235          int ui=(int)mpz_get_si(&u->z);
2236          if ((((ui<<3)>>3)==ui)
2237          && (mpz_cmp_si(&u->z,(long)ui)==0))
2238          {
2239            mpz_clear(&u->z);
2240            FreeSizeOf((ADDRESS)u,rnumber);
2241            return INT_TO_SR(ui);
2242          }
2243        }
2244        return u;
2245      }
2246    }
2247  }
2248  return INT_TO_SR(1);
2249}
Note: See TracBrowser for help on using the repository browser.