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

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