source: git/libpolys/coeffs/longrat.cc @ 8181f5

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