source: git/Singular/longrat.cc @ c4bbf1f

spielwiese
Last change on this file since c4bbf1f was 1b45a8, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: fixed a minor but very nasty bug: nlSetMap was resetting npPrimeM (longrat.cc) git-svn-id: file:///usr/local/Singular/svn/trunk@1392 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 76.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: longrat.cc,v 1.16 1998-04-17 07:53:43 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#ifdef HAVE_GMP
64/*3
65* parameter s in number:
66* 0 (or FALSE): not normalised rational
67* 1 (or TRUE):  normalised rational
68* 3          :  integer with n==NULL
69*/
70/*3
71**  'SR_INT' is the type of those integers small enough to fit into  29  bits.
72**  Therefor the value range of this small integers is: $-2^{28}...2^{28}-1$.
73**
74**  Small integers are represented by an immediate integer handle, containing
75**  the value instead of pointing  to  it,  which  has  the  following  form:
76**
77**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
78**      | guard | sign  | bit   | bit   |       | bit   | tag   | tag   |
79**      | bit   | bit   | 27    | 26    |       | 0     | 0     | 1     |
80**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
81**
82**  Immediate integers handles carry the tag 'SR_INT', i.e. the last bit is 1.
83**  This distuingishes immediate integers from other handles which  point  to
84**  structures aligned on 4 byte boundaries and therefor have last bit  zero.
85**  (The second bit is reserved as tag to allow extensions of  this  scheme.)
86**  Using immediates as pointers and dereferencing them gives address errors.
87**
88**  To aid overflow check the most significant two bits must always be equal,
89**  that is to say that the sign bit of immediate integers has a  guard  bit.
90**
91**  The macros 'INT_TO_SR' and 'SR_TO_INT' should be used to convert  between
92**  a small integer value and its representation as immediate integer handle.
93**
94**  Large integers and rationals are represented by z and n
95**  where n may be undefined (if s==3)
96**  NULL represents only deleted values
97*/
98#define SR_HDL(A) ((long)(A))
99/*#define SR_INT    1*/
100/*#define INT_TO_SR(INT)  ((number) (((long)INT << 2) + SR_INT))*/
101#define SR_TO_INT(SR)   (((long)SR) >> 2)
102
103#define MP_SMALL 1
104//#define mpz_isNeg(A) (mpz_cmp_si(A,(long)0)<0)
105#ifdef HAVE_LIBGMP1
106#define mpz_isNeg(A) ((A)->size<0)
107#define mpz_limb_size(A) ((A)->size)
108#define mpz_limb_d(A) ((A)->d)
109#define MPZ_DIV(A,B,C) mpz_div((A),(B),(C))
110#define MPZ_EXACTDIV(A,B,C) mpz_div((A),(B),(C))
111#else
112#define mpz_isNeg(A) ((A)->_mp_size<0)
113#define mpz_limb_size(A) ((A)->_mp_size)
114#define mpz_limb_d(A) ((A)->_mp_d)
115#define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
116#ifdef HAVE_SMALLGMP
117#define MPZ_EXACTDIV(A,B,C) mpz_tdiv_q((A),(B),(C))
118#else
119#define MPZ_EXACTDIV(A,B,C) mpz_divexact((A),(B),(C))
120#endif
121#define nlGmpSimple(A)
122#endif
123
124
125#ifdef LDEBUG
126#define nlTest(a) nlDBTest(a,__FILE__,__LINE__)
127BOOLEAN nlDBTest(number a, char *f,int l)
128{
129  if (a==NULL)
130  {
131    Print("!!longrat: NULL in %s:%d\n",f,l);
132    return FALSE;
133  }
134  //if ((int)a==1) Print("!! 0x1 as number ? %s %d\n",f,l);
135  if ((((int)a)&3)==3)
136  {
137    Print(" !!longrat:ptr(3) in %s:%d\n",f,l);
138    return FALSE;
139  }
140  if ((((int)a)&3)==1)
141  {
142    if (((((int)a)<<1)>>1)!=((int)a))
143    {
144      Print(" !!longrat:arith:%x in %s:%d\n",(int)a, f,l);
145      return FALSE;
146    }
147    return TRUE;
148  }
149#ifdef MDEBUG
150  mmDBTestBlock(a,sizeof(*a),f,l);
151#endif
152  if (a->debug!=123456)
153  {
154    Print("!!longrat:debug:%d in %s:%d\n",a->debug,f,l);
155    a->debug=123456;
156    return FALSE;
157  }
158  if ((a->s<0)||(a->s>4))
159  {
160    Print("!!longrat:s=%d in %s:%d\n",a->s,f,l);
161    return FALSE;
162  }
163#ifdef MDEBUG
164#ifdef HAVE_LIBGMP2
165  mmDBTestBlock(a->z._mp_d,a->z._mp_alloc*BYTES_PER_MP_LIMB,f,l);
166  if (a->z._mp_alloc==0)
167#else
168  mmDBTestBlock(a->z.d,a->z.alloc*BYTES_PER_MP_LIMB,f,l);
169  if(a->z.alloc==0)
170#endif
171    Print("!!longrat:z->alloc=0 in %s:%l\n",f,l);
172#endif
173  if (a->s<2)
174  {
175#ifdef MDEBUG
176#ifdef HAVE_LIBGMP2
177    mmDBTestBlock(a->n._mp_d,a->n._mp_alloc*BYTES_PER_MP_LIMB,f,-l);
178    if (a->z._mp_alloc==0)
179#else
180    mmDBTestBlock(a->n.d,a->n.alloc*BYTES_PER_MP_LIMB,f,-l);
181    if(a->z.alloc==0)
182#endif
183      Print("!!longrat:n->alloc=0 in %s:%l\n",f,l);
184#endif
185    if (mpz_cmp_si(&a->n,(long)1)==0)
186    {
187      Print("!!longrat:integer as rational in %s:%d\n",f,l);
188      mpz_clear(&a->n);
189      a->s=3;
190      return FALSE;
191    }
192    else if (mpz_isNeg(&a->n))
193    {
194      Print("!!longrat:div. by negative in %s:%d\n",f,l);
195      mpz_neg(&a->z,&a->z);
196      mpz_neg(&a->n,&a->n);
197      return FALSE;
198    }
199    return TRUE;
200  }
201  if (a->s==2)
202  {
203    Print("!!longrat:s=2 in %s:%d\n",f,l);
204    return FALSE;
205  }
206  if (mpz_size1(&a->z)>MP_SMALL) return TRUE;
207  int ui=(int)mpz_get_si(&a->z);
208  if ((((ui<<3)>>3)==ui)
209  && (mpz_cmp_si(&a->z,(long)ui)==0))
210  {
211    Print("!!longrat:im int %d in %s:%d\n",ui,f,l);
212    f=NULL;
213    return FALSE;
214  }
215  return TRUE;
216}
217#endif
218
219void nlNew (number * r)
220{
221  *r=NULL;
222}
223
224/*2
225* z := i
226*/
227#ifdef MDEBUG
228#define nlRInit(A) nlRDBInit(A,__LINE__)
229static inline number nlRDBInit(int i, int n)
230#else
231static inline number nlRInit (int i)
232#endif
233{
234#ifdef MDEBUG
235  number z=(number)mmDBAllocBlock(sizeof(rnumber),"longrat.cc",n);
236#else
237  number z=(number)Alloc(sizeof(rnumber));
238#endif
239#ifdef LDEBUG
240  z->debug=123456;
241#endif
242  mpz_init_set_si(&z->z,(long)i);
243  z->s = 3;
244  return z;
245}
246
247number nlInit (int i)
248{
249  number n;
250  if ( ((i << 3) >> 3) == i ) n=INT_TO_SR(i);
251  else                        n=nlRInit(i);
252#ifdef LDEBUG
253  nlTest(n);
254#endif
255  return n;
256}
257
258
259number nlInit (number u)
260{
261  if (u->s == 3 && mpz_size1(&u->z)<=MP_SMALL)
262  {
263    int ui=(int)mpz_get_si(&u->z);
264    if ((((ui<<3)>>3)==ui)
265        && (mpz_cmp_si(&u->z,(long)ui)==0))
266    {
267      mpz_clear(&u->z);
268      Free((ADDRESS)u,sizeof(rnumber));
269      return INT_TO_SR(ui);
270    }
271  }
272  return u;
273}
274
275int nlSize(number a)
276{
277  if (a==INT_TO_SR(0))
278     return 0; /* rational 0*/
279  if (SR_HDL(a) & SR_INT)
280     return 1; /* immidiate int */
281#ifdef HAVE_LIBGMP2
282  int s=a->z._mp_alloc;
283#else
284  int s=a->z.alloc;
285#endif
286  if (a->s<2)
287  {
288#ifdef HAVE_LIBGMP2
289    s+=a->n._mp_alloc;
290#else
291    s+=a->n.alloc;
292#endif
293  }
294  return s;
295}
296
297/*
298* delete leading 0s from a gmp number
299*/
300#ifdef HAVE_LIBGMP1
301static void nlGmpSimple(MP_INT *z)
302{
303  int k;
304  if (mpz_limb_size(z)<0)
305  {
306    k=-mpz_limb_size(z)-1;
307    while ((k>0) && (mpz_limb_d(z)[k]==0))
308    { k--;printf("X"); }
309    k++;
310    mpz_limb_size(z)=-k;
311  }
312  else
313  {
314    k=mpz_limb_size(z)-1;
315    while ((k>0) && (mpz_limb_d(z)[k]==0))
316    { k--;printf("X"); }
317    k++;
318    mpz_limb_size(z)=k;
319  }
320}
321#endif
322
323/*2
324* convert number to int
325*/
326int nlInt(number &i)
327{
328#ifdef LDEBUG
329  nlTest(i);
330#endif
331  nlNormalize(i);
332  if (SR_HDL(i) &SR_INT) return SR_TO_INT(i);
333  nlGmpSimple(&i->z);
334  if ((i->s!=3)||(mpz_size1(&i->z)>MP_SMALL)) return 0;
335  int ul=(int)mpz_get_si(&i->z);
336  if (mpz_cmp_si(&i->z,(long)ul)!=0) return 0;
337  return ul;
338}
339
340/*2
341* delete a
342*/
343#ifdef LDEBUG
344void nlDBDelete (number * a,char *f, int l)
345#else
346void nlDelete (number * a)
347#endif
348{
349  if (*a!=NULL)
350  {
351#ifdef LDEBUG
352    nlTest(*a);
353#endif
354    if ((SR_HDL(*a) & SR_INT)==0)
355    {
356      switch ((*a)->s)
357      {
358        case 0:
359        case 1:
360          mpz_clear(&(*a)->n);
361        case 3:
362          mpz_clear(&(*a)->z);
363#ifdef LDEBUG
364          (*a)->s=2;
365          (*a)->debug=654321;
366#endif
367      }
368#if defined(MDEBUG) && defined(LDEBUG)
369      mmDBFreeBlock((ADDRESS) *a,sizeof(rnumber),f,l);
370#else
371      Free((ADDRESS) *a,sizeof(rnumber));
372#endif
373    }
374    *a=NULL;
375  }
376}
377
378/*2
379* copy a to b
380*/
381number nlCopy(number a)
382{
383  number b;
384  if ((SR_HDL(a) & SR_INT)||(a==NULL))
385  {
386    return a;
387  }
388#ifdef LDEBUG
389  nlTest(a);
390#endif
391  b=(number)Alloc(sizeof(rnumber));
392#ifdef LDEBUG
393  b->debug=123456;
394#endif
395  switch (a->s)
396  {
397    case 0:
398    case 1:
399            nlGmpSimple(&a->n);
400            mpz_init_set(&b->n,&a->n);
401    case 3:
402            nlGmpSimple(&a->z);
403            mpz_init_set(&b->z,&a->z);
404            break;
405  }
406  b->s = a->s;
407#ifdef LDEBUG
408  nlTest(b);
409#endif
410  return b;
411}
412
413/*2
414* za:= - za
415*/
416number nlNeg (number a)
417{
418#ifdef LDEBUG
419  nlTest(a);
420#endif
421  if(SR_HDL(a) &SR_INT)
422  {
423    int r=SR_TO_INT(a);
424    if (r==(-(1<<28))) a=nlRInit(1<<28);
425    else               a=INT_TO_SR(-r);
426  }
427  else
428  {
429    mpz_neg(&a->z,&a->z);
430    if ((a->s==3) && (mpz_size1(&a->z)<=MP_SMALL))
431    {
432      int ui=(int)mpz_get_si(&a->z);
433      if ((((ui<<3)>>3)==ui)
434      && (mpz_cmp_si(&a->z,(long)ui)==0))
435      {
436        mpz_clear(&a->z);
437        Free((ADDRESS)a,sizeof(rnumber));
438        a=INT_TO_SR(ui);
439      }
440    }
441  }
442#ifdef LDEBUG
443  nlTest(a);
444#endif
445  return a;
446}
447
448/*
449* 1/a
450*/
451number nlInvers(number a)
452{
453#ifdef LDEBUG
454  nlTest(a);
455#endif
456  number n;
457  if (SR_HDL(a) & SR_INT)
458  {
459    if ((a==INT_TO_SR(1)) || (a==INT_TO_SR(-1)))
460    {
461      return a;
462    }
463    if (nlIsZero(a))
464    {
465      WerrorS("div. 1/0");
466      return INT_TO_SR(0);
467    }
468    n=(number)Alloc(sizeof(rnumber));
469#ifdef LDEBUG
470    n->debug=123456;
471#endif
472    n->s=1;
473    if ((int)a>0)
474    {
475      mpz_init_set_si(&n->z,(long)1);
476      mpz_init_set_si(&n->n,(long)SR_TO_INT(a));
477    }
478    else
479    {
480      mpz_init_set_si(&n->z,(long)-1);
481      mpz_init_set_si(&n->n,(long)-SR_TO_INT(a));
482#ifdef LDEBUG
483    nlTest(n);
484#endif
485    }
486#ifdef LDEBUG
487    nlTest(n);
488#endif
489    return n;
490  }
491  n=(number)Alloc(sizeof(rnumber));
492#ifdef LDEBUG
493  n->debug=123456;
494#endif
495  {
496    n->s=a->s;
497    mpz_init_set(&n->n,&a->z);
498    switch (a->s)
499    {
500      case 0:
501      case 1:
502              mpz_init_set(&n->z,&a->n);
503              if (mpz_isNeg(&n->n)) /* && n->s<2*/
504              {
505                mpz_neg(&n->z,&n->z);
506                mpz_neg(&n->n,&n->n);
507              }
508              if (mpz_cmp_si(&n->n,(long)1)==0)
509              {
510                mpz_clear(&n->n);
511                n->s=3;
512                if (mpz_size1(&n->z)<=MP_SMALL)
513                {
514                  int ui=(int)mpz_get_si(&n->z);
515                  if ((((ui<<3)>>3)==ui)
516                  && (mpz_cmp_si(&n->z,(long)ui)==0))
517                  {
518                    mpz_clear(&n->z);
519                    Free((ADDRESS)n,sizeof(rnumber));
520                    n=INT_TO_SR(ui);
521                  }
522                }
523              }
524              break;
525      case 3:
526              n->s=1;
527              if (mpz_isNeg(&n->n)) /* && n->s<2*/
528              {
529                mpz_neg(&n->n,&n->n);
530                mpz_init_set_si(&n->z,(long)-1);
531              }
532              else
533              {
534                mpz_init_set_si(&n->z,(long)1);
535              }
536              break;
537    }
538  }
539#ifdef LDEBUG
540  nlTest(n);
541#endif
542  return n;
543}
544
545/*2
546* u:= a + b
547*/
548number nlAdd (number a, number b)
549{
550  number u;
551  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
552  {
553    int r=SR_HDL(a)+SR_HDL(b)-1;
554    if ( ((r << 1) >> 1) == r )
555    {
556      return (number)r;
557    }
558    u=(number)Alloc(sizeof(rnumber));
559    u->s=3;
560    mpz_init_set_si(&u->z,(long)SR_TO_INT(r));
561#ifdef LDEBUG
562    u->debug=123456;
563    nlTest(u);
564#endif
565    return u;
566  }
567  u=(number)Alloc(sizeof(rnumber));
568#ifdef LDEBUG
569  u->debug=123456;
570#endif
571  mpz_init(&u->z);
572  if (SR_HDL(b) & SR_INT)
573  {
574    number x=a;
575    a=b;
576    b=x;
577  }
578  if (SR_HDL(a) & SR_INT)
579  {
580    switch (b->s)
581    {
582      case 0:
583      case 1:/* a:short, b:1 */
584      {
585        MP_INT x;
586        mpz_init(&x);
587        if ((int)a>0)
588        {
589          mpz_mul_ui(&x,&b->n,SR_TO_INT(a));
590        }
591        else
592        {
593          mpz_mul_ui(&x,&b->n,-SR_TO_INT(a));
594          mpz_neg(&x,&x);
595        }
596        mpz_add(&u->z,&b->z,&x);
597        nlGmpSimple(&u->z);
598        mpz_clear(&x);
599        if (mpz_cmp_ui(&u->z,(long)0)==0)
600        {
601          mpz_clear(&u->z);
602          Free((ADDRESS)u,sizeof(rnumber));
603          return INT_TO_SR(0);
604        }
605        if (mpz_cmp(&u->z,&b->n)==0)
606        {
607          mpz_clear(&u->z);
608          Free((ADDRESS)u,sizeof(rnumber));
609          return INT_TO_SR(1);
610        }
611        mpz_init_set(&u->n,&b->n);
612        u->s = 0;
613        break;
614      }
615      case 3:
616      {
617        if ((int)a>0)
618          mpz_add_ui(&u->z,&b->z,SR_TO_INT(a));
619        else
620          mpz_sub_ui(&u->z,&b->z,-SR_TO_INT(a));
621        nlGmpSimple(&u->z);
622        if (mpz_cmp_ui(&u->z,(long)0)==0)
623        {
624          mpz_clear(&u->z);
625          Free((ADDRESS)u,sizeof(rnumber));
626          return INT_TO_SR(0);
627        }
628        //u->s = 3;
629        if (mpz_size1(&u->z)<=MP_SMALL)
630        {
631          int ui=(int)mpz_get_si(&u->z);
632          if ((((ui<<3)>>3)==ui)
633          && (mpz_cmp_si(&u->z,(long)ui)==0))
634          {
635            mpz_clear(&u->z);
636            Free((ADDRESS)u,sizeof(rnumber));
637            return INT_TO_SR(ui);
638          }
639        }
640        u->s = 3;
641        break;
642      }
643    }
644  }
645  else
646  {
647    switch (a->s)
648    {
649      case 0:
650      case 1:
651      {
652        switch(b->s)
653        {
654          case 0:
655          case 1:
656          {
657            MP_INT x;
658            MP_INT y;
659            mpz_init(&x);
660            mpz_init(&y);
661            mpz_mul(&x,&b->z,&a->n);
662            mpz_mul(&y,&a->z,&b->n);
663            mpz_add(&u->z,&x,&y);
664            nlGmpSimple(&u->z);
665            mpz_clear(&x);
666            mpz_clear(&y);
667            if (mpz_cmp_ui(&u->z,(long)0)==0)
668            {
669              mpz_clear(&u->z);
670              Free((ADDRESS)u,sizeof(rnumber));
671              return INT_TO_SR(0);
672            }
673            mpz_init(&u->n);
674            mpz_mul(&u->n,&a->n,&b->n);
675            nlGmpSimple(&u->n);
676            if (mpz_cmp(&u->z,&u->n)==0)
677            {
678               mpz_clear(&u->z);
679               mpz_clear(&u->n);
680               Free((ADDRESS)u,sizeof(rnumber));
681               return INT_TO_SR(1);
682            }
683            u->s = 0;
684            break;
685          }
686          case 3: /* a:1 b:3 */
687          {
688            MP_INT x;
689            mpz_init(&x);
690            mpz_mul(&x,&b->z,&a->n);
691            mpz_add(&u->z,&a->z,&x);
692            nlGmpSimple(&u->z);
693            mpz_clear(&x);
694            if (mpz_cmp_ui(&u->z,(long)0)==0)
695            {
696              mpz_clear(&u->z);
697              Free((ADDRESS)u,sizeof(rnumber));
698              return INT_TO_SR(0);
699            }
700            if (mpz_cmp(&u->z,&a->n)==0)
701            {
702              mpz_clear(&u->z);
703              Free((ADDRESS)u,sizeof(rnumber));
704              return INT_TO_SR(1);
705            }
706            mpz_init_set(&u->n,&a->n);
707            u->s = 0;
708            break;
709          }
710        } /*switch (b->s) */
711        break;
712      }
713      case 3:
714      {
715        switch(b->s)
716        {
717          case 0:
718          case 1:/* a:3, b:1 */
719          {
720            MP_INT x;
721            mpz_init(&x);
722            mpz_mul(&x,&a->z,&b->n);
723            mpz_add(&u->z,&b->z,&x);
724            nlGmpSimple(&u->z);
725            mpz_clear(&x);
726            if (mpz_cmp_ui(&u->z,(long)0)==0)
727            {
728              mpz_clear(&u->z);
729              Free((ADDRESS)u,sizeof(rnumber));
730              return INT_TO_SR(0);
731            }
732            if (mpz_cmp(&u->z,&b->n)==0)
733            {
734              mpz_clear(&u->z);
735              Free((ADDRESS)u,sizeof(rnumber));
736              return INT_TO_SR(1);
737            }
738            mpz_init_set(&u->n,&b->n);
739            u->s = 0;
740            break;
741          }
742          case 3:
743          {
744            mpz_add(&u->z,&a->z,&b->z);
745            nlGmpSimple(&u->z);
746            if (mpz_cmp_ui(&u->z,(long)0)==0)
747            {
748              mpz_clear(&u->z);
749              Free((ADDRESS)u,sizeof(rnumber));
750              return INT_TO_SR(0);
751            }
752            if (mpz_size1(&u->z)<=MP_SMALL)
753            {
754              int ui=(int)mpz_get_si(&u->z);
755              if ((((ui<<3)>>3)==ui)
756              && (mpz_cmp_si(&u->z,(long)ui)==0))
757              {
758                mpz_clear(&u->z);
759                Free((ADDRESS)u,sizeof(rnumber));
760                return INT_TO_SR(ui);
761              }
762            }
763            u->s = 3;
764            break;
765          }
766        }
767        break;
768      }
769    }
770  }
771#ifdef LDEBUG
772  nlTest(u);
773#endif
774  return u;
775}
776
777/*2
778* u:= a - b
779*/
780number nlSub (number a, number b)
781{
782  number u;
783  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
784  {
785    int r=SR_HDL(a)-SR_HDL(b)+1;
786    if ( ((r << 1) >> 1) == r )
787    {
788      return (number)r;
789    }
790    u=(number)Alloc(sizeof(rnumber));
791    u->s=3;
792    mpz_init_set_si(&u->z,(long)SR_TO_INT(r));
793#ifdef LDEBUG
794    u->debug=123456;
795    nlTest(u);
796#endif
797    return u;
798  }
799  u=(number)Alloc(sizeof(rnumber));
800#ifdef LDEBUG
801  u->debug=123456;
802#endif
803  mpz_init(&u->z);
804  if (SR_HDL(a) & SR_INT)
805  {
806    switch (b->s)
807    {
808      case 0:
809      case 1:/* a:short, b:1 */
810      {
811        MP_INT x;
812        mpz_init(&x);
813        if ((int)a>0)
814        {
815          mpz_mul_ui(&x,&b->n,SR_TO_INT(a));
816        }
817        else
818        {
819          mpz_mul_ui(&x,&b->n,-SR_TO_INT(a));
820          mpz_neg(&x,&x);
821        }
822        mpz_sub(&u->z,&x,&b->z);
823        mpz_clear(&x);
824        if (mpz_cmp_ui(&u->z,(long)0)==0)
825        {
826          mpz_clear(&u->z);
827          Free((ADDRESS)u,sizeof(rnumber));
828          return INT_TO_SR(0);
829        }
830        if (mpz_cmp(&u->z,&b->n)==0)
831        {
832          mpz_clear(&u->z);
833          Free((ADDRESS)u,sizeof(rnumber));
834          return INT_TO_SR(1);
835        }
836        mpz_init_set(&u->n,&b->n);
837        u->s = 0;
838        break;
839      }
840      case 3:
841      {
842        if ((int)a>0)
843        {
844          mpz_sub_ui(&u->z,&b->z,SR_TO_INT(a));
845          mpz_neg(&u->z,&u->z);
846        }
847        else
848        {
849          mpz_add_ui(&u->z,&b->z,-SR_TO_INT(a));
850          mpz_neg(&u->z,&u->z);
851        }
852        nlGmpSimple(&u->z);
853        if (mpz_cmp_ui(&u->z,(long)0)==0)
854        {
855          mpz_clear(&u->z);
856          Free((ADDRESS)u,sizeof(rnumber));
857          return INT_TO_SR(0);
858        }
859        //u->s = 3;
860        if (mpz_size1(&u->z)<=MP_SMALL)
861        {
862          int ui=(int)mpz_get_si(&u->z);
863          if ((((ui<<3)>>3)==ui)
864          && (mpz_cmp_si(&u->z,(long)ui)==0))
865          {
866            mpz_clear(&u->z);
867            Free((ADDRESS)u,sizeof(rnumber));
868            return INT_TO_SR(ui);
869          }
870        }
871        u->s = 3;
872        break;
873      }
874    }
875  }
876  else if (SR_HDL(b) & SR_INT)
877  {
878    switch (a->s)
879    {
880      case 0:
881      case 1:/* b:short, a:1 */
882      {
883        MP_INT x;
884        mpz_init(&x);
885        if ((int)b>0)
886        {
887          mpz_mul_ui(&x,&a->n,SR_TO_INT(b));
888        }
889        else
890        {
891          mpz_mul_ui(&x,&a->n,-SR_TO_INT(b));
892          mpz_neg(&x,&x);
893        }
894        mpz_sub(&u->z,&a->z,&x);
895        mpz_clear(&x);
896        if (mpz_cmp_ui(&u->z,(long)0)==0)
897        {
898          mpz_clear(&u->z);
899          Free((ADDRESS)u,sizeof(rnumber));
900          return INT_TO_SR(0);
901        }
902        if (mpz_cmp(&u->z,&a->n)==0)
903        {
904          mpz_clear(&u->z);
905          Free((ADDRESS)u,sizeof(rnumber));
906          return INT_TO_SR(1);
907        }
908        mpz_init_set(&u->n,&a->n);
909        u->s = 0;
910        break;
911      }
912      case 3:
913      {
914        if ((int)b>0)
915        {
916          mpz_sub_ui(&u->z,&a->z,SR_TO_INT(b));
917        }
918        else
919        {
920          mpz_add_ui(&u->z,&a->z,-SR_TO_INT(b));
921        }
922        if (mpz_cmp_ui(&u->z,(long)0)==0)
923        {
924          mpz_clear(&u->z);
925          Free((ADDRESS)u,sizeof(rnumber));
926          return INT_TO_SR(0);
927        }
928        //u->s = 3;
929        nlGmpSimple(&u->z);
930        if (mpz_size1(&u->z)<=MP_SMALL)
931        {
932          int ui=(int)mpz_get_si(&u->z);
933          if ((((ui<<3)>>3)==ui)
934          && (mpz_cmp_si(&u->z,(long)ui)==0))
935          {
936            mpz_clear(&u->z);
937            Free((ADDRESS)u,sizeof(rnumber));
938            return INT_TO_SR(ui);
939          }
940        }
941        u->s = 3;
942        break;
943      }
944    }
945  }
946  else
947  {
948    switch (a->s)
949    {
950      case 0:
951      case 1:
952      {
953        switch(b->s)
954        {
955          case 0:
956          case 1:
957          {
958            MP_INT x;
959            MP_INT y;
960            mpz_init(&x);
961            mpz_init(&y);
962            mpz_mul(&x,&b->z,&a->n);
963            mpz_mul(&y,&a->z,&b->n);
964            mpz_sub(&u->z,&y,&x);
965            mpz_clear(&x);
966            mpz_clear(&y);
967            if (mpz_cmp_ui(&u->z,(long)0)==0)
968            {
969              mpz_clear(&u->z);
970              Free((ADDRESS)u,sizeof(rnumber));
971              return INT_TO_SR(0);
972            }
973            mpz_init(&u->n);
974            mpz_mul(&u->n,&a->n,&b->n);
975            if (mpz_cmp(&u->z,&u->n)==0)
976            {
977              mpz_clear(&u->z);
978              mpz_clear(&u->n);
979              Free((ADDRESS)u,sizeof(rnumber));
980              return INT_TO_SR(1);
981            }
982            u->s = 0;
983            break;
984          }
985          case 3: /* a:1, b:3 */
986          {
987            MP_INT x;
988            mpz_init(&x);
989            mpz_mul(&x,&b->z,&a->n);
990            mpz_sub(&u->z,&a->z,&x);
991            mpz_clear(&x);
992            if (mpz_cmp_ui(&u->z,(long)0)==0)
993            {
994              mpz_clear(&u->z);
995              Free((ADDRESS)u,sizeof(rnumber));
996              return INT_TO_SR(0);
997            }
998            if (mpz_cmp(&u->z,&a->n)==0)
999            {
1000              mpz_clear(&u->z);
1001              Free((ADDRESS)u,sizeof(rnumber));
1002              return INT_TO_SR(1);
1003            }
1004            mpz_init_set(&u->n,&a->n);
1005            u->s = 0;
1006            break;
1007          }
1008        }
1009        break;
1010      }
1011      case 3:
1012      {
1013        switch(b->s)
1014        {
1015          case 0:
1016          case 1: /* a:3, b:1 */
1017          {
1018            MP_INT x;
1019            mpz_init(&x);
1020            mpz_mul(&x,&a->z,&b->n);
1021            mpz_sub(&u->z,&x,&b->z);
1022            mpz_clear(&x);
1023            if (mpz_cmp_ui(&u->z,(long)0)==0)
1024            {
1025              mpz_clear(&u->z);
1026              Free((ADDRESS)u,sizeof(rnumber));
1027              return INT_TO_SR(0);
1028            }
1029            if (mpz_cmp(&u->z,&b->n)==0)
1030            {
1031              mpz_clear(&u->z);
1032              Free((ADDRESS)u,sizeof(rnumber));
1033              return INT_TO_SR(1);
1034            }
1035            mpz_init_set(&u->n,&b->n);
1036            u->s = 0;
1037            break;
1038          }
1039          case 3: /* a:3 , b:3 */
1040          {
1041            mpz_sub(&u->z,&a->z,&b->z);
1042            nlGmpSimple(&u->z);
1043            if (mpz_cmp_ui(&u->z,(long)0)==0)
1044            {
1045              mpz_clear(&u->z);
1046              Free((ADDRESS)u,sizeof(rnumber));
1047              return INT_TO_SR(0);
1048            }
1049            //u->s = 3;
1050            if (mpz_size1(&u->z)<=MP_SMALL)
1051            {
1052              int ui=(int)mpz_get_si(&u->z);
1053              if ((((ui<<3)>>3)==ui)
1054              && (mpz_cmp_si(&u->z,(long)ui)==0))
1055              {
1056                mpz_clear(&u->z);
1057                Free((ADDRESS)u,sizeof(rnumber));
1058                return INT_TO_SR(ui);
1059              }
1060            }
1061            u->s = 3;
1062            break;
1063          }
1064        }
1065        break;
1066      }
1067    }
1068  }
1069#ifdef LDEBUG
1070  nlTest(u);
1071#endif
1072  return u;
1073}
1074
1075/*2
1076* u := a * b
1077*/
1078number nlMult (number a, number b)
1079{
1080  number u;
1081
1082#ifdef LDEBUG
1083  nlTest(a);
1084  nlTest(b);
1085#endif
1086  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1087  {
1088    if (b==INT_TO_SR(0)) return INT_TO_SR(0);
1089    int r=(SR_HDL(a)-1)*(SR_HDL(b)>>1);
1090    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1))
1091    {
1092      u=((number) ((r>>1)+SR_INT));
1093      if (((SR_HDL(u)<<1)>>1)==SR_HDL(u)) return (u);
1094      return nlRInit(SR_HDL(u)>>2);
1095    }
1096    u=(number)Alloc(sizeof(rnumber));
1097#ifdef LDEBUG
1098    u->debug=123456;
1099#endif
1100    u->s=3;
1101    if ((int)b>0)
1102    {
1103      mpz_init_set_si(&u->z,(long)SR_TO_INT(a));
1104      mpz_mul_ui(&u->z,&u->z,(unsigned long)SR_TO_INT(b));
1105    }
1106//    else if ((int)a>=0)
1107//    {
1108//      mpz_init_set_si(&u->z,(long)SR_TO_INT(b));
1109//      mpz_mul_ui(&u->z,&u->z,(unsigned long)SR_TO_INT(a));
1110//    }
1111    else
1112    {
1113      mpz_init_set_si(&u->z,(long)(-SR_TO_INT(a)));
1114      mpz_mul_ui(&u->z,&u->z,(long)(-SR_TO_INT(b)));
1115    }
1116#ifdef LDEBUG
1117    nlTest(u);
1118#endif
1119    return u;
1120  }
1121  else
1122  {
1123    if ((a==INT_TO_SR(0))
1124    ||(b==INT_TO_SR(0)))
1125      return INT_TO_SR(0);
1126    u=(number)Alloc(sizeof(rnumber));
1127#ifdef LDEBUG
1128    u->debug=123456;
1129#endif
1130    mpz_init(&u->z);
1131    if (SR_HDL(b) & SR_INT)
1132    {
1133      number x=a;
1134      a=b;
1135      b=x;
1136    }
1137    if (SR_HDL(a) & SR_INT)
1138    {
1139      u->s=b->s;
1140      if (u->s==1) u->s=0;
1141      if ((int)a>0)
1142      {
1143        mpz_mul_ui(&u->z,&b->z,(unsigned long)SR_TO_INT(a));
1144      }
1145      else
1146      {
1147        if (a==INT_TO_SR(-1))
1148        {
1149          mpz_set(&u->z,&b->z);
1150          mpz_neg(&u->z,&u->z);
1151          u->s=b->s;
1152        }
1153        else
1154        {
1155          mpz_mul_ui(&u->z,&b->z,(unsigned long)-SR_TO_INT(a));
1156          mpz_neg(&u->z,&u->z);
1157        }
1158      }
1159      nlGmpSimple(&u->z);
1160      if (u->s<2)
1161      {
1162        if (mpz_cmp(&u->z,&b->n)==0)
1163        {
1164          mpz_clear(&u->z);
1165          Free((ADDRESS)u,sizeof(rnumber));
1166          return INT_TO_SR(1);
1167        }
1168        mpz_init_set(&u->n,&b->n);
1169      }
1170      else //u->s==3
1171      {
1172        if (mpz_size1(&u->z)<=MP_SMALL)
1173        {
1174          int ui=(int)mpz_get_si(&u->z);
1175          if ((((ui<<3)>>3)==ui)
1176          && (mpz_cmp_si(&u->z,(long)ui)==0))
1177          {
1178            mpz_clear(&u->z);
1179            Free((ADDRESS)u,sizeof(rnumber));
1180            return INT_TO_SR(ui);
1181          }
1182        }
1183      }
1184    }
1185    else
1186    {
1187      mpz_mul(&u->z,&a->z,&b->z);
1188      u->s = 0;
1189      if(a->s==3)
1190      {
1191        if(b->s==3)
1192        {
1193          u->s = 3;
1194        }
1195        else
1196        {
1197          if (mpz_cmp(&u->z,&b->n)==0)
1198          {
1199            mpz_clear(&u->z);
1200            Free((ADDRESS)u,sizeof(rnumber));
1201            return INT_TO_SR(1);
1202          }
1203          mpz_init_set(&u->n,&b->n);
1204        }
1205      }
1206      else
1207      {
1208        if(b->s==3)
1209        {
1210          if (mpz_cmp(&u->z,&a->n)==0)
1211          {
1212            mpz_clear(&u->z);
1213            Free((ADDRESS)u,sizeof(rnumber));
1214            return INT_TO_SR(1);
1215          }
1216          mpz_init_set(&u->n,&a->n);
1217        }
1218        else
1219        {
1220          mpz_init(&u->n);
1221          mpz_mul(&u->n,&a->n,&b->n);
1222          if (mpz_cmp(&u->z,&u->n)==0)
1223          {
1224            mpz_clear(&u->z);
1225            mpz_clear(&u->n);
1226            Free((ADDRESS)u,sizeof(rnumber));
1227            return INT_TO_SR(1);
1228          }
1229        }
1230      }
1231    }
1232  }
1233#ifdef LDEBUG
1234  nlTest(u);
1235#endif
1236  return u;
1237}
1238
1239/*2
1240* u := a / b in Z, if b | a (else undefined)
1241*/
1242number   nlExactDiv(number a, number b)
1243{
1244  if (b==INT_TO_SR(0))
1245  {
1246    WerrorS("div. by 0");
1247    return INT_TO_SR(0);
1248  }
1249  number u;
1250  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1251  {
1252    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
1253    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
1254    {
1255      return nlRInit(1<<28);
1256    }
1257    int aa=SR_TO_INT(a);
1258    int bb=SR_TO_INT(b);
1259    return INT_TO_SR(aa/bb);
1260  }
1261  number bb=NULL;
1262  if (SR_HDL(b) & SR_INT)
1263  {
1264    bb=nlRInit(SR_TO_INT(b));
1265    b=bb;
1266  }
1267  u=(number)Alloc(sizeof(rnumber));
1268#ifdef LDEBUG
1269  u->debug=123456;
1270#endif
1271  mpz_init(&u->z);
1272  /* u=a/b */
1273  u->s = 3;
1274  MPZ_EXACTDIV(&u->z,&a->z,&b->z);
1275  if (bb!=NULL)
1276  {
1277    mpz_clear(&bb->z);
1278    Free((ADDRESS)bb,sizeof(rnumber));
1279  }
1280  if (mpz_size1(&u->z)<=MP_SMALL)
1281  {
1282    int ui=(int)mpz_get_si(&u->z);
1283    if ((((ui<<3)>>3)==ui)
1284    && (mpz_cmp_si(&u->z,(long)ui)==0))
1285    {
1286      mpz_clear(&u->z);
1287      Free((ADDRESS)u,sizeof(rnumber));
1288      u=INT_TO_SR(ui);
1289    }
1290  }
1291#ifdef LDEBUG
1292  nlTest(u);
1293#endif
1294  return u;
1295}
1296
1297/*2
1298* u := a / b in Z
1299*/
1300number nlIntDiv (number a, number b)
1301{
1302  if (b==INT_TO_SR(0))
1303  {
1304    WerrorS("div. by 0");
1305    return INT_TO_SR(0);
1306  }
1307  number u;
1308  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1309  {
1310    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
1311    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
1312    {
1313      return nlRInit(1<<28);
1314    }
1315    /* consider the signs of a and b:
1316    *  + + -> u=a+b-1
1317    *  + - -> u=a-b-1
1318    *  - + -> u=a-b+1
1319    *  - - -> u=a+b+1
1320    */
1321    int aa=SR_TO_INT(a);
1322    int bb=SR_TO_INT(b);
1323    if (aa>0)
1324    {
1325      if (bb>0)
1326        return INT_TO_SR((aa+bb-1)/bb);
1327      else
1328        return INT_TO_SR((aa-bb-1)/bb);
1329    }
1330    else
1331    {
1332      if (bb>0)
1333      {
1334        return INT_TO_SR((aa-bb+1)/bb);
1335      }
1336      else
1337      {
1338        return INT_TO_SR((aa+bb+1)/bb);
1339      }
1340    }
1341  }
1342  if (SR_HDL(a) & SR_INT)
1343  {
1344    /* the small int -(1<<28) divided by 2^28 is 1   */
1345    if (a==INT_TO_SR(-(1<<28)))
1346    {
1347      if(mpz_cmp_si(&b->z,(long)(1<<28))==0)
1348      {
1349        return INT_TO_SR(-1);
1350      }
1351    }
1352    /* a is a small and b is a large int: -> 0 */
1353    return INT_TO_SR(0);
1354  }
1355  number bb=NULL;
1356  if (SR_HDL(b) & SR_INT)
1357  {
1358    bb=nlRInit(SR_TO_INT(b));
1359    b=bb;
1360  }
1361  u=(number)Alloc(sizeof(rnumber));
1362#ifdef LDEBUG
1363  u->debug=123456;
1364#endif
1365  mpz_init_set(&u->z,&a->z);
1366  /* consider the signs of a and b:
1367  *  + + -> u=a+b-1
1368  *  + - -> u=a-b-1
1369  *  - + -> u=a-b+1
1370  *  - - -> u=a+b+1
1371  */
1372  if (mpz_isNeg(&a->z))
1373  {
1374    if (mpz_isNeg(&b->z))
1375    {
1376      mpz_add(&u->z,&u->z,&b->z);
1377    }
1378    else
1379    {
1380      mpz_sub(&u->z,&u->z,&b->z);
1381    }
1382    mpz_add_ui(&u->z,&u->z,1);
1383  }
1384  else
1385  {
1386    if (mpz_isNeg(&b->z))
1387    {
1388      mpz_sub(&u->z,&u->z,&b->z);
1389    }
1390    else
1391    {
1392      mpz_add(&u->z,&u->z,&b->z);
1393    }
1394    mpz_sub_ui(&u->z,&u->z,1);
1395  }
1396  /* u=u/b */
1397  u->s = 3;
1398  MPZ_DIV(&u->z,&u->z,&b->z);
1399  nlGmpSimple(&u->z);
1400  if (bb!=NULL)
1401  {
1402    mpz_clear(&bb->z);
1403    Free((ADDRESS)bb,sizeof(rnumber));
1404  }
1405  if (mpz_size1(&u->z)<=MP_SMALL)
1406  {
1407    int ui=(int)mpz_get_si(&u->z);
1408    if ((((ui<<3)>>3)==ui)
1409    && (mpz_cmp_si(&u->z,(long)ui)==0))
1410    {
1411      mpz_clear(&u->z);
1412      Free((ADDRESS)u,sizeof(rnumber));
1413      u=INT_TO_SR(ui);
1414    }
1415  }
1416#ifdef LDEBUG
1417  nlTest(u);
1418#endif
1419  return u;
1420}
1421
1422/*2
1423* u := a mod b in Z, u>=0
1424*/
1425number nlIntMod (number a, number b)
1426{
1427  if (b==INT_TO_SR(0))
1428  {
1429    WerrorS("div. by 0");
1430    return INT_TO_SR(0);
1431  }
1432  if (a==INT_TO_SR(0))
1433      return INT_TO_SR(0);
1434  number u;
1435  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1436  {
1437    if ((int)a>0)
1438    {
1439      if ((int)b>0)
1440        return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
1441      else
1442        return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
1443    }
1444    else
1445    {
1446      if ((int)b>0)
1447      {
1448        int i=(-SR_TO_INT(a))%SR_TO_INT(b);
1449        if ( i != 0 ) i = (SR_TO_INT(b))-i;
1450        return INT_TO_SR(i);
1451      }
1452      else
1453      {
1454        int i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
1455        if ( i != 0 ) i = (-SR_TO_INT(b))-i;
1456        return INT_TO_SR(i);
1457      }
1458    }
1459  }
1460  if (SR_HDL(a) & SR_INT)
1461  {
1462    /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
1463    if ((int)a<0)
1464    {
1465      if (mpz_isNeg(&b->z))
1466        return nlSub(a,b);
1467      /*else*/
1468        return nlAdd(a,b);
1469    }
1470    /*else*/
1471      return a;
1472  }
1473  number bb=NULL;
1474  if (SR_HDL(b) & SR_INT)
1475  {
1476    bb=nlRInit(SR_TO_INT(b));
1477    b=bb;
1478  }
1479  u=(number)Alloc(sizeof(rnumber));
1480#ifdef LDEBUG
1481  u->debug=123456;
1482#endif
1483  mpz_init(&u->z);
1484  u->s = 3;
1485  mpz_mod(&u->z,&a->z,&b->z);
1486  if (bb!=NULL)
1487  {
1488    mpz_clear(&bb->z);
1489    Free((ADDRESS)bb,sizeof(rnumber));
1490  }
1491  if (mpz_isNeg(&u->z))
1492  {
1493    if (mpz_isNeg(&b->z))
1494      mpz_sub(&u->z,&u->z,&b->z);
1495    else
1496      mpz_add(&u->z,&u->z,&b->z);
1497  }
1498  nlGmpSimple(&u->z);
1499  if (mpz_size1(&u->z)<=MP_SMALL)
1500  {
1501    int ui=(int)mpz_get_si(&u->z);
1502    if ((((ui<<3)>>3)==ui)
1503    && (mpz_cmp_si(&u->z,(long)ui)==0))
1504    {
1505      mpz_clear(&u->z);
1506      Free((ADDRESS)u,sizeof(rnumber));
1507      u=INT_TO_SR(ui);
1508    }
1509  }
1510#ifdef LDEBUG
1511  nlTest(u);
1512#endif
1513  return u;
1514}
1515
1516/*2
1517* u := a / b
1518*/
1519number nlDiv (number a, number b)
1520{
1521  number u;
1522  if (nlIsZero(b))
1523  {
1524    WerrorS("div. by 0");
1525    return INT_TO_SR(0);
1526  }
1527  u=(number)Alloc(sizeof(rnumber));
1528  u->s=0;
1529#ifdef LDEBUG
1530  u->debug=123456;
1531#endif
1532// ---------- short / short ------------------------------------
1533  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1534  {
1535    int i=SR_TO_INT(a);
1536    int j=SR_TO_INT(b);
1537    int r=i%j;
1538    if (r==0)
1539    {
1540      Free((ADDRESS)u,sizeof(rnumber));
1541      return INT_TO_SR(i/j);
1542    }
1543    mpz_init_set_si(&u->z,(long)i);
1544    mpz_init_set_si(&u->n,(long)j);
1545  }
1546  else
1547  {
1548    mpz_init(&u->z);
1549// ---------- short / long ------------------------------------
1550    if (SR_HDL(a) & SR_INT)
1551    {
1552      // short a / (z/n) -> (a*n)/z
1553      if (b->s<2)
1554      {
1555        if ((int)a>=0)
1556          mpz_mul_ui(&u->z,&b->n,SR_TO_INT(a));
1557        else
1558        {
1559          mpz_mul_ui(&u->z,&b->n,-SR_TO_INT(a));
1560          mpz_neg(&u->z,&u->z);
1561        }
1562      }
1563      else
1564      // short a / long z -> a/z
1565      {
1566        mpz_set_si(&u->z,SR_TO_INT(a));
1567      }
1568      if (mpz_cmp(&u->z,&b->z)==0)
1569      {
1570        mpz_clear(&u->z);
1571        Free((ADDRESS)u,sizeof(rnumber));
1572        return INT_TO_SR(1);
1573      }
1574      mpz_init_set(&u->n,&b->z);
1575    }
1576// ---------- long / short ------------------------------------
1577    else if (SR_HDL(b) & SR_INT)
1578    {
1579      mpz_set(&u->z,&a->z);
1580      // (z/n) / b -> z/(n*b)
1581      if (a->s<2)
1582      {
1583        mpz_init_set(&u->n,&a->n);
1584        if ((int)b>0)
1585          mpz_mul_ui(&u->n,&u->n,SR_TO_INT(b));
1586        else
1587        {
1588          mpz_mul_ui(&u->n,&u->n,-SR_TO_INT(b));
1589          mpz_neg(&u->z,&u->z);
1590        }
1591      }
1592      else
1593      // long z / short b -> z/b
1594      {
1595        //mpz_set(&u->z,&a->z);
1596        mpz_init_set_si(&u->n,SR_TO_INT(b));
1597      }
1598    }
1599// ---------- long / long ------------------------------------
1600    else
1601    {
1602      //u->s=0;
1603      mpz_set(&u->z,&a->z);
1604      mpz_init_set(&u->n,&b->z);
1605      if (a->s<2) mpz_mul(&u->n,&u->n,&a->n);
1606      if (b->s<2) mpz_mul(&u->z,&u->z,&b->n);
1607    }
1608  }
1609  if (mpz_isNeg(&u->n))
1610  {
1611    mpz_neg(&u->z,&u->z);
1612    mpz_neg(&u->n,&u->n);
1613  }
1614  if (mpz_cmp_si(&u->n,(long)1)==0)
1615  {
1616    mpz_clear(&u->n);
1617    if (mpz_size1(&u->z)<=MP_SMALL)
1618    {
1619      int ui=(int)mpz_get_si(&u->z);
1620      if ((((ui<<3)>>3)==ui)
1621      && (mpz_cmp_si(&u->z,(long)ui)==0))
1622      {
1623        mpz_clear(&u->z);
1624        Free((ADDRESS)u,sizeof(rnumber));
1625        return INT_TO_SR(ui);
1626      }
1627    }
1628    u->s=3;
1629  }
1630#ifdef LDEBUG
1631  nlTest(u);
1632#endif
1633  return u;
1634}
1635
1636#if (defined(i386)) && (defined(HAVE_LIBGMP1))
1637/*
1638* compute x^exp for x in Z
1639* there seems to be an bug in mpz_pow_ui for i386
1640*/
1641static inline void nlPow (MP_INT * res,MP_INT * x,int exp)
1642{
1643  if (exp==0)
1644  {
1645    mpz_set_si(res,(long)1);
1646  }
1647  else
1648  {
1649    MP_INT xh;
1650    mpz_init(&xh);
1651    mpz_set(res,x);
1652    exp--;
1653    while (exp!=0)
1654    {
1655      mpz_mul(&xh,res,x);
1656      mpz_set(res,&xh);
1657      exp--;
1658    }
1659    mpz_clear(&xh);
1660  }
1661}
1662#endif
1663
1664/*2
1665* u:= x ^ exp
1666*/
1667void nlPower (number x,int exp,number * u)
1668{
1669  if (!nlIsZero(x))
1670  {
1671#ifdef LDEBUG
1672    nlTest(x);
1673#endif
1674    number aa=NULL;
1675    if (SR_HDL(x) & SR_INT)
1676    {
1677      aa=nlRInit(SR_TO_INT(x));
1678      x=aa;
1679    }
1680    *u=(number)Alloc(sizeof(rnumber));
1681#ifdef LDEBUG
1682    (*u)->debug=123456;
1683#endif
1684    mpz_init(&(*u)->z);
1685    (*u)->s = x->s;
1686#if (!defined(i386)) || (defined(HAVE_LIBGMP2))
1687    mpz_pow_ui(&(*u)->z,&x->z,(unsigned long)exp);
1688#else
1689    nlPow(&(*u)->z,&x->z,exp);
1690#endif
1691    if (x->s<2)
1692    {
1693      mpz_init(&(*u)->n);
1694#if (!defined(i386)) || (defined(HAVE_LIBGMP2))
1695      mpz_pow_ui(&(*u)->n,&x->n,(unsigned long)exp);
1696#else
1697      nlPow(&(*u)->n,&x->n,exp);
1698#endif
1699    }
1700    if (aa!=NULL)
1701    {
1702      mpz_clear(&aa->z);
1703      Free((ADDRESS)aa,sizeof(rnumber));
1704    }
1705    if (((*u)->s<=2) && (mpz_cmp_si(&(*u)->n,(long)1)==0))
1706    {
1707      mpz_clear(&(*u)->n);
1708      (*u)->s=3;
1709    }
1710    if (((*u)->s==3) && (mpz_size1(&(*u)->z)<=MP_SMALL))
1711    {
1712      int ui=(int)mpz_get_si(&(*u)->z);
1713      if ((((ui<<3)>>3)==ui)
1714      && (mpz_cmp_si(&(*u)->z,(long)ui)==0))
1715      {
1716        mpz_clear(&(*u)->z);
1717        Free((ADDRESS)*u,sizeof(rnumber));
1718        *u=INT_TO_SR(ui);
1719      }
1720    }
1721  }
1722  else
1723    *u = INT_TO_SR(0);
1724#ifdef LDEBUG
1725  nlTest(*u);
1726#endif
1727}
1728
1729BOOLEAN nlIsZero (number a)
1730{
1731#ifdef LDEBUG
1732  nlTest(a);
1733#endif
1734  //if (a==INT_TO_SR(0)) return TRUE;
1735  //if (SR_HDL(a) & SR_INT) return FALSE;
1736  //number aa=nlCopy(a);
1737  //nlNormalize(aa);
1738  return (a==INT_TO_SR(0));
1739}
1740
1741/*2
1742* za >= 0 ?
1743*/
1744BOOLEAN nlGreaterZero (number a)
1745{
1746#ifdef LDEBUG
1747  nlTest(a);
1748#endif
1749  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>=0;
1750  return (!mpz_isNeg(&a->z));
1751}
1752
1753/*2
1754* a > b ?
1755*/
1756BOOLEAN nlGreater (number a, number b)
1757{
1758#ifdef LDEBUG
1759  nlTest(a);
1760  nlTest(b);
1761#endif
1762  number r;
1763  BOOLEAN rr;
1764  r=nlSub(a,b);
1765  rr=(!nlIsZero(r)) && (nlGreaterZero(r));
1766  nlDelete(&r);
1767  return rr;
1768}
1769
1770/*2
1771* a = b ?
1772*/
1773BOOLEAN nlEqual (number a, number b)
1774{
1775#ifdef LDEBUG
1776  nlTest(a);
1777  nlTest(b);
1778#endif
1779// short - short
1780  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
1781//  long - short
1782  BOOLEAN bo;
1783  if (SR_HDL(b) & SR_INT)
1784  {
1785    if (a->s!=0) return FALSE;
1786    number n=b; b=a; a=n;
1787  }
1788//  short - long
1789  if (SR_HDL(a) & SR_INT)
1790  {
1791    if (b->s!=0)
1792      return FALSE;
1793    if (((int)a > 0) && (mpz_isNeg(&b->z)))
1794      return FALSE;
1795    if (((int)a < 0) && (!mpz_isNeg(&b->z)))
1796      return FALSE;
1797    MP_INT  bb;
1798    mpz_init_set(&bb,&b->n);
1799    if ((int)a<0)
1800    {
1801      mpz_neg(&bb,&bb);
1802      mpz_mul_ui(&bb,&bb,(long)-SR_TO_INT(a));
1803    }
1804    else
1805    {
1806      mpz_mul_ui(&bb,&bb,(long)SR_TO_INT(a));
1807    }
1808    bo=(mpz_cmp(&bb,&b->z)==0);
1809    mpz_clear(&bb);
1810    return bo;
1811  }
1812// long - long
1813  if (((a->s==1) && (b->s==3))
1814  ||  ((b->s==1) && (a->s==3)))
1815    return FALSE;
1816  if (mpz_isNeg(&a->z)&&(!mpz_isNeg(&b->z)))
1817    return FALSE;
1818  if (mpz_isNeg(&b->z)&&(!mpz_isNeg(&a->z)))
1819    return FALSE;
1820  nlGmpSimple(&a->z);
1821  nlGmpSimple(&b->z);
1822  if (a->s<2)
1823    nlGmpSimple(&a->n);
1824  if (b->s<2)
1825    nlGmpSimple(&b->n);
1826  MP_INT  aa;
1827  MP_INT  bb;
1828  mpz_init_set(&aa,&a->z);
1829  mpz_init_set(&bb,&b->z);
1830  if (a->s<2) mpz_mul(&bb,&bb,&a->n);
1831  if (b->s<2) mpz_mul(&aa,&aa,&b->n);
1832  bo=(mpz_cmp(&aa,&bb)==0);
1833  mpz_clear(&aa);
1834  mpz_clear(&bb);
1835  return bo;
1836}
1837
1838/*2
1839* a == 1 ?
1840*/
1841BOOLEAN nlIsOne (number a)
1842{
1843#ifdef LDEBUG
1844  if (a==NULL) return FALSE;
1845  nlTest(a);
1846#endif
1847  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(1));
1848  return FALSE;
1849}
1850
1851/*2
1852* a == -1 ?
1853*/
1854BOOLEAN nlIsMOne (number a)
1855{
1856#ifdef LDEBUG
1857  if (a==NULL) return FALSE;
1858  nlTest(a);
1859#endif
1860  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1));
1861  return FALSE;
1862}
1863
1864/*2
1865* result =gcd(a,b)
1866*/
1867number nlGcd(number a, number b)
1868{
1869  number result;
1870#ifdef LDEBUG
1871  nlTest(a);
1872  nlTest(b);
1873#endif
1874  //nlNormalize(a);
1875  //nlNormalize(b);
1876  if ((SR_HDL(a)==5)||(a==INT_TO_SR(-1))
1877  ||  (SR_HDL(b)==5)||(b==INT_TO_SR(-1))) return INT_TO_SR(1);
1878  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1879  {
1880    int i=SR_TO_INT(a);
1881    int j=SR_TO_INT(b);
1882    if((i==0)||(j==0))
1883      return INT_TO_SR(1);
1884    int l;
1885    i=ABS(i);
1886    j=ABS(j);
1887    do
1888    {
1889      l=i%j;
1890      i=j;
1891      j=l;
1892    } while (l!=0);
1893    return INT_TO_SR(i);
1894  }
1895  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1896  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1897  number aa=NULL;
1898  if (SR_HDL(a) & SR_INT)
1899  {
1900    aa=nlRInit(SR_TO_INT(a));
1901    a=aa;
1902  }
1903  else
1904  if (SR_HDL(b) & SR_INT)
1905  {
1906    aa=nlRInit(SR_TO_INT(b));
1907    b=aa;
1908  }
1909  result=(number)Alloc(sizeof(rnumber));
1910#ifdef LDEBUG
1911  result->debug=123456;
1912#endif
1913  mpz_init(&result->z);
1914  mpz_gcd(&result->z,&a->z,&b->z);
1915  nlGmpSimple(&result->z);
1916  result->s = 3;
1917  if (mpz_size1(&result->z)<=MP_SMALL)
1918  {
1919    int ui=(int)mpz_get_si(&result->z);
1920    if ((((ui<<3)>>3)==ui)
1921    && (mpz_cmp_si(&result->z,(long)ui)==0))
1922    {
1923      mpz_clear(&result->z);
1924      Free((ADDRESS)result,sizeof(rnumber));
1925      result=INT_TO_SR(ui);
1926    }
1927  }
1928  if (aa!=NULL)
1929  {
1930    mpz_clear(&aa->z);
1931    Free((ADDRESS)aa,sizeof(rnumber));
1932  }
1933#ifdef LDEBUG
1934  nlTest(result);
1935#endif
1936  return result;
1937}
1938
1939/*2
1940* simplify x
1941*/
1942void nlNormalize (number &x)
1943{
1944#ifdef LDEBUG
1945  nlTest(x);
1946#endif
1947  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1948    return;
1949  if (x->s==3)
1950  {
1951    if (mpz_cmp_ui(&x->z,(long)0)==0)
1952    {
1953      nlDelete(&x);
1954      x=nlInit(0);
1955      return;
1956    }
1957    if (mpz_size1(&x->z)<=MP_SMALL)
1958    {
1959      int ui=(int)mpz_get_si(&x->z);
1960      if ((((ui<<3)>>3)==ui)
1961      && (mpz_cmp_si(&x->z,(long)ui)==0))
1962      {
1963        mpz_clear(&x->z);
1964        Free((ADDRESS)x,sizeof(rnumber));
1965        x=INT_TO_SR(ui);
1966        return;
1967      }
1968    }
1969  }
1970  if (x->s!=1) nlGmpSimple(&x->z);
1971  if (x->s==0)
1972  {
1973    BOOLEAN divided=FALSE;
1974    MP_INT gcd;
1975    mpz_init(&gcd);
1976    mpz_gcd(&gcd,&x->z,&x->n);
1977    x->s=1;
1978    if (mpz_cmp_si(&gcd,(long)1)!=0)
1979    {
1980      MP_INT r;
1981      mpz_init(&r);
1982      MPZ_EXACTDIV(&r,&x->z,&gcd);
1983      mpz_set(&x->z,&r);
1984      MPZ_EXACTDIV(&r,&x->n,&gcd);
1985      mpz_set(&x->n,&r);
1986      mpz_clear(&r);
1987      nlGmpSimple(&x->z);
1988      divided=TRUE;
1989    }
1990    mpz_clear(&gcd);
1991    nlGmpSimple(&x->n);
1992    if (mpz_cmp_si(&x->n,(long)1)==0)
1993    {
1994      mpz_clear(&x->n);
1995      if (mpz_size1(&x->z)<=MP_SMALL)
1996      {
1997        int ui=(int)mpz_get_si(&x->z);
1998        if ((((ui<<3)>>3)==ui)
1999        && (mpz_cmp_si(&x->z,(long)ui)==0))
2000        {
2001          mpz_clear(&x->z);
2002#ifdef LDEBUG
2003          x->debug=654324;
2004#endif
2005          Free((ADDRESS)x,sizeof(rnumber));
2006          x=INT_TO_SR(ui);
2007          return;
2008        }
2009      }
2010      x->s=3;
2011    }
2012    else if (divided)
2013    {
2014      _mpz_realloc(&x->n,mpz_size1(&x->n));
2015    }
2016    if (divided) _mpz_realloc(&x->z,mpz_size1(&x->z));
2017  }
2018#ifdef LDEBUG
2019  nlTest(x);
2020#endif
2021}
2022
2023/*2
2024* returns in result->z the lcm(a->z,b->n)
2025*/
2026number nlLcm(number a, number b)
2027{
2028  number result;
2029#ifdef LDEBUG
2030  nlTest(a);
2031  nlTest(b);
2032#endif
2033  if ((SR_HDL(b) & SR_INT)
2034  || (b->s==3))
2035  {
2036    // b is 1/(b->n) => b->n is 1 => result is a
2037    return nlCopy(a);
2038  }
2039  number aa=NULL;
2040  if (SR_HDL(a) & SR_INT)
2041  {
2042    aa=nlRInit(SR_TO_INT(a));
2043    a=aa;
2044  }
2045  result=(number)Alloc(sizeof(rnumber));
2046#ifdef LDEBUG
2047  result->debug=123456;
2048#endif
2049  result->s=3;
2050  MP_INT gcd;
2051  mpz_init(&gcd);
2052  mpz_init(&result->z);
2053  mpz_gcd(&gcd,&a->z,&b->n);
2054  if (mpz_cmp_si(&gcd,(long)1)!=0)
2055  {
2056    MP_INT bt;
2057    mpz_init_set(&bt,&b->n);
2058    MPZ_EXACTDIV(&bt,&bt,&gcd);
2059    mpz_mul(&result->z,&a->z,&bt);
2060    mpz_clear(&bt);
2061  }
2062  else
2063    mpz_mul(&result->z,&a->z,&b->n);
2064  mpz_clear(&gcd);
2065  if (aa!=NULL)
2066  {
2067    mpz_clear(&aa->z);
2068    Free((ADDRESS)aa,sizeof(rnumber));
2069  }
2070  nlGmpSimple(&result->z);
2071  if (mpz_size1(&result->z)<=MP_SMALL)
2072  {
2073    int ui=(int)mpz_get_si(&result->z);
2074    if ((((ui<<3)>>3)==ui)
2075    && (mpz_cmp_si(&result->z,(long)ui)==0))
2076    {
2077      mpz_clear(&result->z);
2078      Free((ADDRESS)result,sizeof(rnumber));
2079      return INT_TO_SR(ui);
2080    }
2081  }
2082#ifdef LDEBUG
2083  nlTest(result);
2084#endif
2085  return result;
2086}
2087
2088int nlModP(number n, int p)
2089{
2090  if (SR_HDL(n) & SR_INT)
2091  {
2092    int i=SR_TO_INT(n);
2093    if (i<0) return (p-((-i)%p));
2094    return i%p;
2095  }
2096  int iz=mpz_mmod_ui(NULL,&n->z,(unsigned long)p);
2097  if (n->s!=3)
2098  {
2099    int in=mpz_mmod_ui(NULL,&n->n,(unsigned long)p);
2100    return (int)npDiv((number)iz,(number)in);
2101  }
2102  return iz;
2103}
2104
2105#else
2106
2107/*0 --------------------------------------------implementation*/
2108
2109typedef unsigned int int32;
2110#define MAX_INT16 0XFFFF
2111#define SIGNUM 0X8000
2112#define BASIS 0X10000
2113#define MOD_BASIS 0XFFFF
2114#define SHIFT_BASIS 0X10
2115
2116/*3
2117*  x%SIGNUM == length(x)
2118*  x/SIGNUM == 1 : x < 0
2119*/
2120#define NLLENGTH(a) ((int)(*(a)&0X7FFF))
2121#define NLNEG(a) (*(a)^SIGNUM)
2122#define NLSIGN(a) (*(a)>SIGNUM)
2123
2124/*3
2125*allocate memory for a long integer
2126*/
2127static lint nlAllocLint (int L)
2128{
2129  lint x;
2130
2131  if (L < SIGNUM)
2132  {
2133    x = (lint)Alloc((L+1)*sizeof(int16));
2134    *x = L;
2135    return x;
2136  }
2137  else
2138  {
2139    WerrorS("number too large");
2140    return NULL;
2141  }
2142}
2143
2144/*3
2145*deallocate Memory of a long integer
2146*/
2147static void nlDeleteLint (lint x)
2148{
2149  if (x)
2150  {
2151    Free((ADDRESS)x,(NLLENGTH(x)+1)*sizeof(int16));
2152  }
2153}
2154
2155/*3
2156*nlShrink the memory used by c
2157*/
2158static void nlShrink (lint * c)
2159{
2160  int j, i;
2161  lint  cc = *c;
2162
2163  i = j = NLLENGTH(cc);
2164  while ((i!=0) && (cc[i]==0))
2165  {
2166    i--;
2167  }
2168  if (i != j)
2169  {
2170    cc = nlAllocLint(i);
2171    memcpy(cc+1,(*c)+1,i*sizeof(int16));
2172    nlDeleteLint(*c);
2173    *c = cc;
2174  }
2175}
2176
2177/*3
2178* x == 1 ?
2179*/
2180static BOOLEAN nlIsOneLint (lint x)
2181{
2182  return ((*x == 1) && (x[1] == 1));
2183}
2184
2185/*3
2186*compares absolute values
2187*abs(x) > abs(y) nlCompAbs(x,y) = 1;
2188*abs(x) < abs(y)               -1;
2189*abs(x) = abs(y)                0;
2190*/
2191static int nlCompAbs (lint x, lint y)
2192{
2193  int i = NLLENGTH(x);
2194
2195  if (i != NLLENGTH(y))
2196  {
2197    if (i > NLLENGTH(y))
2198    {
2199      return 1;
2200    }
2201    else
2202    {
2203      return -1;
2204    }
2205  }
2206  else
2207  {
2208    do
2209    {
2210      if (x[i] > y[i])
2211      {
2212        return 1;
2213      }
2214      if (x[i] < y[i])
2215      {
2216        return -1;
2217      }
2218      i--;
2219    }
2220    while (i!=0);
2221    return 0;
2222  }
2223}
2224
2225/*3
2226*compares long integers
2227*(x) > (y)     nlComp(x,y) = 1;
2228*(x) < (y)                -1;
2229*(x) = (y)                 0;
2230*/
2231static int nlComp (lint x, lint y)
2232{
2233  int i = NLSIGN(x);
2234
2235  if (i != NLSIGN(y))
2236  {
2237    if (i != 0)
2238    {
2239      return -1;
2240    }
2241    else
2242    {
2243      return 1;
2244    }
2245  }
2246  else if (i != 0)
2247  {
2248    return -nlCompAbs(x, y);
2249  }
2250  else
2251  {
2252    return nlCompAbs(x, y);
2253  }
2254}
2255
2256/*3
2257* copy f to t
2258*/
2259static lint nlCopyLint (lint f)
2260{
2261  int i;
2262  lint t;
2263
2264  if (f!=NULL)
2265  {
2266    i = (NLLENGTH(f)+1)*sizeof(int16);
2267    t = (lint)Alloc(i);
2268    memcpy(t, f, i);
2269    return t;
2270  }
2271  else
2272  {
2273    return NULL;
2274  }
2275}
2276
2277/*3
2278* n := 1
2279*/
2280static lint nlInitOne()
2281{
2282  lint n = (lint)Alloc(2*sizeof(int16));
2283  *n = 1;
2284  n[1] = 1;
2285  return n;
2286}
2287
2288/*3
2289* n := i
2290*/
2291static lint nlToLint (int16 i)
2292{
2293  if (i)
2294  {
2295    lint n = (lint)Alloc(2*sizeof(int16));
2296    *n = 1;
2297    n[1] = i;
2298    return n;
2299  }
2300  else
2301  {
2302    return NULL;
2303  }
2304}
2305
2306/*3
2307* Add longs, signs are equal
2308*/
2309static lint nlAddBasic (lint a, lint b)
2310{
2311  int i, max, min;
2312  int32 w;
2313  lint  xx,h;
2314
2315  max = NLLENGTH(a);
2316  min = NLLENGTH(b);
2317  if ((min == 1) && (max == 1))
2318  {
2319    w = (int32) a[1] + (int32) b[1];
2320    if (w < BASIS)
2321    {
2322      xx = (lint)Alloc(2*sizeof(int16));
2323      *xx = 1;
2324      xx[1] = (int16) w;
2325    }
2326    else
2327    {
2328      xx = (lint)Alloc(3 * sizeof(int16));
2329      *xx = 2;
2330      xx[1] = (int16) (MOD_BASIS&w);
2331      xx[2] = 1;
2332    }
2333    if (NLSIGN(a))
2334    {
2335      *xx += SIGNUM;
2336    }
2337    return xx;
2338  }
2339  if(min>max)
2340  {
2341    i = max;
2342    max = min;
2343    min = i;
2344    h = b;
2345  }
2346  else
2347  {
2348    h = a;
2349  }
2350  xx = (lint)Alloc((max+1)*sizeof(int16));
2351  *xx = (int16)max;
2352  w = (int32) a[1] + (int32) b[1];
2353  xx[1] = (int16) (MOD_BASIS&w);
2354  for (i=2; i<=min; i++)
2355  {
2356    w = (w>>SHIFT_BASIS) + (int32) a[i] + (int32) b[i];
2357    xx[i] = (int16) (MOD_BASIS&w);
2358  }
2359  i = min;
2360  loop
2361  {
2362    i++;
2363    if ((w<BASIS) || (i>max))
2364    {
2365      break;
2366    }
2367    w = (w>>SHIFT_BASIS) + (int32) h[i];
2368    xx[i] = (int16) (MOD_BASIS&w);
2369  }
2370  for(; i<=max; i++)
2371  {
2372    xx[i] = h[i];
2373  }
2374  if (w>=BASIS)
2375  {
2376    h = (lint)Alloc((max+2)*sizeof(int16));
2377    memcpy(h+1,xx+1,max*sizeof(int16));
2378    max++;
2379    *h = (int16)max;
2380    h[max] = 1;
2381    nlDeleteLint(xx);
2382    xx = h;
2383  }
2384  if (NLSIGN(a))
2385  {
2386    *xx += SIGNUM;
2387  }
2388  return xx;
2389}
2390
2391/*3
2392* Subtract longs, signs are equal, a > b
2393*/
2394static lint nlSubBasic (lint a, lint b)
2395{
2396  int i,j1,j2;
2397  int16 w;
2398  lint xx;
2399  j1 = NLLENGTH(a);
2400  if (j1 == 1)
2401  {
2402    xx = (lint)Alloc(2*sizeof(int16));
2403    *xx = 1;
2404    xx[1] = a[1]-b[1];
2405    if (NLSIGN(a))
2406    {
2407      *xx += SIGNUM;
2408    }
2409    return xx;
2410  }
2411  j2 = NLLENGTH(b);
2412  xx = (lint)Alloc((j1+1)*sizeof(int16));
2413  *xx = (int16)j1;
2414  if(a[1] < b[1])
2415  {
2416    xx[1] = MAX_INT16-b[1]+a[1]+1;
2417    w = 1;
2418  }
2419  else
2420  {
2421    xx[1] = a[1]-b[1];
2422    w = 0;
2423  }
2424  for (i=2; i<=j2; i++)
2425  {
2426    if(a[i]>b[i])
2427    {
2428      xx[i] = a[i]-b[i]-w;
2429      w = 0;
2430    }
2431    else if(a[i]<b[i])
2432    {
2433      if(w!=0)
2434      {
2435        xx[i] = MAX_INT16 - b[i] + a[i];
2436      }
2437      else
2438      {
2439        xx[i] = MAX_INT16 - b[i] + a[i] + 1;
2440      }
2441      w = 1;
2442    }
2443    else
2444    {
2445      if(w!=0)
2446      {
2447        xx[i] = MAX_INT16;
2448      }
2449      else
2450      {
2451        xx[i] = 0;
2452      }
2453    }
2454  }
2455  i = j2;
2456  while(w!=0)
2457  {
2458    i++;
2459    if(a[i])
2460    {
2461      xx[i] = a[i]-w;
2462      w = 0;
2463    }
2464    else
2465    {
2466      xx[i] = MAX_INT16;
2467    }
2468  }
2469  while(i < j1)
2470  {
2471    i++;
2472    xx[i] = a[i];
2473  }
2474  nlShrink(&xx);
2475  if (NLSIGN(a))
2476  {
2477    *xx += SIGNUM;
2478  }
2479  return xx;
2480}
2481
2482/*3
2483* c:= a + b
2484*/
2485static lint nlAddLong (lint a,lint b)
2486{
2487  if (a==NULL)
2488    return nlCopyLint(b);
2489  else if (b==NULL)
2490    return nlCopyLint(a);
2491  if (NLSIGN(a) == NLSIGN(b))
2492    return nlAddBasic(a,b);
2493  else
2494  {
2495    switch(nlCompAbs(a, b))
2496    {
2497      case -1:
2498        return nlSubBasic(b, a);
2499      case 1:
2500        return nlSubBasic(a, b);
2501      default:
2502        ;
2503    }
2504  }
2505  return NULL;
2506}
2507
2508/*3
2509* c:= a - b
2510*/
2511static lint nlSubLong (lint a,lint b)
2512{
2513  lint c;
2514  if(b==NULL)
2515    return nlCopyLint(a);
2516  if(a==NULL)
2517  {
2518    c = nlCopyLint(b);
2519    *c = NLNEG(c);
2520    return c;
2521  }
2522  if (NLSIGN(a) != NLSIGN(b))
2523    return nlAddBasic(a,b);
2524  else
2525  {
2526    switch(nlCompAbs(a, b))
2527    {
2528      case 1:
2529        return nlSubBasic(a, b);
2530      case -1:
2531        c = nlSubBasic(b, a);
2532        *c = NLNEG(c);
2533        return c;
2534      default:
2535        ;
2536    }
2537  }
2538  return NULL;
2539}
2540
2541/*3
2542* a[i] = a[i] + b[i]*x, i = 1 to length(b)
2543*/
2544static void nlAddMult (lint a, lint  b, int32 x)
2545{
2546  int32 w;
2547  int i, l = NLLENGTH(b);
2548  w = (int32) a[1] + (int32) b[1]*x;
2549  a[1] = (int16) (MOD_BASIS&w);
2550  for (i=2; i<=l; i++)
2551  {
2552    w = (w>>SHIFT_BASIS) + (int32) a[i] + (int32) b[i]*x;
2553    a[i] = (int16) (MOD_BASIS&w);
2554  }
2555  if (w>=BASIS)
2556    a[l+1] = (int16)(w>>SHIFT_BASIS);
2557}
2558
2559/*3
2560* a[i] = b[i]*x, i = 1 to length(b)
2561*/
2562static void nlIniMult (lint a, lint  b, int32 x)
2563{
2564  int i, l = NLLENGTH(b);
2565  int32 w;
2566  w = (int32) b[1]*x;
2567  a[1] = (int16) (MOD_BASIS&w);
2568  for (i=2; i<=l; i++)
2569  {
2570    w = (w>>SHIFT_BASIS) + (int32) b[i]*x;
2571    a[i] = (int16) (MOD_BASIS&w);
2572  }
2573  if (w>=BASIS)
2574    a[l+1] = (int16)(w>>SHIFT_BASIS);
2575}
2576
2577/*3
2578* mult long, c:= a * b
2579*/
2580static lint nlMultLong (lint a,lint b)
2581{
2582  int la, lb, i;
2583  int32 h;
2584  lint xx;
2585  if ((a==NULL) || (b==NULL))
2586    return NULL;
2587  la = NLLENGTH(a);
2588  lb = NLLENGTH(b);
2589  if ((la == 1) && (lb == 1))
2590  {
2591    h = (int32) a[1] * (int32) b[1];
2592    if (h < BASIS)
2593    {
2594      xx = (lint)Alloc(2*sizeof(int16));
2595      *xx = 1;
2596      xx[1] = (int16) h;
2597    }
2598    else
2599    {
2600      xx = (lint)Alloc(3 * sizeof(int16));
2601      *xx = 2;
2602      xx[1] = (int16) (MOD_BASIS&h);
2603      xx[2] = (int16) (h>>SHIFT_BASIS);
2604    }
2605    if ((NLSIGN(a)) != (NLSIGN(b)))
2606      *xx += SIGNUM;
2607    return xx;
2608  }
2609  i = la+lb;
2610  xx = (lint)Alloc((i+1)*sizeof(int16));
2611  *xx = (int16)i;
2612  memset(xx+1,0,i*sizeof(int16));
2613  if (la < lb)
2614  {
2615    nlIniMult(xx, b, (int32) a[1]);
2616    for (i=1; i<la; i++)
2617      nlAddMult(xx+i, b, (int32) a[i+1]);
2618  }
2619  else
2620  {
2621    nlIniMult(xx, a, (int32) b[1]);
2622    for (i=1; i<lb; i++)
2623      nlAddMult(xx+i, a, (int32) b[i+1]);
2624  }
2625  nlShrink(&xx);
2626  if ((NLSIGN(a)) != (NLSIGN(b)))
2627    *xx += SIGNUM;
2628  return xx;
2629}
2630
2631/*3
2632* aq := aq div b, b small
2633*/
2634static void nlQuotSmall (lint * aq, int16 b)
2635{
2636  int al,j;
2637  int16 a1,c1;
2638  int32 db,da,dc;
2639  lint a,cc;
2640  db = (int32) b;
2641  a = *aq;
2642  al = NLLENGTH(a);
2643  a1 = a[al];
2644  if (al == 1)
2645  {
2646    a[1] = a1 / b;
2647    return;
2648  }
2649  if (a1 >= b)
2650  {
2651    cc = a;
2652    c1 = a1 / b;
2653    cc[al] = c1;
2654    a1 -= b*c1;
2655  }
2656  else
2657    cc = nlAllocLint(al-1);
2658  loop
2659  {
2660    al--;
2661    if (al==0)
2662    {
2663      if(a != cc)
2664        nlDeleteLint(a);
2665      *aq = cc;
2666      return;
2667    }
2668    if (a1==0)
2669    {
2670      j = al;
2671      loop
2672      {
2673        a1 = a[j];
2674        if (a1 >= b)
2675        {
2676          c1 = a1 / b;
2677          cc[j] = c1;
2678          a1 -= b*c1;
2679        }
2680        else
2681          cc[j] = 0;
2682        j--;
2683        if (j==0)
2684        {
2685          if(a != cc)
2686            nlDeleteLint(a);
2687          *aq = cc;
2688          return;
2689        }
2690        else if (a1!=0)
2691        {
2692          al = j;
2693          break;
2694        }
2695      }
2696    }
2697    da = BASIS*(int32) a1 + (int32) a[al];
2698    dc = da / db;
2699    cc[al] = (int16) dc;
2700    da -= db*dc;
2701    a1 = (int16) da;
2702  }
2703}
2704
2705/*3
2706* a := a - b, b1...bj
2707*/
2708static void nlSubSimple (lint a, lint b, int j)
2709{
2710  int i;
2711  int16 w;
2712  if (a[1] < b[1])
2713  {
2714    a[1] += MAX_INT16-b[1]+1;
2715    w = 1;
2716  }
2717  else
2718  {
2719    a[1] -= b[1];
2720    w = 0;
2721  }
2722  for (i=2; i<=j; i++)
2723  {
2724    if(a[i] > b[i])
2725    {
2726      a[i] -= b[i]+w;
2727      w = 0;
2728    }
2729    else if(a[i] < b[i])
2730    {
2731      if (w)
2732        a[i] += MAX_INT16 - b[i];
2733      else
2734        a[i] += MAX_INT16 - b[i] + 1;
2735      w = 1;
2736    }
2737    else
2738    {
2739      if(w)
2740        a[i] = MAX_INT16;
2741      else
2742        a[i] = 0;
2743    }
2744  }
2745  if(w!=0)
2746    a[j+1]--;
2747}
2748
2749/*3
2750* a := a - b*x, b1...bj
2751*/
2752static void nlSubMult(lint a, lint b, int32 x, int l)
2753{
2754  int i;
2755  int16 k, w;
2756  int32 y;
2757  y = (int32)b[1]*x;
2758  k = (int16)(MOD_BASIS&y);
2759  if (a[1] < k)
2760  {
2761    a[1] += MAX_INT16-k+1;
2762    w = 1;
2763  }
2764  else
2765  {
2766    a[1] -= k;
2767    w = 0;
2768  }
2769  for(i=2; i<=l; i++)
2770  {
2771    y = (y>>SHIFT_BASIS) + (int32)b[i]*x;
2772    k = (int16)(MOD_BASIS&y);
2773    if(a[i] > k)
2774    {
2775      a[i] -= k+w;
2776      w = 0;
2777    }
2778    else if(a[i] < k)
2779    {
2780      if (w!=0)
2781        a[i] += MAX_INT16 - k;
2782      else
2783        a[i] += MAX_INT16 - k + 1;
2784      w = 1;
2785    }
2786    else
2787    {
2788      if(w!=0)
2789        a[i] = MAX_INT16;
2790      else
2791        a[i] = 0;
2792    }
2793  }
2794  w += (int16)(y>>SHIFT_BASIS);
2795  if(w!=0)
2796    a[l+1] -= w;
2797}
2798
2799/*3
2800* compare xi and yi from i=l to 1
2801*/
2802static BOOLEAN nlCompHigh (lint x, lint y, int l)
2803{
2804  int i;
2805  for (i=l; i!=0; i--)
2806  {
2807    if (x[i] > y[i])
2808      return TRUE;
2809    else if (x[i] < y[i])
2810      return FALSE;
2811  }
2812  return TRUE;
2813}
2814
2815/*3
2816* aq := aq div bn, bn long
2817*/
2818static void nlQuotLong (lint * aq, lint bn)
2819{
2820  int bl,al,si,j;
2821  int16 b1,a1,a2;
2822  int32 db,db1,dc,dc1,dc2;
2823  lint  an,cc,ansi;
2824  BOOLEAN rel;
2825  bl = NLLENGTH(bn);
2826  b1 = bn[bl];
2827  db = BASIS*(int32)b1 + (int32)bn[bl-1];
2828  if ((bl > 2) && ((b1 < MAX_INT16) || (bn[bl-1] < MAX_INT16)))
2829    db++;
2830  db1 = (int32)(b1+1);
2831  an = *aq;
2832  al = NLLENGTH(an);
2833  a1 = an[al];
2834  a2 = an[al-1];
2835  si = al-bl;
2836  ansi = an+si;
2837  rel = nlCompHigh(ansi, bn, bl);
2838  if (rel)
2839  {
2840    cc = nlAllocLint(si+1);
2841    dc = (BASIS*(int32)a1 + (int32)a2) / db;
2842    if (dc > 1)
2843    {
2844      nlSubMult(ansi, bn, dc, bl);
2845      rel = nlCompHigh(ansi, bn, bl);
2846    }
2847    else
2848      dc = 0;
2849    while (rel)
2850    {
2851      dc++;
2852      nlSubSimple(ansi, bn, bl);
2853      rel = nlCompHigh(ansi, bn, bl);
2854    }
2855    cc[si+1] = (int16)dc;
2856    a1 = an[al];
2857    a2 = an[al-1];
2858  }
2859  else
2860    cc = nlAllocLint(si);
2861  loop
2862  {
2863    al--;
2864    if (al < bl)
2865    {
2866      nlDeleteLint(an);
2867      *aq = cc;
2868      return;
2869    }
2870    si--;
2871    ansi--;
2872    if (a1==0)
2873    {
2874      j = al;
2875      loop
2876      {
2877        a1 = a2;
2878        a2 = an[j-1];
2879        rel = nlCompHigh(ansi, bn, bl);
2880        if (rel)
2881        {
2882          dc = (BASIS*(int32)a1 + (int32)a2) / db;
2883          if (dc > 1)
2884          {
2885            nlSubMult(ansi, bn, dc, bl);
2886            rel = nlCompHigh(ansi, bn, bl);
2887          }
2888          else
2889            dc = 0;
2890          while (rel)
2891          {
2892            dc++;
2893            nlSubSimple(ansi, bn, bl);
2894            rel = nlCompHigh(ansi, bn, bl);
2895          }
2896          cc[si+1] = (int16)dc;
2897          a1 = an[j];
2898          a2 = an[j-1];
2899        }
2900        else
2901          cc[si+1] = 0;
2902        si--;
2903        ansi--;
2904        j--;
2905        if (j < bl)
2906        {
2907          nlDeleteLint(an);
2908          *aq = cc;
2909          return;
2910        }
2911        else if(a1!=0)
2912        {
2913          al = j;
2914          break;
2915        }
2916      }
2917    }
2918    dc2 = dc1 = (BASIS*(int32)a1+(int32)a2) / db1;
2919    loop
2920    {
2921      nlSubMult(ansi, bn, dc1, bl);
2922      a1 = an[al+1];
2923      if (a1!=0)
2924      {
2925        dc1 = (BASIS*(int32)a1+(int32)an[al]) / db1;
2926        dc2 += dc1;
2927      }
2928      else
2929        break;
2930    }
2931    a1 = an[al];
2932    a2 = an[al-1];
2933    rel = nlCompHigh(ansi, bn, bl);
2934    if (rel)
2935    {
2936      dc = (BASIS*(int32)a1 + (int32)a2) / db;
2937      if (dc > 1)
2938      {
2939        nlSubMult(ansi, bn, dc, bl);
2940        rel = nlCompHigh(ansi, bn, bl);
2941      }
2942      else
2943        dc = 0;
2944      while (rel)
2945      {
2946        dc++;
2947        nlSubSimple(ansi, bn, bl);
2948        rel = nlCompHigh(ansi, bn, bl);
2949      }
2950      a1 = an[al];
2951      a2 = an[al-1];
2952    }
2953    else
2954      dc = 0;
2955    cc[si+1] = (int16) (dc2+dc);
2956  }
2957}
2958
2959/*3
2960* aq := aq div b
2961*/
2962static void nlQuotient (lint * aq, lint b)
2963{
2964  lint a;
2965  int bl;
2966  BOOLEAN neg;
2967  a = *aq;
2968  neg = (NLSIGN(a) != NLSIGN(b));
2969  bl = NLLENGTH(b);
2970  if(bl == 1)
2971    nlQuotSmall(&a,b[1]);
2972  else
2973    nlQuotLong(&a,b);
2974  if(neg && !NLSIGN(a))
2975    *a += SIGNUM;
2976  else if (!neg && NLSIGN(a))
2977    *a -= SIGNUM;
2978  *aq = a;
2979}
2980
2981/*3
2982* r := a mod b, b small
2983* al = length(a)
2984*/
2985static void nlRemSmall (lint a, int al, int16 b, int16 * r)
2986{
2987  int j;
2988  int16 a1,c1;
2989  int32 db,da,dc;
2990  db = (int32)b;
2991  a1 = a[al];
2992  if (al == 1)
2993  {
2994    c1 = a1 / b;
2995    a1 -= b*c1;
2996    *r = a1;
2997    return;
2998  }
2999  if (a1 >= b)
3000  {
3001    c1 = a1 / b;
3002    a1 -= b*c1;
3003  }
3004  loop
3005  {
3006    al--;
3007    if (al==0)
3008    {
3009      *r = a1;
3010      return;
3011    }
3012    if (a1==0)
3013    {
3014      j = al;
3015      loop
3016      {
3017        a1 = a[j];
3018        if (a1 >= b)
3019        {
3020          c1 = a1 / b;
3021          a1 -= b*c1;
3022        }
3023        j--;
3024        if (j==0)
3025        {
3026          *r = a1;
3027          return;
3028        }
3029        else if (a1!=0)
3030        {
3031          al = j;
3032          break;
3033        }
3034      }
3035    }
3036    da = BASIS*(int32)a1 + (int32)a[al];
3037    dc = da / db;
3038    da -= db*dc;
3039    a1 = (int16)da;
3040  }
3041}
3042
3043static void nlShrink0(lint c, int j, int * i)
3044{
3045  while((j!=0) && (c[j]==0)) j--;
3046  *i = j;
3047}
3048
3049/*3
3050* rl := an mod bn, bn long, but <= an
3051* al = length(an)
3052*/
3053static void nlRemLong(lint an, lint bn, int al, int bl, int * rl)
3054{
3055  int si,j;
3056  int16 b1,a1,a2;
3057  int32 db,db1,dc;
3058  lint ansi;
3059  BOOLEAN rel;
3060  b1 = bn[bl];
3061  db = (int32)bn[bl-1];
3062  db += (int32)b1<<SHIFT_BASIS;
3063  if ((bl > 2) && ((b1 < MAX_INT16) || (bn[bl-1] < MAX_INT16)))
3064    db++;
3065  db1 = (int32)b1+1;
3066  a1 = an[al];
3067  a2 = an[al-1];
3068  si = al-bl;
3069  ansi = an+si;
3070  rel = nlCompHigh(ansi, bn, bl);
3071  if (rel)
3072  {
3073    dc = (BASIS*(int32)a1 + (int32)a2) / db;
3074    if (dc > 1)
3075    {
3076      nlSubMult(ansi, bn, dc, bl);
3077      rel = nlCompHigh(ansi, bn, bl);
3078    }
3079    while (rel)
3080    {
3081      nlSubSimple(ansi, bn, bl);
3082      rel = nlCompHigh(ansi, bn, bl);
3083    }
3084    a1 = an[al];
3085    a2 = an[al-1];
3086  }
3087  loop
3088  {
3089    al--;
3090    if (al < bl)
3091    {
3092      nlShrink0(an,bl,rl);
3093      return;
3094    }
3095    si--;
3096    ansi--;
3097    if(a1==0)
3098    {
3099      j = al;
3100      loop
3101      {
3102        a1 = a2;
3103        a2 = an[j-1];
3104        rel = nlCompHigh(ansi, bn, bl);
3105        if (rel)
3106        {
3107          dc = (BASIS*(int32)a1 + (int32)a2) / db;
3108          if (dc > 1)
3109          {
3110            nlSubMult(ansi, bn, dc, bl);
3111            rel = nlCompHigh(ansi, bn, bl);
3112          }
3113          while (rel)
3114          {
3115            nlSubSimple(ansi, bn, bl);
3116            rel = nlCompHigh(ansi, bn, bl);
3117          }
3118          a1 = an[j];
3119          a2 = an[j-1];
3120        }
3121        si--;
3122        ansi--;
3123        j--;
3124        if (j < bl)
3125        {
3126          nlShrink0(an,bl,rl);
3127          return;
3128        }
3129        else if (a1!=0)
3130        {
3131          al = j;
3132          break;
3133        }
3134      }
3135    }
3136    dc = (BASIS*(int32)a1+(int32)a2) / db1;
3137    loop
3138    {
3139      nlSubMult(ansi, bn, dc, bl);
3140      a1 = an[al+1];
3141      if (a1!=0)
3142        dc = (BASIS*(int32)a1+(int32)an[al]) / db1;
3143      else
3144        break;
3145    }
3146    a1 = an[al];
3147    a2 = an[al-1];
3148    rel = nlCompHigh(ansi, bn, bl);
3149    if (rel)
3150    {
3151      dc = (BASIS*(int32)a1 + (int32)a2) / db;
3152      if (dc > 1)
3153      {
3154        nlSubMult(ansi, bn, dc, bl);
3155        rel = nlCompHigh(ansi, bn, bl);
3156      }
3157      while (rel)
3158      {
3159        nlSubSimple(ansi, bn, bl);
3160        rel = nlCompHigh(ansi, bn, bl);
3161      }
3162      a1 = an[al];
3163      a2 = an[al-1];
3164    }
3165  }
3166}
3167
3168/*3
3169* z = GCD(a, b)
3170*/
3171static void nlGcdLong (lint a, lint b, lint * z)
3172{
3173  lint h, x, y;
3174  int xl,yl,hl;
3175  int16 xc,yc,hc;
3176
3177  x = nlCopyLint(a);
3178  if (NLSIGN(x))
3179    *x %= SIGNUM;
3180  if (nlIsOneLint(x))
3181  {
3182    *z=x;
3183    return;
3184  }
3185  y = nlCopyLint(b);
3186  if (NLSIGN(y))
3187    *y %= SIGNUM;
3188  if (nlIsOneLint(y))
3189  {
3190    *z=y;
3191    nlDeleteLint(x);
3192    return;
3193  }
3194  switch (nlCompAbs(x,y))
3195  {
3196    case -1: h=x;x=y;y=h;h=NULL;
3197            break;
3198    case 0: *z=x;
3199            nlDeleteLint(y);
3200            return;
3201    case 1: break;
3202  }
3203  xl = NLLENGTH(x);
3204  yl = NLLENGTH(y);
3205  /* we have: x > y >=0 as lint, xl>=yl>=0 as int*/
3206  while (yl > 1)
3207  {
3208    nlRemLong(x,y,xl,yl,&hl);
3209    h = x;
3210    x = y;
3211    y = h;
3212    xl = yl;
3213    yl = hl;
3214  }
3215  if(yl==0)
3216  {
3217    /*result is x with xl >1*/
3218    nlDeleteLint(y);
3219    nlShrink(&x);
3220    *z = x;
3221    return;
3222  }
3223  else
3224  {
3225    yc = y[1];
3226    if (xl>1)
3227      nlRemSmall(x,xl,yc,&hc);
3228    else
3229    {
3230      hc = x[1];
3231      if (hc > yc)
3232      {
3233        xc=yc;
3234        yc=hc;
3235        hc=xc;
3236      }
3237    }
3238  }
3239  nlDeleteLint(x);
3240  nlDeleteLint(y);
3241  /*we have now yc >hc*/
3242  while(hc!=0)
3243  {
3244    xc = yc;
3245    yc = hc;
3246    hc = xc % yc;
3247  }
3248  *z = nlToLint(yc);
3249}
3250
3251/*-----------------------------------------------------------------*/
3252
3253/*0
3254* Begin of operations with rational numbers
3255*/
3256
3257/*3
3258* parameter s in number:
3259*  0 :  integer
3260*  1 :  normalised rational
3261* -1 :  not normalised rational
3262*/
3263
3264#define NLISINT(a)  !((a)->s)
3265#define NLISRAT(a)   ((a)->s)
3266
3267#ifdef LDEBUG
3268#define nlTest(a) nlDBTest(a,__FILE__,__LINE__)
3269BOOLEAN nlDBTest(number a, char *f,int l)
3270{
3271  BOOLEAN res;
3272  if (a==NULL)
3273  {
3274    return TRUE;
3275  }
3276  else
3277  {
3278    if (a->z!=NULL)
3279    {
3280      res=mmDBTestBlock(a->z,(NLLENGTH(a->z)+1)*sizeof(int16),f,l);
3281    }
3282    if (NLISRAT(a))
3283    {
3284      res = (res && mmDBTestBlock(a->n,(NLLENGTH(a->n)+1)*sizeof(int16),f,l));
3285    }
3286    return res;
3287  }
3288}
3289#endif
3290
3291/*2
3292* delete a
3293*/
3294#ifdef LDEBUG
3295void nlDBDelete (number *a, char *f, int l)
3296#else
3297void nlDelete (number *a)
3298#endif
3299{
3300  number b=*a;
3301#ifdef LDEBUG
3302  nlTest(b);
3303#endif
3304
3305  if (b!=NULL)
3306  {
3307    nlDeleteLint(b->z);
3308    if (NLISRAT(b))
3309    {
3310      nlDeleteLint(b->n);
3311    }
3312#if defined(LDEBUG) && defined(MDEBUG)
3313    mmDBFreeBlock((ADDRESS)b,sizeof(rnumber),f,l);
3314#else
3315    Free((ADDRESS)b,sizeof(rnumber));
3316#endif
3317    *a=NULL;
3318  }
3319}
3320
3321/*2
3322* r = 0
3323*/
3324void nlNew (number * r)
3325{
3326  *r=NULL;
3327}
3328
3329/*2
3330* result->z =gcd(a->z,b->z)
3331*/
3332number nlGcd(number a, number b)
3333{
3334  number result=(number)Alloc(sizeof(rnumber));
3335
3336  result->s = 0;
3337  nlGcdLong(a->z,b->z,&result->z);
3338#ifdef LDEBUG
3339  nlTest(result);
3340#endif
3341  return result;
3342}
3343
3344/*2
3345* result->z = lcm(a->z,b->n)
3346*/
3347number nlLcm(number a, number b)
3348{
3349  lint t,bt;
3350  number result=(number)Alloc(sizeof(rnumber));
3351
3352  result->s = 0;
3353  if (NLISINT(b))
3354  {
3355    result->z = nlCopyLint(a->z);
3356  }
3357  else
3358  {
3359    nlGcdLong(a->z,b->n,&t);
3360    if (!nlIsOneLint(t))
3361    {
3362      bt = nlCopyLint(b->n);
3363      nlQuotient(&bt,t);
3364      result->z = nlMultLong(a->z,bt);
3365      nlDeleteLint(bt);
3366    }
3367    else
3368    {
3369      result->z = nlMultLong(a->z,b->n);
3370    }
3371    nlDeleteLint(t);
3372  }
3373#ifdef LDEBUG
3374  nlTest(result);
3375#endif
3376  return result;
3377}
3378
3379/*2
3380* a == 0 ?
3381*/
3382BOOLEAN nlIsZero (number a)
3383{
3384#ifdef LDEBUG
3385  nlTest(a);
3386#endif
3387  return (a == NULL);
3388}
3389
3390/*2
3391* a = b ?
3392*/
3393BOOLEAN nlEqual (number a, number b)
3394{
3395  number h;
3396
3397  if(a==NULL)
3398    return (b==NULL);
3399  if(b==NULL)
3400    return(a==NULL);
3401  if(NLISINT(a) && NLISINT(b))
3402  {
3403    return !nlComp(a->z,b->z);
3404  }
3405  h = nlSub(a, b);
3406  if (nlIsZero(h))
3407  {
3408    nlDelete(&h);
3409    return TRUE;
3410  }
3411  else
3412  {
3413    nlDelete(&h);
3414    return FALSE;
3415  }
3416}
3417
3418/*2
3419* a == 1 ?
3420*/
3421BOOLEAN nlIsOne (number a)
3422{
3423#ifdef LDEBUG
3424  nlTest(a);
3425#endif
3426  if(nlIsZero(a) || NLSIGN(a->z))
3427  {
3428    return FALSE;
3429  }
3430  if(NLISINT(a))
3431  {
3432    return ((a->z[1]==1) && (a->z[0]==1));
3433  }
3434  if (nlCompAbs(a->z, a->n))
3435  {
3436    return FALSE;
3437  }
3438  nlDeleteLint(a->z);
3439  nlDeleteLint(a->n);
3440  a->z = nlInitOne();
3441  a->s = 0;
3442  return TRUE;
3443}
3444
3445/*2
3446* a == -1 ?
3447*/
3448BOOLEAN nlIsMOne (number a)
3449{
3450#ifdef LDEBUG
3451  nlTest(a);
3452#endif
3453  if(nlIsZero(a) || !NLSIGN(a->z))
3454  {
3455    return FALSE;
3456  }
3457  if(NLISINT(a))
3458  {
3459    return ((a->z[1]==1) && (a->z[0]==SIGNUM+1));
3460  }
3461  if (nlCompAbs(a->z, a->n))
3462  {
3463    return FALSE;
3464  }
3465  nlDeleteLint(a->z);
3466  nlDeleteLint(a->n);
3467  a->z = nlInitOne();
3468  a->z[0] += SIGNUM;
3469  a->s = 0;
3470  return TRUE;
3471}
3472
3473/*2
3474* z := i
3475*/
3476number nlInit (int i)
3477{
3478  lint nom;
3479  BOOLEAN neg = FALSE;
3480  number z = NULL;
3481
3482  if (i!=0)
3483  {
3484    z = (number)Alloc(sizeof(rnumber));
3485    z->s = 0;
3486    if (i < 0)
3487    {
3488      neg = TRUE;
3489      i = -i;
3490    }
3491    if (i>=BASIS)  // init a 4 byte number
3492    {
3493      nom = nlAllocLint(2);
3494      nom[1] = (int16) (MOD_BASIS&i);
3495      nom[2] = (int16) (i>>SHIFT_BASIS);
3496    }
3497    else           // init a 2 byte number
3498    {
3499      nom = nlAllocLint(1);
3500      nom[1] = (int16) i;
3501    }
3502    if (neg)
3503    {
3504      *nom += SIGNUM;
3505    }
3506    z->z = nom;
3507#ifdef LDEBUG
3508    nlTest(z);
3509#endif
3510  }
3511  return z;
3512}
3513
3514/*2
3515* convert number to int
3516*/
3517int nlInt(number &a)
3518{
3519  lint x;
3520  int res, i;
3521
3522  if (nlIsZero(a))
3523  {
3524    return 0;
3525  }
3526  nlNormalize(a);
3527  if (NLISRAT(a))
3528  {
3529    return 0;
3530  }
3531  x = a->z;
3532  i = NLLENGTH(x);
3533  if (i > 2)
3534  {
3535    return 0;
3536  }
3537  res = x[1];
3538  if (i == 2)
3539  {
3540    if (x[2] < SIGNUM)
3541    {
3542      res += BASIS*x[2];
3543    }
3544    else
3545    {
3546      return 0;
3547    }
3548  }
3549  if (NLSIGN(x))
3550  {
3551    return -res;
3552  }
3553  else
3554  {
3555    return res;
3556  }
3557}
3558
3559/*2
3560* simplify x
3561*/
3562void nlNormalize (number &x)
3563{
3564  lint zh, nh, d;
3565
3566  if ((x==NULL)||(x->s>=0))
3567  {
3568    return;
3569  }
3570  x->s = 1;
3571  zh = x->z;
3572  nh = x->n;
3573  nlGcdLong(zh, nh, &d);
3574  if (!nlIsOneLint(d))
3575  {
3576    nlQuotient(&zh,d);
3577    x->z = zh;
3578    nlQuotient(&nh,d);
3579    if (nlIsOneLint(nh))
3580    {
3581      nlDeleteLint(nh);
3582      x->s = 0;
3583    }
3584    else
3585    {
3586      x->n = nh;
3587    }
3588  }
3589  nlDeleteLint(d);
3590#ifdef LDEBUG
3591  nlTest(x);
3592#endif
3593}
3594
3595/*2
3596* za >= 0 ?
3597*/
3598BOOLEAN nlGreaterZero (number a)
3599{
3600  if (nlIsZero(a))
3601  {
3602    return TRUE;
3603  }
3604  return !NLSIGN(a->z);
3605}
3606
3607/*2
3608* a := -a
3609*/
3610number nlNeg (number a)
3611{
3612  if (!nlIsZero(a))
3613  {
3614    a->z[0] = NLNEG(a->z);
3615  }
3616  return a;
3617}
3618
3619/*
3620* 1/a
3621*/
3622number nlInvers(number a)
3623{
3624  number n=(number)Alloc(sizeof(rnumber));
3625
3626  if (nlIsZero(a))
3627  {
3628    WerrorS("div. 1/0");
3629    Free((ADDRESS)n,sizeof(rnumber));
3630    return NULL;
3631  }
3632  nlNormalize(a);
3633  if(NLISRAT(a))
3634  {
3635    n->z = nlCopyLint(a->n);
3636  }
3637  else
3638  {
3639    n->z = nlInitOne();
3640  }
3641  if((a->z[1]!=1) || (NLLENGTH(a->z)!=1))
3642  {
3643    n->n = nlCopyLint(a->z);
3644    n->s = 1;
3645    if(n->n[0]>SIGNUM)
3646    {
3647      n->n[0] -= SIGNUM;
3648      n->z[0] += SIGNUM;
3649    }
3650  }
3651  else
3652  {
3653    if (NLSIGN(a->z))
3654    {
3655      n->z[0] += SIGNUM;
3656    }
3657    n->s = 0;
3658  }
3659#ifdef LDEBUG
3660  nlTest(n);
3661#endif
3662  return n;
3663}
3664
3665/*2
3666* copy a to b
3667*/
3668number nlCopy(number a)
3669{
3670  number b=(number)Alloc0(sizeof(rnumber));
3671
3672  if (!nlIsZero(a))
3673  {
3674    b->s = a->s;
3675    b->z = nlCopyLint(a->z);
3676    if (NLISRAT(a))
3677    {
3678      b->n = nlCopyLint(a->n);
3679    }
3680#ifdef LDEBUG
3681    nlTest(b);
3682#endif
3683    return b;
3684  }
3685  else
3686  {
3687    Free((ADDRESS)b,sizeof(rnumber));
3688    return NULL;
3689  }
3690}
3691
3692/*2
3693* lu:= la + li
3694*/
3695number nlAdd (number la, number li)
3696{
3697  number lu;
3698  lint x, y;
3699  BOOLEAN a0 = nlIsZero(la), i0 = nlIsZero(li);
3700
3701  if (a0 || i0)
3702  {
3703    if (a0)
3704    {
3705      return nlCopy(li);
3706    }
3707    else
3708    {
3709      return nlCopy(la);
3710    }
3711  }
3712  lu = (number)Alloc(sizeof(rnumber));
3713  if (NLISINT(la))
3714  {
3715    if (NLISINT(li))
3716    {
3717      lu->z = nlAddLong(la->z, li->z);
3718      if (lu->z==NULL)
3719      {
3720        Free((ADDRESS)lu,sizeof(rnumber));
3721        return NULL;
3722      }
3723      lu->s = 0;
3724#ifdef LDEBUG
3725      nlTest(lu);
3726#endif
3727      return lu;
3728    }
3729    else
3730    {
3731      x = nlMultLong(la->z, li->n);
3732      lu->z = nlAddLong(x, li->z);
3733      nlDeleteLint(x);
3734      lu->n = nlCopyLint(li->n);
3735    }
3736  }
3737  else
3738  {
3739    if (NLISINT(li))
3740    {
3741      y = nlMultLong(li->z, la->n);
3742      lu->z = nlAddLong(la->z, y);
3743      nlDeleteLint(y);
3744      lu->n = nlCopyLint(la->n);
3745    }
3746    else
3747    {
3748      x = nlMultLong(la->z, li->n);
3749      y = nlMultLong(li->z, la->n);
3750      lu->z = nlAddLong(x, y);
3751      nlDeleteLint(x);
3752      nlDeleteLint(y);
3753      lu->n = nlMultLong(la->n, li->n);
3754    }
3755  }
3756  if (lu->z==NULL)
3757  {
3758    nlDeleteLint(lu->n);
3759    Free((ADDRESS)lu,sizeof(rnumber));
3760    return NULL;
3761  }
3762  lu->s = -1;
3763#ifdef LDEBUG
3764  nlTest(lu);
3765#endif
3766  return lu;
3767}
3768
3769/*2
3770* lu:= la - li
3771*/
3772number nlSub (number la, number li)
3773{
3774  number lu;
3775  lint x, y;
3776  BOOLEAN a0 = nlIsZero(la), i0 = nlIsZero(li);
3777
3778  if (a0 || i0)
3779  {
3780    if (a0)
3781    {
3782      lu = nlCopy(li);
3783      lu = nlNeg(lu);
3784      return lu;
3785    }
3786    else
3787    {
3788      return nlCopy(la);
3789    }
3790  }
3791  lu = (number)Alloc(sizeof(rnumber));
3792  if (NLISINT(la))
3793  {
3794    if (NLISINT(li))
3795    {
3796      lu->z = nlSubLong(la->z, li->z);
3797      if (lu->z==NULL)
3798      {
3799        Free((ADDRESS)lu,sizeof(rnumber));
3800        return NULL;
3801      }
3802      lu->s = 0;
3803#ifdef LDEBUG
3804      nlTest(lu);
3805#endif
3806      return lu;
3807    }
3808    else
3809    {
3810      x = nlMultLong(la->z, li->n);
3811      lu->z = nlSubLong(x, li->z);
3812      nlDeleteLint(x);
3813      lu->n = nlCopyLint(li->n);
3814    }
3815  }
3816  else
3817  {
3818    if (NLISINT(li))
3819    {
3820      y = nlMultLong(li->z, la->n);
3821      lu->z = nlSubLong(la->z, y);
3822      nlDeleteLint(y);
3823      lu->n = nlCopyLint(la->n);
3824    }
3825    else
3826    {
3827      x = nlMultLong(la->z, li->n);
3828      y = nlMultLong(li->z, la->n);
3829      lu->z = nlSubLong(x, y);
3830      nlDeleteLint(x);
3831      nlDeleteLint(y);
3832      lu->n = nlMultLong(la->n, li->n);
3833    }
3834  }
3835  if (lu->z==NULL)
3836  {
3837    nlDeleteLint(lu->n);
3838    Free((ADDRESS)lu,sizeof(rnumber));
3839    return NULL;
3840  }
3841  lu->s = -1;
3842#ifdef LDEBUG
3843  nlTest(lu);
3844#endif
3845  return lu;
3846}
3847
3848/*2
3849* lo := la * li
3850*/
3851number nlMult (number la, number li)
3852{
3853  number lo=(number)Alloc(sizeof(rnumber));
3854
3855  if (nlIsZero(la) || nlIsZero(li))
3856  {
3857    Free((ADDRESS)lo,sizeof(rnumber));
3858    return NULL;
3859  }
3860  lo->z = nlMultLong(la->z,li->z);
3861  if(NLISINT(la))
3862  {
3863    if(NLISINT(li))
3864    {
3865      lo->s = 0;
3866#ifdef LDEBUG
3867      nlTest(lo);
3868#endif
3869      return lo;
3870    }
3871    else
3872    {
3873      lo->n = nlCopyLint(li->n);
3874    }
3875  }
3876  else
3877  {
3878    if(NLISINT(li))
3879    {
3880      lo->n = nlCopyLint(la->n);
3881    }
3882    else
3883    {
3884      lo->n = nlMultLong(la->n,li->n);
3885    }
3886  }
3887  lo->s = -1;
3888#ifdef LDEBUG
3889  nlTest(lo);
3890#endif
3891  return lo;
3892}
3893
3894/*2
3895* lo := la / li in Z
3896* with |li| > la - (lo * li) >= 0
3897*/
3898number nlIntDiv (number la, number li)
3899{
3900  int c;
3901  number lo=(number)Alloc(sizeof(rnumber));
3902
3903  if (nlIsZero(la) || nlIsZero(li))
3904  {
3905    if (nlIsZero(li))
3906    {
3907      WerrorS("div. by 0");
3908    }
3909    Free((ADDRESS)lo,sizeof(rnumber));
3910    return NULL;
3911  }
3912  lo->s = 0;
3913  c = nlCompAbs(li->z, la->z);
3914  if (c >= 0)
3915  {
3916    if (c == 0)
3917    {
3918      lo->z = nlInitOne();
3919      if (NLSIGN(la->z) != NLSIGN(li->z))
3920      {
3921        lo->z[0] += SIGNUM;
3922      }
3923    }
3924    else if (NLSIGN(la->z))
3925    {
3926      lo->z = nlInitOne();
3927      if (!NLSIGN(li->z))
3928      {
3929        lo->z[0] += SIGNUM;
3930      }
3931    }
3932    else
3933    {
3934      Free((ADDRESS)lo,sizeof(rnumber));
3935      return NULL;
3936    }
3937#ifdef LDEBUG
3938    nlTest(lo);
3939#endif
3940    return lo;
3941  }
3942  lo->z = nlCopyLint(la->z);
3943  if (NLLENGTH(la->z)==1)
3944  {
3945    lo->z[1] /= li->z[1];
3946    if (NLSIGN(la->z))
3947    {
3948      int16 h = la->z[1] - (lo->z[1]*li->z[1]);
3949      if (h!=0)
3950      {
3951        lo->z[1]++;
3952      }
3953    }
3954    if (NLSIGN(li->z))
3955    {
3956      lo->z[0] = NLNEG(lo->z);
3957    }
3958#ifdef LDEBUG
3959    nlTest(lo);
3960#endif
3961    return lo;
3962  }
3963  if ((li->z[1]!=1) || (NLLENGTH(li->z)>1))
3964  {
3965    nlQuotient(&lo->z, li->z);
3966    if (NLSIGN(la->z))
3967    {
3968      lint o1, t, h = nlMultLong(lo->z, li->z);
3969      if (nlCompAbs(h, la->z))
3970      {
3971        o1 = nlInitOne();
3972        if (NLSIGN(li->z))
3973        {
3974          t = nlAddLong(lo->z,o1);
3975        }
3976        else
3977        {
3978          t = nlSubLong(lo->z,o1);
3979        }
3980        nlDeleteLint(o1);
3981        nlDeleteLint(lo->z);
3982        lo->z = t;
3983      }
3984      nlDeleteLint(h);
3985    }
3986  }
3987  else if (li->z[0] == (SIGNUM+1))
3988  {
3989    lo->z[0] = NLNEG(lo->z);
3990  }
3991#ifdef LDEBUG
3992  nlTest(lo);
3993#endif
3994  return lo;
3995}
3996
3997/*2
3998* lo := la mod li == la - ((la / li) * li)
3999* with |li| > lo >= 0
4000*/
4001number nlIntMod (number la, number li)
4002{
4003  int16 xl, yl;
4004  int c;
4005  lint r, t;
4006  number lo=(number)Alloc(sizeof(rnumber));
4007
4008  if (nlIsZero(la) || nlIsZero(li))
4009  {
4010    if (nlIsZero(li))
4011    {
4012      WerrorS("mod. by 0");
4013    }
4014    Free((ADDRESS)lo,sizeof(rnumber));
4015    return NULL;
4016  }
4017  c = nlCompAbs(li->z, la->z);
4018  if (c == 0)
4019  {
4020    Free((ADDRESS)lo,sizeof(rnumber));
4021    return NULL;
4022  }
4023  if (c > 0)
4024  {
4025    if (NLSIGN(la->z))
4026    {
4027      if (NLSIGN(li->z))
4028      {
4029        lo->z = nlSubLong(la->z,li->z);
4030      }
4031      else
4032      {
4033        lo->z = nlAddLong(la->z,li->z);
4034      }
4035    }
4036    else
4037    {
4038      lo->z = nlCopyLint(la->z);
4039    }
4040#ifdef LDEBUG
4041    nlTest(lo);
4042#endif
4043    return lo;
4044  }
4045  r = nlCopyLint(la->z);
4046  xl = NLLENGTH(r);
4047  if (xl == 1)
4048  {
4049    yl = r[1] / li->z[1];
4050    r[1] -= yl * li->z[1];
4051    if (r[1] != 0)
4052    {
4053      if (NLSIGN(la->z))
4054      {
4055        r[1] = li->z[1]-r[1];
4056      }
4057      lo->z = r;
4058#ifdef LDEBUG
4059      nlTest(lo);
4060#endif
4061      return lo;
4062    }
4063    else
4064    {
4065      nlDeleteLint(r);
4066      Free((ADDRESS)lo,sizeof(rnumber));
4067      return NULL;
4068    }
4069  }
4070  yl = NLLENGTH(li->z);
4071  if (yl > 1)
4072  {
4073    nlRemLong(r,li->z,xl,yl,&c);
4074    xl = c;
4075    if (xl != 0)
4076    {
4077      nlAllocLint(xl);
4078      for (c=xl; c>0; c--)
4079      {
4080        t[c] = r[c];
4081      }
4082    }
4083  }
4084  else if (li->z[1]!=1)
4085  {
4086    nlRemSmall(r,xl,li->z[1],&xl);
4087    if (xl != 0)
4088    {
4089      t = nlAllocLint(1);
4090      t[1] = xl;
4091    }
4092  }
4093  else
4094  {
4095    xl = 0;
4096  }
4097  nlDeleteLint(r);
4098  if (xl == 0)
4099  {
4100    Free((ADDRESS)lo,sizeof(rnumber));
4101    return NULL;
4102  }
4103  if (NLSIGN(la->z))
4104  {
4105    if (NLSIGN(li->z))
4106    {
4107      r = nlSubLong(t,li->z);
4108    }
4109    else
4110    {
4111      r = nlSubLong(li->z,t);
4112    }
4113    nlDeleteLint(t);
4114    t = r;
4115  }
4116  lo->z = t;
4117#ifdef LDEBUG
4118  nlTest(lo);
4119#endif
4120  return lo;
4121}
4122
4123/*2
4124* lo := la / li
4125*/
4126number nlDiv (number la, number li)
4127{
4128  number lo=(number)Alloc(sizeof(rnumber));
4129
4130  if (nlIsZero(la) || nlIsZero(li))
4131  {
4132    if (nlIsZero(li))
4133    {
4134      WerrorS("div. by 0");
4135    }
4136    Free((ADDRESS)lo,sizeof(rnumber));
4137    return NULL;
4138  }
4139  if (NLISINT(la))
4140  {
4141    lo->n = nlCopyLint(li->z);
4142  }
4143  else
4144  {
4145    lo->n = nlMultLong(la->n, li->z);
4146  }
4147  if (NLISINT(li))
4148  {
4149    lo->z = nlCopyLint(la->z);
4150  }
4151  else
4152  {
4153    lo->z = nlMultLong(la->z, li->n);
4154  }
4155  if(NLSIGN(lo->n))
4156  {
4157    lo->n[0] -= SIGNUM;
4158    lo->z[0] = NLNEG(lo->z);
4159  }
4160  lo->s = -1;
4161  if(nlIsOneLint(lo->n))
4162  {
4163    nlDeleteLint(lo->n);
4164    lo->s = 0;
4165  }
4166#ifdef LDEBUG
4167  nlTest(lo);
4168#endif
4169  return lo;
4170}
4171
4172/*3
4173* Power of long; lu:= x^exp
4174*/
4175static void nlPow (lint x,int exp,lint *lu)
4176{
4177  lint yh, y, w;
4178  y = nlInitOne();
4179  w = nlCopyLint(x);
4180  while (exp)
4181  {
4182    if (exp%2 != 0)
4183    {
4184      yh = nlMultLong(y,w);
4185      nlDeleteLint(y);
4186      y = yh;
4187    }
4188    if (exp > 1)
4189    {
4190      yh = nlMultLong(w,w);
4191      nlDeleteLint(w);
4192      w = yh;
4193    }
4194    exp = exp / 2;
4195  }
4196  *lu = y;
4197}
4198
4199/*2
4200* lu:= x ^ exp
4201*/
4202void nlPower (number x,int exp,number * lu)
4203{
4204  number b=NULL;
4205
4206  if (!nlIsZero(x))
4207  {
4208    b = (number)Alloc(sizeof(rnumber));
4209    if (exp == 0)
4210    {
4211      b->z = nlInitOne();
4212      *lu = b;
4213#ifdef LDEBUG
4214      nlTest(*lu);
4215#endif
4216      return;
4217    }
4218    nlNormalize(x);
4219    if (exp > 0)
4220    {
4221      b->s = x->s;
4222      if ((x->z[1]!=1) || (NLLENGTH(x->z)>1))
4223      {
4224        nlPow(x->z,exp,&b->z);
4225      }
4226      else
4227      {
4228        b->z = nlInitOne();
4229      }
4230      if (NLISRAT(x))
4231      {
4232        nlPow(x->n,exp,&b->n);
4233      }
4234    }
4235    else
4236    {
4237      if ((x->z[1]!=1) || (NLLENGTH(x->z)>1))
4238      {
4239        b->s = 1;
4240        nlPow(x->z,exp,&b->n);
4241        if (NLSIGN(b->n))
4242        {
4243          b->n[0] -= SIGNUM;
4244        }
4245      }
4246      else
4247      {
4248        b->s = 0;
4249      }
4250      if (NLISRAT(x))
4251      {
4252        nlPow(x->n,exp,&b->z);
4253      }
4254      else
4255      {
4256        b->z = nlInitOne();
4257      }
4258    }
4259    if (exp%2 != 0) // make sign
4260    {
4261      if(NLSIGN(b->z) != NLSIGN(x->z))
4262      {
4263        b->z[0] = NLNEG(b->z);
4264      }
4265    }
4266    else
4267    {
4268      if(NLSIGN(b->z))
4269      {
4270        b->z[0] -= SIGNUM;
4271      }
4272    }
4273  }
4274  *lu = b;
4275#ifdef LDEBUG
4276  nlTest(*lu);
4277#endif
4278}
4279
4280/*2
4281* a > b ?
4282*/
4283BOOLEAN nlGreater (number a, number b)
4284{
4285  number r;
4286  BOOLEAN rr;
4287
4288  r=nlSub(a,b);
4289  rr=(!nlIsZero(r)) && (nlGreaterZero(r));
4290  nlDelete(&r);
4291  return rr;
4292}
4293
4294/*3
4295* a := b * a + r
4296*/
4297void nlDivMod(lint a, int16 al, int16 b, int16 * r)
4298{
4299  int32 da;
4300  const int32 db = (int32)b;
4301  int16 a1 = a[al];
4302
4303  a[al] /= b;
4304  a1 -= b*a[al];
4305  loop
4306  {
4307    al--;
4308    if (al==0)
4309    {
4310      *r = a1;
4311      return;
4312    }
4313    while (a1==0)
4314    {
4315      a1 = a[al];
4316      a[al] /= b;
4317      a1 -= b*a[al];
4318      al--;
4319      if (al==0)
4320      {
4321        *r = a1;
4322        return;
4323      }
4324    }
4325    da = BASIS*(int32) a1 + (int32) a[al];
4326    a[al] = da / db;
4327    da -= db*a[al];
4328    a1 = (int16)da;
4329  }
4330}
4331
4332/*2
4333* returns (n->z mod p) / (n->n mod p)
4334*/
4335int nlModP(number n, int p)
4336{
4337  number a;
4338  int16 iz, in;
4339  int jz, jn, al;
4340
4341  if (nlIsZero(n))
4342  {
4343    return 0;
4344  }
4345  a = nCopy(n);
4346  al = NLLENGTH(a->z);
4347  if (al == 1)
4348  {
4349    iz = (a->z[1])%p;
4350  }
4351  else
4352  {
4353    nlRemSmall(a->z,al,p,&iz);
4354  }
4355  jz = iz;
4356  if (NLSIGN(a->z))
4357  {
4358    jz = p-jz;
4359  }
4360  if (NLISRAT(a))
4361  {
4362    al = *(a->n);
4363    if (al == 1)
4364    {
4365      in = (a->n[1])%p;
4366    }
4367    else
4368    {
4369      nlRemSmall(a->n,al,p,&in);
4370    }
4371    jn = in;
4372    jz = (int)npDiv((number)jz,(number)jn);
4373  }
4374  nlDelete(&a);
4375  return jz;
4376}
4377
4378int nlSize(number a)
4379{
4380  /*to be implemtented*/ return 0;
4381}
4382#endif
4383
Note: See TracBrowser for help on using the repository browser.