source: git/coeffs/longrat.cc @ c3d175

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