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

jengelh-datetimespielwiese
Last change on this file since a44bcf was a44bcf, checked in by Hans Schoenemann <hannes@…>, 10 years ago
add: checks for rational coeff. in rational functions
  • Property mode set to 100644
File size: 59.0 KB
RevLine 
[35aab3]1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
6*/
[d663f2]7
[22a09d]8#include "config.h"
9#include <misc/auxiliary.h>
10
11#ifdef HAVE_FACTORY
12#include <factory/factory.h>
13#endif
[35aab3]14
[ea11dc]15#include <coeffs/longrat.h>
[246d15]16
[9eb0f9]17
[246d15]18// 64 bit version:
[7aff627]19//#if SIZEOF_LONG == 8
20#if 0
[246d15]21#define MAX_NUM_SIZE 60
[3feb05]22#define POW_2_28 (1L<<60)
[246d15]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);
[cc3240c]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)
[246d15]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);
[896561]46      FREE_RNUMBER(x);
[246d15]47      return INT_TO_SR(ui);
48    }
49  }
50  return x;
51}
52
[35aab3]53#ifndef LONGRAT_CC
54#define LONGRAT_CC
55
56#include <string.h>
57#include <float.h>
[9eb0f9]58
[2d805a]59#include <coeffs/coeffs.h>
[31213a4]60#include <reporter/reporter.h>
61#include <omalloc/omalloc.h>
[9eb0f9]62
[2d805a]63#include <coeffs/numbers.h>
64#include <coeffs/modulop.h>
65#include <coeffs/shortfl.h>
66#include <coeffs/mpr_complex.h>
[35aab3]67
68#ifndef BYTES_PER_MP_LIMB
69#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
70#endif
71
[8181f5]72//#define SR_HDL(A) ((long)(A))
[46f6af6]73/*#define SR_INT    1L*/
[35aab3]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
[73a9ffb]85/// Our Type!
86static const n_coeffType ID = n_Q;
[35aab3]87
[73a9ffb]88void    _nlDelete_NoImm(number *a);
[fa31d2]89
[35aab3]90/***************************************************************
91 *
92 * Routines which are never inlined by p_Numbers.h
93 *
94 *******************************************************************/
95#ifndef P_NUMBERS_H
96
[246d15]97number nlShort3_noinline(number x) // assume x->s==3
98{
99  return nlShort3(x);
100}
101
[35aab3]102
[eb059c5]103number nlOne=INT_TO_SR(1);
[35aab3]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
[fa31d2]118static number nlMapP(number from, const coeffs src, const coeffs dst)
[35aab3]119{
[fa31d2]120  assume( getCoeffType(dst) == ID );
121  assume( getCoeffType(src) == n_Zp );
[9bb5457]122
[35aab3]123  number to;
[fa31d2]124  to = nlInit(npInt(from,src), dst);
[35aab3]125  return to;
126}
127
[fa31d2]128static number nlMapLongR(number from, const coeffs src, const coeffs dst);
129static number nlMapR(number from, const coeffs src, const coeffs dst);
[35aab3]130
[d663f2]131
[894f5b1]132#ifdef HAVE_RINGS
133/*2
134* convert from a GMP integer
135*/
[9bb5457]136number nlMapGMP(number from, const coeffs /*src*/, const coeffs /*dst*/)
[894f5b1]137{
[896561]138  number z=ALLOC_RNUMBER();
[894f5b1]139#if defined(LDEBUG)
140  z->debug=123456;
141#endif
[a604c3]142  mpz_init_set(z->z,(mpz_ptr) from);
[20704f]143  //mpz_init_set_ui(&z->n,1);
[894f5b1]144  z->s = 3;
[246d15]145  z=nlShort3(z);
[894f5b1]146  return z;
147}
148
149/*2
150* convert from an machine long
151*/
[9bb5457]152number nlMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
[894f5b1]153{
[896561]154  number z=ALLOC_RNUMBER();
[894f5b1]155#if defined(LDEBUG)
156  z->debug=123456;
157#endif
[a604c3]158  mpz_init_set_ui(z->z,(unsigned long) from);
[894f5b1]159  z->s = 3;
[246d15]160  z=nlShort3(z);
[894f5b1]161  return z;
162}
163#endif
164
[35aab3]165
166#ifdef LDEBUG
[9bb5457]167BOOLEAN nlDBTest(number a, const char *f,const int l, const coeffs /*r*/)
[35aab3]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);
[7447d8]175  if ((((long)a)&3L)==3L)
[35aab3]176  {
177    Print(" !!longrat:ptr(3) in %s:%d\n",f,l);
178    return FALSE;
179  }
[7447d8]180  if ((((long)a)&3L)==1L)
[35aab3]181  {
[246d15]182    if (((((LONG)(long)a)<<1)>>1)!=((LONG)(long)a))
[35aab3]183    {
[2132395]184      Print(" !!longrat:arith:%lx in %s:%d\n",(long)a, f,l);
[35aab3]185      return FALSE;
186    }
187    return TRUE;
188  }
[4a2260e]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. */
[d12f186]193  omCheckIf(omCheckAddrSize(a,sizeof(*a)), return FALSE);
[35aab3]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  }
[c28ecf]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. */
[896561]209  //omCheckAddrSize(a->z[0]._mp_d,a->z[0]._mp_alloc*BYTES_PER_MP_LIMB);
[a604c3]210  if (a->z[0]._mp_alloc==0)
[2132395]211    Print("!!longrat:z->alloc=0 in %s:%d\n",f,l);
[35aab3]212
213  if (a->s<2)
214  {
[a122f4]215    if ((a->n[0]._mp_d[0]==0)&&(a->n[0]._mp_alloc<=1))
[9d68fd]216    {
217      Print("!!longrat: n==0 in %s:%d\n",f,l);
218      return FALSE;
219    }
[c28ecf]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. */
[896561]224    //omCheckIf(omCheckAddrSize(a->n[0]._mp_d,a->n[0]._mp_alloc*BYTES_PER_MP_LIMB), return FALSE);
[a604c3]225    if (a->z[0]._mp_alloc==0)
[2132395]226      Print("!!longrat:n->alloc=0 in %s:%d\n",f,l);
[a604c3]227    if ((mpz_size1(a->n) ==1) && (mpz_cmp_si(a->n,(long)1)==0))
[35aab3]228    {
229      Print("!!longrat:integer as rational in %s:%d\n",f,l);
[8fd3bd5]230      mpz_clear(a->n); a->s=3;
[35aab3]231      return FALSE;
232    }
[a604c3]233    else if (mpz_isNeg(a->n))
[35aab3]234    {
235      Print("!!longrat:div. by negative in %s:%d\n",f,l);
[a604c3]236      mpz_neg(a->z,a->z);
237      mpz_neg(a->n,a->n);
[35aab3]238      return FALSE;
239    }
240    return TRUE;
241  }
[9257fe]242  //if (a->s==2)
243  //{
244  //  Print("!!longrat:s=2 in %s:%d\n",f,l);
245  //  return FALSE;
246  //}
[a604c3]247  if (mpz_size1(a->z)>MP_SMALL) return TRUE;
[246d15]248  LONG ui=(int)mpz_get_si(a->z);
[35aab3]249  if ((((ui<<3)>>3)==ui)
[a604c3]250  && (mpz_cmp_si(a->z,(long)ui)==0))
[35aab3]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
[a934fc8]260#ifdef HAVE_FACTORY
[9bb5457]261CanonicalForm nlConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs /*r*/ )
[a934fc8]262{
[abb4787]263  if (setChar) setCharacteristic( 0 );
264
[a934fc8]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    {
[8b46459]274      mpz_t dummy;
[8710ff0]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      }
[a934fc8]282    }
283    else
284    {
285      // assume s==0 or s==1
[8b46459]286      mpz_t num, den;
[a934fc8]287      On(SW_RATIONAL);
[8b46459]288      mpz_init_set( num, n->z );
289      mpz_init_set( den, n->n );
[a934fc8]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  {
[8710ff0]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);
[a934fc8]307    return nlInit(n.intval(),r);
308  }
309  else
310  {
[d663f2]311    number z = ALLOC_RNUMBER(); //Q!? // (number)omAllocBin(rnumber_bin);
[a934fc8]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
[246d15]328number nlRInit (long i);
[35aab3]329
[371bc5f]330static number nlMapR(number from, const coeffs src, const coeffs dst)
[35aab3]331{
[fa31d2]332  assume( getCoeffType(dst) == ID );
333  assume( getCoeffType(src) == n_R );
[9bb5457]334
[35aab3]335  double f=nrFloat(from);
[eb059c5]336  if (f==0.0) return INT_TO_SR(0);
[35aab3]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  }
[a642c64]352  number re=nlRInit(1);
[bbf7a8]353  mpz_set_d(re->z,f);
[a642c64]354  memcpy(&(re->n),&h1,sizeof(h1));
355  re->s=0; /* not normalized */
[371bc5f]356  if(f_sign==-1) re=nlNeg(re,dst);
357  nlNormalize(re,dst);
[a642c64]358  return re;
[35aab3]359}
360
[9bb5457]361static number nlMapLongR(number from, const coeffs src, const coeffs dst)
[35aab3]362{
[fa31d2]363  assume( getCoeffType(dst) == ID );
364  assume( getCoeffType(src) == n_long_R );
[9bb5457]365
[35aab3]366  gmp_float *ff=(gmp_float*)from;
367  mpf_t *f=ff->_mpfp();
368  number res;
[a604c3]369  mpz_ptr dest,ndest;
[35aab3]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)
[eb059c5]376    return INT_TO_SR(0);
[35aab3]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;
[896561]393  res = ALLOC_RNUMBER();
[35aab3]394#if defined(LDEBUG)
395  res->debug=123456;
396#endif
[a604c3]397  dest = res->z;
[35aab3]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;
[a604c3]409    ndest = res->n;
[35aab3]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)
[371bc5f]429    nlNormalize(res,dst);
[bbf7a8]430  else if (mpz_size1(res->z)<=MP_SMALL)
[35aab3]431  {
432    // res is new, res->ref is 1
[246d15]433    res=nlShort3(res);
[35aab3]434  }
[2b957a]435  nlTest(res, dst);
[35aab3]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);
[611871]464//  memcpy(&(r->n),&h1,sizeof(h1));
[35aab3]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
[611871]479//    memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
[35aab3]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
[611871]489//      memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
[35aab3]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
[611871]495//      memcpy(mpz_limb_d(&r->z)-f_shift,((*f)[0]._mp_d),
[35aab3]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
[121fd9]504int nlSize(number a, const coeffs)
[35aab3]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 */
[a604c3]510  int s=a->z[0]._mp_alloc;
[9bb5457]511//  while ((s>0) &&(a->z._mp_d[s]==0L)) s--;
[e08bce]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
[790096]516//  s++;
[35aab3]517  if (a->s<2)
518  {
[a604c3]519    int d=a->n[0]._mp_alloc;
[9bb5457]520//    while ((d>0) && (a->n._mp_d[d]==0L)) d--;
[e08bce]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
[1da0b5]525    s+=d;
[35aab3]526  }
527  return s;
528}
529
530/*2
531* convert number to int
532*/
[a642c64]533int nlInt(number &i, const coeffs r)
[35aab3]534{
[2b957a]535  nlTest(i, r);
[a642c64]536  nlNormalize(i,r);
[35aab3]537  if (SR_HDL(i) &SR_INT) return SR_TO_INT(i);
538  if (i->s==3)
539  {
[a604c3]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;
[35aab3]543    return ul;
544  }
[a604c3]545  mpz_t tmp;
[35aab3]546  int ul;
[a604c3]547  mpz_init(tmp);
548  MPZ_DIV(tmp,i->z,i->n);
549  if(mpz_size1(tmp)>MP_SMALL) ul=0;
[35aab3]550  else
551  {
[a604c3]552    ul=(int)mpz_get_si(tmp);
553    if (mpz_cmp_si(tmp,(long)ul)!=0) ul=0;
[35aab3]554  }
[a604c3]555  mpz_clear(tmp);
[35aab3]556  return ul;
557}
558
[7f7573]559/*2
560* convert number to bigint
561*/
[a642c64]562number nlBigInt(number &i, const coeffs r)
[7f7573]563{
[2b957a]564  nlTest(i, r);
[a642c64]565  nlNormalize(i,r);
[7f7573]566  if (SR_HDL(i) &SR_INT) return (i);
567  if (i->s==3)
568  {
[a642c64]569    return nlCopy(i,r);
[7f7573]570  }
571  number tmp=nlRInit(1);
[a604c3]572  MPZ_DIV(tmp->z,i->z,i->n);
[246d15]573  tmp=nlShort3(tmp);
[7f7573]574  return tmp;
575}
576
[35aab3]577/*
578* 1/a
579*/
[a642c64]580number nlInvers(number a, const coeffs r)
[35aab3]581{
[2b957a]582  nlTest(a, r);
[35aab3]583  number n;
584  if (SR_HDL(a) & SR_INT)
585  {
[246d15]586    if ((a==INT_TO_SR(1L)) || (a==INT_TO_SR(-1L)))
[35aab3]587    {
588      return a;
589    }
[a642c64]590    if (nlIsZero(a,r))
[35aab3]591    {
[b7e838]592      WerrorS(nDivBy0);
[35aab3]593      return INT_TO_SR(0);
594    }
[896561]595    n=ALLOC_RNUMBER();
[35aab3]596#if defined(LDEBUG)
597    n->debug=123456;
598#endif
599    n->s=1;
[7447d8]600    if ((long)a>0L)
[35aab3]601    {
[a604c3]602      mpz_init_set_si(n->z,(long)1);
603      mpz_init_set_si(n->n,(long)SR_TO_INT(a));
[35aab3]604    }
605    else
606    {
[a604c3]607      mpz_init_set_si(n->z,(long)-1);
608      mpz_init_set_si(n->n,(long)-SR_TO_INT(a));
[35aab3]609    }
[2b957a]610    nlTest(n, r);
[35aab3]611    return n;
612  }
[896561]613  n=ALLOC_RNUMBER();
[35aab3]614#if defined(LDEBUG)
615  n->debug=123456;
616#endif
617  {
618    n->s=a->s;
[a604c3]619    mpz_init_set(n->n,a->z);
[35aab3]620    switch (a->s)
621    {
622      case 0:
623      case 1:
[a604c3]624              mpz_init_set(n->z,a->n);
625              if (mpz_isNeg(n->n)) /* && n->s<2*/
[35aab3]626              {
[a604c3]627                mpz_neg(n->z,n->z);
628                mpz_neg(n->n,n->n);
[35aab3]629              }
[a604c3]630              if (mpz_cmp_si(n->n,(long)1)==0)
[35aab3]631              {
[a604c3]632                mpz_clear(n->n);
[35aab3]633                n->s=3;
[246d15]634                n=nlShort3(n);
[35aab3]635              }
636              break;
637      case 3:
638              n->s=1;
[a604c3]639              if (mpz_isNeg(n->n)) /* && n->s<2*/
[35aab3]640              {
[a604c3]641                mpz_neg(n->n,n->n);
642                mpz_init_set_si(n->z,(long)-1);
[35aab3]643              }
644              else
645              {
[a604c3]646                mpz_init_set_si(n->z,(long)1);
[35aab3]647              }
648              break;
649    }
650  }
[2b957a]651  nlTest(n, r);
[35aab3]652  return n;
653}
654
655
656/*2
657* u := a / b in Z, if b | a (else undefined)
658*/
[2b957a]659number   nlExactDiv(number a, number b, const coeffs r)
[35aab3]660{
661  if (b==INT_TO_SR(0))
662  {
[b7e838]663    WerrorS(nDivBy0);
[35aab3]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)   */
[246d15]672    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
[35aab3]673    {
[246d15]674      return nlRInit(POW_2_28);
[35aab3]675    }
[246d15]676    long aa=SR_TO_INT(a);
677    long bb=SR_TO_INT(b);
[35aab3]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  }
[896561]686  u=ALLOC_RNUMBER();
[35aab3]687#if defined(LDEBUG)
688  u->debug=123456;
689#endif
[a604c3]690  mpz_init(u->z);
[35aab3]691  /* u=a/b */
692  u->s = 3;
[a604c3]693  MPZ_EXACTDIV(u->z,a->z,b->z);
[35aab3]694  if (bb!=NULL)
695  {
[a604c3]696    mpz_clear(bb->z);
[35aab3]697#if defined(LDEBUG)
698    bb->debug=654324;
699#endif
[d663f2]700    FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
[35aab3]701  }
[246d15]702  u=nlShort3(u);
[2b957a]703  nlTest(u, r);
[35aab3]704  return u;
705}
706
707/*2
708* u := a / b in Z
709*/
[2b957a]710number nlIntDiv (number a, number b, const coeffs r)
[35aab3]711{
712  if (b==INT_TO_SR(0))
713  {
[b7e838]714    WerrorS(nDivBy0);
[35aab3]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)   */
[246d15]723    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
[35aab3]724    {
[246d15]725      return nlRInit(POW_2_28);
[35aab3]726    }
[246d15]727    long aa=SR_TO_INT(a);
728    long bb=SR_TO_INT(b);
[7326063]729    return INT_TO_SR(aa/bb);
[35aab3]730  }
731  if (SR_HDL(a) & SR_INT)
732  {
733    /* the small int -(1<<28) divided by 2^28 is 1   */
[246d15]734    if (a==INT_TO_SR(-(POW_2_28)))
[35aab3]735    {
[246d15]736      if(mpz_cmp_si(b->z,(POW_2_28))==0)
[35aab3]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  }
[896561]750  u=ALLOC_RNUMBER();
[35aab3]751#if defined(LDEBUG)
752  u->debug=123456;
753#endif
[98bda9]754  assume(a->s==3);
755  assume(b->s==3);
[a604c3]756  mpz_init_set(u->z,a->z);
[35aab3]757  /* u=u/b */
758  u->s = 3;
[a604c3]759  MPZ_DIV(u->z,u->z,b->z);
[35aab3]760  if (bb!=NULL)
761  {
[a604c3]762    mpz_clear(bb->z);
[35aab3]763#if defined(LDEBUG)
764    bb->debug=654324;
765#endif
[de27d8]766    FREE_RNUMBER(bb);
[35aab3]767  }
[246d15]768  u=nlShort3(u);
[de27d8]769  nlTest(u,r);
[35aab3]770  return u;
771}
772
773/*2
774* u := a mod b in Z, u>=0
775*/
[a642c64]776number nlIntMod (number a, number b, const coeffs r)
[35aab3]777{
778  if (b==INT_TO_SR(0))
779  {
[b7e838]780    WerrorS(nDivBy0);
[35aab3]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  {
[de27d8]788    LONG bb=SR_TO_INT(b);
789    LONG c=SR_TO_INT(a)%bb;
790    return INT_TO_SR(c);
[35aab3]791  }
792  if (SR_HDL(a) & SR_INT)
793  {
[824c47]794    /* a is a small and b is a large int: -> a */
[de27d8]795    return a;
[35aab3]796  }
797  number bb=NULL;
798  if (SR_HDL(b) & SR_INT)
799  {
800    bb=nlRInit(SR_TO_INT(b));
801    b=bb;
802  }
[896561]803  u=ALLOC_RNUMBER();
[35aab3]804#if defined(LDEBUG)
805  u->debug=123456;
806#endif
[a604c3]807  mpz_init(u->z);
[35aab3]808  u->s = 3;
[a604c3]809  mpz_mod(u->z,a->z,b->z);
[35aab3]810  if (bb!=NULL)
811  {
[a604c3]812    mpz_clear(bb->z);
[35aab3]813#if defined(LDEBUG)
814    bb->debug=654324;
815#endif
[de27d8]816    FREE_RNUMBER(bb);
[35aab3]817  }
[a604c3]818  if (mpz_isNeg(u->z))
[35aab3]819  {
[a604c3]820    if (mpz_isNeg(b->z))
821      mpz_sub(u->z,u->z,b->z);
[35aab3]822    else
[a604c3]823      mpz_add(u->z,u->z,b->z);
[35aab3]824  }
[246d15]825  u=nlShort3(u);
[de27d8]826  nlTest(u,r);
[35aab3]827  return u;
828}
829
830/*2
831* u := a / b
832*/
[a642c64]833number nlDiv (number a, number b, const coeffs r)
[35aab3]834{
835  number u;
[a642c64]836  if (nlIsZero(b,r))
[35aab3]837  {
[b7e838]838    WerrorS(nDivBy0);
[35aab3]839    return INT_TO_SR(0);
840  }
[896561]841  u=ALLOC_RNUMBER();
[35aab3]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  {
[135e78]849    LONG i=SR_TO_INT(a);
850    LONG j=SR_TO_INT(b);
[1a97da]851    if ((i==-POW_2_28) && (j== -1L))
[e11ed0]852    {
[896561]853      FREE_RNUMBER(u);
[1a97da]854      return nlRInit(POW_2_28);
855    }
[135e78]856    LONG r=i%j;
[35aab3]857    if (r==0)
858    {
[d663f2]859      FREE_RNUMBER(u); // omFreeBin((void *)u, rnumber_bin);
[35aab3]860      return INT_TO_SR(i/j);
861    }
[a604c3]862    mpz_init_set_si(u->z,(long)i);
863    mpz_init_set_si(u->n,(long)j);
[35aab3]864  }
865  else
866  {
[a604c3]867    mpz_init(u->z);
[35aab3]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      {
[a604c3]874        mpz_mul_si(u->z,b->n,SR_TO_INT(a));
[35aab3]875      }
876      else
877      // short a / long z -> a/z
878      {
[a604c3]879        mpz_set_si(u->z,SR_TO_INT(a));
[35aab3]880      }
[a604c3]881      if (mpz_cmp(u->z,b->z)==0)
[35aab3]882      {
[a604c3]883        mpz_clear(u->z);
[896561]884        FREE_RNUMBER(u);
[35aab3]885        return INT_TO_SR(1);
886      }
[a604c3]887      mpz_init_set(u->n,b->z);
[35aab3]888    }
889// ---------- long / short ------------------------------------
890    else if (SR_HDL(b) & SR_INT)
891    {
[a604c3]892      mpz_set(u->z,a->z);
[35aab3]893      // (z/n) / b -> z/(n*b)
894      if (a->s<2)
895      {
[a604c3]896        mpz_init_set(u->n,a->n);
[7447d8]897        if ((long)b>0L)
[a604c3]898          mpz_mul_ui(u->n,u->n,SR_TO_INT(b));
[35aab3]899        else
900        {
[a604c3]901          mpz_mul_ui(u->n,u->n,-SR_TO_INT(b));
902          mpz_neg(u->z,u->z);
[35aab3]903        }
904      }
905      else
906      // long z / short b -> z/b
907      {
[a604c3]908        //mpz_set(u->z,a->z);
909        mpz_init_set_si(u->n,SR_TO_INT(b));
[35aab3]910      }
911    }
912// ---------- long / long ------------------------------------
913    else
914    {
[a604c3]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);
[35aab3]919    }
920  }
[a604c3]921  if (mpz_isNeg(u->n))
[35aab3]922  {
[a604c3]923    mpz_neg(u->z,u->z);
924    mpz_neg(u->n,u->n);
[35aab3]925  }
[a604c3]926  if (mpz_cmp_si(u->n,(long)1)==0)
[35aab3]927  {
[a604c3]928    mpz_clear(u->n);
[35aab3]929    u->s=3;
[246d15]930    u=nlShort3(u);
[35aab3]931  }
[2b957a]932  nlTest(u, r);
[35aab3]933  return u;
934}
935
936/*2
937* u:= x ^ exp
938*/
[a642c64]939void nlPower (number x,int exp,number * u, const coeffs r)
[35aab3]940{
[6242bb3]941  *u = INT_TO_SR(0); // 0^e, e!=0
[ef7324]942  if (exp==0)
943    *u= INT_TO_SR(1);
944  else if (!nlIsZero(x,r))
[35aab3]945  {
[2b957a]946    nlTest(x, r);
[35aab3]947    number aa=NULL;
948    if (SR_HDL(x) & SR_INT)
949    {
950      aa=nlRInit(SR_TO_INT(x));
951      x=aa;
952    }
[ee7dcf]953    else if (x->s==0)
[a642c64]954      nlNormalize(x,r);
[896561]955    *u=ALLOC_RNUMBER();
[35aab3]956#if defined(LDEBUG)
957    (*u)->debug=123456;
958#endif
[a604c3]959    mpz_init((*u)->z);
960    mpz_pow_ui((*u)->z,x->z,(unsigned long)exp);
[35aab3]961    if (x->s<2)
962    {
[a604c3]963      if (mpz_cmp_si(x->n,(long)1)==0)
[35aab3]964      {
965        x->s=3;
[a604c3]966        mpz_clear(x->n);
[35aab3]967      }
968      else
969      {
[a604c3]970        mpz_init((*u)->n);
971        mpz_pow_ui((*u)->n,x->n,(unsigned long)exp);
[35aab3]972      }
973    }
974    (*u)->s = x->s;
[246d15]975    if ((*u)->s==3) *u=nlShort3(*u);
[35aab3]976    if (aa!=NULL)
977    {
[a604c3]978      mpz_clear(aa->z);
[896561]979      FREE_RNUMBER(aa);
[35aab3]980    }
981  }
982#ifdef LDEBUG
[deca7f6]983  if (exp<0) Print("nlPower: neg. exp. %d\n",exp);
[2b957a]984  nlTest(*u, r);
[35aab3]985#endif
986}
987
988
989/*2
990* za >= 0 ?
991*/
[2b957a]992BOOLEAN nlGreaterZero (number a, const coeffs r)
[35aab3]993{
[2b957a]994  nlTest(a, r);
[246d15]995  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>1L /* represents number(0) */;
[a604c3]996  return (!mpz_isNeg(a->z));
[35aab3]997}
998
999/*2
1000* a > b ?
1001*/
[a642c64]1002BOOLEAN nlGreater (number a, number b, const coeffs r)
[35aab3]1003{
[2b957a]1004  nlTest(a, r);
1005  nlTest(b, r);
[a642c64]1006  number re;
[35aab3]1007  BOOLEAN rr;
[a642c64]1008  re=nlSub(a,b,r);
1009  rr=(!nlIsZero(re,r)) && (nlGreaterZero(re,r));
1010  nlDelete(&re,r);
[35aab3]1011  return rr;
1012}
1013
1014/*2
1015* a == -1 ?
1016*/
[2b957a]1017BOOLEAN nlIsMOne (number a, const coeffs r)
[35aab3]1018{
1019#ifdef LDEBUG
1020  if (a==NULL) return FALSE;
[2b957a]1021  nlTest(a, r);
[35aab3]1022#endif
[246d15]1023  //if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1L));
[5f4ae4]1024  //return FALSE;
[246d15]1025  return  (a==INT_TO_SR(-1L));
[35aab3]1026}
1027
1028/*2
1029* result =gcd(a,b)
1030*/
[a642c64]1031number nlGcd(number a, number b, const coeffs r)
[35aab3]1032{
1033  number result;
[2b957a]1034  nlTest(a, r);
1035  nlTest(b, r);
[35aab3]1036  //nlNormalize(a);
1037  //nlNormalize(b);
[246d15]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);
[f553f3]1041  if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
[a642c64]1042    return nlCopy(b,r);
[f553f3]1043  if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
[a642c64]1044    return nlCopy(a,r);
[35aab3]1045  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1046  {
[246d15]1047    long i=SR_TO_INT(a);
1048    long j=SR_TO_INT(b);
1049    if((i==0L)||(j==0L))
[35aab3]1050      return INT_TO_SR(1);
[246d15]1051    long l;
[35aab3]1052    i=ABS(i);
1053    j=ABS(j);
1054    do
1055    {
1056      l=i%j;
1057      i=j;
1058      j=l;
[246d15]1059    } while (l!=0L);
[794463]1060    if (i==POW_2_28)
1061      result=nlRInit(POW_2_28);
1062    else
1063     result=INT_TO_SR(i);
[2b957a]1064    nlTest(result,r);
[794463]1065    return result;
[35aab3]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  {
[794463]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);
[35aab3]1077  }
1078  else
1079  if (SR_HDL(b) & SR_INT)
1080  {
[794463]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);
[35aab3]1087  }
1088  else
1089  {
[a44bcf]1090    result=ALLOC0_RNUMBER();
[a604c3]1091    mpz_init(result->z);
1092    mpz_gcd(result->z,a->z,b->z);
[ea66a97]1093    result->s = 3;
1094  #ifdef LDEBUG
1095    result->debug=123456;
1096  #endif
[246d15]1097    result=nlShort3(result);
[35aab3]1098  }
[2b957a]1099  nlTest(result, r);
[35aab3]1100  return result;
1101}
1102
[18057e]1103number nlShort1(number x) // assume x->s==0/1
[e8a0aa]1104{
[18057e]1105  assume(x->s<2);
[a604c3]1106  if (mpz_cmp_ui(x->z,(long)0)==0)
[e8a0aa]1107  {
[106406]1108    _nlDelete_NoImm(&x);
[e8a0aa]1109    return INT_TO_SR(0);
1110  }
[18057e]1111  if (x->s<2)
1112  {
[a604c3]1113    if (mpz_cmp(x->z,x->n)==0)
[18057e]1114    {
[106406]1115      _nlDelete_NoImm(&x);
[18057e]1116      return INT_TO_SR(1);
1117    }
1118  }
1119  return x;
1120}
[35aab3]1121/*2
1122* simplify x
1123*/
[a642c64]1124void nlNormalize (number &x, const coeffs r)
[35aab3]1125{
1126  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1127    return;
1128  if (x->s==3)
1129  {
[246d15]1130    x=nlShort3_noinline(x);
[2b957a]1131    nlTest(x,r);
[18057e]1132    return;
[35aab3]1133  }
[3ccafbf]1134  else if (x->s==0)
[35aab3]1135  {
[a604c3]1136    if (mpz_cmp_si(x->n,(long)1)==0)
[35aab3]1137    {
[a604c3]1138      mpz_clear(x->n);
[3ccafbf]1139      x->s=3;
[246d15]1140      x=nlShort3(x);
[eb059c5]1141    }
[3ccafbf]1142    else
[eb059c5]1143    {
[a604c3]1144      mpz_t gcd;
1145      mpz_init(gcd);
1146      mpz_gcd(gcd,x->z,x->n);
[3ccafbf]1147      x->s=1;
[a604c3]1148      if (mpz_cmp_si(gcd,(long)1)!=0)
[eb059c5]1149      {
[a604c3]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)
[3ccafbf]1158        {
[a604c3]1159          mpz_clear(x->n);
[3ccafbf]1160          x->s=3;
[246d15]1161          x=nlShort3_noinline(x);
[3ccafbf]1162        }
[eb059c5]1163      }
[a604c3]1164      mpz_clear(gcd);
[35aab3]1165    }
1166  }
[2b957a]1167  nlTest(x, r);
[35aab3]1168}
1169
1170/*2
1171* returns in result->z the lcm(a->z,b->n)
1172*/
[a642c64]1173number nlLcm(number a, number b, const coeffs r)
[35aab3]1174{
1175  number result;
[2b957a]1176  nlTest(a, r);
1177  nlTest(b, r);
[35aab3]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
[a642c64]1182    return nlCopy(a,r);
[35aab3]1183  }
[896561]1184  result=ALLOC_RNUMBER();
[35aab3]1185#if defined(LDEBUG)
1186  result->debug=123456;
1187#endif
1188  result->s=3;
[a604c3]1189  mpz_t gcd;
1190  mpz_init(gcd);
1191  mpz_init(result->z);
[35aab3]1192  if (SR_HDL(a) & SR_INT)
[a604c3]1193    mpz_gcd_ui(gcd,b->n,ABS(SR_TO_INT(a)));
[35aab3]1194  else
[a604c3]1195    mpz_gcd(gcd,a->z,b->n);
1196  if (mpz_cmp_si(gcd,(long)1)!=0)
[35aab3]1197  {
[a604c3]1198    mpz_t bt;
1199    mpz_init_set(bt,b->n);
1200    MPZ_EXACTDIV(bt,bt,gcd);
[35aab3]1201    if (SR_HDL(a) & SR_INT)
[a604c3]1202      mpz_mul_si(result->z,bt,SR_TO_INT(a));
[35aab3]1203    else
[a604c3]1204      mpz_mul(result->z,bt,a->z);
1205    mpz_clear(bt);
[35aab3]1206  }
1207  else
1208    if (SR_HDL(a) & SR_INT)
[a604c3]1209      mpz_mul_si(result->z,b->n,SR_TO_INT(a));
[35aab3]1210    else
[a604c3]1211      mpz_mul(result->z,b->n,a->z);
1212  mpz_clear(gcd);
[246d15]1213  result=nlShort3(result);
[2b957a]1214  nlTest(result, r);
[35aab3]1215  return result;
1216}
1217
[79020f]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)
[35aab3]1221{
[79020f]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)
[35aab3]1233  {
[79020f]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 );
[35aab3]1239  }
[79020f]1240
1241  const unsigned long PP = p;
1242
1243  // numerator modulo char. should fit into int
[e11ed0]1244  number z = n_Init( static_cast<int>(mpz_fdiv_ui(q->z, PP)), Zp );
[79020f]1245
1246  // denominator != 1?
1247  if (q->s!=3)
[35aab3]1248  {
[79020f]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;
[35aab3]1258  }
[79020f]1259
1260  return z;
[35aab3]1261}
1262
[894f5b1]1263#ifdef HAVE_RINGS
1264/*2
[1919112]1265* convert number i (from Q) to GMP and warn if denom != 1
[894f5b1]1266*/
[2518db]1267void nlGMP(number &i, number n, const coeffs r)
[894f5b1]1268{
1269  // Hier brauche ich einfach die GMP Zahl
[2b957a]1270  nlTest(i, r);
[2518db]1271  nlNormalize(i, r);
[894f5b1]1272  if (SR_HDL(i) & SR_INT)
1273  {
[246d15]1274    mpz_set_si((mpz_ptr) n, SR_TO_INT(i));
[0959dc]1275    return;
[894f5b1]1276  }
1277  if (i->s!=3)
1278  {
1279     WarnS("Omitted denominator during coefficient mapping !");
1280  }
[a604c3]1281  mpz_set((mpz_ptr) n, i->z);
[894f5b1]1282}
1283#endif
1284
[35aab3]1285/*2
1286* acces to denominator, other 1 for integers
1287*/
[a642c64]1288number   nlGetDenom(number &n, const coeffs r)
[35aab3]1289{
1290  if (!(SR_HDL(n) & SR_INT))
1291  {
1292    if (n->s==0)
1293    {
[a642c64]1294      nlNormalize(n,r);
[35aab3]1295    }
1296    if (!(SR_HDL(n) & SR_INT))
1297    {
1298      if (n->s!=3)
1299      {
[896561]1300        number u=ALLOC_RNUMBER();
[35aab3]1301        u->s=3;
1302#if defined(LDEBUG)
1303        u->debug=123456;
1304#endif
[a604c3]1305        mpz_init_set(u->z,n->n);
[246d15]1306        u=nlShort3_noinline(u);
[35aab3]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*/
[a642c64]1317number   nlGetNumerator(number &n, const coeffs r)
[35aab3]1318{
1319  if (!(SR_HDL(n) & SR_INT))
1320  {
1321    if (n->s==0)
1322    {
[a642c64]1323      nlNormalize(n,r);
[35aab3]1324    }
1325    if (!(SR_HDL(n) & SR_INT))
1326    {
[896561]1327      number u=ALLOC_RNUMBER();
[35aab3]1328#if defined(LDEBUG)
1329      u->debug=123456;
1330#endif
1331      u->s=3;
[a604c3]1332      mpz_init_set(u->z,n->z);
[35aab3]1333      if (n->s!=3)
1334      {
[246d15]1335        u=nlShort3_noinline(u);
[35aab3]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;
[a604c3]1363    if (((long)a > 0L) && (mpz_isNeg(b->z)))
[35aab3]1364      return FALSE;
[a604c3]1365    if (((long)a < 0L) && (!mpz_isNeg(b->z)))
[35aab3]1366      return FALSE;
[a604c3]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);
[35aab3]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;
[a604c3]1378  if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
[35aab3]1379    return FALSE;
[a604c3]1380  if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
[35aab3]1381    return FALSE;
[a604c3]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);
[35aab3]1391  return bo;
1392}
1393
1394// copy not immediate number a
[e69b7c]1395number _nlCopy_NoImm(number a)
[35aab3]1396{
1397  assume(!((SR_HDL(a) & SR_INT)||(a==NULL)));
[e69b7c]1398  //nlTest(a, r);
[896561]1399  number b=ALLOC_RNUMBER();
[35aab3]1400#if defined(LDEBUG)
1401  b->debug=123456;
1402#endif
1403  switch (a->s)
1404  {
1405    case 0:
1406    case 1:
[a604c3]1407            mpz_init_set(b->n,a->n);
[35aab3]1408    case 3:
[a604c3]1409            mpz_init_set(b->z,a->z);
[35aab3]1410            break;
1411  }
1412  b->s = a->s;
1413  return b;
1414}
1415
[fee4c2]1416void _nlDelete_NoImm(number *a)
[35aab3]1417{
1418  {
1419    switch ((*a)->s)
1420    {
1421      case 0:
1422      case 1:
[a604c3]1423        mpz_clear((*a)->n);
[35aab3]1424      case 3:
[a604c3]1425        mpz_clear((*a)->z);
[35aab3]1426#ifdef LDEBUG
1427        (*a)->s=2;
1428#endif
1429    }
[d663f2]1430    FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
[35aab3]1431  }
1432}
1433
[e69b7c]1434number _nlNeg_NoImm(number a)
[35aab3]1435{
1436  {
[a604c3]1437    mpz_neg(a->z,a->z);
[246d15]1438    if (a->s==3)
[35aab3]1439    {
[246d15]1440      a=nlShort3(a);
[35aab3]1441    }
1442  }
1443  return a;
1444}
1445
[e69b7c]1446number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
[35aab3]1447{
[896561]1448  number u=ALLOC_RNUMBER();
[35aab3]1449#if defined(LDEBUG)
1450  u->debug=123456;
1451#endif
[a604c3]1452  mpz_init(u->z);
[35aab3]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      {
[a604c3]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)
[35aab3]1472        {
[a604c3]1473          mpz_clear(u->z);
[896561]1474          FREE_RNUMBER(u);
[35aab3]1475          return INT_TO_SR(0);
1476        }
[a604c3]1477        if (mpz_cmp(u->z,b->n)==0)
[35aab3]1478        {
[a604c3]1479          mpz_clear(u->z);
[896561]1480          FREE_RNUMBER(u);
[35aab3]1481          return INT_TO_SR(1);
1482        }
[a604c3]1483        mpz_init_set(u->n,b->n);
[35aab3]1484        u->s = 0;
1485        break;
1486      }
1487      case 3:
1488      {
[7447d8]1489        if ((long)a>0L)
[a604c3]1490          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
[35aab3]1491        else
[a604c3]1492          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1493        if (mpz_cmp_ui(u->z,(long)0)==0)
[35aab3]1494        {
[a604c3]1495          mpz_clear(u->z);
[896561]1496          FREE_RNUMBER(u);
[35aab3]1497          return INT_TO_SR(0);
1498        }
1499        u->s = 3;
[246d15]1500        u=nlShort3(u);
[35aab3]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          {
[a604c3]1517            mpz_t x;
1518            mpz_init(x);
[e8a0aa]1519
[a604c3]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);
[e8a0aa]1524
[a604c3]1525            if (mpz_cmp_ui(u->z,(long)0)==0)
[35aab3]1526            {
[a604c3]1527              mpz_clear(u->z);
[896561]1528              FREE_RNUMBER(u);
[35aab3]1529              return INT_TO_SR(0);
1530            }
[a604c3]1531            mpz_init(u->n);
1532            mpz_mul(u->n,a->n,b->n);
1533            if (mpz_cmp(u->z,u->n)==0)
[35aab3]1534            {
[a604c3]1535               mpz_clear(u->z);
1536               mpz_clear(u->n);
[896561]1537               FREE_RNUMBER(u);
[35aab3]1538               return INT_TO_SR(1);
1539            }
1540            u->s = 0;
1541            break;
1542          }
1543          case 3: /* a:1 b:3 */
1544          {
[a604c3]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)
[35aab3]1548            {
[a604c3]1549              mpz_clear(u->z);
[896561]1550              FREE_RNUMBER(u);
[35aab3]1551              return INT_TO_SR(0);
1552            }
[a604c3]1553            if (mpz_cmp(u->z,a->n)==0)
[35aab3]1554            {
[a604c3]1555              mpz_clear(u->z);
[896561]1556              FREE_RNUMBER(u);
[35aab3]1557              return INT_TO_SR(1);
1558            }
[a604c3]1559            mpz_init_set(u->n,a->n);
[35aab3]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          {
[a604c3]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)
[35aab3]1576            {
[a604c3]1577              mpz_clear(u->z);
[896561]1578              FREE_RNUMBER(u);
[35aab3]1579              return INT_TO_SR(0);
1580            }
[a604c3]1581            if (mpz_cmp(u->z,b->n)==0)
[35aab3]1582            {
[a604c3]1583              mpz_clear(u->z);
[896561]1584              FREE_RNUMBER(u);
[35aab3]1585              return INT_TO_SR(1);
1586            }
[a604c3]1587            mpz_init_set(u->n,b->n);
[35aab3]1588            u->s = 0;
1589            break;
1590          }
1591          case 3:
1592          {
[a604c3]1593            mpz_add(u->z,a->z,b->z);
1594            if (mpz_cmp_ui(u->z,(long)0)==0)
[35aab3]1595            {
[a604c3]1596              mpz_clear(u->z);
[896561]1597              FREE_RNUMBER(u);
[35aab3]1598              return INT_TO_SR(0);
1599            }
1600            u->s = 3;
[246d15]1601            u=nlShort3(u);
[35aab3]1602            break;
1603          }
1604        }
1605        break;
1606      }
1607    }
1608  }
1609  return u;
1610}
1611
[65c5f5]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
[e69b7c]1751number _nlSub_aNoImm_OR_bNoImm(number a, number b)
[35aab3]1752{
[896561]1753  number u=ALLOC_RNUMBER();
[35aab3]1754#if defined(LDEBUG)
1755  u->debug=123456;
1756#endif
[a604c3]1757  mpz_init(u->z);
[35aab3]1758  if (SR_HDL(a) & SR_INT)
1759  {
1760    switch (b->s)
1761    {
1762      case 0:
1763      case 1:/* a:short, b:1 */
1764      {
[a604c3]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)
[35aab3]1771        {
[a604c3]1772          mpz_clear(u->z);
[896561]1773          FREE_RNUMBER(u);
[35aab3]1774          return INT_TO_SR(0);
1775        }
[a604c3]1776        if (mpz_cmp(u->z,b->n)==0)
[35aab3]1777        {
[a604c3]1778          mpz_clear(u->z);
[896561]1779          FREE_RNUMBER(u);
[35aab3]1780          return INT_TO_SR(1);
1781        }
[a604c3]1782        mpz_init_set(u->n,b->n);
[35aab3]1783        u->s = 0;
1784        break;
1785      }
1786      case 3:
1787      {
[7447d8]1788        if ((long)a>0L)
[35aab3]1789        {
[a604c3]1790          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
1791          mpz_neg(u->z,u->z);
[35aab3]1792        }
1793        else
1794        {
[a604c3]1795          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
1796          mpz_neg(u->z,u->z);
[35aab3]1797        }
[a604c3]1798        if (mpz_cmp_ui(u->z,(long)0)==0)
[35aab3]1799        {
[a604c3]1800          mpz_clear(u->z);
[896561]1801          FREE_RNUMBER(u);
[35aab3]1802          return INT_TO_SR(0);
1803        }
1804        u->s = 3;
[246d15]1805        u=nlShort3(u);
[35aab3]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      {
[a604c3]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)
[35aab3]1823        {
[a604c3]1824          mpz_clear(u->z);
[896561]1825          FREE_RNUMBER(u);
[35aab3]1826          return INT_TO_SR(0);
1827        }
[a604c3]1828        if (mpz_cmp(u->z,a->n)==0)
[35aab3]1829        {
[a604c3]1830          mpz_clear(u->z);
[896561]1831          FREE_RNUMBER(u);
[35aab3]1832          return INT_TO_SR(1);
1833        }
[a604c3]1834        mpz_init_set(u->n,a->n);
[35aab3]1835        u->s = 0;
1836        break;
1837      }
1838      case 3:
1839      {
[7447d8]1840        if ((long)b>0L)
[35aab3]1841        {
[a604c3]1842          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
[35aab3]1843        }
1844        else
1845        {
[a604c3]1846          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
[35aab3]1847        }
[a604c3]1848        if (mpz_cmp_ui(u->z,(long)0)==0)
[35aab3]1849        {
[a604c3]1850          mpz_clear(u->z);
[896561]1851          FREE_RNUMBER(u);
[35aab3]1852          return INT_TO_SR(0);
1853        }
1854        u->s = 3;
[246d15]1855        u=nlShort3(u);
[35aab3]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          {
[a604c3]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)
[35aab3]1882            {
[a604c3]1883              mpz_clear(u->z);
[896561]1884              FREE_RNUMBER(u);
[35aab3]1885              return INT_TO_SR(0);
1886            }
[a604c3]1887            mpz_init(u->n);
1888            mpz_mul(u->n,a->n,b->n);
1889            if (mpz_cmp(u->z,u->n)==0)
[35aab3]1890            {
[a604c3]1891              mpz_clear(u->z);
1892              mpz_clear(u->n);
[896561]1893              FREE_RNUMBER(u);
[35aab3]1894              return INT_TO_SR(1);
1895            }
1896            u->s = 0;
1897            break;
1898          }
1899          case 3: /* a:1, b:3 */
1900          {
[a604c3]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)
[35aab3]1907            {
[a604c3]1908              mpz_clear(u->z);
[896561]1909              FREE_RNUMBER(u);
[35aab3]1910              return INT_TO_SR(0);
1911            }
[a604c3]1912            if (mpz_cmp(u->z,a->n)==0)
[35aab3]1913            {
[a604c3]1914              mpz_clear(u->z);
[896561]1915              FREE_RNUMBER(u);
[35aab3]1916              return INT_TO_SR(1);
1917            }
[a604c3]1918            mpz_init_set(u->n,a->n);
[35aab3]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          {
[a604c3]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)
[35aab3]1938            {
[a604c3]1939              mpz_clear(u->z);
[896561]1940              FREE_RNUMBER(u);
[35aab3]1941              return INT_TO_SR(0);
1942            }
[a604c3]1943            if (mpz_cmp(u->z,b->n)==0)
[35aab3]1944            {
[a604c3]1945              mpz_clear(u->z);
[896561]1946              FREE_RNUMBER(u);
[35aab3]1947              return INT_TO_SR(1);
1948            }
[a604c3]1949            mpz_init_set(u->n,b->n);
[35aab3]1950            u->s = 0;
1951            break;
1952          }
1953          case 3: /* a:3 , b:3 */
1954          {
[a604c3]1955            mpz_sub(u->z,a->z,b->z);
1956            if (mpz_cmp_ui(u->z,(long)0)==0)
[35aab3]1957            {
[a604c3]1958              mpz_clear(u->z);
[896561]1959              FREE_RNUMBER(u);
[35aab3]1960              return INT_TO_SR(0);
1961            }
1962            u->s = 3;
[246d15]1963            u=nlShort3(u);
[35aab3]1964            break;
1965          }
1966        }
1967        break;
1968      }
1969    }
1970  }
1971  return u;
1972}
1973
1974// a and b are intermediate, but a*b not
[e69b7c]1975number _nlMult_aImm_bImm_rNoImm(number a, number b)
[35aab3]1976{
[896561]1977  number u=ALLOC_RNUMBER();
[35aab3]1978#if defined(LDEBUG)
1979  u->debug=123456;
1980#endif
1981  u->s=3;
[246d15]1982  mpz_init_set_si(u->z,SR_TO_INT(a));
1983  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
[35aab3]1984  return u;
1985}
1986
1987// a or b are not immediate
[e69b7c]1988number _nlMult_aNoImm_OR_bNoImm(number a, number b)
[35aab3]1989{
1990  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
[896561]1991  number u=ALLOC_RNUMBER();
[35aab3]1992#if defined(LDEBUG)
1993  u->debug=123456;
1994#endif
[a604c3]1995  mpz_init(u->z);
[35aab3]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;
[7447d8]2006    if ((long)a>0L)
[35aab3]2007    {
[a604c3]2008      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
[35aab3]2009    }
2010    else
2011    {
2012      if (a==INT_TO_SR(-1))
2013      {
[a604c3]2014        mpz_set(u->z,b->z);
2015        mpz_neg(u->z,u->z);
[35aab3]2016        u->s=b->s;
2017      }
2018      else
2019      {
[a604c3]2020        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2021        mpz_neg(u->z,u->z);
[35aab3]2022      }
2023    }
2024    if (u->s<2)
2025    {
[a604c3]2026      if (mpz_cmp(u->z,b->n)==0)
[35aab3]2027      {
[a604c3]2028        mpz_clear(u->z);
[896561]2029        FREE_RNUMBER(u);
[35aab3]2030        return INT_TO_SR(1);
2031      }
[a604c3]2032      mpz_init_set(u->n,b->n);
[35aab3]2033    }
2034    else //u->s==3
2035    {
[246d15]2036      u=nlShort3(u);
[35aab3]2037    }
2038  }
2039  else
2040  {
[a604c3]2041    mpz_mul(u->z,a->z,b->z);
[35aab3]2042    u->s = 0;
2043    if(a->s==3)
2044    {
2045      if(b->s==3)
2046      {
2047        u->s = 3;
2048      }
2049      else
2050      {
[a604c3]2051        if (mpz_cmp(u->z,b->n)==0)
[35aab3]2052        {
[a604c3]2053          mpz_clear(u->z);
[896561]2054          FREE_RNUMBER(u);
[35aab3]2055          return INT_TO_SR(1);
2056        }
[a604c3]2057        mpz_init_set(u->n,b->n);
[35aab3]2058      }
2059    }
2060    else
2061    {
2062      if(b->s==3)
2063      {
[a604c3]2064        if (mpz_cmp(u->z,a->n)==0)
[35aab3]2065        {
[a604c3]2066          mpz_clear(u->z);
[896561]2067          FREE_RNUMBER(u);
[35aab3]2068          return INT_TO_SR(1);
2069        }
[a604c3]2070        mpz_init_set(u->n,a->n);
[35aab3]2071      }
2072      else
2073      {
[a604c3]2074        mpz_init(u->n);
2075        mpz_mul(u->n,a->n,b->n);
2076        if (mpz_cmp(u->z,u->n)==0)
[35aab3]2077        {
[a604c3]2078          mpz_clear(u->z);
2079          mpz_clear(u->n);
[896561]2080          FREE_RNUMBER(u);
[35aab3]2081          return INT_TO_SR(1);
2082        }
2083      }
2084    }
2085  }
2086  return u;
2087}
2088
[371bc5f]2089/*2
[fa31d2]2090* copy a to b for mapping
[371bc5f]2091*/
[fa31d2]2092number nlCopyMap(number a, const coeffs src, const coeffs dst)
[371bc5f]2093{
[fa31d2]2094  assume( getCoeffType(dst) == ID );
2095  assume( getCoeffType(src) == ID );
[9bb5457]2096
[371bc5f]2097  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2098  {
2099    return a;
2100  }
[e69b7c]2101  return _nlCopy_NoImm(a);
[371bc5f]2102}
[2b957a]2103
[371bc5f]2104nMapFunc nlSetMap(const coeffs src, const coeffs dst)
2105{
[fa31d2]2106  assume( getCoeffType(dst) == ID );
[2b957a]2107//  assume( getCoeffType(src) == ID );
[9bb5457]2108
[1cce47]2109  if (nCoeff_is_Q(src))
[371bc5f]2110  {
[c14846c]2111    return ndCopyMap;
[371bc5f]2112  }
[1cce47]2113  if (nCoeff_is_Zp(src))
[371bc5f]2114  {
2115    return nlMapP;
2116  }
[1cce47]2117  if (nCoeff_is_R(src))
[371bc5f]2118  {
2119    return nlMapR;
2120  }
[1cce47]2121  if (nCoeff_is_long_R(src))
[371bc5f]2122  {
2123    return nlMapLongR; /* long R -> Q */
2124  }
2125#ifdef HAVE_RINGS
[1cce47]2126  if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
[371bc5f]2127  {
2128    return nlMapGMP;
2129  }
[1cce47]2130  if (nCoeff_is_Ring_2toM(src))
[371bc5f]2131  {
2132    return nlMapMachineInt;
2133  }
2134#endif
2135  return NULL;
2136}
[35aab3]2137/*2
2138* z := i
2139*/
[246d15]2140number nlRInit (long i)
[35aab3]2141{
[896561]2142  number z=ALLOC_RNUMBER();
[35aab3]2143#if defined(LDEBUG)
2144  z->debug=123456;
2145#endif
[246d15]2146  mpz_init_set_si(z->z,i);
[35aab3]2147  z->s = 3;
2148  return z;
2149}
2150
2151/*2
2152* z := i/j
2153*/
[a642c64]2154number nlInit2 (int i, int j, const coeffs r)
[35aab3]2155{
[896561]2156  number z=ALLOC_RNUMBER();
[35aab3]2157#if defined(LDEBUG)
2158  z->debug=123456;
2159#endif
[a604c3]2160  mpz_init_set_si(z->z,(long)i);
2161  mpz_init_set_si(z->n,(long)j);
[35aab3]2162  z->s = 0;
[a642c64]2163  nlNormalize(z,r);
[35aab3]2164  return z;
2165}
2166
[a642c64]2167number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
[a1ab2a]2168{
[896561]2169  number z=ALLOC_RNUMBER();
[a1ab2a]2170#if defined(LDEBUG)
2171  z->debug=123456;
2172#endif
[a604c3]2173  mpz_init_set(z->z,i);
2174  mpz_init_set(z->n,j);
[a1ab2a]2175  z->s = 0;
[a642c64]2176  nlNormalize(z,r);
[a1ab2a]2177  return z;
2178}
2179
[35aab3]2180#else // DO_LINLINE
2181
2182// declare immedate routines
[2f3764]2183number nlRInit (long i);
[35aab3]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);
[65c5f5]2188void    _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b);
[35aab3]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*/
[2b957a]2207LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
[35aab3]2208{
[2b957a]2209  nlTest(a, r);
2210  nlTest(b, r);
[35aab3]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
[2f3764]2216LINLINE number nlInit (long i, const coeffs r)
[35aab3]2217{
2218  number n;
[246d15]2219  LONG ii=(LONG)i;
2220  if ( ((ii << 3) >> 3) == ii ) n=INT_TO_SR(ii);
2221  else                          n=nlRInit(ii);
[2b957a]2222  nlTest(n, r);
[35aab3]2223  return n;
2224}
2225
2226/*2
2227* a == 1 ?
2228*/
[2b957a]2229LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
[35aab3]2230{
2231#ifdef LDEBUG
2232  if (a==NULL) return FALSE;
[2b957a]2233  nlTest(a, r);
[35aab3]2234#endif
[89ca47]2235  return (a==INT_TO_SR(1));
[35aab3]2236}
2237
[121fd9]2238LINLINE BOOLEAN nlIsZero (number a, const coeffs)
[35aab3]2239{
[cc3240c]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
[35aab3]2251  return (a==INT_TO_SR(0));
[cc3240c]2252  #endif
[35aab3]2253}
2254
2255/*2
2256* copy a to b
2257*/
[121fd9]2258LINLINE number nlCopy(number a, const coeffs)
[35aab3]2259{
2260  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2261  {
2262    return a;
2263  }
[e69b7c]2264  return _nlCopy_NoImm(a);
[35aab3]2265}
2266
2267
2268/*2
2269* delete a
2270*/
[a642c64]2271LINLINE void nlDelete (number * a, const coeffs r)
[35aab3]2272{
2273  if (*a!=NULL)
2274  {
[2b957a]2275    nlTest(*a, r);
[35aab3]2276    if ((SR_HDL(*a) & SR_INT)==0)
2277    {
[fee4c2]2278      _nlDelete_NoImm(a);
[35aab3]2279    }
2280    *a=NULL;
2281  }
2282}
2283
2284/*2
2285* za:= - za
2286*/
[e69b7c]2287LINLINE number nlNeg (number a, const coeffs R)
[35aab3]2288{
[e69b7c]2289  nlTest(a, R);
[35aab3]2290  if(SR_HDL(a) &SR_INT)
2291  {
[135e78]2292    LONG r=SR_TO_INT(a);
[246d15]2293    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
[35aab3]2294    else               a=INT_TO_SR(-r);
2295    return a;
2296  }
[630ca9]2297  a = _nlNeg_NoImm(a);
2298  nlTest(a, R);
2299  return a;
2300
[35aab3]2301}
2302
2303/*2
2304* u:= a + b
2305*/
[e69b7c]2306LINLINE number nlAdd (number a, number b, const coeffs R)
[35aab3]2307{
2308  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2309  {
[246d15]2310    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
[35aab3]2311    if ( ((r << 1) >> 1) == r )
[37a003]2312      return (number)(long)r;
[35aab3]2313    else
2314      return nlRInit(SR_TO_INT(r));
2315  }
[630ca9]2316  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2317  nlTest(u, R);
2318  return u;
[35aab3]2319}
2320
[18057e]2321number nlShort1(number a);
[246d15]2322number nlShort3_noinline(number x);
[e8a0aa]2323
[cc3240c]2324LINLINE void nlInpAdd(number &a, number b, const coeffs r)
[e8a0aa]2325{
[65c5f5]2326  // a=a+b
[e8a0aa]2327  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2328  {
[246d15]2329    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
[e8a0aa]2330    if ( ((r << 1) >> 1) == r )
[cc3240c]2331      a=(number)(long)r;
[e8a0aa]2332    else
[cc3240c]2333      a=nlRInit(SR_TO_INT(r));
[e8a0aa]2334  }
2335  else
2336  {
[65c5f5]2337    _nlInpAdd_aNoImm_OR_bNoImm(a,b);
2338    nlTest(a,r);
[e8a0aa]2339  }
2340}
2341
[7dc3e1c]2342LINLINE number nlMult (number a, number b, const coeffs R)
[35aab3]2343{
[2b957a]2344  nlTest(a, R);
2345  nlTest(b, R);
[35aab3]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  {
[12777c]2350    LONG r=(LONG)((unsigned LONG)(SR_HDL(a)-1L))*((unsigned LONG)(SR_HDL(b)>>1));
[246d15]2351    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
[35aab3]2352    {
2353      number u=((number) ((r>>1)+SR_INT));
[246d15]2354      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
[35aab3]2355      return nlRInit(SR_HDL(u)>>2);
2356    }
[630ca9]2357    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2358    nlTest(u, R);
2359    return u;
[9bb5457]2360
[35aab3]2361  }
[630ca9]2362  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2363  nlTest(u, R);
2364  return u;
[9bb5457]2365
[35aab3]2366}
2367
2368
2369/*2
2370* u:= a - b
2371*/
[a642c64]2372LINLINE number nlSub (number a, number b, const coeffs r)
[35aab3]2373{
2374  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2375  {
[246d15]2376    LONG r=SR_HDL(a)-SR_HDL(b)+1;
[35aab3]2377    if ( ((r << 1) >> 1) == r )
2378    {
[246d15]2379      return (number)(long)r;
[35aab3]2380    }
2381    else
2382      return nlRInit(SR_TO_INT(r));
2383  }
[630ca9]2384  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2385  nlTest(u, r);
2386  return u;
[9bb5457]2387
[35aab3]2388}
2389
[a642c64]2390LINLINE void nlInpMult(number &a, number b, const coeffs r)
[4e35a89]2391{
[896561]2392  number aa=a;
2393  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
[4e35a89]2394  {
[a642c64]2395    number n=nlMult(aa,b,r);
[4e35a89]2396    nlDelete(&a,r);
2397    a=n;
2398  }
2399  else
2400  {
[896561]2401    mpz_mul(aa->z,a->z,b->z);
2402    if (aa->s==3)
[4e35a89]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}
[35aab3]2420#endif // DO_LINLINE
2421
2422#ifndef P_NUMBERS_H
2423
[d663f2]2424static void nlMPZ(mpz_t m, number &n, const coeffs r)
2425{
2426  assume( getCoeffType(r) == ID );
[e11ed0]2427
[d663f2]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
[2d2e40]2435static number nlInitMPZ(mpz_t m, const coeffs)
[d663f2]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;
[e11ed0]2441  return z;
[d663f2]2442}
2443
2444
[a642c64]2445void nlInpGcd(number &a, number b, const coeffs r)
[35aab3]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  {
[a604c3]2455    mpz_gcd(a->z,a->z,b->z);
[246d15]2456    a=nlShort3_noinline(a);
[35aab3]2457  }
2458}
[f59722]2459
[a642c64]2460void nlInpIntDiv(number &a, number b, const coeffs r)
[35aab3]2461{
2462  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2463  {
[2b957a]2464    number n=nlIntDiv(a,b, r);
[35aab3]2465    nlDelete(&a,r);
2466    a=n;
2467  }
2468  else
2469  {
[a604c3]2470    if (mpz_isNeg(a->z))
[35aab3]2471    {
[a604c3]2472      if (mpz_isNeg(b->z))
[35aab3]2473      {
[a604c3]2474        mpz_add(a->z,a->z,b->z);
[35aab3]2475      }
2476      else
2477      {
[a604c3]2478        mpz_sub(a->z,a->z,b->z);
[35aab3]2479      }
[a604c3]2480      mpz_add_ui(a->z,a->z,1);
[35aab3]2481    }
2482    else
2483    {
[a604c3]2484      if (mpz_isNeg(b->z))
[35aab3]2485      {
[a604c3]2486        mpz_sub(a->z,a->z,b->z);
[35aab3]2487      }
2488      else
2489      {
[a604c3]2490        mpz_add(a->z,a->z,b->z);
[35aab3]2491      }
[a604c3]2492      mpz_sub_ui(a->z,a->z,1);
[35aab3]2493    }
[a604c3]2494    MPZ_DIV(a->z,a->z,b->z);
[246d15]2495    a=nlShort3_noinline(a);
[35aab3]2496  }
2497}
2498
[a642c64]2499number nlFarey(number nN, number nP, const coeffs r)
[3149a5]2500{
[a604c3]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);
[3149a5]2514  number z=INT_TO_SR(0);
[a604c3]2515  while(mpz_cmp_si(N,(long)0)!=0)
[3149a5]2516  {
[a604c3]2517    mpz_mul(tmp,N,N);
2518    mpz_add(tmp,tmp,tmp);
2519    if (mpz_cmp(tmp,P)<0)
[3149a5]2520    {
[a604c3]2521       if (mpz_isNeg(B))
[26300d]2522       {
[a604c3]2523         mpz_neg(B,B);
2524         mpz_neg(N,N);
[26300d]2525       }
[9f4478]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       }
[3149a5]2545       break;
2546    }
[7453b4]2547    //mpz_mod(D,E,N);
2548    //mpz_div(tmp,E,N);
2549    mpz_divmod(tmp,D,E,N);
[a604c3]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);
[3149a5]2565  return z;
2566}
[bbf7a8]2567
[3216ec]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
[03f7b5]2602void    nlCoeffWrite  (const coeffs, BOOLEAN /*details*/)
[8167af]2603{
2604  PrintS("//   characteristic : 0\n");
2605}
2606
[de27d8]2607number   nlChineseRemainderSym(number *x, number *q,int rl, BOOLEAN sym, const coeffs CF)
[e8c8d5]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
[de27d8]2612  Off(SW_RATIONAL);
[e8c8d5]2613  CFArray X(rl), Q(rl);
2614  int i;
2615  for(i=rl-1;i>=0;i--)
2616  {
[de27d8]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
[e8c8d5]2619  }
2620  CanonicalForm xnew,qnew;
2621  chineseRemainder(X,Q,xnew,qnew);
[de27d8]2622  number n=CF->convFactoryNSingN(xnew,CF);
2623  if (sym)
[e8c8d5]2624  {
[de27d8]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);
[e8c8d5]2635  }
2636  return n;
2637#else
2638  WerrorS("not implemented");
[de27d8]2639  return nlInit(0,CF);
[e8c8d5]2640#endif
2641}
[de27d8]2642number   nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
2643{
2644  return nlChineseRemainderSym(x,q,rl,TRUE,C);
2645}
[e8c8d5]2646
[558f3cc]2647static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
[de88371]2648{
2649  assume(cf != NULL);
2650  assume(getCoeffType(cf) == ID);
[9952bd]2651
2652  numberCollectionEnumerator.Reset();
2653
2654  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2655  {
2656    c = n_Init(1, cf);
2657    return;
2658  }
2659 
[de88371]2660  // all coeffs are given by integers!!!
2661
[4f9652]2662  // part 1, find a small candidate for gcd
2663  number cand1,cand;
2664  int s1,s;
2665  s=2147483647; // max. int
[9952bd]2666
[dc79bd]2667  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
[9952bd]2668
2669  int normalcount = 0;
[4f9652]2670  do
2671  {
[9952bd]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
[4f9652]2678    s1=mpz_size1(cand1->z);
2679    if (s>s1)
2680    {
2681      cand=cand1;
2682      s=s1;
2683    }
[e11ed0]2684  } while (numberCollectionEnumerator.MoveNext() );
2685
[9952bd]2686//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
2687
[d1d27f]2688  cand=nlCopy(cand,cf);
[4f9652]2689  // part 2: compute gcd(cand,all coeffs)
[9952bd]2690
[4f9652]2691  numberCollectionEnumerator.Reset();
[9952bd]2692 
2693  while (numberCollectionEnumerator.MoveNext() )
[4f9652]2694  {
[9952bd]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   
[39d731]2704    if(nlIsOne(cand,cf))
2705    {
[9952bd]2706      c = cand;
2707
[39d731]2708      if(!lc_is_pos)
2709      {
2710        // make the leading coeff positive
[9952bd]2711        c = nlNeg(c, cf);
[39d731]2712        numberCollectionEnumerator.Reset();
[9952bd]2713       
2714        while (numberCollectionEnumerator.MoveNext() )
[39d731]2715        {
[9952bd]2716          number& nn = numberCollectionEnumerator.Current();
2717          nn = nlNeg(nn, cf);
2718        }
[39d731]2719      }
2720      return;
2721    }
[9952bd]2722  } 
[e11ed0]2723
[4f9652]2724  // part3: all coeffs = all coeffs / cand
[9952bd]2725  if (!lc_is_pos)
2726    cand = nlNeg(cand,cf);
2727 
2728  c = cand;
[4f9652]2729  numberCollectionEnumerator.Reset();
[e11ed0]2730
[9952bd]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  }
[de88371]2738}
2739
[558f3cc]2740static void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
[de88371]2741{
2742  assume(cf != NULL);
2743  assume(getCoeffType(cf) == ID);
[9952bd]2744
2745  numberCollectionEnumerator.Reset();
2746 
2747  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2748  {
2749    c = n_Init(1, cf);
[dc79bd]2750//    assume( n_GreaterZero(c, cf) );
[9952bd]2751    return;
2752  }
2753
[d1d27f]2754  // all coeffs are given by integers after returning from this routine
2755
2756  // part 1, collect product of all denominators /gcds
[e11ed0]2757  number cand;
[d1d27f]2758  cand=ALLOC_RNUMBER();
[9952bd]2759#if defined(LDEBUG)
[d1d27f]2760  cand->debug=123456;
[9952bd]2761#endif
[d1d27f]2762  cand->s=3;
2763
2764  int s=0;
[9952bd]2765 
[dc79bd]2766  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
[9952bd]2767
[d1d27f]2768  do
2769  {
[560a3d]2770    number& cand1 = numberCollectionEnumerator.Current();
[e11ed0]2771
[d1d27f]2772    if (!(SR_HDL(cand1)&SR_INT))
2773    {
[e11ed0]2774      nlNormalize(cand1, cf);
[d1d27f]2775      if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
[9952bd]2776      && (cand1->s==1))             // and is a normalised rational
[d1d27f]2777      {
2778        if (s==0) // first denom, we meet
2779        {
[9952bd]2780          mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
[e11ed0]2781          s=1;
[d1d27f]2782        }
2783        else // we have already something
2784        {
[9952bd]2785          mpz_lcm(cand->z, cand->z, cand1->n);
[d1d27f]2786        }
2787      }
2788    }
[9952bd]2789  }
2790  while (numberCollectionEnumerator.MoveNext() );
2791 
[de88371]2792
[d1d27f]2793  if (s==0) // nothing to do, all coeffs are already integers
2794  {
[9952bd]2795//    mpz_clear(tmp);
[e11ed0]2796    FREE_RNUMBER(cand);
[39d731]2797    if (lc_is_pos)
2798      c=nlInit(1,cf);
2799    else
2800    {
2801      // make the leading coeff positive
2802      c=nlInit(-1,cf);
[9952bd]2803
2804      // TODO: incorporate the following into the loop below?
[39d731]2805      numberCollectionEnumerator.Reset();
[9952bd]2806      while (numberCollectionEnumerator.MoveNext() )
[39d731]2807      {
[9952bd]2808        number& n = numberCollectionEnumerator.Current();
2809        n = nlNeg(n, cf);
2810      } 
[39d731]2811    }
[dc79bd]2812//    assume( n_GreaterZero(c, cf) );
[d1d27f]2813    return;
2814  }
[9952bd]2815
2816  cand = nlShort3(cand);
[e11ed0]2817
[d1d27f]2818  // part2: all coeffs = all coeffs * cand
[e11ed0]2819  // make the lead coeff positive
[d1d27f]2820  numberCollectionEnumerator.Reset();
[9952bd]2821 
2822  if (!lc_is_pos)
2823    cand = nlNeg(cand, cf);
2824 
[e11ed0]2825  c = cand;
[9952bd]2826 
2827  while (numberCollectionEnumerator.MoveNext() )
[d1d27f]2828  {
[560a3d]2829    number &n = numberCollectionEnumerator.Current();
2830    n_InpMult(n, cand, cf);
[9952bd]2831  }
[e11ed0]2832
[de88371]2833}
2834
[9bb5457]2835BOOLEAN nlInitChar(coeffs r, void*)
[e3c02d2]2836{
[73a9ffb]2837  assume( getCoeffType(r) == ID );
2838  //const int ch = (int)(long)(p);
[9bb5457]2839
[e3c02d2]2840  r->cfKillChar=NULL;
[121fd9]2841  r->nCoeffIsEqual=ndCoeffIsEqual;
[577888]2842  r->cfKillChar = ndKillChar; /* dummy */
[e3c02d2]2843
[d663f2]2844  r->cfInitMPZ = nlInitMPZ;
2845  r->cfMPZ  = nlMPZ;
[e11ed0]2846
[7bbbef]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;
[bbf7a8]2854  r->cfInit = nlInit;
[7bbbef]2855  r->cfSize  = nlSize;
2856  r->cfInt  = nlInt;
[e11ed0]2857
[de27d8]2858  r->cfChineseRemainder=nlChineseRemainderSym;
[e8c8d5]2859  r->cfFarey=nlFarey;
[bbf7a8]2860  #ifdef HAVE_RINGS
[7bbbef]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
[bbf7a8]2865  #endif
[7bbbef]2866  r->cfNeg   = nlNeg;
2867  r->cfInvers= nlInvers;
[2b957a]2868  r->cfCopy  = nlCopy;
2869  r->cfRePart = nlCopy;
[7bbbef]2870  //r->cfImPart = ndReturn0;
[ce1f78]2871  r->cfWriteLong = nlWrite;
[7bbbef]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;
[bbf7a8]2881  r->cfGetDenom = nlGetDenom;
2882  r->cfGetNumerator = nlGetNumerator;
[7bbbef]2883  r->cfGcd  = nlGcd;
[3216ec]2884  r->cfExtGcd = nlExtGcd; // only for ring stuff and Z
[7bbbef]2885  r->cfLcm  = nlLcm;
[bbf7a8]2886  r->cfDelete= nlDelete;
[e3c02d2]2887  r->cfSetMap = nlSetMap;
[7bbbef]2888  //r->cfName = ndName;
2889  r->cfInpMult=nlInpMult;
2890  r->cfInit_bigint=nlCopyMap;
[8167af]2891  r->cfCoeffWrite=nlCoeffWrite;
[de88371]2892
2893  r->cfClearContent = nlClearContent;
2894  r->cfClearDenominators = nlClearDenominators;
[e11ed0]2895
[e3c02d2]2896#ifdef LDEBUG
[bbf7a8]2897  // debug stuff
[7bbbef]2898  r->cfDBTest=nlDBTest;
[e3c02d2]2899#endif
[e5d267]2900#ifdef HAVE_FACTORY
2901  r->convSingNFactoryN=nlConvSingNFactoryN;
2902  r->convFactoryNSingN=nlConvFactoryNSingN;
2903#endif
[039016]2904
[8835e2b]2905  // the variables: general stuff (required)
[bbf7a8]2906  r->nNULL = INT_TO_SR(0);
[aff5ae]2907  //r->type = n_Q;
[bbf7a8]2908  r->ch = 0;
2909  r->has_simple_Alloc=FALSE;
2910  r->has_simple_Inverse=FALSE;
[9bb5457]2911
[8835e2b]2912  // variables for this type of coeffs:
2913  // (none)
[5d594a9]2914  return FALSE;
[e3c02d2]2915}
[35aab3]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  {
[246d15]2944    LONG ui=(int)mpz_get_si(&r->z);
[35aab3]2945    if ((((ui<<3)>>3)==ui)
[a604c3]2946    && (mpz_cmp_si(x->z,(long)ui)==0))
[35aab3]2947    {
2948      mpz_clear(&r->z);
[d663f2]2949      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
[35aab3]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.