source: git/libpolys/coeffs/longrat.cc @ bd795d

spielwiese
Last change on this file since bd795d was fea2af, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
SW cannot be built without factory FIX: fixed libpolys and Singular/eigenval_ip.cc (quick and dirty) TODO: make check in libpolys/tests -> all dynamic tests are missing SearchPATH :( TODO: undefined reference to strat_fac_debug & strat_nr? TODO: kernel? Singular?
  • Property mode set to 100644
File size: 55.2 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 "config.h"
10#include <misc/auxiliary.h>
11
12#ifdef HAVE_FACTORY
13#include <factory/factory.h>
14#endif
15
16#include <coeffs/longrat.h>
17
18
19// 64 bit version:
20#if 0
21#define MAX_NUM_SIZE 60
22#define POW_2_28 0x10000000000000L
23#define LONG long
24#else
25#define MAX_NUM_SIZE 28
26#define POW_2_28 (1L<<28)
27#define LONG int
28#endif
29
30static inline number nlShort3(number x) // assume x->s==3
31{
32  assume(x->s==3);
33  if ((mpz_cmp_ui(x->z,(long)0)==0)
34  || (mpz_size1(x->z)<=MP_SMALL))
35  {
36    LONG ui=mpz_get_si(x->z);
37    if ((((ui<<3)>>3)==ui)
38    && (mpz_cmp_si(x->z,(long)ui)==0))
39    {
40      mpz_clear(x->z);
41      FREE_RNUMBER(x);
42      return INT_TO_SR(ui);
43    }
44  }
45  return x;
46}
47
48#ifndef LONGRAT_CC
49#define LONGRAT_CC
50
51#include <string.h>
52#include <float.h>
53
54#include <coeffs/coeffs.h>
55#include <reporter/reporter.h>
56#include <omalloc/omalloc.h>
57
58#include <coeffs/numbers.h>
59#include <coeffs/modulop.h>
60#include <coeffs/shortfl.h>
61#include <coeffs/mpr_complex.h>
62
63#ifndef BYTES_PER_MP_LIMB
64#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
65#endif
66
67/*-----------------------------------------------------------------*/
68/*3
69* parameter s in number:
70* 0 (or FALSE): not normalised rational
71* 1 (or TRUE):  normalised rational
72* 3          :  integer with n==NULL
73*/
74/*3
75**  'SR_INT' is the type of those integers small enough to fit into  29  bits.
76**  Therefor the value range of this small integers is: $-2^{28}...2^{28}-1$.
77**
78**  Small integers are represented by an immediate integer handle, containing
79**  the value instead of pointing  to  it,  which  has  the  following  form:
80**
81**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
82**      | guard | sign  | bit   | bit   |       | bit   | tag   | tag   |
83**      | bit   | bit   | 27    | 26    |       | 0     | 0     | 1     |
84**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
85**
86**  Immediate integers handles carry the tag 'SR_INT', i.e. the last bit is 1.
87**  This distuingishes immediate integers from other handles which  point  to
88**  structures aligned on 4 byte boundaries and therefor have last bit  zero.
89**  (The second bit is reserved as tag to allow extensions of  this  scheme.)
90**  Using immediates as pointers and dereferencing them gives address errors.
91**
92**  To aid overflow check the most significant two bits must always be equal,
93**  that is to say that the sign bit of immediate integers has a  guard  bit.
94**
95**  The macros 'INT_TO_SR' and 'SR_TO_INT' should be used to convert  between
96**  a small integer value and its representation as immediate integer handle.
97**
98**  Large integers and rationals are represented by z and n
99**  where n may be undefined (if s==3)
100**  NULL represents only deleted values
101*/
102#define SR_HDL(A) ((long)(A))
103/*#define SR_INT    1L*/
104/*#define INT_TO_SR(INT)  ((number) (((long)INT << 2) + SR_INT))*/
105// #define SR_TO_INT(SR)   (((long)SR) >> 2)
106
107#define MP_SMALL 1
108//#define mpz_isNeg(A) (mpz_cmp_si(A,(long)0)<0)
109#define mpz_isNeg(A) ((A)->_mp_size<0)
110#define mpz_limb_size(A) ((A)->_mp_size)
111#define mpz_limb_d(A) ((A)->_mp_d)
112#define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
113#define MPZ_EXACTDIV(A,B,C) mpz_divexact((A),(B),(C))
114
115/// Our Type!
116static const n_coeffType ID = n_Q;
117
118void    _nlDelete_NoImm(number *a);
119
120/***************************************************************
121 *
122 * Routines which are never inlined by p_Numbers.h
123 *
124 *******************************************************************/
125#ifndef P_NUMBERS_H
126
127number nlShort3_noinline(number x) // assume x->s==3
128{
129  return nlShort3(x);
130}
131
132
133number nlOne=INT_TO_SR(1);
134
135#if (__GNU_MP_VERSION*10+__GNU_MP_VERSION_MINOR < 31)
136void mpz_mul_si (mpz_ptr r, mpz_srcptr s, long int si)
137{
138  if (si>=0)
139    mpz_mul_ui(r,s,si);
140  else
141  {
142    mpz_mul_ui(r,s,-si);
143    mpz_neg(r,r);
144  }
145}
146#endif
147
148static number nlMapP(number from, const coeffs src, const coeffs dst)
149{
150  assume( getCoeffType(dst) == ID );
151  assume( getCoeffType(src) == n_Zp );
152
153  number to;
154  to = nlInit(npInt(from,src), dst);
155  return to;
156}
157
158static number nlMapLongR(number from, const coeffs src, const coeffs dst);
159static number nlMapR(number from, const coeffs src, const coeffs dst);
160
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      mpz_t 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      mpz_t 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 = ALLOC_RNUMBER(); //Q!? // (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    FREE_RNUMBER(bb); // 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    FREE_RNUMBER(bb); // 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    FREE_RNUMBER(bb); // 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      FREE_RNUMBER(u); // 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
1258// Map q \in QQ \to Zp
1259// src = Q, dst = Zp (or an extension of Zp?)
1260number nlModP(number q, const coeffs Q, const coeffs Zp)
1261{
1262  assume( getCoeffType(Q) == ID );
1263
1264  const int p = n_GetChar(Zp);
1265  assume( p > 0 );
1266
1267  const long P = p;
1268  assume( P > 0 );
1269
1270  // embedded long within q => only long numerator has to be converted
1271  // to int (modulo char.)
1272  if (SR_HDL(q) & SR_INT)
1273  {
1274    long i = SR_TO_INT(q);
1275    if (i<0L)
1276      return n_Init( static_cast<int>( P - ((-i)%P) ), Zp);
1277    else
1278      return n_Init( static_cast<int>( i % P ), Zp );
1279  }
1280
1281  const unsigned long PP = p;
1282
1283  // numerator modulo char. should fit into int
1284  number z = n_Init( static_cast<int>(mpz_fdiv_ui(q->z, PP)), Zp ); 
1285
1286  // denominator != 1?
1287  if (q->s!=3)
1288  {
1289    // denominator modulo char. should fit into int
1290    number n = n_Init( static_cast<int>(mpz_fdiv_ui(q->n, PP)), Zp );
1291
1292    number res = n_Div( z, n, Zp );
1293
1294    n_Delete(&z, Zp);
1295    n_Delete(&n, Zp);
1296
1297    return res;
1298  }
1299
1300  return z;
1301}
1302
1303#ifdef HAVE_RINGS
1304/*2
1305* convert number i (from Q) to GMP and warn if denom != 1
1306*/
1307void nlGMP(number &i, number n, const coeffs r)
1308{
1309  // Hier brauche ich einfach die GMP Zahl
1310  nlTest(i, r);
1311  nlNormalize(i, r);
1312  if (SR_HDL(i) & SR_INT)
1313  {
1314    mpz_set_si((mpz_ptr) n, SR_TO_INT(i));
1315    return;
1316  }
1317  if (i->s!=3)
1318  {
1319     WarnS("Omitted denominator during coefficient mapping !");
1320  }
1321  mpz_set((mpz_ptr) n, i->z);
1322}
1323#endif
1324
1325/*2
1326* acces to denominator, other 1 for integers
1327*/
1328number   nlGetDenom(number &n, const coeffs r)
1329{
1330  if (!(SR_HDL(n) & SR_INT))
1331  {
1332    if (n->s==0)
1333    {
1334      nlNormalize(n,r);
1335    }
1336    if (!(SR_HDL(n) & SR_INT))
1337    {
1338      if (n->s!=3)
1339      {
1340        number u=ALLOC_RNUMBER();
1341        u->s=3;
1342#if defined(LDEBUG)
1343        u->debug=123456;
1344#endif
1345        mpz_init_set(u->z,n->n);
1346        u=nlShort3_noinline(u);
1347        return u;
1348      }
1349    }
1350  }
1351  return INT_TO_SR(1);
1352}
1353
1354/*2
1355* acces to Nominator, nlCopy(n) for integers
1356*/
1357number   nlGetNumerator(number &n, const coeffs r)
1358{
1359  if (!(SR_HDL(n) & SR_INT))
1360  {
1361    if (n->s==0)
1362    {
1363      nlNormalize(n,r);
1364    }
1365    if (!(SR_HDL(n) & SR_INT))
1366    {
1367      number u=ALLOC_RNUMBER();
1368#if defined(LDEBUG)
1369      u->debug=123456;
1370#endif
1371      u->s=3;
1372      mpz_init_set(u->z,n->z);
1373      if (n->s!=3)
1374      {
1375        u=nlShort3_noinline(u);
1376      }
1377      return u;
1378    }
1379  }
1380  return n; // imm. int
1381}
1382
1383/***************************************************************
1384 *
1385 * routines which are needed by Inline(d) routines
1386 *
1387 *******************************************************************/
1388BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
1389{
1390  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1391//  long - short
1392  BOOLEAN bo;
1393  if (SR_HDL(b) & SR_INT)
1394  {
1395    if (a->s!=0) return FALSE;
1396    number n=b; b=a; a=n;
1397  }
1398//  short - long
1399  if (SR_HDL(a) & SR_INT)
1400  {
1401    if (b->s!=0)
1402      return FALSE;
1403    if (((long)a > 0L) && (mpz_isNeg(b->z)))
1404      return FALSE;
1405    if (((long)a < 0L) && (!mpz_isNeg(b->z)))
1406      return FALSE;
1407    mpz_t  bb;
1408    mpz_init_set(bb,b->n);
1409    mpz_mul_si(bb,bb,(long)SR_TO_INT(a));
1410    bo=(mpz_cmp(bb,b->z)==0);
1411    mpz_clear(bb);
1412    return bo;
1413  }
1414// long - long
1415  if (((a->s==1) && (b->s==3))
1416  ||  ((b->s==1) && (a->s==3)))
1417    return FALSE;
1418  if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
1419    return FALSE;
1420  if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
1421    return FALSE;
1422  mpz_t  aa;
1423  mpz_t  bb;
1424  mpz_init_set(aa,a->z);
1425  mpz_init_set(bb,b->z);
1426  if (a->s<2) mpz_mul(bb,bb,a->n);
1427  if (b->s<2) mpz_mul(aa,aa,b->n);
1428  bo=(mpz_cmp(aa,bb)==0);
1429  mpz_clear(aa);
1430  mpz_clear(bb);
1431  return bo;
1432}
1433
1434// copy not immediate number a
1435number _nlCopy_NoImm(number a)
1436{
1437  assume(!((SR_HDL(a) & SR_INT)||(a==NULL)));
1438  //nlTest(a, r);
1439  number b=ALLOC_RNUMBER();
1440#if defined(LDEBUG)
1441  b->debug=123456;
1442#endif
1443  switch (a->s)
1444  {
1445    case 0:
1446    case 1:
1447            mpz_init_set(b->n,a->n);
1448    case 3:
1449            mpz_init_set(b->z,a->z);
1450            break;
1451  }
1452  b->s = a->s;
1453  return b;
1454}
1455
1456void _nlDelete_NoImm(number *a)
1457{
1458  {
1459    switch ((*a)->s)
1460    {
1461      case 0:
1462      case 1:
1463        mpz_clear((*a)->n);
1464      case 3:
1465        mpz_clear((*a)->z);
1466#ifdef LDEBUG
1467        (*a)->s=2;
1468#endif
1469    }
1470    FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
1471  }
1472}
1473
1474number _nlNeg_NoImm(number a)
1475{
1476  {
1477    mpz_neg(a->z,a->z);
1478    if (a->s==3)
1479    {
1480      a=nlShort3(a);
1481    }
1482  }
1483  return a;
1484}
1485
1486number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1487{
1488  number u=ALLOC_RNUMBER();
1489#if defined(LDEBUG)
1490  u->debug=123456;
1491#endif
1492  mpz_init(u->z);
1493  if (SR_HDL(b) & SR_INT)
1494  {
1495    number x=a;
1496    a=b;
1497    b=x;
1498  }
1499  if (SR_HDL(a) & SR_INT)
1500  {
1501    switch (b->s)
1502    {
1503      case 0:
1504      case 1:/* a:short, b:1 */
1505      {
1506        mpz_t x;
1507        mpz_init(x);
1508        mpz_mul_si(x,b->n,SR_TO_INT(a));
1509        mpz_add(u->z,b->z,x);
1510        mpz_clear(x);
1511        if (mpz_cmp_ui(u->z,(long)0)==0)
1512        {
1513          mpz_clear(u->z);
1514          FREE_RNUMBER(u);
1515          return INT_TO_SR(0);
1516        }
1517        if (mpz_cmp(u->z,b->n)==0)
1518        {
1519          mpz_clear(u->z);
1520          FREE_RNUMBER(u);
1521          return INT_TO_SR(1);
1522        }
1523        mpz_init_set(u->n,b->n);
1524        u->s = 0;
1525        break;
1526      }
1527      case 3:
1528      {
1529        if ((long)a>0L)
1530          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1531        else
1532          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1533        if (mpz_cmp_ui(u->z,(long)0)==0)
1534        {
1535          mpz_clear(u->z);
1536          FREE_RNUMBER(u);
1537          return INT_TO_SR(0);
1538        }
1539        u->s = 3;
1540        u=nlShort3(u);
1541        break;
1542      }
1543    }
1544  }
1545  else
1546  {
1547    switch (a->s)
1548    {
1549      case 0:
1550      case 1:
1551      {
1552        switch(b->s)
1553        {
1554          case 0:
1555          case 1:
1556          {
1557            mpz_t x;
1558            mpz_init(x);
1559
1560            mpz_mul(x,b->z,a->n);
1561            mpz_mul(u->z,a->z,b->n);
1562            mpz_add(u->z,u->z,x);
1563            mpz_clear(x);
1564
1565            if (mpz_cmp_ui(u->z,(long)0)==0)
1566            {
1567              mpz_clear(u->z);
1568              FREE_RNUMBER(u);
1569              return INT_TO_SR(0);
1570            }
1571            mpz_init(u->n);
1572            mpz_mul(u->n,a->n,b->n);
1573            if (mpz_cmp(u->z,u->n)==0)
1574            {
1575               mpz_clear(u->z);
1576               mpz_clear(u->n);
1577               FREE_RNUMBER(u);
1578               return INT_TO_SR(1);
1579            }
1580            u->s = 0;
1581            break;
1582          }
1583          case 3: /* a:1 b:3 */
1584          {
1585            mpz_mul(u->z,b->z,a->n);
1586            mpz_add(u->z,u->z,a->z);
1587            if (mpz_cmp_ui(u->z,(long)0)==0)
1588            {
1589              mpz_clear(u->z);
1590              FREE_RNUMBER(u);
1591              return INT_TO_SR(0);
1592            }
1593            if (mpz_cmp(u->z,a->n)==0)
1594            {
1595              mpz_clear(u->z);
1596              FREE_RNUMBER(u);
1597              return INT_TO_SR(1);
1598            }
1599            mpz_init_set(u->n,a->n);
1600            u->s = 0;
1601            break;
1602          }
1603        } /*switch (b->s) */
1604        break;
1605      }
1606      case 3:
1607      {
1608        switch(b->s)
1609        {
1610          case 0:
1611          case 1:/* a:3, b:1 */
1612          {
1613            mpz_mul(u->z,a->z,b->n);
1614            mpz_add(u->z,u->z,b->z);
1615            if (mpz_cmp_ui(u->z,(long)0)==0)
1616            {
1617              mpz_clear(u->z);
1618              FREE_RNUMBER(u);
1619              return INT_TO_SR(0);
1620            }
1621            if (mpz_cmp(u->z,b->n)==0)
1622            {
1623              mpz_clear(u->z);
1624              FREE_RNUMBER(u);
1625              return INT_TO_SR(1);
1626            }
1627            mpz_init_set(u->n,b->n);
1628            u->s = 0;
1629            break;
1630          }
1631          case 3:
1632          {
1633            mpz_add(u->z,a->z,b->z);
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            u->s = 3;
1641            u=nlShort3(u);
1642            break;
1643          }
1644        }
1645        break;
1646      }
1647    }
1648  }
1649  return u;
1650}
1651
1652number _nlSub_aNoImm_OR_bNoImm(number a, number b)
1653{
1654  number u=ALLOC_RNUMBER();
1655#if defined(LDEBUG)
1656  u->debug=123456;
1657#endif
1658  mpz_init(u->z);
1659  if (SR_HDL(a) & SR_INT)
1660  {
1661    switch (b->s)
1662    {
1663      case 0:
1664      case 1:/* a:short, b:1 */
1665      {
1666        mpz_t x;
1667        mpz_init(x);
1668        mpz_mul_si(x,b->n,SR_TO_INT(a));
1669        mpz_sub(u->z,x,b->z);
1670        mpz_clear(x);
1671        if (mpz_cmp_ui(u->z,(long)0)==0)
1672        {
1673          mpz_clear(u->z);
1674          FREE_RNUMBER(u);
1675          return INT_TO_SR(0);
1676        }
1677        if (mpz_cmp(u->z,b->n)==0)
1678        {
1679          mpz_clear(u->z);
1680          FREE_RNUMBER(u);
1681          return INT_TO_SR(1);
1682        }
1683        mpz_init_set(u->n,b->n);
1684        u->s = 0;
1685        break;
1686      }
1687      case 3:
1688      {
1689        if ((long)a>0L)
1690        {
1691          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
1692          mpz_neg(u->z,u->z);
1693        }
1694        else
1695        {
1696          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
1697          mpz_neg(u->z,u->z);
1698        }
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        u->s = 3;
1706        u=nlShort3(u);
1707        break;
1708      }
1709    }
1710  }
1711  else if (SR_HDL(b) & SR_INT)
1712  {
1713    switch (a->s)
1714    {
1715      case 0:
1716      case 1:/* b:short, a:1 */
1717      {
1718        mpz_t x;
1719        mpz_init(x);
1720        mpz_mul_si(x,a->n,SR_TO_INT(b));
1721        mpz_sub(u->z,a->z,x);
1722        mpz_clear(x);
1723        if (mpz_cmp_ui(u->z,(long)0)==0)
1724        {
1725          mpz_clear(u->z);
1726          FREE_RNUMBER(u);
1727          return INT_TO_SR(0);
1728        }
1729        if (mpz_cmp(u->z,a->n)==0)
1730        {
1731          mpz_clear(u->z);
1732          FREE_RNUMBER(u);
1733          return INT_TO_SR(1);
1734        }
1735        mpz_init_set(u->n,a->n);
1736        u->s = 0;
1737        break;
1738      }
1739      case 3:
1740      {
1741        if ((long)b>0L)
1742        {
1743          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
1744        }
1745        else
1746        {
1747          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
1748        }
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        u->s = 3;
1756        u=nlShort3(u);
1757        break;
1758      }
1759    }
1760  }
1761  else
1762  {
1763    switch (a->s)
1764    {
1765      case 0:
1766      case 1:
1767      {
1768        switch(b->s)
1769        {
1770          case 0:
1771          case 1:
1772          {
1773            mpz_t x;
1774            mpz_t y;
1775            mpz_init(x);
1776            mpz_init(y);
1777            mpz_mul(x,b->z,a->n);
1778            mpz_mul(y,a->z,b->n);
1779            mpz_sub(u->z,y,x);
1780            mpz_clear(x);
1781            mpz_clear(y);
1782            if (mpz_cmp_ui(u->z,(long)0)==0)
1783            {
1784              mpz_clear(u->z);
1785              FREE_RNUMBER(u);
1786              return INT_TO_SR(0);
1787            }
1788            mpz_init(u->n);
1789            mpz_mul(u->n,a->n,b->n);
1790            if (mpz_cmp(u->z,u->n)==0)
1791            {
1792              mpz_clear(u->z);
1793              mpz_clear(u->n);
1794              FREE_RNUMBER(u);
1795              return INT_TO_SR(1);
1796            }
1797            u->s = 0;
1798            break;
1799          }
1800          case 3: /* a:1, b:3 */
1801          {
1802            mpz_t x;
1803            mpz_init(x);
1804            mpz_mul(x,b->z,a->n);
1805            mpz_sub(u->z,a->z,x);
1806            mpz_clear(x);
1807            if (mpz_cmp_ui(u->z,(long)0)==0)
1808            {
1809              mpz_clear(u->z);
1810              FREE_RNUMBER(u);
1811              return INT_TO_SR(0);
1812            }
1813            if (mpz_cmp(u->z,a->n)==0)
1814            {
1815              mpz_clear(u->z);
1816              FREE_RNUMBER(u);
1817              return INT_TO_SR(1);
1818            }
1819            mpz_init_set(u->n,a->n);
1820            u->s = 0;
1821            break;
1822          }
1823        }
1824        break;
1825      }
1826      case 3:
1827      {
1828        switch(b->s)
1829        {
1830          case 0:
1831          case 1: /* a:3, b:1 */
1832          {
1833            mpz_t x;
1834            mpz_init(x);
1835            mpz_mul(x,a->z,b->n);
1836            mpz_sub(u->z,x,b->z);
1837            mpz_clear(x);
1838            if (mpz_cmp_ui(u->z,(long)0)==0)
1839            {
1840              mpz_clear(u->z);
1841              FREE_RNUMBER(u);
1842              return INT_TO_SR(0);
1843            }
1844            if (mpz_cmp(u->z,b->n)==0)
1845            {
1846              mpz_clear(u->z);
1847              FREE_RNUMBER(u);
1848              return INT_TO_SR(1);
1849            }
1850            mpz_init_set(u->n,b->n);
1851            u->s = 0;
1852            break;
1853          }
1854          case 3: /* a:3 , b:3 */
1855          {
1856            mpz_sub(u->z,a->z,b->z);
1857            if (mpz_cmp_ui(u->z,(long)0)==0)
1858            {
1859              mpz_clear(u->z);
1860              FREE_RNUMBER(u);
1861              return INT_TO_SR(0);
1862            }
1863            u->s = 3;
1864            u=nlShort3(u);
1865            break;
1866          }
1867        }
1868        break;
1869      }
1870    }
1871  }
1872  return u;
1873}
1874
1875// a and b are intermediate, but a*b not
1876number _nlMult_aImm_bImm_rNoImm(number a, number b)
1877{
1878  number u=ALLOC_RNUMBER();
1879#if defined(LDEBUG)
1880  u->debug=123456;
1881#endif
1882  u->s=3;
1883  mpz_init_set_si(u->z,SR_TO_INT(a));
1884  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
1885  return u;
1886}
1887
1888// a or b are not immediate
1889number _nlMult_aNoImm_OR_bNoImm(number a, number b)
1890{
1891  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1892  number u=ALLOC_RNUMBER();
1893#if defined(LDEBUG)
1894  u->debug=123456;
1895#endif
1896  mpz_init(u->z);
1897  if (SR_HDL(b) & SR_INT)
1898  {
1899    number x=a;
1900    a=b;
1901    b=x;
1902  }
1903  if (SR_HDL(a) & SR_INT)
1904  {
1905    u->s=b->s;
1906    if (u->s==1) u->s=0;
1907    if ((long)a>0L)
1908    {
1909      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
1910    }
1911    else
1912    {
1913      if (a==INT_TO_SR(-1))
1914      {
1915        mpz_set(u->z,b->z);
1916        mpz_neg(u->z,u->z);
1917        u->s=b->s;
1918      }
1919      else
1920      {
1921        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
1922        mpz_neg(u->z,u->z);
1923      }
1924    }
1925    if (u->s<2)
1926    {
1927      if (mpz_cmp(u->z,b->n)==0)
1928      {
1929        mpz_clear(u->z);
1930        FREE_RNUMBER(u);
1931        return INT_TO_SR(1);
1932      }
1933      mpz_init_set(u->n,b->n);
1934    }
1935    else //u->s==3
1936    {
1937      u=nlShort3(u);
1938    }
1939  }
1940  else
1941  {
1942    mpz_mul(u->z,a->z,b->z);
1943    u->s = 0;
1944    if(a->s==3)
1945    {
1946      if(b->s==3)
1947      {
1948        u->s = 3;
1949      }
1950      else
1951      {
1952        if (mpz_cmp(u->z,b->n)==0)
1953        {
1954          mpz_clear(u->z);
1955          FREE_RNUMBER(u);
1956          return INT_TO_SR(1);
1957        }
1958        mpz_init_set(u->n,b->n);
1959      }
1960    }
1961    else
1962    {
1963      if(b->s==3)
1964      {
1965        if (mpz_cmp(u->z,a->n)==0)
1966        {
1967          mpz_clear(u->z);
1968          FREE_RNUMBER(u);
1969          return INT_TO_SR(1);
1970        }
1971        mpz_init_set(u->n,a->n);
1972      }
1973      else
1974      {
1975        mpz_init(u->n);
1976        mpz_mul(u->n,a->n,b->n);
1977        if (mpz_cmp(u->z,u->n)==0)
1978        {
1979          mpz_clear(u->z);
1980          mpz_clear(u->n);
1981          FREE_RNUMBER(u);
1982          return INT_TO_SR(1);
1983        }
1984      }
1985    }
1986  }
1987  return u;
1988}
1989
1990/*2
1991* copy a to b for mapping
1992*/
1993number nlCopyMap(number a, const coeffs src, const coeffs dst)
1994{
1995  assume( getCoeffType(dst) == ID );
1996  assume( getCoeffType(src) == ID );
1997
1998  if ((SR_HDL(a) & SR_INT)||(a==NULL))
1999  {
2000    return a;
2001  }
2002  return _nlCopy_NoImm(a);
2003}
2004
2005nMapFunc nlSetMap(const coeffs src, const coeffs dst)
2006{
2007  assume( getCoeffType(dst) == ID );
2008//  assume( getCoeffType(src) == ID );
2009
2010  if (nCoeff_is_Q(src))
2011  {
2012    return ndCopyMap;
2013  }
2014  if (nCoeff_is_Zp(src))
2015  {
2016    return nlMapP;
2017  }
2018  if (nCoeff_is_R(src))
2019  {
2020    return nlMapR;
2021  }
2022  if (nCoeff_is_long_R(src))
2023  {
2024    return nlMapLongR; /* long R -> Q */
2025  }
2026#ifdef HAVE_RINGS
2027  if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
2028  {
2029    return nlMapGMP;
2030  }
2031  if (nCoeff_is_Ring_2toM(src))
2032  {
2033    return nlMapMachineInt;
2034  }
2035#endif
2036  return NULL;
2037}
2038/*2
2039* z := i
2040*/
2041number nlRInit (long i)
2042{
2043  number z=ALLOC_RNUMBER();
2044#if defined(LDEBUG)
2045  z->debug=123456;
2046#endif
2047  mpz_init_set_si(z->z,i);
2048  z->s = 3;
2049  return z;
2050}
2051
2052/*2
2053* z := i/j
2054*/
2055number nlInit2 (int i, int j, const coeffs r)
2056{
2057  number z=ALLOC_RNUMBER();
2058#if defined(LDEBUG)
2059  z->debug=123456;
2060#endif
2061  mpz_init_set_si(z->z,(long)i);
2062  mpz_init_set_si(z->n,(long)j);
2063  z->s = 0;
2064  nlNormalize(z,r);
2065  return z;
2066}
2067
2068number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
2069{
2070  number z=ALLOC_RNUMBER();
2071#if defined(LDEBUG)
2072  z->debug=123456;
2073#endif
2074  mpz_init_set(z->z,i);
2075  mpz_init_set(z->n,j);
2076  z->s = 0;
2077  nlNormalize(z,r);
2078  return z;
2079}
2080
2081#else // DO_LINLINE
2082
2083// declare immedate routines
2084number nlRInit (int i);
2085BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2086number  _nlCopy_NoImm(number a);
2087number  _nlNeg_NoImm(number a);
2088number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2089number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2090number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2091number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2092
2093#endif
2094
2095/***************************************************************
2096 *
2097 * Routines which might be inlined by p_Numbers.h
2098 *
2099 *******************************************************************/
2100#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2101
2102// routines which are always inlined/static
2103
2104/*2
2105* a = b ?
2106*/
2107LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
2108{
2109  nlTest(a, r);
2110  nlTest(b, r);
2111// short - short
2112  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2113  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2114}
2115
2116LINLINE number nlInit (int i, const coeffs r)
2117{
2118  number n;
2119  LONG ii=(LONG)i;
2120  if ( ((ii << 3) >> 3) == ii ) n=INT_TO_SR(ii);
2121  else                          n=nlRInit(ii);
2122  nlTest(n, r);
2123  return n;
2124}
2125
2126/*2
2127* a == 1 ?
2128*/
2129LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
2130{
2131#ifdef LDEBUG
2132  if (a==NULL) return FALSE;
2133  nlTest(a, r);
2134#endif
2135  return (a==INT_TO_SR(1));
2136}
2137
2138LINLINE BOOLEAN nlIsZero (number a, const coeffs)
2139{
2140  return (a==INT_TO_SR(0));
2141  //return (mpz_cmp_si(a->z,(long)0)==0);
2142}
2143
2144/*2
2145* copy a to b
2146*/
2147LINLINE number nlCopy(number a, const coeffs)
2148{
2149  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2150  {
2151    return a;
2152  }
2153  return _nlCopy_NoImm(a);
2154}
2155
2156
2157/*2
2158* delete a
2159*/
2160LINLINE void nlDelete (number * a, const coeffs r)
2161{
2162  if (*a!=NULL)
2163  {
2164    nlTest(*a, r);
2165    if ((SR_HDL(*a) & SR_INT)==0)
2166    {
2167      _nlDelete_NoImm(a);
2168    }
2169    *a=NULL;
2170  }
2171}
2172
2173/*2
2174* za:= - za
2175*/
2176LINLINE number nlNeg (number a, const coeffs R)
2177{
2178  nlTest(a, R);
2179  if(SR_HDL(a) &SR_INT)
2180  {
2181    int r=SR_TO_INT(a);
2182    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2183    else               a=INT_TO_SR(-r);
2184    return a;
2185  }
2186  a = _nlNeg_NoImm(a);
2187  nlTest(a, R);
2188  return a;
2189
2190}
2191
2192/*2
2193* u:= a + b
2194*/
2195LINLINE number nlAdd (number a, number b, const coeffs R)
2196{
2197  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2198  {
2199    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2200    if ( ((r << 1) >> 1) == r )
2201      return (number)(long)r;
2202    else
2203      return nlRInit(SR_TO_INT(r));
2204  }
2205  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2206  nlTest(u, R);
2207  return u;
2208}
2209
2210number nlShort1(number a);
2211number nlShort3_noinline(number x);
2212
2213static inline number nlInpAdd_(number a, number b, const coeffs r)
2214{
2215  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2216  {
2217    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2218    if ( ((r << 1) >> 1) == r )
2219      return (number)(long)r;
2220    else
2221      return nlRInit(SR_TO_INT(r));
2222  }
2223  // a=a+b
2224  if (SR_HDL(b) & SR_INT)
2225  {
2226    switch (a->s)
2227    {
2228      case 0:
2229      case 1:/* b:short, a:1 */
2230      {
2231        mpz_t x;
2232        mpz_init(x);
2233        mpz_mul_si(x,a->n,SR_TO_INT(b));
2234        mpz_add(a->z,a->z,x);
2235        mpz_clear(x);
2236        a->s = 0;
2237        a=nlShort1(a);
2238        break;
2239      }
2240      case 3:
2241      {
2242        if ((long)b>0L)
2243          mpz_add_ui(a->z,a->z,SR_TO_INT(b));
2244        else
2245          mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
2246        a->s = 3;
2247        a=nlShort3_noinline(a);
2248        break;
2249      }
2250    }
2251    return a;
2252  }
2253  if (SR_HDL(a) & SR_INT)
2254  {
2255    number u=ALLOC_RNUMBER();
2256    #if defined(LDEBUG)
2257    u->debug=123456;
2258    #endif
2259    mpz_init(u->z);
2260    switch (b->s)
2261    {
2262      case 0:
2263      case 1:/* a:short, b:1 */
2264      {
2265        mpz_t x;
2266        mpz_init(x);
2267
2268        mpz_mul_si(x,b->n,SR_TO_INT(a));
2269        mpz_add(u->z,b->z,x);
2270        mpz_clear(x);
2271        // result cannot be 0, if coeffs are normalized
2272        mpz_init_set(u->n,b->n);
2273        u->s = 0;
2274        u=nlShort1(u);
2275        break;
2276      }
2277      case 3:
2278      {
2279        if ((long)a>0L)
2280          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
2281        else
2282          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
2283        // result cannot be 0, if coeffs are normalized
2284        u->s = 3;
2285        u=nlShort3_noinline(u);
2286        break;
2287      }
2288    }
2289    nlTest(u, r);
2290    return u;
2291  }
2292  else
2293  {
2294    switch (a->s)
2295    {
2296      case 0:
2297      case 1:
2298      {
2299        switch(b->s)
2300        {
2301          case 0:
2302          case 1: /* a:1 b:1 */
2303          {
2304            mpz_t x;
2305            mpz_t y;
2306            mpz_init(x);
2307            mpz_init(y);
2308            mpz_mul(x,b->z,a->n);
2309            mpz_mul(y,a->z,b->n);
2310            mpz_add(a->z,x,y);
2311            mpz_clear(x);
2312            mpz_clear(y);
2313            mpz_mul(a->n,a->n,b->n);
2314            a->s = 0;
2315            break;
2316          }
2317          case 3: /* a:1 b:3 */
2318          {
2319            mpz_t x;
2320            mpz_init(x);
2321            mpz_mul(x,b->z,a->n);
2322            mpz_add(a->z,a->z,x);
2323            mpz_clear(x);
2324            a->s = 0;
2325            break;
2326          }
2327        } /*switch (b->s) */
2328        a=nlShort1(a);
2329        break;
2330      }
2331      case 3:
2332      {
2333        switch(b->s)
2334        {
2335          case 0:
2336          case 1:/* a:3, b:1 */
2337          {
2338            mpz_t x;
2339            mpz_init(x);
2340            mpz_mul(x,a->z,b->n);
2341            mpz_add(a->z,b->z,x);
2342            mpz_clear(x);
2343            mpz_init_set(a->n,b->n);
2344            a->s = 0;
2345            a=nlShort1(a);
2346            break;
2347          }
2348          case 3:
2349          {
2350            mpz_add(a->z,a->z,b->z);
2351            a->s = 3;
2352            a=nlShort3_noinline(a);
2353            break;
2354          }
2355        }
2356        break;
2357      }
2358    }
2359    nlTest(a, r);
2360    return a;
2361  }
2362}
2363
2364LINLINE void nlInpAdd(number &a, number b, const coeffs r)
2365{
2366  a = nlInpAdd_(a, b, r);
2367}
2368
2369
2370LINLINE number nlMult (number a, number b, const coeffs R)
2371{
2372  nlTest(a, R);
2373  nlTest(b, R);
2374  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2375  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2376  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2377  {
2378    LONG r=(SR_HDL(a)-1L)*(SR_HDL(b)>>1);
2379    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2380    {
2381      number u=((number) ((r>>1)+SR_INT));
2382      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2383      return nlRInit(SR_HDL(u)>>2);
2384    }
2385    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2386    nlTest(u, R);
2387    return u;
2388
2389  }
2390  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2391  nlTest(u, R);
2392  return u;
2393
2394}
2395
2396
2397/*2
2398* u:= a - b
2399*/
2400LINLINE number nlSub (number a, number b, const coeffs r)
2401{
2402  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2403  {
2404    LONG r=SR_HDL(a)-SR_HDL(b)+1;
2405    if ( ((r << 1) >> 1) == r )
2406    {
2407      return (number)(long)r;
2408    }
2409    else
2410      return nlRInit(SR_TO_INT(r));
2411  }
2412  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2413  nlTest(u, r);
2414  return u;
2415
2416}
2417
2418LINLINE void nlInpMult(number &a, number b, const coeffs r)
2419{
2420  number aa=a;
2421  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2422  {
2423    number n=nlMult(aa,b,r);
2424    nlDelete(&a,r);
2425    a=n;
2426  }
2427  else
2428  {
2429    mpz_mul(aa->z,a->z,b->z);
2430    if (aa->s==3)
2431    {
2432      if(b->s!=3)
2433      {
2434        mpz_init_set(a->n,b->n);
2435        a->s=0;
2436      }
2437    }
2438    else
2439    {
2440      if(b->s!=3)
2441      {
2442        mpz_mul(a->n,a->n,b->n);
2443      }
2444      a->s=0;
2445    }
2446  }
2447}
2448#endif // DO_LINLINE
2449
2450#ifndef P_NUMBERS_H
2451
2452static void nlMPZ(mpz_t m, number &n, const coeffs r)
2453{
2454  assume( getCoeffType(r) == ID );
2455   
2456  nlTest(n, r);
2457  nlNormalize(n, r);
2458  if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n));     /* n fits in an int */
2459  else             mpz_init_set(m, (mpz_ptr)n->z);
2460}
2461
2462
2463static number nlInitMPZ(mpz_t m, const coeffs)
2464{
2465  number z = ALLOC_RNUMBER();
2466  mpz_init_set(z->z, m);
2467  mpz_init_set_ui(z->n, 1);
2468  z->s = 3;
2469  return z;   
2470}
2471
2472
2473void nlInpGcd(number &a, number b, const coeffs r)
2474{
2475  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2476  {
2477    number n=nlGcd(a,b,r);
2478    nlDelete(&a,r);
2479    a=n;
2480  }
2481  else
2482  {
2483    mpz_gcd(a->z,a->z,b->z);
2484    a=nlShort3_noinline(a);
2485  }
2486}
2487
2488void nlInpIntDiv(number &a, number b, const coeffs r)
2489{
2490  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2491  {
2492    number n=nlIntDiv(a,b, r);
2493    nlDelete(&a,r);
2494    a=n;
2495  }
2496  else
2497  {
2498    if (mpz_isNeg(a->z))
2499    {
2500      if (mpz_isNeg(b->z))
2501      {
2502        mpz_add(a->z,a->z,b->z);
2503      }
2504      else
2505      {
2506        mpz_sub(a->z,a->z,b->z);
2507      }
2508      mpz_add_ui(a->z,a->z,1);
2509    }
2510    else
2511    {
2512      if (mpz_isNeg(b->z))
2513      {
2514        mpz_sub(a->z,a->z,b->z);
2515      }
2516      else
2517      {
2518        mpz_add(a->z,a->z,b->z);
2519      }
2520      mpz_sub_ui(a->z,a->z,1);
2521    }
2522    MPZ_DIV(a->z,a->z,b->z);
2523    a=nlShort3_noinline(a);
2524  }
2525}
2526
2527number nlFarey(number nN, number nP, const coeffs r)
2528{
2529  mpz_t tmp; mpz_init(tmp);
2530  mpz_t A,B,C,D,E,N,P;
2531  if (SR_HDL(nN) & SR_INT) mpz_init_set_si(N,SR_TO_INT(nN));
2532  else                     mpz_init_set(N,nN->z);
2533  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2534  else                     mpz_init_set(P,nP->z);
2535  assume(!mpz_isNeg(P));
2536  if (mpz_isNeg(N))  mpz_add(N,N,P);
2537  mpz_init_set_si(A,(long)0);
2538  mpz_init_set_ui(B,(unsigned long)1);
2539  mpz_init_set_si(C,(long)0);
2540  mpz_init(D);
2541  mpz_init_set(E,P);
2542  number z=INT_TO_SR(0);
2543  while(mpz_cmp_si(N,(long)0)!=0)
2544  {
2545    mpz_mul(tmp,N,N);
2546    mpz_add(tmp,tmp,tmp);
2547    if (mpz_cmp(tmp,P)<0)
2548    {
2549       // return N/B
2550       z=ALLOC_RNUMBER();
2551       #ifdef LDEBUG
2552       z->debug=123456;
2553       #endif
2554       if (mpz_isNeg(B))
2555       {
2556         mpz_neg(B,B);
2557         mpz_neg(N,N);
2558       }
2559       mpz_init_set(z->z,N);
2560       mpz_init_set(z->n,B);
2561       z->s = 0;
2562       nlNormalize(z,r);
2563       break;
2564    }
2565    //mpz_mod(D,E,N);
2566    //mpz_div(tmp,E,N);
2567    mpz_divmod(tmp,D,E,N);
2568    mpz_mul(tmp,tmp,B);
2569    mpz_sub(C,A,tmp);
2570    mpz_set(E,N);
2571    mpz_set(N,D);
2572    mpz_set(A,B);
2573    mpz_set(B,C);
2574  }
2575  mpz_clear(tmp);
2576  mpz_clear(A);
2577  mpz_clear(B);
2578  mpz_clear(C);
2579  mpz_clear(D);
2580  mpz_clear(E);
2581  mpz_clear(N);
2582  mpz_clear(P);
2583  return z;
2584}
2585
2586void    nlCoeffWrite  (const coeffs, BOOLEAN /*details*/)
2587{
2588  PrintS("//   characteristic : 0\n");
2589}
2590
2591number   nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
2592// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
2593{
2594#ifdef HAVE_FACTORY
2595  setCharacteristic( 0 ); // only in char 0
2596  CFArray X(rl), Q(rl);
2597  int i;
2598  for(i=rl-1;i>=0;i--)
2599  {
2600    X[i]=C->convSingNFactoryN(x[i],FALSE,C); // may be larger MAX_INT
2601    Q[i]=C->convSingNFactoryN(q[i],FALSE,C); // may be larger MAX_INT
2602  }
2603  CanonicalForm xnew,qnew;
2604  chineseRemainder(X,Q,xnew,qnew);
2605  number n=C->convFactoryNSingN(xnew,C);
2606  number p=C->convFactoryNSingN(qnew,C);
2607  number p2=nlIntDiv(p,nlInit(2, C),C);
2608  if (nlGreater(n,p2,C))
2609  {
2610     number n2=nlSub(n,p,C);
2611     nlDelete(&n,C);
2612     n=n2;
2613  }
2614  nlDelete(&p,C);
2615  nlDelete(&p2,C);
2616  return n;
2617#else
2618  WerrorS("not implemented");
2619  return nlInit(0, C);
2620#endif
2621}
2622
2623BOOLEAN nlInitChar(coeffs r, void*)
2624{
2625  assume( getCoeffType(r) == ID );
2626  //const int ch = (int)(long)(p);
2627
2628  r->cfKillChar=NULL;
2629  r->nCoeffIsEqual=ndCoeffIsEqual;
2630  r->cfKillChar = ndKillChar; /* dummy */
2631
2632  r->cfInitMPZ = nlInitMPZ;
2633  r->cfMPZ  = nlMPZ;
2634   
2635  r->cfMult  = nlMult;
2636  r->cfSub   = nlSub;
2637  r->cfAdd   = nlAdd;
2638  r->cfDiv   = nlDiv;
2639  r->cfIntDiv= nlIntDiv;
2640  r->cfIntMod= nlIntMod;
2641  r->cfExactDiv= nlExactDiv;
2642  r->cfInit = nlInit;
2643  r->cfSize  = nlSize;
2644  r->cfInt  = nlInt;
2645   
2646  r->cfChineseRemainder=nlChineseRemainder;
2647  r->cfFarey=nlFarey;
2648  #ifdef HAVE_RINGS
2649  //r->cfDivComp = NULL; // only for ring stuff
2650  //r->cfIsUnit = NULL; // only for ring stuff
2651  //r->cfGetUnit = NULL; // only for ring stuff
2652  //r->cfExtGcd = NULL; // only for ring stuff
2653  //r->cfDivBy = NULL; // only for ring stuff
2654  #endif
2655  r->cfNeg   = nlNeg;
2656  r->cfInvers= nlInvers;
2657  r->cfCopy  = nlCopy;
2658  r->cfRePart = nlCopy;
2659  //r->cfImPart = ndReturn0;
2660  r->cfWrite = nlWrite;
2661  r->cfRead = nlRead;
2662  r->cfNormalize=nlNormalize;
2663  r->cfGreater = nlGreater;
2664  r->cfEqual = nlEqual;
2665  r->cfIsZero = nlIsZero;
2666  r->cfIsOne = nlIsOne;
2667  r->cfIsMOne = nlIsMOne;
2668  r->cfGreaterZero = nlGreaterZero;
2669  r->cfPower = nlPower;
2670  r->cfGetDenom = nlGetDenom;
2671  r->cfGetNumerator = nlGetNumerator;
2672  r->cfGcd  = nlGcd;
2673  r->cfLcm  = nlLcm;
2674  r->cfDelete= nlDelete;
2675  r->cfSetMap = nlSetMap;
2676  //r->cfName = ndName;
2677  r->cfInpMult=nlInpMult;
2678  r->cfInit_bigint=nlCopyMap;
2679  r->cfCoeffWrite=nlCoeffWrite;
2680#ifdef LDEBUG
2681  // debug stuff
2682  r->cfDBTest=nlDBTest;
2683#endif
2684#ifdef HAVE_FACTORY
2685  r->convSingNFactoryN=nlConvSingNFactoryN;
2686  r->convFactoryNSingN=nlConvFactoryNSingN;
2687#endif
2688
2689  // the variables: general stuff (required)
2690  r->nNULL = INT_TO_SR(0);
2691  //r->type = n_Q;
2692  r->ch = 0;
2693  r->has_simple_Alloc=FALSE;
2694  r->has_simple_Inverse=FALSE;
2695
2696  // variables for this type of coeffs:
2697  // (none)
2698  return FALSE;
2699}
2700#if 0
2701number nlMod(number a, number b)
2702{
2703  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
2704  {
2705    int bi=SR_TO_INT(b);
2706    int ai=SR_TO_INT(a);
2707    int bb=ABS(bi);
2708    int c=ai%bb;
2709    if (c<0)  c+=bb;
2710    return (INT_TO_SR(c));
2711  }
2712  number al;
2713  number bl;
2714  if (SR_HDL(a))&SR_INT)
2715    al=nlRInit(SR_TO_INT(a));
2716  else
2717    al=nlCopy(a);
2718  if (SR_HDL(b))&SR_INT)
2719    bl=nlRInit(SR_TO_INT(b));
2720  else
2721    bl=nlCopy(b);
2722  number r=nlRInit(0);
2723  mpz_mod(r->z,al->z,bl->z);
2724  nlDelete(&al);
2725  nlDelete(&bl);
2726  if (mpz_size1(&r->z)<=MP_SMALL)
2727  {
2728    LONG ui=(int)mpz_get_si(&r->z);
2729    if ((((ui<<3)>>3)==ui)
2730    && (mpz_cmp_si(x->z,(long)ui)==0))
2731    {
2732      mpz_clear(&r->z);
2733      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
2734      r=INT_TO_SR(ui);
2735    }
2736  }
2737  return r;
2738}
2739#endif
2740#endif // not P_NUMBERS_H
2741#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.