source: git/libpolys/coeffs/longrat.cc @ 0acf3e

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