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

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