source: git/libpolys/coeffs/longrat.cc @ 73a9ffb

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