source: git/Singular/longrat.cc @ b719a3

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