source: git/Singular/longrat.cc @ d4373f

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