source: git/libpolys/coeffs/longrat.cc @ 70305c5

spielwiese
Last change on this file since 70305c5 was 5cd658f, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: avoid spurious messages in debug output (unless debugging mpr stuff)
  • 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#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  return u;
1972}
1973
1974// a and b are intermediate, but a*b not
1975number _nlMult_aImm_bImm_rNoImm(number a, number b)
1976{
1977  number u=ALLOC_RNUMBER();
1978#if defined(LDEBUG)
1979  u->debug=123456;
1980#endif
1981  u->s=3;
1982  mpz_init_set_si(u->z,SR_TO_INT(a));
1983  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
1984  return u;
1985}
1986
1987// a or b are not immediate
1988number _nlMult_aNoImm_OR_bNoImm(number a, number b)
1989{
1990  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1991  number u=ALLOC_RNUMBER();
1992#if defined(LDEBUG)
1993  u->debug=123456;
1994#endif
1995  mpz_init(u->z);
1996  if (SR_HDL(b) & SR_INT)
1997  {
1998    number x=a;
1999    a=b;
2000    b=x;
2001  }
2002  if (SR_HDL(a) & SR_INT)
2003  {
2004    u->s=b->s;
2005    if (u->s==1) u->s=0;
2006    if ((long)a>0L)
2007    {
2008      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2009    }
2010    else
2011    {
2012      if (a==INT_TO_SR(-1))
2013      {
2014        mpz_set(u->z,b->z);
2015        mpz_neg(u->z,u->z);
2016        u->s=b->s;
2017      }
2018      else
2019      {
2020        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2021        mpz_neg(u->z,u->z);
2022      }
2023    }
2024    if (u->s<2)
2025    {
2026      if (mpz_cmp(u->z,b->n)==0)
2027      {
2028        mpz_clear(u->z);
2029        FREE_RNUMBER(u);
2030        return INT_TO_SR(1);
2031      }
2032      mpz_init_set(u->n,b->n);
2033    }
2034    else //u->s==3
2035    {
2036      u=nlShort3(u);
2037    }
2038  }
2039  else
2040  {
2041    mpz_mul(u->z,a->z,b->z);
2042    u->s = 0;
2043    if(a->s==3)
2044    {
2045      if(b->s==3)
2046      {
2047        u->s = 3;
2048      }
2049      else
2050      {
2051        if (mpz_cmp(u->z,b->n)==0)
2052        {
2053          mpz_clear(u->z);
2054          FREE_RNUMBER(u);
2055          return INT_TO_SR(1);
2056        }
2057        mpz_init_set(u->n,b->n);
2058      }
2059    }
2060    else
2061    {
2062      if(b->s==3)
2063      {
2064        if (mpz_cmp(u->z,a->n)==0)
2065        {
2066          mpz_clear(u->z);
2067          FREE_RNUMBER(u);
2068          return INT_TO_SR(1);
2069        }
2070        mpz_init_set(u->n,a->n);
2071      }
2072      else
2073      {
2074        mpz_init(u->n);
2075        mpz_mul(u->n,a->n,b->n);
2076        if (mpz_cmp(u->z,u->n)==0)
2077        {
2078          mpz_clear(u->z);
2079          mpz_clear(u->n);
2080          FREE_RNUMBER(u);
2081          return INT_TO_SR(1);
2082        }
2083      }
2084    }
2085  }
2086  return u;
2087}
2088
2089/*2
2090* copy a to b for mapping
2091*/
2092number nlCopyMap(number a, const coeffs src, const coeffs dst)
2093{
2094  assume( getCoeffType(dst) == ID );
2095  assume( getCoeffType(src) == ID );
2096
2097  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2098  {
2099    return a;
2100  }
2101  return _nlCopy_NoImm(a);
2102}
2103
2104nMapFunc nlSetMap(const coeffs src, const coeffs dst)
2105{
2106  assume( getCoeffType(dst) == ID );
2107//  assume( getCoeffType(src) == ID );
2108
2109  if (nCoeff_is_Q(src))
2110  {
2111    return ndCopyMap;
2112  }
2113  if (nCoeff_is_Zp(src))
2114  {
2115    return nlMapP;
2116  }
2117  if (nCoeff_is_R(src))
2118  {
2119    return nlMapR;
2120  }
2121  if (nCoeff_is_long_R(src))
2122  {
2123    return nlMapLongR; /* long R -> Q */
2124  }
2125#ifdef HAVE_RINGS
2126  if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
2127  {
2128    return nlMapGMP;
2129  }
2130  if (nCoeff_is_Ring_2toM(src))
2131  {
2132    return nlMapMachineInt;
2133  }
2134#endif
2135  return NULL;
2136}
2137/*2
2138* z := i
2139*/
2140number nlRInit (long i)
2141{
2142  number z=ALLOC_RNUMBER();
2143#if defined(LDEBUG)
2144  z->debug=123456;
2145#endif
2146  mpz_init_set_si(z->z,i);
2147  z->s = 3;
2148  return z;
2149}
2150
2151/*2
2152* z := i/j
2153*/
2154number nlInit2 (int i, int j, const coeffs r)
2155{
2156  number z=ALLOC_RNUMBER();
2157#if defined(LDEBUG)
2158  z->debug=123456;
2159#endif
2160  mpz_init_set_si(z->z,(long)i);
2161  mpz_init_set_si(z->n,(long)j);
2162  z->s = 0;
2163  nlNormalize(z,r);
2164  return z;
2165}
2166
2167number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
2168{
2169  number z=ALLOC_RNUMBER();
2170#if defined(LDEBUG)
2171  z->debug=123456;
2172#endif
2173  mpz_init_set(z->z,i);
2174  mpz_init_set(z->n,j);
2175  z->s = 0;
2176  nlNormalize(z,r);
2177  return z;
2178}
2179
2180#else // DO_LINLINE
2181
2182// declare immedate routines
2183number nlRInit (long i);
2184BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2185number  _nlCopy_NoImm(number a);
2186number  _nlNeg_NoImm(number a);
2187number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2188void    _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b);
2189number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2190number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2191number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2192
2193#endif
2194
2195/***************************************************************
2196 *
2197 * Routines which might be inlined by p_Numbers.h
2198 *
2199 *******************************************************************/
2200#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2201
2202// routines which are always inlined/static
2203
2204/*2
2205* a = b ?
2206*/
2207LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
2208{
2209  nlTest(a, r);
2210  nlTest(b, r);
2211// short - short
2212  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2213  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2214}
2215
2216LINLINE number nlInit (long i, const coeffs r)
2217{
2218  number n;
2219  LONG ii=(LONG)i;
2220  if ( ((ii << 3) >> 3) == ii ) n=INT_TO_SR(ii);
2221  else                          n=nlRInit(ii);
2222  nlTest(n, r);
2223  return n;
2224}
2225
2226/*2
2227* a == 1 ?
2228*/
2229LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
2230{
2231#ifdef LDEBUG
2232  if (a==NULL) return FALSE;
2233  nlTest(a, r);
2234#endif
2235  return (a==INT_TO_SR(1));
2236}
2237
2238LINLINE BOOLEAN nlIsZero (number a, const coeffs)
2239{
2240  #if 0
2241  if (a==INT_TO_SR(0)) return TRUE;
2242  if ((SR_HDL(a) & SR_INT)||(a==NULL)) return FALSE;
2243  if (mpz_cmp_si(a->z,(long)0)==0)
2244  {
2245    printf("gmp-0 in nlIsZero\n");
2246    dErrorBreak();
2247    return TRUE;
2248  }
2249  return FALSE;
2250  #else
2251  return (a==INT_TO_SR(0));
2252  #endif
2253}
2254
2255/*2
2256* copy a to b
2257*/
2258LINLINE number nlCopy(number a, const coeffs)
2259{
2260  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2261  {
2262    return a;
2263  }
2264  return _nlCopy_NoImm(a);
2265}
2266
2267
2268/*2
2269* delete a
2270*/
2271LINLINE void nlDelete (number * a, const coeffs r)
2272{
2273  if (*a!=NULL)
2274  {
2275    nlTest(*a, r);
2276    if ((SR_HDL(*a) & SR_INT)==0)
2277    {
2278      _nlDelete_NoImm(a);
2279    }
2280    *a=NULL;
2281  }
2282}
2283
2284/*2
2285* za:= - za
2286*/
2287LINLINE number nlNeg (number a, const coeffs R)
2288{
2289  nlTest(a, R);
2290  if(SR_HDL(a) &SR_INT)
2291  {
2292    LONG r=SR_TO_INT(a);
2293    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2294    else               a=INT_TO_SR(-r);
2295    return a;
2296  }
2297  a = _nlNeg_NoImm(a);
2298  nlTest(a, R);
2299  return a;
2300
2301}
2302
2303/*2
2304* u:= a + b
2305*/
2306LINLINE number nlAdd (number a, number b, const coeffs R)
2307{
2308  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2309  {
2310    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2311    if ( ((r << 1) >> 1) == r )
2312      return (number)(long)r;
2313    else
2314      return nlRInit(SR_TO_INT(r));
2315  }
2316  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2317  nlTest(u, R);
2318  return u;
2319}
2320
2321number nlShort1(number a);
2322number nlShort3_noinline(number x);
2323
2324LINLINE void nlInpAdd(number &a, number b, const coeffs r)
2325{
2326  // a=a+b
2327  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2328  {
2329    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2330    if ( ((r << 1) >> 1) == r )
2331      a=(number)(long)r;
2332    else
2333      a=nlRInit(SR_TO_INT(r));
2334  }
2335  else
2336  {
2337    _nlInpAdd_aNoImm_OR_bNoImm(a,b);
2338    nlTest(a,r);
2339  }
2340}
2341
2342LINLINE number nlMult (number a, number b, const coeffs R)
2343{
2344  nlTest(a, R);
2345  nlTest(b, R);
2346  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2347  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2348  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2349  {
2350    LONG r=(LONG)((unsigned LONG)(SR_HDL(a)-1L))*((unsigned LONG)(SR_HDL(b)>>1));
2351    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2352    {
2353      number u=((number) ((r>>1)+SR_INT));
2354      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2355      return nlRInit(SR_HDL(u)>>2);
2356    }
2357    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2358    nlTest(u, R);
2359    return u;
2360
2361  }
2362  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2363  nlTest(u, R);
2364  return u;
2365
2366}
2367
2368
2369/*2
2370* u:= a - b
2371*/
2372LINLINE number nlSub (number a, number b, const coeffs r)
2373{
2374  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2375  {
2376    LONG r=SR_HDL(a)-SR_HDL(b)+1;
2377    if ( ((r << 1) >> 1) == r )
2378    {
2379      return (number)(long)r;
2380    }
2381    else
2382      return nlRInit(SR_TO_INT(r));
2383  }
2384  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2385  nlTest(u, r);
2386  return u;
2387
2388}
2389
2390LINLINE void nlInpMult(number &a, number b, const coeffs r)
2391{
2392  number aa=a;
2393  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2394  {
2395    number n=nlMult(aa,b,r);
2396    nlDelete(&a,r);
2397    a=n;
2398  }
2399  else
2400  {
2401    mpz_mul(aa->z,a->z,b->z);
2402    if (aa->s==3)
2403    {
2404      if(b->s!=3)
2405      {
2406        mpz_init_set(a->n,b->n);
2407        a->s=0;
2408      }
2409    }
2410    else
2411    {
2412      if(b->s!=3)
2413      {
2414        mpz_mul(a->n,a->n,b->n);
2415      }
2416      a->s=0;
2417    }
2418  }
2419}
2420#endif // DO_LINLINE
2421
2422#ifndef P_NUMBERS_H
2423
2424static void nlMPZ(mpz_t m, number &n, const coeffs r)
2425{
2426  assume( getCoeffType(r) == ID );
2427
2428  nlTest(n, r);
2429  nlNormalize(n, r);
2430  if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n));     /* n fits in an int */
2431  else             mpz_init_set(m, (mpz_ptr)n->z);
2432}
2433
2434
2435static number nlInitMPZ(mpz_t m, const coeffs)
2436{
2437  number z = ALLOC_RNUMBER();
2438  mpz_init_set(z->z, m);
2439  mpz_init_set_ui(z->n, 1);
2440  z->s = 3;
2441  return z;
2442}
2443
2444
2445void nlInpGcd(number &a, number b, const coeffs r)
2446{
2447  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2448  {
2449    number n=nlGcd(a,b,r);
2450    nlDelete(&a,r);
2451    a=n;
2452  }
2453  else
2454  {
2455    mpz_gcd(a->z,a->z,b->z);
2456    a=nlShort3_noinline(a);
2457  }
2458}
2459
2460void nlInpIntDiv(number &a, number b, const coeffs r)
2461{
2462  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2463  {
2464    number n=nlIntDiv(a,b, r);
2465    nlDelete(&a,r);
2466    a=n;
2467  }
2468  else
2469  {
2470    if (mpz_isNeg(a->z))
2471    {
2472      if (mpz_isNeg(b->z))
2473      {
2474        mpz_add(a->z,a->z,b->z);
2475      }
2476      else
2477      {
2478        mpz_sub(a->z,a->z,b->z);
2479      }
2480      mpz_add_ui(a->z,a->z,1);
2481    }
2482    else
2483    {
2484      if (mpz_isNeg(b->z))
2485      {
2486        mpz_sub(a->z,a->z,b->z);
2487      }
2488      else
2489      {
2490        mpz_add(a->z,a->z,b->z);
2491      }
2492      mpz_sub_ui(a->z,a->z,1);
2493    }
2494    MPZ_DIV(a->z,a->z,b->z);
2495    a=nlShort3_noinline(a);
2496  }
2497}
2498
2499number nlFarey(number nN, number nP, const coeffs r)
2500{
2501  mpz_t tmp; mpz_init(tmp);
2502  mpz_t A,B,C,D,E,N,P;
2503  if (SR_HDL(nN) & SR_INT) mpz_init_set_si(N,SR_TO_INT(nN));
2504  else                     mpz_init_set(N,nN->z);
2505  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2506  else                     mpz_init_set(P,nP->z);
2507  assume(!mpz_isNeg(P));
2508  if (mpz_isNeg(N))  mpz_add(N,N,P);
2509  mpz_init_set_si(A,(long)0);
2510  mpz_init_set_ui(B,(unsigned long)1);
2511  mpz_init_set_si(C,(long)0);
2512  mpz_init(D);
2513  mpz_init_set(E,P);
2514  number z=INT_TO_SR(0);
2515  while(mpz_cmp_si(N,(long)0)!=0)
2516  {
2517    mpz_mul(tmp,N,N);
2518    mpz_add(tmp,tmp,tmp);
2519    if (mpz_cmp(tmp,P)<0)
2520    {
2521       if (mpz_isNeg(B))
2522       {
2523         mpz_neg(B,B);
2524         mpz_neg(N,N);
2525       }
2526       // check for gcd(N,B)==1
2527       mpz_gcd(tmp,N,B);
2528       if (mpz_cmp_ui(tmp,1)==0)
2529       {
2530         // return N/B
2531         z=ALLOC_RNUMBER();
2532         #ifdef LDEBUG
2533         z->debug=123456;
2534         #endif
2535         mpz_init_set(z->z,N);
2536         mpz_init_set(z->n,B);
2537         z->s = 0;
2538         nlNormalize(z,r);
2539       }
2540       else
2541       {
2542         // return nN (the input) instead of "fail"
2543         z=nlCopy(nN,r);
2544       }
2545       break;
2546    }
2547    //mpz_mod(D,E,N);
2548    //mpz_div(tmp,E,N);
2549    mpz_divmod(tmp,D,E,N);
2550    mpz_mul(tmp,tmp,B);
2551    mpz_sub(C,A,tmp);
2552    mpz_set(E,N);
2553    mpz_set(N,D);
2554    mpz_set(A,B);
2555    mpz_set(B,C);
2556  }
2557  mpz_clear(tmp);
2558  mpz_clear(A);
2559  mpz_clear(B);
2560  mpz_clear(C);
2561  mpz_clear(D);
2562  mpz_clear(E);
2563  mpz_clear(N);
2564  mpz_clear(P);
2565  return z;
2566}
2567
2568number nlExtGcd(number a, number b, number *s, number *t, const coeffs)
2569{
2570  mpz_t aa,bb;
2571  *s=ALLOC_RNUMBER();
2572  mpz_init((*s)->z); (*s)->s=3;
2573  (*t)=ALLOC_RNUMBER();
2574  mpz_init((*t)->z); (*t)->s=3;
2575  number g=ALLOC_RNUMBER();
2576  mpz_init(g->z); g->s=3;
2577  if (SR_HDL(a) & SR_INT)
2578  {
2579    mpz_init_set_si(aa,SR_TO_INT(a));
2580  }
2581  else
2582  {
2583    mpz_init_set(aa,a->z);
2584  }
2585  if (SR_HDL(b) & SR_INT)
2586  {
2587    mpz_init_set_si(bb,SR_TO_INT(b));
2588  }
2589  else
2590  {
2591    mpz_init_set(bb,b->z);
2592  }
2593  mpz_gcdext(g->z,(*s)->z,(*t)->z,aa,bb);
2594  mpz_clear(aa);
2595  mpz_clear(bb);
2596  (*s)=nlShort3((*s));
2597  (*t)=nlShort3((*t));
2598  g=nlShort3(g);
2599  return g;
2600}
2601
2602void    nlCoeffWrite  (const coeffs, BOOLEAN /*details*/)
2603{
2604  PrintS("//   characteristic : 0\n");
2605}
2606
2607number   nlChineseRemainderSym(number *x, number *q,int rl, BOOLEAN sym, const coeffs CF)
2608// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
2609{
2610#ifdef HAVE_FACTORY
2611  setCharacteristic( 0 ); // only in char 0
2612  Off(SW_RATIONAL);
2613  CFArray X(rl), Q(rl);
2614  int i;
2615  for(i=rl-1;i>=0;i--)
2616  {
2617    X[i]=CF->convSingNFactoryN(x[i],FALSE,CF); // may be larger MAX_INT
2618    Q[i]=CF->convSingNFactoryN(q[i],FALSE,CF); // may be larger MAX_INT
2619  }
2620  CanonicalForm xnew,qnew;
2621  chineseRemainder(X,Q,xnew,qnew);
2622  number n=CF->convFactoryNSingN(xnew,CF);
2623  if (sym)
2624  {
2625    number p=CF->convFactoryNSingN(qnew,CF);
2626    number p2=nlIntDiv(p,nlInit(2, CF),CF);
2627    if (nlGreater(n,p2,CF))
2628    {
2629       number n2=nlSub(n,p,CF);
2630       nlDelete(&n,CF);
2631       n=n2;
2632    }
2633    nlDelete(&p2,CF);
2634    nlDelete(&p,CF);
2635  }
2636  return n;
2637#else
2638  WerrorS("not implemented");
2639  return nlInit(0,CF);
2640#endif
2641}
2642number   nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
2643{
2644  return nlChineseRemainderSym(x,q,rl,TRUE,C);
2645}
2646
2647static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2648{
2649  assume(cf != NULL);
2650  assume(getCoeffType(cf) == ID);
2651
2652  numberCollectionEnumerator.Reset();
2653
2654  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2655  {
2656    c = n_Init(1, cf);
2657    return;
2658  }
2659 
2660  // all coeffs are given by integers!!!
2661
2662  // part 1, find a small candidate for gcd
2663  number cand1,cand;
2664  int s1,s;
2665  s=2147483647; // max. int
2666
2667  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
2668
2669  int normalcount = 0;
2670  do
2671  {
2672    number& n = numberCollectionEnumerator.Current();
2673    nlNormalize(n, cf); ++normalcount;
2674    cand1 = n;
2675   
2676    if (SR_HDL(cand1)&SR_INT) { cand=cand1; break; }
2677    assume(cand1->s==3); // all coeffs should be integers // ==0?!! after printing
2678    s1=mpz_size1(cand1->z);
2679    if (s>s1)
2680    {
2681      cand=cand1;
2682      s=s1;
2683    }
2684  } while (numberCollectionEnumerator.MoveNext() );
2685
2686//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
2687
2688  cand=nlCopy(cand,cf);
2689  // part 2: compute gcd(cand,all coeffs)
2690
2691  numberCollectionEnumerator.Reset();
2692 
2693  while (numberCollectionEnumerator.MoveNext() )
2694  {
2695    number& n = numberCollectionEnumerator.Current();
2696
2697    if( (--normalcount) <= 0)
2698      nlNormalize(n, cf);
2699
2700    nlInpGcd(cand, n, cf);
2701
2702    assume( nlGreaterZero(cand,cf) );
2703   
2704    if(nlIsOne(cand,cf))
2705    {
2706      c = cand;
2707
2708      if(!lc_is_pos)
2709      {
2710        // make the leading coeff positive
2711        c = nlNeg(c, cf);
2712        numberCollectionEnumerator.Reset();
2713       
2714        while (numberCollectionEnumerator.MoveNext() )
2715        {
2716          number& nn = numberCollectionEnumerator.Current();
2717          nn = nlNeg(nn, cf);
2718        }
2719      }
2720      return;
2721    }
2722  } 
2723
2724  // part3: all coeffs = all coeffs / cand
2725  if (!lc_is_pos)
2726    cand = nlNeg(cand,cf);
2727 
2728  c = cand;
2729  numberCollectionEnumerator.Reset();
2730
2731  while (numberCollectionEnumerator.MoveNext() )
2732  {
2733    number& n = numberCollectionEnumerator.Current();
2734    number t=nlIntDiv(n, cand, cf); // simple integer exact division, no ratios to remain
2735    nlDelete(&n, cf);
2736    n = t;
2737  }
2738}
2739
2740static void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2741{
2742  assume(cf != NULL);
2743  assume(getCoeffType(cf) == ID);
2744
2745  numberCollectionEnumerator.Reset();
2746 
2747  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2748  {
2749    c = n_Init(1, cf);
2750//    assume( n_GreaterZero(c, cf) );
2751    return;
2752  }
2753
2754  // all coeffs are given by integers after returning from this routine
2755
2756  // part 1, collect product of all denominators /gcds
2757  number cand;
2758  cand=ALLOC_RNUMBER();
2759#if defined(LDEBUG)
2760  cand->debug=123456;
2761#endif
2762  cand->s=3;
2763
2764  int s=0;
2765 
2766  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
2767
2768  do
2769  {
2770    number& cand1 = numberCollectionEnumerator.Current();
2771
2772    if (!(SR_HDL(cand1)&SR_INT))
2773    {
2774      nlNormalize(cand1, cf);
2775      if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
2776      && (cand1->s==1))             // and is a normalised rational
2777      {
2778        if (s==0) // first denom, we meet
2779        {
2780          mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
2781          s=1;
2782        }
2783        else // we have already something
2784        {
2785          mpz_lcm(cand->z, cand->z, cand1->n);
2786        }
2787      }
2788    }
2789  }
2790  while (numberCollectionEnumerator.MoveNext() );
2791 
2792
2793  if (s==0) // nothing to do, all coeffs are already integers
2794  {
2795//    mpz_clear(tmp);
2796    FREE_RNUMBER(cand);
2797    if (lc_is_pos)
2798      c=nlInit(1,cf);
2799    else
2800    {
2801      // make the leading coeff positive
2802      c=nlInit(-1,cf);
2803
2804      // TODO: incorporate the following into the loop below?
2805      numberCollectionEnumerator.Reset();
2806      while (numberCollectionEnumerator.MoveNext() )
2807      {
2808        number& n = numberCollectionEnumerator.Current();
2809        n = nlNeg(n, cf);
2810      } 
2811    }
2812//    assume( n_GreaterZero(c, cf) );
2813    return;
2814  }
2815
2816  cand = nlShort3(cand);
2817
2818  // part2: all coeffs = all coeffs * cand
2819  // make the lead coeff positive
2820  numberCollectionEnumerator.Reset();
2821 
2822  if (!lc_is_pos)
2823    cand = nlNeg(cand, cf);
2824 
2825  c = cand;
2826 
2827  while (numberCollectionEnumerator.MoveNext() )
2828  {
2829    number &n = numberCollectionEnumerator.Current();
2830    n_InpMult(n, cand, cf);
2831  }
2832
2833}
2834
2835BOOLEAN nlInitChar(coeffs r, void*)
2836{
2837  assume( getCoeffType(r) == ID );
2838  //const int ch = (int)(long)(p);
2839
2840  r->cfKillChar=NULL;
2841  r->nCoeffIsEqual=ndCoeffIsEqual;
2842  r->cfKillChar = ndKillChar; /* dummy */
2843
2844  r->cfInitMPZ = nlInitMPZ;
2845  r->cfMPZ  = nlMPZ;
2846
2847  r->cfMult  = nlMult;
2848  r->cfSub   = nlSub;
2849  r->cfAdd   = nlAdd;
2850  r->cfDiv   = nlDiv;
2851  r->cfIntDiv= nlIntDiv;
2852  r->cfIntMod= nlIntMod;
2853  r->cfExactDiv= nlExactDiv;
2854  r->cfInit = nlInit;
2855  r->cfSize  = nlSize;
2856  r->cfInt  = nlInt;
2857
2858  r->cfChineseRemainder=nlChineseRemainderSym;
2859  r->cfFarey=nlFarey;
2860  #ifdef HAVE_RINGS
2861  //r->cfDivComp = NULL; // only for ring stuff
2862  //r->cfIsUnit = NULL; // only for ring stuff
2863  //r->cfGetUnit = NULL; // only for ring stuff
2864  //r->cfDivBy = NULL; // only for ring stuff
2865  #endif
2866  r->cfNeg   = nlNeg;
2867  r->cfInvers= nlInvers;
2868  r->cfCopy  = nlCopy;
2869  r->cfRePart = nlCopy;
2870  //r->cfImPart = ndReturn0;
2871  r->cfWriteLong = nlWrite;
2872  r->cfRead = nlRead;
2873  r->cfNormalize=nlNormalize;
2874  r->cfGreater = nlGreater;
2875  r->cfEqual = nlEqual;
2876  r->cfIsZero = nlIsZero;
2877  r->cfIsOne = nlIsOne;
2878  r->cfIsMOne = nlIsMOne;
2879  r->cfGreaterZero = nlGreaterZero;
2880  r->cfPower = nlPower;
2881  r->cfGetDenom = nlGetDenom;
2882  r->cfGetNumerator = nlGetNumerator;
2883  r->cfGcd  = nlGcd;
2884  r->cfExtGcd = nlExtGcd; // only for ring stuff and Z
2885  r->cfLcm  = nlLcm;
2886  r->cfDelete= nlDelete;
2887  r->cfSetMap = nlSetMap;
2888  //r->cfName = ndName;
2889  r->cfInpMult=nlInpMult;
2890  r->cfInit_bigint=nlCopyMap;
2891  r->cfCoeffWrite=nlCoeffWrite;
2892
2893  r->cfClearContent = nlClearContent;
2894  r->cfClearDenominators = nlClearDenominators;
2895
2896#ifdef LDEBUG
2897  // debug stuff
2898  r->cfDBTest=nlDBTest;
2899#endif
2900#ifdef HAVE_FACTORY
2901  r->convSingNFactoryN=nlConvSingNFactoryN;
2902  r->convFactoryNSingN=nlConvFactoryNSingN;
2903#endif
2904
2905  // the variables: general stuff (required)
2906  r->nNULL = INT_TO_SR(0);
2907  //r->type = n_Q;
2908  r->ch = 0;
2909  r->has_simple_Alloc=FALSE;
2910  r->has_simple_Inverse=FALSE;
2911
2912  // variables for this type of coeffs:
2913  // (none)
2914  return FALSE;
2915}
2916#if 0
2917number nlMod(number a, number b)
2918{
2919  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
2920  {
2921    int bi=SR_TO_INT(b);
2922    int ai=SR_TO_INT(a);
2923    int bb=ABS(bi);
2924    int c=ai%bb;
2925    if (c<0)  c+=bb;
2926    return (INT_TO_SR(c));
2927  }
2928  number al;
2929  number bl;
2930  if (SR_HDL(a))&SR_INT)
2931    al=nlRInit(SR_TO_INT(a));
2932  else
2933    al=nlCopy(a);
2934  if (SR_HDL(b))&SR_INT)
2935    bl=nlRInit(SR_TO_INT(b));
2936  else
2937    bl=nlCopy(b);
2938  number r=nlRInit(0);
2939  mpz_mod(r->z,al->z,bl->z);
2940  nlDelete(&al);
2941  nlDelete(&bl);
2942  if (mpz_size1(&r->z)<=MP_SMALL)
2943  {
2944    LONG ui=(int)mpz_get_si(&r->z);
2945    if ((((ui<<3)>>3)==ui)
2946    && (mpz_cmp_si(x->z,(long)ui)==0))
2947    {
2948      mpz_clear(&r->z);
2949      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
2950      r=INT_TO_SR(ui);
2951    }
2952  }
2953  return r;
2954}
2955#endif
2956#endif // not P_NUMBERS_H
2957#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.