source: git/libpolys/coeffs/longrat.cc @ 868d375

spielwiese
Last change on this file since 868d375 was 6fb433, checked in by Hans Schoenemann <hannes@…>, 2 years ago
simplify nlCopy (no NULL)
  • Property mode set to 100644
File size: 73.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
6*/
7
8#include "misc/auxiliary.h"
9
10#include "factory/factory.h"
11
12#include "misc/sirandom.h"
13#include "misc/prime.h"
14#include "reporter/reporter.h"
15
16#include "coeffs/coeffs.h"
17#include "coeffs/numbers.h"
18#include "coeffs/rmodulon.h" // ZnmInfo
19#include "coeffs/longrat.h"
20#include "coeffs/shortfl.h"
21#include "coeffs/modulop.h"
22#include "coeffs/mpr_complex.h"
23
24#include <string.h>
25#include <float.h>
26
27// allow inlining only from p_Numbers.h and if ! LDEBUG
28#if defined(DO_LINLINE) && defined(P_NUMBERS_H) && !defined(LDEBUG)
29#define LINLINE static FORCE_INLINE
30#else
31#define LINLINE
32#undef DO_LINLINE
33#endif // DO_LINLINE
34
35LINLINE BOOLEAN  nlEqual(number a, number b, const coeffs r);
36LINLINE number   nlInit(long i, const coeffs r);
37LINLINE BOOLEAN  nlIsOne(number a, const coeffs r);
38LINLINE BOOLEAN  nlIsZero(number za, const coeffs r);
39LINLINE number   nlCopy(number a, const coeffs r);
40LINLINE number   nl_Copy(number a, const coeffs r);
41LINLINE void     nlDelete(number *a, const coeffs r);
42LINLINE number   nlNeg(number za, const coeffs r);
43LINLINE number   nlAdd(number la, number li, const coeffs r);
44LINLINE number   nlSub(number la, number li, const coeffs r);
45LINLINE number   nlMult(number a, number b, const coeffs r);
46LINLINE void nlInpAdd(number &a, number b, const coeffs r);
47LINLINE void nlInpMult(number &a, number b, const coeffs r);
48
49number nlRInit (long i);
50
51
52// number nlInitMPZ(mpz_t m, const coeffs r);
53// void nlMPZ(mpz_t m, number &n, const coeffs r);
54
55void     nlNormalize(number &x, const coeffs r);
56
57number   nlGcd(number a, number b, const coeffs r);
58number nlExtGcd(number a, number b, number *s, number *t, const coeffs);
59number   nlNormalizeHelper(number a, number b, const coeffs r);   /*special routine !*/
60BOOLEAN  nlGreater(number a, number b, const coeffs r);
61BOOLEAN  nlIsMOne(number a, const coeffs r);
62long     nlInt(number &n, const coeffs r);
63number   nlBigInt(number &n);
64
65BOOLEAN  nlGreaterZero(number za, const coeffs r);
66number   nlInvers(number a, const coeffs r);
67number   nlDiv(number a, number b, const coeffs r);
68number   nlExactDiv(number a, number b, const coeffs r);
69number   nlIntDiv(number a, number b, const coeffs r);
70number   nlIntMod(number a, number b, const coeffs r);
71void     nlPower(number x, int exp, number *lu, const coeffs r);
72const char *   nlRead (const char *s, number *a, const coeffs r);
73void     nlWrite(number a, const coeffs r);
74
75number   nlFarey(number nN, number nP, const coeffs CF);
76
77#ifdef LDEBUG
78BOOLEAN  nlDBTest(number a, const char *f, const int l);
79#endif
80
81nMapFunc nlSetMap(const coeffs src, const coeffs dst);
82
83// in-place operations
84void nlInpIntDiv(number &a, number b, const coeffs r);
85
86#ifdef LDEBUG
87#define nlTest(a, r) nlDBTest(a,__FILE__,__LINE__, r)
88BOOLEAN nlDBTest(number a, const char *f,int l, const coeffs r);
89#else
90#define nlTest(a, r) do {} while (0)
91#endif
92
93
94// 64 bit version:
95//#if SIZEOF_LONG == 8
96#if 0
97#define MAX_NUM_SIZE 60
98#define POW_2_28 (1L<<60)
99#define POW_2_28_32 (1L<<28)
100#define LONG long
101#else
102#define MAX_NUM_SIZE 28
103#define POW_2_28 (1L<<28)
104#define POW_2_28_32 (1L<<28)
105#define LONG int
106#endif
107
108
109static inline number nlShort3(number x) // assume x->s==3
110{
111  assume(x->s==3);
112  if (mpz_sgn1(x->z)==0)
113  {
114    mpz_clear(x->z);
115    FREE_RNUMBER(x);
116    return INT_TO_SR(0);
117  }
118  if (mpz_size1(x->z)<=MP_SMALL)
119  {
120    LONG ui=mpz_get_si(x->z);
121    if ((((ui<<3)>>3)==ui)
122    && (mpz_cmp_si(x->z,(long)ui)==0))
123    {
124      mpz_clear(x->z);
125      FREE_RNUMBER(x);
126      return INT_TO_SR(ui);
127    }
128  }
129  return x;
130}
131
132#ifndef LONGRAT_CC
133#define LONGRAT_CC
134
135#ifndef BYTES_PER_MP_LIMB
136#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
137#endif
138
139//#define SR_HDL(A) ((long)(A))
140/*#define SR_INT    1L*/
141/*#define INT_TO_SR(INT)  ((number) (((long)INT << 2) + SR_INT))*/
142// #define SR_TO_INT(SR)   (((long)SR) >> 2)
143
144#define MP_SMALL 1
145//#define mpz_isNeg(A) (mpz_sgn1(A)<0)
146#define mpz_isNeg(A) ((A)->_mp_size<0)
147#define mpz_limb_size(A) ((A)->_mp_size)
148#define mpz_limb_d(A) ((A)->_mp_d)
149
150void    _nlDelete_NoImm(number *a);
151
152/***************************************************************
153 *
154 * Routines which are never inlined by p_Numbers.h
155 *
156 *******************************************************************/
157#ifndef P_NUMBERS_H
158
159number nlShort3_noinline(number x) // assume x->s==3
160{
161  return nlShort3(x);
162}
163
164static number nlInitMPZ(mpz_t m, const coeffs)
165{
166  number z = ALLOC_RNUMBER();
167  z->s = 3;
168  #ifdef LDEBUG
169  z->debug=123456;
170  #endif
171  mpz_init_set(z->z, m);
172  z=nlShort3(z);
173  return z;
174}
175
176#if (__GNU_MP_VERSION*10+__GNU_MP_VERSION_MINOR < 31)
177void mpz_mul_si (mpz_ptr r, mpz_srcptr s, long int si)
178{
179  if (si>=0)
180    mpz_mul_ui(r,s,si);
181  else
182  {
183    mpz_mul_ui(r,s,-si);
184    mpz_neg(r,r);
185  }
186}
187#endif
188
189static number nlMapP(number from, const coeffs src, const coeffs dst)
190{
191  assume( getCoeffType(src) == n_Zp );
192
193  number to = nlInit(npInt(from,src), dst); // FIXME? TODO? // extern long     npInt         (number &n, const coeffs r);
194
195  return to;
196}
197
198static number nlMapLongR(number from, const coeffs src, const coeffs dst);
199static number nlMapR(number from, const coeffs src, const coeffs dst);
200
201
202#ifdef HAVE_RINGS
203/*2
204* convert from a GMP integer
205*/
206static inline number nlMapGMP(number from, const coeffs /*src*/, const coeffs dst)
207{
208  return nlInitMPZ((mpz_ptr)from,dst);
209}
210
211number nlMapZ(number from, const coeffs src, const coeffs dst)
212{
213  if (SR_HDL(from) & SR_INT)
214  {
215    return from;
216  }
217  return nlInitMPZ((mpz_ptr)from,dst);
218}
219
220/*2
221* convert from an machine long
222*/
223number nlMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
224{
225  number z=ALLOC_RNUMBER();
226#if defined(LDEBUG)
227  z->debug=123456;
228#endif
229  mpz_init_set_ui(z->z,(unsigned long) from);
230  z->s = 3;
231  z=nlShort3(z);
232  return z;
233}
234#endif
235
236
237#ifdef LDEBUG
238BOOLEAN nlDBTest(number a, const char *f,const int l, const coeffs /*r*/)
239{
240  if (a==NULL)
241  {
242    Print("!!longrat: NULL in %s:%d\n",f,l);
243    return FALSE;
244  }
245  //if ((int)a==1) Print("!! 0x1 as number ? %s %d\n",f,l);
246  if ((((long)a)&3L)==3L)
247  {
248    Print(" !!longrat:ptr(3) in %s:%d\n",f,l);
249    return FALSE;
250  }
251  if ((((long)a)&3L)==1L)
252  {
253    if (((((LONG)(long)a)<<1)>>1)!=((LONG)(long)a))
254    {
255      Print(" !!longrat:arith:%lx in %s:%d\n",(long)a, f,l);
256      return FALSE;
257    }
258    return TRUE;
259  }
260  /* TODO: If next line is active, then computations in algebraic field
261           extensions over Q will throw a lot of assume violations although
262           everything is computed correctly and no seg fault appears.
263           Maybe the test is not appropriate in this case. */
264  omCheckIf(omCheckAddrSize(a,sizeof(*a)), return FALSE);
265  if (a->debug!=123456)
266  {
267    Print("!!longrat:debug:%d in %s:%d\n",a->debug,f,l);
268    a->debug=123456;
269    return FALSE;
270  }
271  if ((a->s<0)||(a->s>4))
272  {
273    Print("!!longrat:s=%d in %s:%d\n",a->s,f,l);
274    return FALSE;
275  }
276  /* TODO: If next line is active, then computations in algebraic field
277           extensions over Q will throw a lot of assume violations although
278           everything is computed correctly and no seg fault appears.
279           Maybe the test is not appropriate in this case. */
280  //omCheckAddrSize(a->z[0]._mp_d,a->z[0]._mp_alloc*BYTES_PER_MP_LIMB);
281  if (a->z[0]._mp_alloc==0)
282    Print("!!longrat:z->alloc=0 in %s:%d\n",f,l);
283
284  if (a->s<2)
285  {
286    if ((a->n[0]._mp_d[0]==0)&&(a->n[0]._mp_alloc<=1))
287    {
288      Print("!!longrat: n==0 in %s:%d\n",f,l);
289      return FALSE;
290    }
291    /* TODO: If next line is active, then computations in algebraic field
292             extensions over Q will throw a lot of assume violations although
293             everything is computed correctly and no seg fault appears.
294             Maybe the test is not appropriate in this case. */
295    //omCheckIf(omCheckAddrSize(a->n[0]._mp_d,a->n[0]._mp_alloc*BYTES_PER_MP_LIMB), return FALSE);
296    if (a->z[0]._mp_alloc==0)
297      Print("!!longrat:n->alloc=0 in %s:%d\n",f,l);
298    if ((mpz_size1(a->n) ==1) && (mpz_cmp_si(a->n,1L)==0))
299    {
300      Print("!!longrat:integer as rational in %s:%d\n",f,l);
301      mpz_clear(a->n); a->s=3;
302      return FALSE;
303    }
304    else if (mpz_isNeg(a->n))
305    {
306      Print("!!longrat:div. by negative in %s:%d\n",f,l);
307      mpz_neg(a->z,a->z);
308      mpz_neg(a->n,a->n);
309      return FALSE;
310    }
311    return TRUE;
312  }
313  //if (a->s==2)
314  //{
315  //  Print("!!longrat:s=2 in %s:%d\n",f,l);
316  //  return FALSE;
317  //}
318  if (mpz_size1(a->z)>MP_SMALL) return TRUE;
319  LONG ui=(LONG)mpz_get_si(a->z);
320  if ((((ui<<3)>>3)==ui)
321  && (mpz_cmp_si(a->z,(long)ui)==0))
322  {
323    Print("!!longrat:im int %d in %s:%d\n",ui,f,l);
324    return FALSE;
325  }
326  return TRUE;
327}
328#endif
329
330static CanonicalForm nlConvSingNFactoryN( number n, const BOOLEAN setChar, const coeffs /*r*/ )
331{
332  if (setChar) setCharacteristic( 0 );
333
334  CanonicalForm term;
335  if ( SR_HDL(n) & SR_INT )
336  {
337    long nn=SR_TO_INT(n);
338    term = nn;
339  }
340  else
341  {
342    if ( n->s == 3 )
343    {
344      mpz_t dummy;
345      long lz=mpz_get_si(n->z);
346      if (mpz_cmp_si(n->z,lz)==0) term=lz;
347      else
348      {
349        mpz_init_set( dummy,n->z );
350        term = make_cf( dummy );
351      }
352    }
353    else
354    {
355      // assume s==0 or s==1
356      mpz_t num, den;
357      On(SW_RATIONAL);
358      mpz_init_set( num, n->z );
359      mpz_init_set( den, n->n );
360      term = make_cf( num, den, ( n->s != 1 ));
361    }
362  }
363  return term;
364}
365
366number nlRInit (long i);
367
368static number nlConvFactoryNSingN( const CanonicalForm f, const coeffs r)
369{
370  if (f.isImm())
371  {
372    return nlInit(f.intval(),r);
373  }
374  else
375  {
376    number z = ALLOC_RNUMBER();
377#if defined(LDEBUG)
378    z->debug=123456;
379#endif
380    gmp_numerator( f, z->z );
381    if ( f.den().isOne() )
382    {
383      z->s = 3;
384      z=nlShort3(z);
385    }
386    else
387    {
388      gmp_denominator( f, z->n );
389      z->s = 1;
390    }
391    return z;
392  }
393}
394
395static number nlMapR(number from, const coeffs src, const coeffs dst)
396{
397  assume( getCoeffType(src) == n_R );
398
399  double f=nrFloat(from); // FIXME? TODO? // extern float   nrFloat(number n);
400  if (f==0.0) return INT_TO_SR(0);
401  int f_sign=1;
402  if (f<0.0)
403  {
404    f_sign=-1;
405    f=-f;
406  }
407  int i=0;
408  mpz_t h1;
409  mpz_init_set_ui(h1,1);
410  while((FLT_RADIX*f) < DBL_MAX && i<DBL_MANT_DIG)
411  {
412    f*=FLT_RADIX;
413    mpz_mul_ui(h1,h1,FLT_RADIX);
414    i++;
415  }
416  number re=nlRInit(1);
417  mpz_set_d(re->z,f);
418  memcpy(&(re->n),&h1,sizeof(h1));
419  re->s=0; /* not normalized */
420  if(f_sign==-1) re=nlNeg(re,dst);
421  nlNormalize(re,dst);
422  return re;
423}
424
425static number nlMapR_BI(number from, const coeffs src, const coeffs dst)
426{
427  assume( getCoeffType(src) == n_R );
428
429  double f=nrFloat(from); // FIXME? TODO? // extern float   nrFloat(number n);
430  if (f==0.0) return INT_TO_SR(0);
431  long l=long(f);
432  return nlInit(l,dst);
433}
434
435static number nlMapLongR(number from, const coeffs src, const coeffs dst)
436{
437  assume( getCoeffType(src) == n_long_R );
438
439  gmp_float *ff=(gmp_float*)from;
440  mpf_t *f=ff->_mpfp();
441  number res;
442  mpz_ptr dest,ndest;
443  int size, i,negative;
444  int e,al,bl;
445  mp_ptr qp,dd,nn;
446
447  size = (*f)[0]._mp_size;
448  if (size == 0)
449    return INT_TO_SR(0);
450  if(size<0)
451  {
452    negative = 1;
453    size = -size;
454  }
455  else
456    negative = 0;
457
458  qp = (*f)[0]._mp_d;
459  while(qp[0]==0)
460  {
461    qp++;
462    size--;
463  }
464
465  e=(*f)[0]._mp_exp-size;
466  res = ALLOC_RNUMBER();
467#if defined(LDEBUG)
468  res->debug=123456;
469#endif
470  dest = res->z;
471
472  void* (*allocfunc) (size_t);
473  mp_get_memory_functions (&allocfunc,NULL, NULL);
474  if (e<0)
475  {
476    al = dest->_mp_size = size;
477    if (al<2) al = 2;
478    dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
479    for (i=0;i<size;i++) dd[i] = qp[i];
480    bl = 1-e;
481    nn = (mp_ptr)allocfunc(sizeof(mp_limb_t)*bl);
482    memset(nn,0,sizeof(mp_limb_t)*bl);
483    nn[bl-1] = 1;
484    ndest = res->n;
485    ndest->_mp_d = nn;
486    ndest->_mp_alloc = ndest->_mp_size = bl;
487    res->s = 0;
488  }
489  else
490  {
491    al = dest->_mp_size = size+e;
492    if (al<2) al = 2;
493    dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
494    memset(dd,0,sizeof(mp_limb_t)*al);
495    for (i=0;i<size;i++) dd[i+e] = qp[i];
496    for (i=0;i<e;i++) dd[i] = 0;
497    res->s = 3;
498  }
499
500  dest->_mp_d = dd;
501  dest->_mp_alloc = al;
502  if (negative) mpz_neg(dest,dest);
503
504  if (res->s==0)
505    nlNormalize(res,dst);
506  else if (mpz_size1(res->z)<=MP_SMALL)
507  {
508    // res is new, res->ref is 1
509    res=nlShort3(res);
510  }
511  nlTest(res, dst);
512  return res;
513}
514
515static number nlMapLongR_BI(number from, const coeffs src, const coeffs dst)
516{
517  assume( getCoeffType(src) == n_long_R );
518
519  gmp_float *ff=(gmp_float*)from;
520  if (mpf_fits_slong_p(ff->t))
521  {
522    long l=mpf_get_si(ff->t);
523    return nlInit(l,dst);
524  }
525  char *out=floatToStr(*(gmp_float*)from, src->float_len);
526  char *p=strchr(out,'.');
527  *p='\0';
528  number res;
529  res = ALLOC_RNUMBER();
530#if defined(LDEBUG)
531  res->debug=123456;
532#endif
533  res->s=3;
534  mpz_init(res->z);
535  if (out[0]=='-')
536  {
537    mpz_set_str(res->z,out+1,10);
538    res=nlNeg(res,dst);
539  }
540  else
541  {
542    mpz_set_str(res->z,out,10);
543  }
544  omFree( (void *)out );
545  return res;
546}
547
548static number nlMapC(number from, const coeffs src, const coeffs dst)
549{
550  assume( getCoeffType(src) == n_long_C );
551  if ( ! ((gmp_complex*)from)->imag().isZero() )
552    return INT_TO_SR(0);
553
554  if (dst->is_field==FALSE) /* ->ZZ */
555  {
556    char *s=floatToStr(((gmp_complex*)from)->real(),src->float_len);
557    mpz_t z;
558    mpz_init(z);
559    char *ss=nEatLong(s,z);
560    if (*ss=='\0')
561    {
562      omFree(s);
563      number n=nlInitMPZ(z,dst);
564      mpz_clear(z);
565      return n;
566    }
567    omFree(s);
568    mpz_clear(z);
569    WarnS("conversion problem in CC -> ZZ mapping");
570    return INT_TO_SR(0);
571  }
572     
573  mpf_t *f = ((gmp_complex*)from)->real()._mpfp();
574
575  number res;
576  mpz_ptr dest,ndest;
577  int size, i,negative;
578  int e,al,bl;
579  mp_ptr qp,dd,nn;
580
581  size = (*f)[0]._mp_size;
582  if (size == 0)
583    return INT_TO_SR(0);
584  if(size<0)
585  {
586    negative = 1;
587    size = -size;
588  }
589  else
590    negative = 0;
591
592  qp = (*f)[0]._mp_d;
593  while(qp[0]==0)
594  {
595    qp++;
596    size--;
597  }
598
599  e=(*f)[0]._mp_exp-size;
600  res = ALLOC_RNUMBER();
601#if defined(LDEBUG)
602  res->debug=123456;
603#endif
604  dest = res->z;
605
606  void* (*allocfunc) (size_t);
607  mp_get_memory_functions (&allocfunc,NULL, NULL);
608  if (e<0)
609  {
610    al = dest->_mp_size = size;
611    if (al<2) al = 2;
612    dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
613    for (i=0;i<size;i++) dd[i] = qp[i];
614    bl = 1-e;
615    nn = (mp_ptr)allocfunc(sizeof(mp_limb_t)*bl);
616    memset(nn,0,sizeof(mp_limb_t)*bl);
617    nn[bl-1] = 1;
618    ndest = res->n;
619    ndest->_mp_d = nn;
620    ndest->_mp_alloc = ndest->_mp_size = bl;
621    res->s = 0;
622  }
623  else
624  {
625    al = dest->_mp_size = size+e;
626    if (al<2) al = 2;
627    dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
628    memset(dd,0,sizeof(mp_limb_t)*al);
629    for (i=0;i<size;i++) dd[i+e] = qp[i];
630    for (i=0;i<e;i++) dd[i] = 0;
631    res->s = 3;
632  }
633
634  dest->_mp_d = dd;
635  dest->_mp_alloc = al;
636  if (negative) mpz_neg(dest,dest);
637
638  if (res->s==0)
639    nlNormalize(res,dst);
640  else if (mpz_size1(res->z)<=MP_SMALL)
641  {
642    // res is new, res->ref is 1
643    res=nlShort3(res);
644  }
645  nlTest(res, dst);
646  return res;
647}
648
649//static number nlMapLongR(number from)
650//{
651//  gmp_float *ff=(gmp_float*)from;
652//  const mpf_t *f=ff->mpfp();
653//  int f_size=ABS((*f)[0]._mp_size);
654//  if (f_size==0)
655//    return nlInit(0);
656//  int f_sign=1;
657//  number work=ngcCopy(from);
658//  if (!ngcGreaterZero(work))
659//  {
660//    f_sign=-1;
661//    work=ngcNeg(work);
662//  }
663//  int i=0;
664//  mpz_t h1;
665//  mpz_init_set_ui(h1,1);
666//  while((FLT_RADIX*f) < DBL_MAX && i<DBL_MANT_DIG)
667//  {
668//    f*=FLT_RADIX;
669//    mpz_mul_ui(h1,h1,FLT_RADIX);
670//    i++;
671//  }
672//  number r=nlRInit(1);
673//  mpz_set_d(&(r->z),f);
674//  memcpy(&(r->n),&h1,sizeof(h1));
675//  r->s=0; /* not normalized */
676//  nlNormalize(r);
677//  return r;
678//
679//
680//  number r=nlRInit(1);
681//  int f_shift=f_size+(*f)[0]._mp_exp;
682//  if ( f_shift > 0)
683//  {
684//    r->s=0;
685//    mpz_init(&r->n);
686//    mpz_setbit(&r->n,f_shift*BYTES_PER_MP_LIMB*8);
687//    mpz_setbit(&r->z,f_size*BYTES_PER_MP_LIMB*8-1);
688//    // now r->z has enough space
689//    memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
690//    nlNormalize(r);
691//  }
692//  else
693//  {
694//    r->s=3;
695//    if (f_shift==0)
696//    {
697//      mpz_setbit(&r->z,f_size*BYTES_PER_MP_LIMB*8-1);
698//      // now r->z has enough space
699//      memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
700//    }
701//    else /* f_shift < 0 */
702//    {
703//      mpz_setbit(&r->z,(f_size-f_shift)*BYTES_PER_MP_LIMB*8-1);
704//      // now r->z has enough space
705//      memcpy(mpz_limb_d(&r->z)-f_shift,((*f)[0]._mp_d),
706//        f_size*BYTES_PER_MP_LIMB);
707//    }
708//  }
709//  if ((*f)[0]._mp_size<0);
710//    r=nlNeg(r);
711//  return r;
712//}
713
714int nlSize(number a, const coeffs)
715{
716  if (a==INT_TO_SR(0))
717     return 0; /* rational 0*/
718  if (SR_HDL(a) & SR_INT)
719     return 1; /* immediate int */
720  int s=a->z[0]._mp_alloc;
721//  while ((s>0) &&(a->z._mp_d[s]==0L)) s--;
722//#if SIZEOF_LONG == 8
723//  if (a->z._mp_d[s] < (unsigned long)0x100000000L) s=s*2-1;
724//  else s *=2;
725//#endif
726//  s++;
727  if (a->s<2)
728  {
729    int d=a->n[0]._mp_alloc;
730//    while ((d>0) && (a->n._mp_d[d]==0L)) d--;
731//#if SIZEOF_LONG == 8
732//    if (a->n._mp_d[d] < (unsigned long)0x100000000L) d=d*2-1;
733//    else d *=2;
734//#endif
735    s+=d;
736  }
737  return s;
738}
739
740/*2
741* convert number to int
742*/
743long nlInt(number &i, const coeffs r)
744{
745  nlTest(i, r);
746  nlNormalize(i,r);
747  if (SR_HDL(i) & SR_INT)
748  {
749    return SR_TO_INT(i);
750  }
751  if (i->s==3)
752  {
753    if(mpz_size1(i->z)>MP_SMALL) return 0;
754    long ul=mpz_get_si(i->z);
755    if (mpz_cmp_si(i->z,ul)!=0) return 0;
756    return ul;
757  }
758  mpz_t tmp;
759  long ul;
760  mpz_init(tmp);
761  mpz_tdiv_q(tmp,i->z,i->n);
762  if(mpz_size1(tmp)>MP_SMALL) ul=0;
763  else
764  {
765    ul=mpz_get_si(tmp);
766    if (mpz_cmp_si(tmp,ul)!=0) ul=0;
767  }
768  mpz_clear(tmp);
769  return ul;
770}
771
772/*2
773* convert number to bigint
774*/
775number nlBigInt(number &i, const coeffs r)
776{
777  nlTest(i, r);
778  nlNormalize(i,r);
779  if (SR_HDL(i) & SR_INT) return (i);
780  if (i->s==3)
781  {
782    return nlCopy(i,r);
783  }
784  number tmp=nlRInit(1);
785  mpz_tdiv_q(tmp->z,i->z,i->n);
786  tmp=nlShort3(tmp);
787  return tmp;
788}
789
790/*
791* 1/a
792*/
793number nlInvers(number a, const coeffs r)
794{
795  nlTest(a, r);
796  number n;
797  if (SR_HDL(a) & SR_INT)
798  {
799    if ((a==INT_TO_SR(1L)) || (a==INT_TO_SR(-1L)))
800    {
801      return a;
802    }
803    if (nlIsZero(a,r))
804    {
805      WerrorS(nDivBy0);
806      return INT_TO_SR(0);
807    }
808    n=ALLOC_RNUMBER();
809#if defined(LDEBUG)
810    n->debug=123456;
811#endif
812    n->s=1;
813    if (((long)a)>0L)
814    {
815      mpz_init_set_ui(n->z,1L);
816      mpz_init_set_si(n->n,(long)SR_TO_INT(a));
817    }
818    else
819    {
820      mpz_init_set_si(n->z,-1L);
821      mpz_init_set_si(n->n,(long)-SR_TO_INT(a));
822    }
823    nlTest(n, r);
824    return n;
825  }
826  n=ALLOC_RNUMBER();
827#if defined(LDEBUG)
828  n->debug=123456;
829#endif
830  {
831    mpz_init_set(n->n,a->z);
832    switch (a->s)
833    {
834      case 0:
835      case 1:
836              n->s=a->s;
837              mpz_init_set(n->z,a->n);
838              if (mpz_isNeg(n->n)) /* && n->s<2*/
839              {
840                mpz_neg(n->z,n->z);
841                mpz_neg(n->n,n->n);
842              }
843              if (mpz_cmp_ui(n->n,1L)==0)
844              {
845                mpz_clear(n->n);
846                n->s=3;
847                n=nlShort3(n);
848              }
849              break;
850      case 3:
851              // i.e. |a| > 2^...
852              n->s=1;
853              if (mpz_isNeg(n->n)) /* && n->s<2*/
854              {
855                mpz_neg(n->n,n->n);
856                mpz_init_set_si(n->z,-1L);
857              }
858              else
859              {
860                mpz_init_set_ui(n->z,1L);
861              }
862              break;
863    }
864  }
865  nlTest(n, r);
866  return n;
867}
868
869
870/*2
871* u := a / b in Z, if b | a (else undefined)
872*/
873number   nlExactDiv(number a, number b, const coeffs r)
874{
875  if (b==INT_TO_SR(0))
876  {
877    WerrorS(nDivBy0);
878    return INT_TO_SR(0);
879  }
880  number u;
881  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
882  {
883    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
884    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
885    {
886      return nlRInit(POW_2_28);
887    }
888    long aa=SR_TO_INT(a);
889    long bb=SR_TO_INT(b);
890    return INT_TO_SR(aa/bb);
891  }
892  number aa=NULL;
893  number bb=NULL;
894  if (SR_HDL(a) & SR_INT)
895  {
896    aa=nlRInit(SR_TO_INT(a));
897    a=aa;
898  }
899  if (SR_HDL(b) & SR_INT)
900  {
901    bb=nlRInit(SR_TO_INT(b));
902    b=bb;
903  }
904  u=ALLOC_RNUMBER();
905#if defined(LDEBUG)
906  u->debug=123456;
907#endif
908  mpz_init(u->z);
909  /* u=a/b */
910  u->s = 3;
911  assume(a->s==3);
912  assume(b->s==3);
913  mpz_divexact(u->z,a->z,b->z);
914  if (aa!=NULL)
915  {
916    mpz_clear(aa->z);
917#if defined(LDEBUG)
918    aa->debug=654324;
919#endif
920    FREE_RNUMBER(aa); // omFreeBin((void *)aa, rnumber_bin);
921  }
922  if (bb!=NULL)
923  {
924    mpz_clear(bb->z);
925#if defined(LDEBUG)
926    bb->debug=654324;
927#endif
928    FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
929  }
930  u=nlShort3(u);
931  nlTest(u, r);
932  return u;
933}
934
935/*2
936* u := a / b in Z
937*/
938number nlIntDiv (number a, number b, const coeffs r)
939{
940  if (b==INT_TO_SR(0))
941  {
942    WerrorS(nDivBy0);
943    return INT_TO_SR(0);
944  }
945  number u;
946  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
947  {
948    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
949    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
950    {
951      return nlRInit(POW_2_28);
952    }
953    LONG aa=SR_TO_INT(a);
954    LONG bb=SR_TO_INT(b);
955    LONG rr=aa%bb;
956    if (rr<0) rr+=ABS(bb);
957    LONG cc=(aa-rr)/bb;
958    return INT_TO_SR(cc);
959  }
960  number aa=NULL;
961  if (SR_HDL(a) & SR_INT)
962  {
963    /* the small int -(1<<28) divided by 2^28 is 1   */
964    if (a==INT_TO_SR(-(POW_2_28)))
965    {
966      if(mpz_cmp_si(b->z,(POW_2_28))==0)
967      {
968        return INT_TO_SR(-1);
969      }
970    }
971    aa=nlRInit(SR_TO_INT(a));
972    a=aa;
973  }
974  number bb=NULL;
975  if (SR_HDL(b) & SR_INT)
976  {
977    bb=nlRInit(SR_TO_INT(b));
978    b=bb;
979  }
980  u=ALLOC_RNUMBER();
981#if defined(LDEBUG)
982  u->debug=123456;
983#endif
984  assume(a->s==3);
985  assume(b->s==3);
986  /* u=u/b */
987  mpz_t rr;
988  mpz_init(rr);
989  mpz_mod(rr,a->z,b->z);
990  u->s = 3;
991  mpz_init(u->z);
992  mpz_sub(u->z,a->z,rr);
993  mpz_clear(rr);
994  mpz_divexact(u->z,u->z,b->z);
995  if (aa!=NULL)
996  {
997    mpz_clear(aa->z);
998#if defined(LDEBUG)
999    aa->debug=654324;
1000#endif
1001    FREE_RNUMBER(aa);
1002  }
1003  if (bb!=NULL)
1004  {
1005    mpz_clear(bb->z);
1006#if defined(LDEBUG)
1007    bb->debug=654324;
1008#endif
1009    FREE_RNUMBER(bb);
1010  }
1011  u=nlShort3(u);
1012  nlTest(u,r);
1013  return u;
1014}
1015
1016/*2
1017* u := a mod b in Z, u>=0
1018*/
1019number nlIntMod (number a, number b, const coeffs r)
1020{
1021  if (b==INT_TO_SR(0))
1022  {
1023    WerrorS(nDivBy0);
1024    return INT_TO_SR(0);
1025  }
1026  if (a==INT_TO_SR(0))
1027    return INT_TO_SR(0);
1028  number u;
1029  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1030  {
1031    LONG aa=SR_TO_INT(a);
1032    LONG bb=SR_TO_INT(b);
1033    LONG c=aa % bb;
1034    if (c<0) c+=ABS(bb);
1035    return INT_TO_SR(c);
1036  }
1037  if (SR_HDL(a) & SR_INT)
1038  {
1039    LONG ai=SR_TO_INT(a);
1040    mpz_t aa;
1041    mpz_init_set_si(aa, ai);
1042    u=ALLOC_RNUMBER();
1043#if defined(LDEBUG)
1044    u->debug=123456;
1045#endif
1046    u->s = 3;
1047    mpz_init(u->z);
1048    mpz_mod(u->z, aa, b->z);
1049    mpz_clear(aa);
1050    u=nlShort3(u);
1051    nlTest(u,r);
1052    return u;
1053  }
1054  number bb=NULL;
1055  if (SR_HDL(b) & SR_INT)
1056  {
1057    bb=nlRInit(SR_TO_INT(b));
1058    b=bb;
1059  }
1060  u=ALLOC_RNUMBER();
1061#if defined(LDEBUG)
1062  u->debug=123456;
1063#endif
1064  mpz_init(u->z);
1065  u->s = 3;
1066  mpz_mod(u->z, a->z, b->z);
1067  if (bb!=NULL)
1068  {
1069    mpz_clear(bb->z);
1070#if defined(LDEBUG)
1071    bb->debug=654324;
1072#endif
1073    FREE_RNUMBER(bb);
1074  }
1075  u=nlShort3(u);
1076  nlTest(u,r);
1077  return u;
1078}
1079
1080BOOLEAN nlDivBy (number a,number b, const coeffs)
1081{
1082  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1083  {
1084    return ((SR_TO_INT(a) % SR_TO_INT(b))==0);
1085  }
1086  if (SR_HDL(b) & SR_INT)
1087  {
1088    return (mpz_divisible_ui_p(a->z,SR_TO_INT(b))!=0);
1089  }
1090  if (SR_HDL(a) & SR_INT) return FALSE;
1091  return mpz_divisible_p(a->z, b->z) != 0;
1092}
1093
1094int nlDivComp(number a, number b, const coeffs r)
1095{
1096  if (nlDivBy(a, b, r))
1097  {
1098    if (nlDivBy(b, a, r)) return 2;
1099    return -1;
1100  }
1101  if (nlDivBy(b, a, r)) return 1;
1102  return 0;
1103}
1104
1105number  nlGetUnit (number n, const coeffs cf)
1106{
1107  if (nlGreaterZero(n,cf)) return INT_TO_SR(1);
1108  else                     return INT_TO_SR(-1);
1109}
1110
1111coeffs nlQuot1(number c, const coeffs r)
1112{
1113  long ch = r->cfInt(c, r);
1114  int p=IsPrime(ch);
1115  coeffs rr=NULL;
1116  if (((long)p)==ch)
1117  {
1118    rr = nInitChar(n_Zp,(void*)ch);
1119  }
1120  #ifdef HAVE_RINGS
1121  else
1122  {
1123    mpz_t dummy;
1124    mpz_init_set_ui(dummy, ch);
1125    ZnmInfo info;
1126    info.base = dummy;
1127    info.exp = (unsigned long) 1;
1128    rr = nInitChar(n_Zn, (void*)&info);
1129    mpz_clear(dummy);
1130  }
1131  #endif
1132  return(rr);
1133}
1134
1135
1136BOOLEAN nlIsUnit (number a, const coeffs)
1137{
1138  return ((SR_HDL(a) & SR_INT) && (ABS(SR_TO_INT(a))==1));
1139}
1140
1141
1142/*2
1143* u := a / b
1144*/
1145number nlDiv (number a, number b, const coeffs r)
1146{
1147  if (nlIsZero(b,r))
1148  {
1149    WerrorS(nDivBy0);
1150    return INT_TO_SR(0);
1151  }
1152  number u;
1153// ---------- short / short ------------------------------------
1154  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1155  {
1156    LONG i=SR_TO_INT(a);
1157    LONG j=SR_TO_INT(b);
1158    if (j==1L) return a;
1159    if ((i==-POW_2_28) && (j== -1L))
1160    {
1161      return nlRInit(POW_2_28);
1162    }
1163    LONG r=i%j;
1164    if (r==0)
1165    {
1166      return INT_TO_SR(i/j);
1167    }
1168    u=ALLOC_RNUMBER();
1169    u->s=0;
1170    #if defined(LDEBUG)
1171    u->debug=123456;
1172    #endif
1173    mpz_init_set_si(u->z,(long)i);
1174    mpz_init_set_si(u->n,(long)j);
1175  }
1176  else
1177  {
1178    u=ALLOC_RNUMBER();
1179    u->s=0;
1180    #if defined(LDEBUG)
1181    u->debug=123456;
1182    #endif
1183    mpz_init(u->z);
1184// ---------- short / long ------------------------------------
1185    if (SR_HDL(a) & SR_INT)
1186    {
1187      // short a / (z/n) -> (a*n)/z
1188      if (b->s<2)
1189      {
1190        mpz_mul_si(u->z,b->n,SR_TO_INT(a));
1191      }
1192      else
1193      // short a / long z -> a/z
1194      {
1195        mpz_set_si(u->z,SR_TO_INT(a));
1196      }
1197      if (mpz_cmp(u->z,b->z)==0)
1198      {
1199        mpz_clear(u->z);
1200        FREE_RNUMBER(u);
1201        return INT_TO_SR(1);
1202      }
1203      mpz_init_set(u->n,b->z);
1204    }
1205// ---------- long / short ------------------------------------
1206    else if (SR_HDL(b) & SR_INT)
1207    {
1208      mpz_set(u->z,a->z);
1209      // (z/n) / b -> z/(n*b)
1210      if (a->s<2)
1211      {
1212        mpz_init_set(u->n,a->n);
1213        if (((long)b)>0L)
1214          mpz_mul_ui(u->n,u->n,SR_TO_INT(b));
1215        else
1216        {
1217          mpz_mul_ui(u->n,u->n,-SR_TO_INT(b));
1218          mpz_neg(u->z,u->z);
1219        }
1220      }
1221      else
1222      // long z / short b -> z/b
1223      {
1224        //mpz_set(u->z,a->z);
1225        mpz_init_set_si(u->n,SR_TO_INT(b));
1226      }
1227    }
1228// ---------- long / long ------------------------------------
1229    else
1230    {
1231      mpz_set(u->z,a->z);
1232      mpz_init_set(u->n,b->z);
1233      if (a->s<2) mpz_mul(u->n,u->n,a->n);
1234      if (b->s<2) mpz_mul(u->z,u->z,b->n);
1235    }
1236  }
1237  if (mpz_isNeg(u->n))
1238  {
1239    mpz_neg(u->z,u->z);
1240    mpz_neg(u->n,u->n);
1241  }
1242  if (mpz_cmp_si(u->n,1L)==0)
1243  {
1244    mpz_clear(u->n);
1245    u->s=3;
1246    u=nlShort3(u);
1247  }
1248  nlTest(u, r);
1249  return u;
1250}
1251
1252/*2
1253* u:= x ^ exp
1254*/
1255void nlPower (number x,int exp,number * u, const coeffs r)
1256{
1257  *u = INT_TO_SR(0); // 0^e, e!=0
1258  if (exp==0)
1259    *u= INT_TO_SR(1);
1260  else if (!nlIsZero(x,r))
1261  {
1262    nlTest(x, r);
1263    number aa=NULL;
1264    if (SR_HDL(x) & SR_INT)
1265    {
1266      aa=nlRInit(SR_TO_INT(x));
1267      x=aa;
1268    }
1269    else if (x->s==0)
1270      nlNormalize(x,r);
1271    *u=ALLOC_RNUMBER();
1272#if defined(LDEBUG)
1273    (*u)->debug=123456;
1274#endif
1275    mpz_init((*u)->z);
1276    mpz_pow_ui((*u)->z,x->z,(unsigned long)exp);
1277    if (x->s<2)
1278    {
1279      if (mpz_cmp_si(x->n,1L)==0)
1280      {
1281        x->s=3;
1282        mpz_clear(x->n);
1283      }
1284      else
1285      {
1286        mpz_init((*u)->n);
1287        mpz_pow_ui((*u)->n,x->n,(unsigned long)exp);
1288      }
1289    }
1290    (*u)->s = x->s;
1291    if ((*u)->s==3) *u=nlShort3(*u);
1292    if (aa!=NULL)
1293    {
1294      mpz_clear(aa->z);
1295      FREE_RNUMBER(aa);
1296    }
1297  }
1298#ifdef LDEBUG
1299  if (exp<0) Print("nlPower: neg. exp. %d\n",exp);
1300  nlTest(*u, r);
1301#endif
1302}
1303
1304
1305/*2
1306* za >= 0 ?
1307*/
1308BOOLEAN nlGreaterZero (number a, const coeffs r)
1309{
1310  nlTest(a, r);
1311  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>1L /* represents number(0) */;
1312  return (!mpz_isNeg(a->z));
1313}
1314
1315/*2
1316* a > b ?
1317*/
1318BOOLEAN nlGreater (number a, number b, const coeffs r)
1319{
1320  nlTest(a, r);
1321  nlTest(b, r);
1322  number re;
1323  BOOLEAN rr;
1324  re=nlSub(a,b,r);
1325  rr=(!nlIsZero(re,r)) && (nlGreaterZero(re,r));
1326  nlDelete(&re,r);
1327  return rr;
1328}
1329
1330/*2
1331* a == -1 ?
1332*/
1333BOOLEAN nlIsMOne (number a, const coeffs r)
1334{
1335#ifdef LDEBUG
1336  if (a==NULL) return FALSE;
1337  nlTest(a, r);
1338#endif
1339  return  (a==INT_TO_SR(-1L));
1340}
1341
1342/*2
1343* result =gcd(a,b)
1344*/
1345number nlGcd(number a, number b, const coeffs r)
1346{
1347  number result;
1348  nlTest(a, r);
1349  nlTest(b, r);
1350  //nlNormalize(a);
1351  //nlNormalize(b);
1352  if ((a==INT_TO_SR(1L))||(a==INT_TO_SR(-1L))
1353  ||  (b==INT_TO_SR(1L))||(b==INT_TO_SR(-1L)))
1354    return INT_TO_SR(1L);
1355  if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
1356    return nlCopy(b,r);
1357  if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
1358    return nlCopy(a,r);
1359  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1360  {
1361    long i=SR_TO_INT(a);
1362    long j=SR_TO_INT(b);
1363    long l;
1364    i=ABS(i);
1365    j=ABS(j);
1366    do
1367    {
1368      l=i%j;
1369      i=j;
1370      j=l;
1371    } while (l!=0L);
1372    if (i==POW_2_28)
1373      result=nlRInit(POW_2_28);
1374    else
1375     result=INT_TO_SR(i);
1376    nlTest(result,r);
1377    return result;
1378  }
1379  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1380  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1381  if (SR_HDL(a) & SR_INT)
1382  {
1383    LONG aa=ABS(SR_TO_INT(a));
1384    unsigned long t=mpz_gcd_ui(NULL,b->z,(long)aa);
1385    if (t==POW_2_28)
1386      result=nlRInit(POW_2_28);
1387    else
1388     result=INT_TO_SR(t);
1389  }
1390  else
1391  if (SR_HDL(b) & SR_INT)
1392  {
1393    LONG bb=ABS(SR_TO_INT(b));
1394    unsigned long t=mpz_gcd_ui(NULL,a->z,(long)bb);
1395    if (t==POW_2_28)
1396      result=nlRInit(POW_2_28);
1397    else
1398     result=INT_TO_SR(t);
1399  }
1400  else
1401  {
1402    result=ALLOC0_RNUMBER();
1403    result->s = 3;
1404  #ifdef LDEBUG
1405    result->debug=123456;
1406  #endif
1407    mpz_init(result->z);
1408    mpz_gcd(result->z,a->z,b->z);
1409    result=nlShort3(result);
1410  }
1411  nlTest(result, r);
1412  return result;
1413}
1414
1415static int int_extgcd(int a, int b, int * u, int* x, int * v, int* y)
1416{
1417  int q, r;
1418  if (a==0)
1419  {
1420    *u = 0;
1421    *v = 1;
1422    *x = -1;
1423    *y = 0;
1424    return b;
1425  }
1426  if (b==0)
1427  {
1428    *u = 1;
1429    *v = 0;
1430    *x = 0;
1431    *y = 1;
1432    return a;
1433  }
1434  *u=1;
1435  *v=0;
1436  *x=0;
1437  *y=1;
1438  do
1439  {
1440    q = a/b;
1441    r = a%b;
1442    assume (q*b+r == a);
1443    a = b;
1444    b = r;
1445
1446    r = -(*v)*q+(*u);
1447    (*u) =(*v);
1448    (*v) = r;
1449
1450    r = -(*y)*q+(*x);
1451    (*x) = (*y);
1452    (*y) = r;
1453  } while (b);
1454
1455  return a;
1456}
1457
1458//number nlGcd_dummy(number a, number b, const coeffs r)
1459//{
1460//  extern char my_yylinebuf[80];
1461//  Print("nlGcd in >>%s<<\n",my_yylinebuf);
1462//  return nlGcd(a,b,r);;
1463//}
1464
1465number nlShort1(number x) // assume x->s==0/1
1466{
1467  assume(x->s<2);
1468  if (mpz_sgn1(x->z)==0)
1469  {
1470    _nlDelete_NoImm(&x);
1471    return INT_TO_SR(0);
1472  }
1473  if (x->s<2)
1474  {
1475    if (mpz_cmp(x->z,x->n)==0)
1476    {
1477      _nlDelete_NoImm(&x);
1478      return INT_TO_SR(1);
1479    }
1480  }
1481  return x;
1482}
1483/*2
1484* simplify x
1485*/
1486void nlNormalize (number &x, const coeffs r)
1487{
1488  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1489    return;
1490  if (x->s==3)
1491  {
1492    x=nlShort3_noinline(x);
1493    nlTest(x,r);
1494    return;
1495  }
1496  else if (x->s==0)
1497  {
1498    if (mpz_cmp_si(x->n,1L)==0)
1499    {
1500      mpz_clear(x->n);
1501      x->s=3;
1502      x=nlShort3(x);
1503    }
1504    else
1505    {
1506      mpz_t gcd;
1507      mpz_init(gcd);
1508      mpz_gcd(gcd,x->z,x->n);
1509      x->s=1;
1510      if (mpz_cmp_si(gcd,1L)!=0)
1511      {
1512        mpz_divexact(x->z,x->z,gcd);
1513        mpz_divexact(x->n,x->n,gcd);
1514        if (mpz_cmp_si(x->n,1L)==0)
1515        {
1516          mpz_clear(x->n);
1517          x->s=3;
1518          x=nlShort3_noinline(x);
1519        }
1520      }
1521      mpz_clear(gcd);
1522    }
1523  }
1524  nlTest(x, r);
1525}
1526
1527/*2
1528* returns in result->z the lcm(a->z,b->n)
1529*/
1530number nlNormalizeHelper(number a, number b, const coeffs r)
1531{
1532  number result;
1533  nlTest(a, r);
1534  nlTest(b, r);
1535  if ((SR_HDL(b) & SR_INT)
1536  || (b->s==3))
1537  {
1538    // b is 1/(b->n) => b->n is 1 => result is a
1539    return nlCopy(a,r);
1540  }
1541  result=ALLOC_RNUMBER();
1542#if defined(LDEBUG)
1543  result->debug=123456;
1544#endif
1545  result->s=3;
1546  mpz_t gcd;
1547  mpz_init(gcd);
1548  mpz_init(result->z);
1549  if (SR_HDL(a) & SR_INT)
1550    mpz_gcd_ui(gcd,b->n,ABS(SR_TO_INT(a)));
1551  else
1552    mpz_gcd(gcd,a->z,b->n);
1553  if (mpz_cmp_si(gcd,1L)!=0)
1554  {
1555    mpz_t bt;
1556    mpz_init(bt);
1557    mpz_divexact(bt,b->n,gcd);
1558    if (SR_HDL(a) & SR_INT)
1559      mpz_mul_si(result->z,bt,SR_TO_INT(a));
1560    else
1561      mpz_mul(result->z,bt,a->z);
1562    mpz_clear(bt);
1563  }
1564  else
1565    if (SR_HDL(a) & SR_INT)
1566      mpz_mul_si(result->z,b->n,SR_TO_INT(a));
1567    else
1568      mpz_mul(result->z,b->n,a->z);
1569  mpz_clear(gcd);
1570  result=nlShort3(result);
1571  nlTest(result, r);
1572  return result;
1573}
1574
1575// Map q \in QQ or ZZ \to Zp or an extension of it
1576// src = Q or Z, dst = Zp (or an extension of Zp)
1577number nlModP(number q, const coeffs /*Q*/, const coeffs Zp)
1578{
1579  const int p = n_GetChar(Zp);
1580  assume( p > 0 );
1581
1582  const long P = p;
1583  assume( P > 0 );
1584
1585  // embedded long within q => only long numerator has to be converted
1586  // to int (modulo char.)
1587  if (SR_HDL(q) & SR_INT)
1588  {
1589    long i = SR_TO_INT(q);
1590    return n_Init( i, Zp );
1591  }
1592
1593  const unsigned long PP = p;
1594
1595  // numerator modulo char. should fit into int
1596  number z = n_Init( static_cast<long>(mpz_fdiv_ui(q->z, PP)), Zp );
1597
1598  // denominator != 1?
1599  if (q->s!=3)
1600  {
1601    // denominator modulo char. should fit into int
1602    number n = n_Init( static_cast<long>(mpz_fdiv_ui(q->n, PP)), Zp );
1603
1604    number res = n_Div( z, n, Zp );
1605
1606    n_Delete(&z, Zp);
1607    n_Delete(&n, Zp);
1608
1609    return res;
1610  }
1611
1612  return z;
1613}
1614
1615#ifdef HAVE_RINGS
1616/*2
1617* convert number i (from Q) to GMP and warn if denom != 1
1618*/
1619void nlGMP(number &i, mpz_t n, const coeffs r)
1620{
1621  // Hier brauche ich einfach die GMP Zahl
1622  nlTest(i, r);
1623  nlNormalize(i, r);
1624  if (SR_HDL(i) & SR_INT)
1625  {
1626    mpz_set_si(n, SR_TO_INT(i));
1627    return;
1628  }
1629  if (i->s!=3)
1630  {
1631    WarnS("Omitted denominator during coefficient mapping !");
1632  }
1633  mpz_set(n, i->z);
1634}
1635#endif
1636
1637/*2
1638* acces to denominator, other 1 for integers
1639*/
1640number nlGetDenom(number &n, const coeffs r)
1641{
1642  if (!(SR_HDL(n) & SR_INT))
1643  {
1644    if (n->s==0)
1645    {
1646      nlNormalize(n,r);
1647    }
1648    if (!(SR_HDL(n) & SR_INT))
1649    {
1650      if (n->s!=3)
1651      {
1652        number u=ALLOC_RNUMBER();
1653        u->s=3;
1654#if defined(LDEBUG)
1655        u->debug=123456;
1656#endif
1657        mpz_init_set(u->z,n->n);
1658        u=nlShort3_noinline(u);
1659        return u;
1660      }
1661    }
1662  }
1663  return INT_TO_SR(1);
1664}
1665
1666/*2
1667* acces to Nominator, nlCopy(n) for integers
1668*/
1669number nlGetNumerator(number &n, const coeffs r)
1670{
1671  if (!(SR_HDL(n) & SR_INT))
1672  {
1673    if (n->s==0)
1674    {
1675      nlNormalize(n,r);
1676    }
1677    if (!(SR_HDL(n) & SR_INT))
1678    {
1679      number u=ALLOC_RNUMBER();
1680#if defined(LDEBUG)
1681      u->debug=123456;
1682#endif
1683      u->s=3;
1684      mpz_init_set(u->z,n->z);
1685      if (n->s!=3)
1686      {
1687        u=nlShort3_noinline(u);
1688      }
1689      return u;
1690    }
1691  }
1692  return n; // imm. int
1693}
1694
1695/***************************************************************
1696 *
1697 * routines which are needed by Inline(d) routines
1698 *
1699 *******************************************************************/
1700BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
1701{
1702  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1703//  long - short
1704  BOOLEAN bo;
1705  if (SR_HDL(b) & SR_INT)
1706  {
1707    if (a->s!=0) return FALSE;
1708    number n=b; b=a; a=n;
1709  }
1710//  short - long
1711  if (SR_HDL(a) & SR_INT)
1712  {
1713    if (b->s!=0)
1714      return FALSE;
1715    if ((((long)a) > 0L) && (mpz_isNeg(b->z)))
1716      return FALSE;
1717    if ((((long)a) < 0L) && (!mpz_isNeg(b->z)))
1718      return FALSE;
1719    mpz_t  bb;
1720    mpz_init(bb);
1721    mpz_mul_si(bb,b->n,(long)SR_TO_INT(a));
1722    bo=(mpz_cmp(bb,b->z)==0);
1723    mpz_clear(bb);
1724    return bo;
1725  }
1726// long - long
1727  if (((a->s==1) && (b->s==3))
1728  ||  ((b->s==1) && (a->s==3)))
1729    return FALSE;
1730  if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
1731    return FALSE;
1732  if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
1733    return FALSE;
1734  mpz_t  aa;
1735  mpz_t  bb;
1736  mpz_init_set(aa,a->z);
1737  mpz_init_set(bb,b->z);
1738  if (a->s<2) mpz_mul(bb,bb,a->n);
1739  if (b->s<2) mpz_mul(aa,aa,b->n);
1740  bo=(mpz_cmp(aa,bb)==0);
1741  mpz_clear(aa);
1742  mpz_clear(bb);
1743  return bo;
1744}
1745
1746// copy not immediate number a
1747number _nlCopy_NoImm(number a)
1748{
1749  assume(!(SR_HDL(a) & SR_INT));
1750  //nlTest(a, r);
1751  number b=ALLOC_RNUMBER();
1752#if defined(LDEBUG)
1753  b->debug=123456;
1754#endif
1755  switch (a->s)
1756  {
1757    case 0:
1758    case 1:
1759            mpz_init_set(b->n,a->n);
1760    case 3:
1761            mpz_init_set(b->z,a->z);
1762            break;
1763  }
1764  b->s = a->s;
1765  return b;
1766}
1767
1768void _nlDelete_NoImm(number *a)
1769{
1770  {
1771    switch ((*a)->s)
1772    {
1773      case 0:
1774      case 1:
1775        mpz_clear((*a)->n);
1776      case 3:
1777        mpz_clear((*a)->z);
1778    }
1779    #ifdef LDEBUG
1780    memset(*a,0,sizeof(**a));
1781    #endif
1782    FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
1783  }
1784}
1785
1786number _nlNeg_NoImm(number a)
1787{
1788  mpz_neg(a->z,a->z);
1789  if (a->s==3)
1790  {
1791    a=nlShort3(a);
1792  }
1793  return a;
1794}
1795
1796// conditio to use nlNormalize_Gcd in intermediate computations:
1797#define GCD_NORM_COND(OLD,NEW) (mpz_size1(NEW->z)>mpz_size1(OLD->z))
1798
1799static void nlNormalize_Gcd(number &x)
1800{
1801  mpz_t gcd;
1802  mpz_init(gcd);
1803  mpz_gcd(gcd,x->z,x->n);
1804  x->s=1;
1805  if (mpz_cmp_si(gcd,1L)!=0)
1806  {
1807    mpz_divexact(x->z,x->z,gcd);
1808    mpz_divexact(x->n,x->n,gcd);
1809    if (mpz_cmp_si(x->n,1L)==0)
1810    {
1811      mpz_clear(x->n);
1812      x->s=3;
1813      x=nlShort3_noinline(x);
1814    }
1815  }
1816  mpz_clear(gcd);
1817}
1818
1819number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
1820{
1821  number u=ALLOC_RNUMBER();
1822#if defined(LDEBUG)
1823  u->debug=123456;
1824#endif
1825  mpz_init(u->z);
1826  if (SR_HDL(b) & SR_INT)
1827  {
1828    number x=a;
1829    a=b;
1830    b=x;
1831  }
1832  if (SR_HDL(a) & SR_INT)
1833  {
1834    switch (b->s)
1835    {
1836      case 0:
1837      case 1:/* a:short, b:1 */
1838      {
1839        mpz_t x;
1840        mpz_init(x);
1841        mpz_mul_si(x,b->n,SR_TO_INT(a));
1842        mpz_add(u->z,b->z,x);
1843        mpz_clear(x);
1844        if (mpz_sgn1(u->z)==0)
1845        {
1846          mpz_clear(u->z);
1847          FREE_RNUMBER(u);
1848          return INT_TO_SR(0);
1849        }
1850        if (mpz_cmp(u->z,b->n)==0)
1851        {
1852          mpz_clear(u->z);
1853          FREE_RNUMBER(u);
1854          return INT_TO_SR(1);
1855        }
1856        mpz_init_set(u->n,b->n);
1857        u->s = 0;
1858        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1859        break;
1860      }
1861      case 3:
1862      {
1863        if (((long)a)>0L)
1864          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1865        else
1866          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1867        u->s = 3;
1868        u=nlShort3(u);
1869        break;
1870      }
1871    }
1872  }
1873  else
1874  {
1875    switch (a->s)
1876    {
1877      case 0:
1878      case 1:
1879      {
1880        switch(b->s)
1881        {
1882          case 0:
1883          case 1:
1884          {
1885            mpz_t x;
1886            mpz_init(x);
1887
1888            mpz_mul(x,b->z,a->n);
1889            mpz_mul(u->z,a->z,b->n);
1890            mpz_add(u->z,u->z,x);
1891            mpz_clear(x);
1892
1893            if (mpz_sgn1(u->z)==0)
1894            {
1895              mpz_clear(u->z);
1896              FREE_RNUMBER(u);
1897              return INT_TO_SR(0);
1898            }
1899            mpz_init(u->n);
1900            mpz_mul(u->n,a->n,b->n);
1901            if (mpz_cmp(u->z,u->n)==0)
1902            {
1903               mpz_clear(u->z);
1904               mpz_clear(u->n);
1905               FREE_RNUMBER(u);
1906               return INT_TO_SR(1);
1907            }
1908            u->s = 0;
1909            if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1910            break;
1911          }
1912          case 3: /* a:1 b:3 */
1913          {
1914            mpz_mul(u->z,b->z,a->n);
1915            mpz_add(u->z,u->z,a->z);
1916            if (mpz_sgn1(u->z)==0)
1917            {
1918              mpz_clear(u->z);
1919              FREE_RNUMBER(u);
1920              return INT_TO_SR(0);
1921            }
1922            if (mpz_cmp(u->z,a->n)==0)
1923            {
1924              mpz_clear(u->z);
1925              FREE_RNUMBER(u);
1926              return INT_TO_SR(1);
1927            }
1928            mpz_init_set(u->n,a->n);
1929            u->s = 0;
1930            if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
1931            break;
1932          }
1933        } /*switch (b->s) */
1934        break;
1935      }
1936      case 3:
1937      {
1938        switch(b->s)
1939        {
1940          case 0:
1941          case 1:/* a:3, b:1 */
1942          {
1943            mpz_mul(u->z,a->z,b->n);
1944            mpz_add(u->z,u->z,b->z);
1945            if (mpz_sgn1(u->z)==0)
1946            {
1947              mpz_clear(u->z);
1948              FREE_RNUMBER(u);
1949              return INT_TO_SR(0);
1950            }
1951            if (mpz_cmp(u->z,b->n)==0)
1952            {
1953              mpz_clear(u->z);
1954              FREE_RNUMBER(u);
1955              return INT_TO_SR(1);
1956            }
1957            mpz_init_set(u->n,b->n);
1958            u->s = 0;
1959            if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1960            break;
1961          }
1962          case 3:
1963          {
1964            mpz_add(u->z,a->z,b->z);
1965            u->s = 3;
1966            u=nlShort3(u);
1967            break;
1968          }
1969        }
1970        break;
1971      }
1972    }
1973  }
1974  return u;
1975}
1976
1977void _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b)
1978{
1979  if (SR_HDL(b) & SR_INT)
1980  {
1981    switch (a->s)
1982    {
1983      case 0:
1984      case 1:/* b:short, a:1 */
1985      {
1986        mpz_t x;
1987        mpz_init(x);
1988        mpz_mul_si(x,a->n,SR_TO_INT(b));
1989        mpz_add(a->z,a->z,x);
1990        mpz_clear(x);
1991        nlNormalize_Gcd(a);
1992        break;
1993      }
1994      case 3:
1995      {
1996        if (((long)b)>0L)
1997          mpz_add_ui(a->z,a->z,SR_TO_INT(b));
1998        else
1999          mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
2000        a->s = 3;
2001        a=nlShort3_noinline(a);
2002        break;
2003      }
2004    }
2005    return;
2006  }
2007  else if (SR_HDL(a) & SR_INT)
2008  {
2009    number u=ALLOC_RNUMBER();
2010    #if defined(LDEBUG)
2011    u->debug=123456;
2012    #endif
2013    mpz_init(u->z);
2014    switch (b->s)
2015    {
2016      case 0:
2017      case 1:/* a:short, b:1 */
2018      {
2019        mpz_t x;
2020        mpz_init(x);
2021
2022        mpz_mul_si(x,b->n,SR_TO_INT(a));
2023        mpz_add(u->z,b->z,x);
2024        mpz_clear(x);
2025        // result cannot be 0, if coeffs are normalized
2026        mpz_init_set(u->n,b->n);
2027        u->s=0;
2028        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2029        else { u=nlShort1(u); }
2030        break;
2031      }
2032      case 3:
2033      {
2034        if (((long)a)>0L)
2035          mpz_add_ui(u->z,b->z,SR_TO_INT(a));
2036        else
2037          mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
2038        // result cannot be 0, if coeffs are normalized
2039        u->s = 3;
2040        u=nlShort3_noinline(u);
2041        break;
2042      }
2043    }
2044    a=u;
2045  }
2046  else
2047  {
2048    switch (a->s)
2049    {
2050      case 0:
2051      case 1:
2052      {
2053        switch(b->s)
2054        {
2055          case 0:
2056          case 1: /* a:1 b:1 */
2057          {
2058            mpz_t x;
2059            mpz_t y;
2060            mpz_init(x);
2061            mpz_init(y);
2062            mpz_mul(x,b->z,a->n);
2063            mpz_mul(y,a->z,b->n);
2064            mpz_add(a->z,x,y);
2065            mpz_clear(x);
2066            mpz_clear(y);
2067            mpz_mul(a->n,a->n,b->n);
2068            a->s=0;
2069            if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
2070            else { a=nlShort1(a);}
2071            break;
2072          }
2073          case 3: /* a:1 b:3 */
2074          {
2075            mpz_t x;
2076            mpz_init(x);
2077            mpz_mul(x,b->z,a->n);
2078            mpz_add(a->z,a->z,x);
2079            mpz_clear(x);
2080            a->s=0;
2081            if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
2082            else { a=nlShort1(a);}
2083            break;
2084          }
2085        } /*switch (b->s) */
2086        break;
2087      }
2088      case 3:
2089      {
2090        switch(b->s)
2091        {
2092          case 0:
2093          case 1:/* a:3, b:1 */
2094          {
2095            mpz_t x;
2096            mpz_init(x);
2097            mpz_mul(x,a->z,b->n);
2098            mpz_add(a->z,b->z,x);
2099            mpz_clear(x);
2100            mpz_init_set(a->n,b->n);
2101            a->s=0;
2102            if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
2103            else { a=nlShort1(a);}
2104            break;
2105          }
2106          case 3:
2107          {
2108            mpz_add(a->z,a->z,b->z);
2109            a->s = 3;
2110            a=nlShort3_noinline(a);
2111            break;
2112          }
2113        }
2114        break;
2115      }
2116    }
2117  }
2118}
2119
2120number _nlSub_aNoImm_OR_bNoImm(number a, number b)
2121{
2122  number u=ALLOC_RNUMBER();
2123#if defined(LDEBUG)
2124  u->debug=123456;
2125#endif
2126  mpz_init(u->z);
2127  if (SR_HDL(a) & SR_INT)
2128  {
2129    switch (b->s)
2130    {
2131      case 0:
2132      case 1:/* a:short, b:1 */
2133      {
2134        mpz_t x;
2135        mpz_init(x);
2136        mpz_mul_si(x,b->n,SR_TO_INT(a));
2137        mpz_sub(u->z,x,b->z);
2138        mpz_clear(x);
2139        if (mpz_sgn1(u->z)==0)
2140        {
2141          mpz_clear(u->z);
2142          FREE_RNUMBER(u);
2143          return INT_TO_SR(0);
2144        }
2145        if (mpz_cmp(u->z,b->n)==0)
2146        {
2147          mpz_clear(u->z);
2148          FREE_RNUMBER(u);
2149          return INT_TO_SR(1);
2150        }
2151        mpz_init_set(u->n,b->n);
2152        u->s=0;
2153        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2154        break;
2155      }
2156      case 3:
2157      {
2158        if (((long)a)>0L)
2159        {
2160          mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
2161          mpz_neg(u->z,u->z);
2162        }
2163        else
2164        {
2165          mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
2166          mpz_neg(u->z,u->z);
2167        }
2168        u->s = 3;
2169        u=nlShort3(u);
2170        break;
2171      }
2172    }
2173  }
2174  else if (SR_HDL(b) & SR_INT)
2175  {
2176    switch (a->s)
2177    {
2178      case 0:
2179      case 1:/* b:short, a:1 */
2180      {
2181        mpz_t x;
2182        mpz_init(x);
2183        mpz_mul_si(x,a->n,SR_TO_INT(b));
2184        mpz_sub(u->z,a->z,x);
2185        mpz_clear(x);
2186        if (mpz_sgn1(u->z)==0)
2187        {
2188          mpz_clear(u->z);
2189          FREE_RNUMBER(u);
2190          return INT_TO_SR(0);
2191        }
2192        if (mpz_cmp(u->z,a->n)==0)
2193        {
2194          mpz_clear(u->z);
2195          FREE_RNUMBER(u);
2196          return INT_TO_SR(1);
2197        }
2198        mpz_init_set(u->n,a->n);
2199        u->s=0;
2200        if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2201        break;
2202      }
2203      case 3:
2204      {
2205        if (((long)b)>0L)
2206        {
2207          mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
2208        }
2209        else
2210        {
2211          mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
2212        }
2213        u->s = 3;
2214        u=nlShort3(u);
2215        break;
2216      }
2217    }
2218  }
2219  else
2220  {
2221    switch (a->s)
2222    {
2223      case 0:
2224      case 1:
2225      {
2226        switch(b->s)
2227        {
2228          case 0:
2229          case 1:
2230          {
2231            mpz_t x;
2232            mpz_t y;
2233            mpz_init(x);
2234            mpz_init(y);
2235            mpz_mul(x,b->z,a->n);
2236            mpz_mul(y,a->z,b->n);
2237            mpz_sub(u->z,y,x);
2238            mpz_clear(x);
2239            mpz_clear(y);
2240            if (mpz_sgn1(u->z)==0)
2241            {
2242              mpz_clear(u->z);
2243              FREE_RNUMBER(u);
2244              return INT_TO_SR(0);
2245            }
2246            mpz_init(u->n);
2247            mpz_mul(u->n,a->n,b->n);
2248            if (mpz_cmp(u->z,u->n)==0)
2249            {
2250              mpz_clear(u->z);
2251              mpz_clear(u->n);
2252              FREE_RNUMBER(u);
2253              return INT_TO_SR(1);
2254            }
2255            u->s=0;
2256            if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2257            break;
2258          }
2259          case 3: /* a:1, b:3 */
2260          {
2261            mpz_t x;
2262            mpz_init(x);
2263            mpz_mul(x,b->z,a->n);
2264            mpz_sub(u->z,a->z,x);
2265            mpz_clear(x);
2266            if (mpz_sgn1(u->z)==0)
2267            {
2268              mpz_clear(u->z);
2269              FREE_RNUMBER(u);
2270              return INT_TO_SR(0);
2271            }
2272            if (mpz_cmp(u->z,a->n)==0)
2273            {
2274              mpz_clear(u->z);
2275              FREE_RNUMBER(u);
2276              return INT_TO_SR(1);
2277            }
2278            mpz_init_set(u->n,a->n);
2279            u->s=0;
2280            if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2281            break;
2282          }
2283        }
2284        break;
2285      }
2286      case 3:
2287      {
2288        switch(b->s)
2289        {
2290          case 0:
2291          case 1: /* a:3, b:1 */
2292          {
2293            mpz_t x;
2294            mpz_init(x);
2295            mpz_mul(x,a->z,b->n);
2296            mpz_sub(u->z,x,b->z);
2297            mpz_clear(x);
2298            if (mpz_sgn1(u->z)==0)
2299            {
2300              mpz_clear(u->z);
2301              FREE_RNUMBER(u);
2302              return INT_TO_SR(0);
2303            }
2304            if (mpz_cmp(u->z,b->n)==0)
2305            {
2306              mpz_clear(u->z);
2307              FREE_RNUMBER(u);
2308              return INT_TO_SR(1);
2309            }
2310            mpz_init_set(u->n,b->n);
2311            u->s=0;
2312            if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2313            break;
2314          }
2315          case 3: /* a:3 , b:3 */
2316          {
2317            mpz_sub(u->z,a->z,b->z);
2318            u->s = 3;
2319            u=nlShort3(u);
2320            break;
2321          }
2322        }
2323        break;
2324      }
2325    }
2326  }
2327  return u;
2328}
2329
2330// a and b are intermediate, but a*b not
2331number _nlMult_aImm_bImm_rNoImm(number a, number b)
2332{
2333  number u=ALLOC_RNUMBER();
2334#if defined(LDEBUG)
2335  u->debug=123456;
2336#endif
2337  u->s=3;
2338  mpz_init_set_si(u->z,SR_TO_INT(a));
2339  mpz_mul_si(u->z,u->z,SR_TO_INT(b));
2340  return u;
2341}
2342
2343// a or b are not immediate
2344number _nlMult_aNoImm_OR_bNoImm(number a, number b)
2345{
2346  assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2347  number u=ALLOC_RNUMBER();
2348#if defined(LDEBUG)
2349  u->debug=123456;
2350#endif
2351  mpz_init(u->z);
2352  if (SR_HDL(b) & SR_INT)
2353  {
2354    number x=a;
2355    a=b;
2356    b=x;
2357  }
2358  if (SR_HDL(a) & SR_INT)
2359  {
2360    u->s=b->s;
2361    if (u->s==1) u->s=0;
2362    if (((long)a)>0L)
2363    {
2364      mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2365    }
2366    else
2367    {
2368      if (a==INT_TO_SR(-1))
2369      {
2370        mpz_set(u->z,b->z);
2371        mpz_neg(u->z,u->z);
2372        u->s=b->s;
2373      }
2374      else
2375      {
2376        mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2377        mpz_neg(u->z,u->z);
2378      }
2379    }
2380    if (u->s<2)
2381    {
2382      if (mpz_cmp(u->z,b->n)==0)
2383      {
2384        mpz_clear(u->z);
2385        FREE_RNUMBER(u);
2386        return INT_TO_SR(1);
2387      }
2388      mpz_init_set(u->n,b->n);
2389      if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2390    }
2391    else //u->s==3
2392    {
2393      u=nlShort3(u);
2394    }
2395  }
2396  else
2397  {
2398    mpz_mul(u->z,a->z,b->z);
2399    u->s = 0;
2400    if(a->s==3)
2401    {
2402      if(b->s==3)
2403      {
2404        u->s = 3;
2405      }
2406      else
2407      {
2408        if (mpz_cmp(u->z,b->n)==0)
2409        {
2410          mpz_clear(u->z);
2411          FREE_RNUMBER(u);
2412          return INT_TO_SR(1);
2413        }
2414        mpz_init_set(u->n,b->n);
2415        if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2416      }
2417    }
2418    else
2419    {
2420      if(b->s==3)
2421      {
2422        if (mpz_cmp(u->z,a->n)==0)
2423        {
2424          mpz_clear(u->z);
2425          FREE_RNUMBER(u);
2426          return INT_TO_SR(1);
2427        }
2428        mpz_init_set(u->n,a->n);
2429        if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2430      }
2431      else
2432      {
2433        mpz_init(u->n);
2434        mpz_mul(u->n,a->n,b->n);
2435        if (mpz_cmp(u->z,u->n)==0)
2436        {
2437          mpz_clear(u->z);
2438          mpz_clear(u->n);
2439          FREE_RNUMBER(u);
2440          return INT_TO_SR(1);
2441        }
2442        if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2443      }
2444    }
2445  }
2446  return u;
2447}
2448
2449/*2
2450* copy a to b for mapping
2451*/
2452number nlCopyMap(number a, const coeffs /*src*/, const coeffs /*dst*/)
2453{
2454  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2455  {
2456    return a;
2457  }
2458  return _nlCopy_NoImm(a);
2459}
2460
2461number nlMapQtoZ(number a, const coeffs src, const coeffs dst)
2462{
2463  if ((SR_HDL(a) & SR_INT)||(a==NULL))
2464  {
2465    return a;
2466  }
2467  if (a->s==3) return _nlCopy_NoImm(a);
2468  number a0=a;
2469  BOOLEAN a1=FALSE;
2470  if (a->s==0) { a0=_nlCopy_NoImm(a); a1=TRUE; }
2471  number b1=nlGetNumerator(a0,src);
2472  number b2=nlGetDenom(a0,src);
2473  number b=nlIntDiv(b1,b2,dst);
2474  nlDelete(&b1,src);
2475  nlDelete(&b2,src);
2476  if (a1)  _nlDelete_NoImm(&a0);
2477  return b;
2478}
2479
2480nMapFunc nlSetMap(const coeffs src, const coeffs dst)
2481{
2482  if (src->rep==n_rep_gap_rat)  /*Q, coeffs_BIGINT */
2483  {
2484    if ((src->is_field==dst->is_field) /* Q->Q, Z->Z*/
2485    || (src->is_field==FALSE))         /* Z->Q */
2486      return nlCopyMap;
2487    return nlMapQtoZ;        /* Q->Z */
2488  }
2489  if ((src->rep==n_rep_int) && nCoeff_is_Zp(src))
2490  {
2491    return nlMapP;
2492  }
2493  if ((src->rep==n_rep_float) && nCoeff_is_R(src))
2494  {
2495    if (dst->is_field) /* R -> Q */
2496      return nlMapR;
2497    else
2498      return nlMapR_BI; /* R -> bigint */
2499  }
2500  if ((src->rep==n_rep_gmp_float) && nCoeff_is_long_R(src))
2501  {
2502    if (dst->is_field)
2503      return nlMapLongR; /* long R -> Q */
2504    else
2505      return nlMapLongR_BI;
2506  }
2507  if (nCoeff_is_long_C(src))
2508  {
2509    return nlMapC; /* C -> Q */
2510  }
2511#ifdef HAVE_RINGS
2512  if (src->rep==n_rep_gmp) // nCoeff_is_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Zn(src))
2513  {
2514    return nlMapGMP;
2515  }
2516  if (src->rep==n_rep_gap_gmp)
2517  {
2518    return nlMapZ;
2519  }
2520  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src))
2521  {
2522    return nlMapMachineInt;
2523  }
2524#endif
2525  return NULL;
2526}
2527/*2
2528* z := i
2529*/
2530number nlRInit (long i)
2531{
2532  number z=ALLOC_RNUMBER();
2533#if defined(LDEBUG)
2534  z->debug=123456;
2535#endif
2536  mpz_init_set_si(z->z,i);
2537  z->s = 3;
2538  return z;
2539}
2540
2541/*2
2542* z := i/j
2543*/
2544number nlInit2 (int i, int j, const coeffs r)
2545{
2546  number z=ALLOC_RNUMBER();
2547#if defined(LDEBUG)
2548  z->debug=123456;
2549#endif
2550  mpz_init_set_si(z->z,(long)i);
2551  mpz_init_set_si(z->n,(long)j);
2552  z->s = 0;
2553  nlNormalize(z,r);
2554  return z;
2555}
2556
2557number nlInit2gmp (mpz_t i, mpz_t j, const coeffs r)
2558{
2559  number z=ALLOC_RNUMBER();
2560#if defined(LDEBUG)
2561  z->debug=123456;
2562#endif
2563  mpz_init_set(z->z,i);
2564  mpz_init_set(z->n,j);
2565  z->s = 0;
2566  nlNormalize(z,r);
2567  return z;
2568}
2569
2570#else // DO_LINLINE
2571
2572// declare immedate routines
2573number nlRInit (long i);
2574BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b);
2575number  _nlCopy_NoImm(number a);
2576number  _nlNeg_NoImm(number a);
2577number  _nlAdd_aNoImm_OR_bNoImm(number a, number b);
2578void    _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b);
2579number  _nlSub_aNoImm_OR_bNoImm(number a, number b);
2580number  _nlMult_aNoImm_OR_bNoImm(number a, number b);
2581number  _nlMult_aImm_bImm_rNoImm(number a, number b);
2582
2583#endif
2584
2585/***************************************************************
2586 *
2587 * Routines which might be inlined by p_Numbers.h
2588 *
2589 *******************************************************************/
2590#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2591
2592// routines which are always inlined/static
2593
2594/*2
2595* a = b ?
2596*/
2597LINLINE BOOLEAN nlEqual (number a, number b, const coeffs r)
2598{
2599  nlTest(a, r);
2600  nlTest(b, r);
2601// short - short
2602  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2603  return _nlEqual_aNoImm_OR_bNoImm(a, b);
2604}
2605
2606LINLINE number nlInit (long i, const coeffs r)
2607{
2608  number n;
2609  #if MAX_NUM_SIZE == 60
2610  if (((i << 3) >> 3) == i) n=INT_TO_SR(i);
2611  else                      n=nlRInit(i);
2612  #else
2613  LONG ii=(LONG)i;
2614  if ( ((((long)ii)==i) && ((ii << 3) >> 3) == ii )) n=INT_TO_SR(ii);
2615  else                                             n=nlRInit(i);
2616  #endif
2617  nlTest(n, r);
2618  return n;
2619}
2620
2621/*2
2622* a == 1 ?
2623*/
2624LINLINE BOOLEAN nlIsOne (number a, const coeffs r)
2625{
2626#ifdef LDEBUG
2627  if (a==NULL) return FALSE;
2628  nlTest(a, r);
2629#endif
2630  return (a==INT_TO_SR(1));
2631}
2632
2633LINLINE BOOLEAN nlIsZero (number a, const coeffs)
2634{
2635  #if 0
2636  if (a==INT_TO_SR(0)) return TRUE;
2637  if ((SR_HDL(a) & SR_INT)||(a==NULL)) return FALSE;
2638  if (mpz_cmp_si(a->z,0L)==0)
2639  {
2640    printf("gmp-0 in nlIsZero\n");
2641    dErrorBreak();
2642    return TRUE;
2643  }
2644  return FALSE;
2645  #else
2646  return (a==NULL)|| (a==INT_TO_SR(0));
2647  #endif
2648}
2649
2650/*2
2651* copy a to b
2652*/
2653LINLINE number nlCopy(number a, const coeffs)
2654{
2655  if (SR_HDL(a) & SR_INT)
2656  {
2657    return a;
2658  }
2659  return _nlCopy_NoImm(a);
2660}
2661
2662
2663/*2
2664* delete a
2665*/
2666LINLINE void nlDelete (number * a, const coeffs r)
2667{
2668  if (*a!=NULL)
2669  {
2670    nlTest(*a, r);
2671    if ((SR_HDL(*a) & SR_INT)==0)
2672    {
2673      _nlDelete_NoImm(a);
2674    }
2675    *a=NULL;
2676  }
2677}
2678
2679/*2
2680* za:= - za
2681*/
2682LINLINE number nlNeg (number a, const coeffs R)
2683{
2684  nlTest(a, R);
2685  if(SR_HDL(a) &SR_INT)
2686  {
2687    LONG r=SR_TO_INT(a);
2688    if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2689    else               a=INT_TO_SR(-r);
2690    return a;
2691  }
2692  a = _nlNeg_NoImm(a);
2693  nlTest(a, R);
2694  return a;
2695
2696}
2697
2698/*2
2699* u:= a + b
2700*/
2701LINLINE number nlAdd (number a, number b, const coeffs R)
2702{
2703  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2704  {
2705    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2706    if ( ((r << 1) >> 1) == r )
2707      return (number)(long)r;
2708    else
2709      return nlRInit(SR_TO_INT(r));
2710  }
2711  number u =  _nlAdd_aNoImm_OR_bNoImm(a, b);
2712  nlTest(u, R);
2713  return u;
2714}
2715
2716number nlShort1(number a);
2717number nlShort3_noinline(number x);
2718
2719LINLINE void nlInpAdd(number &a, number b, const coeffs r)
2720{
2721  // a=a+b
2722  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2723  {
2724    LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2725    if ( ((r << 1) >> 1) == r )
2726      a=(number)(long)r;
2727    else
2728      a=nlRInit(SR_TO_INT(r));
2729  }
2730  else
2731  {
2732    _nlInpAdd_aNoImm_OR_bNoImm(a,b);
2733    nlTest(a,r);
2734  }
2735}
2736
2737LINLINE number nlMult (number a, number b, const coeffs R)
2738{
2739  nlTest(a, R);
2740  nlTest(b, R);
2741  if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2742  if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2743  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2744  {
2745    LONG r=(LONG)((unsigned LONG)(SR_HDL(a)-1L))*((unsigned LONG)(SR_HDL(b)>>1));
2746    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2747    {
2748      number u=((number) ((r>>1)+SR_INT));
2749      if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2750      return nlRInit(SR_HDL(u)>>2);
2751    }
2752    number u = _nlMult_aImm_bImm_rNoImm(a, b);
2753    nlTest(u, R);
2754    return u;
2755
2756  }
2757  number u = _nlMult_aNoImm_OR_bNoImm(a, b);
2758  nlTest(u, R);
2759  return u;
2760
2761}
2762
2763
2764/*2
2765* u:= a - b
2766*/
2767LINLINE number nlSub (number a, number b, const coeffs r)
2768{
2769  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2770  {
2771    LONG r=SR_HDL(a)-SR_HDL(b)+1;
2772    if ( ((r << 1) >> 1) == r )
2773    {
2774      return (number)(long)r;
2775    }
2776    else
2777      return nlRInit(SR_TO_INT(r));
2778  }
2779  number u = _nlSub_aNoImm_OR_bNoImm(a, b);
2780  nlTest(u, r);
2781  return u;
2782
2783}
2784
2785LINLINE void nlInpMult(number &a, number b, const coeffs r)
2786{
2787  number aa=a;
2788  if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2789  {
2790    number n=nlMult(aa,b,r);
2791    nlDelete(&a,r);
2792    a=n;
2793  }
2794  else
2795  {
2796    mpz_mul(aa->z,a->z,b->z);
2797    if (aa->s==3)
2798    {
2799      if(b->s!=3)
2800      {
2801        mpz_init_set(a->n,b->n);
2802        a->s=0;
2803      }
2804    }
2805    else
2806    {
2807      if(b->s!=3)
2808      {
2809        mpz_mul(a->n,a->n,b->n);
2810      }
2811      a->s=0;
2812    }
2813  }
2814}
2815#endif // DO_LINLINE
2816
2817#ifndef P_NUMBERS_H
2818
2819void nlMPZ(mpz_t m, number &n, const coeffs r)
2820{
2821  nlTest(n, r);
2822  nlNormalize(n, r);
2823  if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n));     /* n fits in an int */
2824  else             mpz_init_set(m, (mpz_ptr)n->z);
2825}
2826
2827
2828number  nlXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
2829{
2830  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2831  {
2832    int uu, vv, x, y;
2833    int g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &uu, &vv, &x, &y);
2834    *s = INT_TO_SR(uu);
2835    *t = INT_TO_SR(vv);
2836    *u = INT_TO_SR(x);
2837    *v = INT_TO_SR(y);
2838    return INT_TO_SR(g);
2839  }
2840  else
2841  {
2842    mpz_t aa, bb;
2843    if (SR_HDL(a) & SR_INT)
2844    {
2845      mpz_init_set_si(aa, SR_TO_INT(a));
2846    }
2847    else
2848    {
2849      mpz_init_set(aa, a->z);
2850    }
2851    if (SR_HDL(b) & SR_INT)
2852    {
2853      mpz_init_set_si(bb, SR_TO_INT(b));
2854    }
2855    else
2856    {
2857      mpz_init_set(bb, b->z);
2858    }
2859    mpz_t erg; mpz_t bs; mpz_t bt;
2860    mpz_init(erg);
2861    mpz_init(bs);
2862    mpz_init(bt);
2863
2864    mpz_gcdext(erg, bs, bt, aa, bb);
2865
2866    mpz_div(aa, aa, erg);
2867    *u=nlInitMPZ(bb,r);
2868    *u=nlNeg(*u,r);
2869    *v=nlInitMPZ(aa,r);
2870
2871    mpz_clear(aa);
2872    mpz_clear(bb);
2873
2874    *s = nlInitMPZ(bs,r);
2875    *t = nlInitMPZ(bt,r);
2876    return nlInitMPZ(erg,r);
2877  }
2878}
2879
2880number nlQuotRem (number a, number b, number * r, const coeffs R)
2881{
2882  assume(SR_TO_INT(b)!=0);
2883  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2884  {
2885    if (r!=NULL)
2886      *r = INT_TO_SR(SR_TO_INT(a) % SR_TO_INT(b));
2887    return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
2888  }
2889  else if (SR_HDL(a) & SR_INT)
2890  {
2891    // -2^xx / 2^xx
2892    if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
2893    {
2894      if (r!=NULL) *r=INT_TO_SR(0);
2895      return nlRInit(POW_2_28);
2896    }
2897    //a is small, b is not, so q=0, r=a
2898    if (r!=NULL)
2899      *r = a;
2900    return INT_TO_SR(0);
2901  }
2902  else if (SR_HDL(b) & SR_INT)
2903  {
2904    unsigned long rr;
2905    mpz_t qq;
2906    mpz_init(qq);
2907    mpz_t rrr;
2908    mpz_init(rrr);
2909    rr = mpz_divmod_ui(qq, rrr, a->z, (unsigned long)ABS(SR_TO_INT(b)));
2910    mpz_clear(rrr);
2911
2912    if (r!=NULL)
2913      *r = INT_TO_SR(rr);
2914    if (SR_TO_INT(b)<0)
2915    {
2916      mpz_neg(qq, qq);
2917    }
2918    return nlInitMPZ(qq,R);
2919  }
2920  mpz_t qq,rr;
2921  mpz_init(qq);
2922  mpz_init(rr);
2923  mpz_divmod(qq, rr, a->z, b->z);
2924  if (r!=NULL)
2925    *r = nlInitMPZ(rr,R);
2926  else
2927  {
2928    mpz_clear(rr);
2929  }
2930  return nlInitMPZ(qq,R);
2931}
2932
2933void nlInpGcd(number &a, number b, const coeffs r)
2934{
2935  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2936  {
2937    number n=nlGcd(a,b,r);
2938    nlDelete(&a,r);
2939    a=n;
2940  }
2941  else
2942  {
2943    mpz_gcd(a->z,a->z,b->z);
2944    a=nlShort3_noinline(a);
2945  }
2946}
2947
2948void nlInpIntDiv(number &a, number b, const coeffs r)
2949{
2950  if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2951  {
2952    number n=nlIntDiv(a,b, r);
2953    nlDelete(&a,r);
2954    a=n;
2955  }
2956  else
2957  {
2958    mpz_t rr;
2959    mpz_init(rr);
2960    mpz_mod(rr,a->z,b->z);
2961    mpz_sub(a->z,a->z,rr);
2962    mpz_clear(rr);
2963    mpz_divexact(a->z,a->z,b->z);
2964    a=nlShort3_noinline(a);
2965  }
2966}
2967
2968number nlFarey(number nN, number nP, const coeffs r)
2969{
2970  mpz_t A,B,C,D,E,N,P,tmp;
2971  if (SR_HDL(nP) & SR_INT) mpz_init_set_si(P,SR_TO_INT(nP));
2972  else                     mpz_init_set(P,nP->z);
2973  const mp_bitcnt_t bits=2*(mpz_size1(P)+1)*GMP_LIMB_BITS;
2974  mpz_init2(N,bits);
2975  if (SR_HDL(nN) & SR_INT) mpz_set_si(N,SR_TO_INT(nN));
2976  else                     mpz_set(N,nN->z);
2977  assume(!mpz_isNeg(P));
2978  if (mpz_isNeg(N))  mpz_add(N,N,P);
2979  mpz_init2(A,bits); mpz_set_ui(A,0L);
2980  mpz_init2(B,bits); mpz_set_ui(B,1L);
2981  mpz_init2(C,bits); mpz_set_ui(C,0L);
2982  mpz_init2(D,bits);
2983  mpz_init2(E,bits); mpz_set(E,P);
2984  mpz_init2(tmp,bits);
2985  number z=INT_TO_SR(0);
2986  while(mpz_sgn1(N)!=0)
2987  {
2988    mpz_mul(tmp,N,N);
2989    mpz_add(tmp,tmp,tmp);
2990    if (mpz_cmp(tmp,P)<0)
2991    {
2992       if (mpz_isNeg(B))
2993       {
2994         mpz_neg(B,B);
2995         mpz_neg(N,N);
2996       }
2997       // check for gcd(N,B)==1
2998       mpz_gcd(tmp,N,B);
2999       if (mpz_cmp_ui(tmp,1)==0)
3000       {
3001         // return N/B
3002         z=ALLOC_RNUMBER();
3003         #ifdef LDEBUG
3004         z->debug=123456;
3005         #endif
3006         memcpy(z->z,N,sizeof(mpz_t));
3007         memcpy(z->n,B,sizeof(mpz_t));
3008         z->s = 0;
3009         nlNormalize(z,r);
3010       }
3011       else
3012       {
3013         // return nN (the input) instead of "fail"
3014         z=nlCopy(nN,r);
3015         mpz_clear(B);
3016         mpz_clear(N);
3017       }
3018       break;
3019    }
3020    //mpz_mod(D,E,N);
3021    //mpz_div(tmp,E,N);
3022    mpz_divmod(tmp,D,E,N);
3023    mpz_mul(tmp,tmp,B);
3024    mpz_sub(C,A,tmp);
3025    mpz_set(E,N);
3026    mpz_set(N,D);
3027    mpz_set(A,B);
3028    mpz_set(B,C);
3029  }
3030  mpz_clear(tmp);
3031  mpz_clear(A);
3032  mpz_clear(C);
3033  mpz_clear(D);
3034  mpz_clear(E);
3035  mpz_clear(P);
3036  return z;
3037}
3038
3039number nlExtGcd(number a, number b, number *s, number *t, const coeffs)
3040{
3041  mpz_ptr aa,bb;
3042  *s=ALLOC_RNUMBER();
3043  mpz_init((*s)->z); (*s)->s=3;
3044  (*t)=ALLOC_RNUMBER();
3045  mpz_init((*t)->z); (*t)->s=3;
3046  number g=ALLOC_RNUMBER();
3047  mpz_init(g->z); g->s=3;
3048  #ifdef LDEBUG
3049  g->debug=123456;
3050  (*s)->debug=123456;
3051  (*t)->debug=123456;
3052  #endif
3053  if (SR_HDL(a) & SR_INT)
3054  {
3055    aa=(mpz_ptr)omAlloc(sizeof(mpz_t));
3056    mpz_init_set_si(aa,SR_TO_INT(a));
3057  }
3058  else
3059  {
3060    aa=a->z;
3061  }
3062  if (SR_HDL(b) & SR_INT)
3063  {
3064    bb=(mpz_ptr)omAlloc(sizeof(mpz_t));
3065    mpz_init_set_si(bb,SR_TO_INT(b));
3066  }
3067  else
3068  {
3069    bb=b->z;
3070  }
3071  mpz_gcdext(g->z,(*s)->z,(*t)->z,aa,bb);
3072  g=nlShort3(g);
3073  (*s)=nlShort3((*s));
3074  (*t)=nlShort3((*t));
3075  if (SR_HDL(a) & SR_INT)
3076  {
3077    mpz_clear(aa);
3078    omFreeSize(aa, sizeof(mpz_t));
3079  }
3080  if (SR_HDL(b) & SR_INT)
3081  {
3082    mpz_clear(bb);
3083    omFreeSize(bb, sizeof(mpz_t));
3084  }
3085  return g;
3086}
3087
3088//void    nlCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
3089//{
3090//  if (r->is_field)  PrintS("QQ");
3091//  else  PrintS("ZZ");
3092//}
3093
3094VAR int n_SwitchChinRem=0;
3095number   nlChineseRemainderSym(number *x, number *q,int rl, BOOLEAN sym, CFArray &inv_cache,const coeffs CF)
3096// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
3097{
3098  setCharacteristic( 0 ); // only in char 0
3099  Off(SW_RATIONAL);
3100  CFArray X(rl), Q(rl);
3101  int i;
3102  for(i=rl-1;i>=0;i--)
3103  {
3104    X[i]=CF->convSingNFactoryN(x[i],FALSE,CF); // may be larger MAX_INT
3105    Q[i]=CF->convSingNFactoryN(q[i],FALSE,CF); // may be larger MAX_INT
3106  }
3107  CanonicalForm xnew,qnew;
3108  if (n_SwitchChinRem)
3109    chineseRemainder(X,Q,xnew,qnew);
3110  else
3111    chineseRemainderCached(X,Q,xnew,qnew,inv_cache);
3112  number n=CF->convFactoryNSingN(xnew,CF);
3113  if (sym)
3114  {
3115    number p=CF->convFactoryNSingN(qnew,CF);
3116    number p2;
3117    if (getCoeffType(CF) == n_Q) p2=nlIntDiv(p,nlInit(2, CF),CF);
3118    else                         p2=CF->cfDiv(p,CF->cfInit(2, CF),CF);
3119    if (CF->cfGreater(n,p2,CF))
3120    {
3121       number n2=CF->cfSub(n,p,CF);
3122       CF->cfDelete(&n,CF);
3123       n=n2;
3124    }
3125    CF->cfDelete(&p2,CF);
3126    CF->cfDelete(&p,CF);
3127  }
3128  CF->cfNormalize(n,CF);
3129  return n;
3130}
3131#if 0
3132number nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
3133{
3134  CFArray inv(rl);
3135  return nlChineseRemainderSym(x,q,rl,TRUE,inv,C);
3136}
3137#endif
3138
3139static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
3140{
3141  assume(cf != NULL);
3142
3143  numberCollectionEnumerator.Reset();
3144
3145  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
3146  {
3147    c = nlInit(1, cf);
3148    return;
3149  }
3150
3151  // all coeffs are given by integers!!!
3152
3153  // part 1, find a small candidate for gcd
3154  number cand1,cand;
3155  int s1,s;
3156  s=2147483647; // max. int
3157
3158  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
3159
3160  int normalcount = 0;
3161  do
3162  {
3163    number& n = numberCollectionEnumerator.Current();
3164    nlNormalize(n, cf); ++normalcount;
3165    cand1 = n;
3166
3167    if (SR_HDL(cand1)&SR_INT) { cand=cand1; break; }
3168    assume(cand1->s==3); // all coeffs should be integers // ==0?!! after printing
3169    s1=mpz_size1(cand1->z);
3170    if (s>s1)
3171    {
3172      cand=cand1;
3173      s=s1;
3174    }
3175  } while (numberCollectionEnumerator.MoveNext() );
3176
3177//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
3178
3179  cand=nlCopy(cand,cf);
3180  // part 2: compute gcd(cand,all coeffs)
3181
3182  numberCollectionEnumerator.Reset();
3183
3184  while (numberCollectionEnumerator.MoveNext() )
3185  {
3186    number& n = numberCollectionEnumerator.Current();
3187
3188    if( (--normalcount) <= 0)
3189      nlNormalize(n, cf);
3190
3191    nlInpGcd(cand, n, cf);
3192    assume( nlGreaterZero(cand,cf) );
3193
3194    if(nlIsOne(cand,cf))
3195    {
3196      c = cand;
3197
3198      if(!lc_is_pos)
3199      {
3200        // make the leading coeff positive
3201        c = nlNeg(c, cf);
3202        numberCollectionEnumerator.Reset();
3203
3204        while (numberCollectionEnumerator.MoveNext() )
3205        {
3206          number& nn = numberCollectionEnumerator.Current();
3207          nn = nlNeg(nn, cf);
3208        }
3209      }
3210      return;
3211    }
3212  }
3213
3214  // part3: all coeffs = all coeffs / cand
3215  if (!lc_is_pos)
3216    cand = nlNeg(cand,cf);
3217
3218  c = cand;
3219  numberCollectionEnumerator.Reset();
3220
3221  while (numberCollectionEnumerator.MoveNext() )
3222  {
3223    number& n = numberCollectionEnumerator.Current();
3224    number t=nlExactDiv(n, cand, cf); // simple integer exact division, no ratios to remain
3225    nlDelete(&n, cf);
3226    n = t;
3227  }
3228}
3229
3230static void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
3231{
3232  assume(cf != NULL);
3233
3234  numberCollectionEnumerator.Reset();
3235
3236  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
3237  {
3238    c = nlInit(1, cf);
3239//    assume( n_GreaterZero(c, cf) );
3240    return;
3241  }
3242
3243  // all coeffs are given by integers after returning from this routine
3244
3245  // part 1, collect product of all denominators /gcds
3246  number cand;
3247  cand=ALLOC_RNUMBER();
3248#if defined(LDEBUG)
3249  cand->debug=123456;
3250#endif
3251  cand->s=3;
3252
3253  int s=0;
3254
3255  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
3256
3257  do
3258  {
3259    number& cand1 = numberCollectionEnumerator.Current();
3260
3261    if (!(SR_HDL(cand1)&SR_INT))
3262    {
3263      nlNormalize(cand1, cf);
3264      if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
3265      && (cand1->s==1))             // and is a normalised rational
3266      {
3267        if (s==0) // first denom, we meet
3268        {
3269          mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
3270          s=1;
3271        }
3272        else // we have already something
3273        {
3274          mpz_lcm(cand->z, cand->z, cand1->n);
3275        }
3276      }
3277    }
3278  }
3279  while (numberCollectionEnumerator.MoveNext() );
3280
3281
3282  if (s==0) // nothing to do, all coeffs are already integers
3283  {
3284//    mpz_clear(tmp);
3285    FREE_RNUMBER(cand);
3286    if (lc_is_pos)
3287      c=nlInit(1,cf);
3288    else
3289    {
3290      // make the leading coeff positive
3291      c=nlInit(-1,cf);
3292
3293      // TODO: incorporate the following into the loop below?
3294      numberCollectionEnumerator.Reset();
3295      while (numberCollectionEnumerator.MoveNext() )
3296      {
3297        number& n = numberCollectionEnumerator.Current();
3298        n = nlNeg(n, cf);
3299      }
3300    }
3301//    assume( n_GreaterZero(c, cf) );
3302    return;
3303  }
3304
3305  cand = nlShort3(cand);
3306
3307  // part2: all coeffs = all coeffs * cand
3308  // make the lead coeff positive
3309  numberCollectionEnumerator.Reset();
3310
3311  if (!lc_is_pos)
3312    cand = nlNeg(cand, cf);
3313
3314  c = cand;
3315
3316  while (numberCollectionEnumerator.MoveNext() )
3317  {
3318    number &n = numberCollectionEnumerator.Current();
3319    nlInpMult(n, cand, cf);
3320  }
3321
3322}
3323
3324char * nlCoeffName(const coeffs r)
3325{
3326  if (r->cfDiv==nlDiv) return (char*)"QQ";
3327  else                 return (char*)"ZZ";
3328}
3329
3330void nlWriteFd(number n, const ssiInfo* d, const coeffs)
3331{
3332  if(SR_HDL(n) & SR_INT)
3333  {
3334    #if SIZEOF_LONG == 4
3335    fprintf(d->f_write,"4 %ld ",SR_TO_INT(n));
3336    #else
3337    long nn=SR_TO_INT(n);
3338    if ((nn<POW_2_28_32)&&(nn>= -POW_2_28_32))
3339    {
3340      int nnn=(int)nn;
3341      fprintf(d->f_write,"4 %d ",nnn);
3342    }
3343    else
3344    {
3345      mpz_t tmp;
3346      mpz_init_set_si(tmp,nn);
3347      fputs("8 ",d->f_write);
3348      mpz_out_str (d->f_write,SSI_BASE, tmp);
3349      fputc(' ',d->f_write);
3350      mpz_clear(tmp);
3351    }
3352    #endif
3353  }
3354  else if (n->s<2)
3355  {
3356    //gmp_fprintf(f,"%d %Zd %Zd ",n->s,n->z,n->n);
3357    fprintf(d->f_write,"%d ",n->s+5);
3358    mpz_out_str (d->f_write,SSI_BASE, n->z);
3359    fputc(' ',d->f_write);
3360    mpz_out_str (d->f_write,SSI_BASE, n->n);
3361    fputc(' ',d->f_write);
3362
3363    //if (d->f_debug!=NULL) gmp_fprintf(d->f_debug,"number: s=%d gmp/gmp \"%Zd %Zd\" ",n->s,n->z,n->n);
3364  }
3365  else /*n->s==3*/
3366  {
3367    //gmp_fprintf(d->f_write,"3 %Zd ",n->z);
3368    fputs("8 ",d->f_write);
3369    mpz_out_str (d->f_write,SSI_BASE, n->z);
3370    fputc(' ',d->f_write);
3371
3372    //if (d->f_debug!=NULL) gmp_fprintf(d->f_debug,"number: gmp \"%Zd\" ",n->z);
3373  }
3374}
3375
3376number nlReadFd(const ssiInfo *d, const coeffs)
3377{
3378  int sub_type=-1;
3379  sub_type=s_readint(d->f_read);
3380  switch(sub_type)
3381  {
3382     case 0:
3383     case 1:
3384       {// read mpz_t, mpz_t
3385         number n=nlRInit(0);
3386         mpz_init(n->n);
3387         s_readmpz(d->f_read,n->z);
3388         s_readmpz(d->f_read,n->n);
3389         n->s=sub_type;
3390         return n;
3391       }
3392
3393     case 3:
3394       {// read mpz_t
3395         number n=nlRInit(0);
3396         s_readmpz(d->f_read,n->z);
3397         n->s=3; /*sub_type*/
3398         #if SIZEOF_LONG == 8
3399         n=nlShort3(n);
3400         #endif
3401         return n;
3402       }
3403     case 4:
3404       {
3405         LONG dd=s_readlong(d->f_read);
3406         //#if SIZEOF_LONG == 8
3407         return INT_TO_SR(dd);
3408         //#else
3409         //return nlInit(dd,NULL);
3410         //#endif
3411       }
3412     case 5:
3413     case 6:
3414       {// read raw mpz_t, mpz_t
3415         number n=nlRInit(0);
3416         mpz_init(n->n);
3417         s_readmpz_base (d->f_read,n->z, SSI_BASE);
3418         s_readmpz_base (d->f_read,n->n, SSI_BASE);
3419         n->s=sub_type-5;
3420         return n;
3421       }
3422     case 8:
3423       {// read raw mpz_t
3424         number n=nlRInit(0);
3425         s_readmpz_base (d->f_read,n->z, SSI_BASE);
3426         n->s=sub_type=3; /*subtype-5*/
3427         #if SIZEOF_LONG == 8
3428         n=nlShort3(n);
3429         #endif
3430         return n;
3431       }
3432
3433     default: Werror("error in reading number: invalid subtype %d",sub_type);
3434              return NULL;
3435  }
3436  return NULL;
3437}
3438
3439BOOLEAN nlCoeffIsEqual(const coeffs r, n_coeffType n, void *p)
3440{
3441  /* test, if r is an instance of nInitCoeffs(n,parameter) */
3442  /* if parameter is not needed */
3443  if (n==r->type)
3444  {
3445    if ((p==NULL)&&(r->cfDiv==nlDiv)) return TRUE;
3446    if ((p!=NULL)&&(r->cfDiv!=nlDiv)) return TRUE;
3447  }
3448  return FALSE;
3449}
3450
3451static number nlLcm(number a,number b,const coeffs r)
3452{
3453  number g=nlGcd(a,b,r);
3454  number n1=nlMult(a,b,r);
3455  number n2=nlExactDiv(n1,g,r);
3456  nlDelete(&g,r);
3457  nlDelete(&n1,r);
3458  return n2;
3459}
3460
3461static number nlRandom(siRandProc p, number v2, number, const coeffs cf)
3462{
3463  number a=nlInit(p(),cf);
3464  if (v2!=NULL)
3465  {
3466    number b=nlInit(p(),cf);
3467    number c=nlDiv(a,b,cf);
3468    nlDelete(&b,cf);
3469    nlDelete(&a,cf);
3470    a=c;
3471  }
3472  return a;
3473}
3474
3475BOOLEAN nlInitChar(coeffs r, void*p)
3476{
3477  r->is_domain=TRUE;
3478  r->rep=n_rep_gap_rat;
3479
3480  r->nCoeffIsEqual=nlCoeffIsEqual;
3481  //r->cfKillChar = ndKillChar; /* dummy */
3482  //r->cfCoeffString=nlCoeffString;
3483  r->cfCoeffName=nlCoeffName;
3484
3485  r->cfInitMPZ = nlInitMPZ;
3486  r->cfMPZ  = nlMPZ;
3487
3488  r->cfMult  = nlMult;
3489  r->cfSub   = nlSub;
3490  r->cfAdd   = nlAdd;
3491  r->cfExactDiv= nlExactDiv;
3492  if (p==NULL) /* Q */
3493  {
3494    r->is_field=TRUE;
3495    r->cfDiv   = nlDiv;
3496    //r->cfGcd  = ndGcd_dummy;
3497    r->cfSubringGcd  = nlGcd;
3498  }
3499  else /* Z: coeffs_BIGINT */
3500  {
3501    r->is_field=FALSE;
3502    r->cfDiv   = nlIntDiv;
3503    r->cfIntMod= nlIntMod;
3504    r->cfGcd  = nlGcd;
3505    r->cfDivBy=nlDivBy;
3506    r->cfDivComp = nlDivComp;
3507    r->cfIsUnit = nlIsUnit;
3508    r->cfGetUnit = nlGetUnit;
3509    r->cfQuot1 = nlQuot1;
3510    r->cfLcm = nlLcm;
3511    r->cfXExtGcd=nlXExtGcd;
3512    r->cfQuotRem=nlQuotRem;
3513  }
3514  r->cfInit = nlInit;
3515  r->cfSize  = nlSize;
3516  r->cfInt  = nlInt;
3517
3518  r->cfChineseRemainder=nlChineseRemainderSym;
3519  r->cfFarey=nlFarey;
3520  r->cfInpNeg   = nlNeg;
3521  r->cfInvers= nlInvers;
3522  r->cfCopy  = nlCopy;
3523  r->cfRePart = nlCopy;
3524  //r->cfImPart = ndReturn0;
3525  r->cfWriteLong = nlWrite;
3526  r->cfRead = nlRead;
3527  r->cfNormalize=nlNormalize;
3528  r->cfGreater = nlGreater;
3529  r->cfEqual = nlEqual;
3530  r->cfIsZero = nlIsZero;
3531  r->cfIsOne = nlIsOne;
3532  r->cfIsMOne = nlIsMOne;
3533  r->cfGreaterZero = nlGreaterZero;
3534  r->cfPower = nlPower;
3535  r->cfGetDenom = nlGetDenom;
3536  r->cfGetNumerator = nlGetNumerator;
3537  r->cfExtGcd = nlExtGcd; // only for ring stuff and Z
3538  r->cfNormalizeHelper  = nlNormalizeHelper;
3539  r->cfDelete= nlDelete;
3540  r->cfSetMap = nlSetMap;
3541  //r->cfName = ndName;
3542  r->cfInpMult=nlInpMult;
3543  r->cfInpAdd=nlInpAdd;
3544  //r->cfCoeffWrite=nlCoeffWrite;
3545
3546  r->cfClearContent = nlClearContent;
3547  r->cfClearDenominators = nlClearDenominators;
3548
3549#ifdef LDEBUG
3550  // debug stuff
3551  r->cfDBTest=nlDBTest;
3552#endif
3553  r->convSingNFactoryN=nlConvSingNFactoryN;
3554  r->convFactoryNSingN=nlConvFactoryNSingN;
3555
3556  r->cfRandom=nlRandom;
3557
3558  // io via ssi
3559  r->cfWriteFd=nlWriteFd;
3560  r->cfReadFd=nlReadFd;
3561
3562  //r->type = n_Q;
3563  r->ch = 0;
3564  r->has_simple_Alloc=FALSE;
3565  r->has_simple_Inverse=FALSE;
3566
3567  // variables for this type of coeffs:
3568  // (none)
3569  return FALSE;
3570}
3571#if 0
3572number nlMod(number a, number b)
3573{
3574  if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
3575  {
3576    int bi=SR_TO_INT(b);
3577    int ai=SR_TO_INT(a);
3578    int bb=ABS(bi);
3579    int c=ai%bb;
3580    if (c<0)  c+=bb;
3581    return (INT_TO_SR(c));
3582  }
3583  number al;
3584  number bl;
3585  if (SR_HDL(a))&SR_INT)
3586    al=nlRInit(SR_TO_INT(a));
3587  else
3588    al=nlCopy(a);
3589  if (SR_HDL(b))&SR_INT)
3590    bl=nlRInit(SR_TO_INT(b));
3591  else
3592    bl=nlCopy(b);
3593  number r=nlRInit(0);
3594  mpz_mod(r->z,al->z,bl->z);
3595  nlDelete(&al);
3596  nlDelete(&bl);
3597  if (mpz_size1(&r->z)<=MP_SMALL)
3598  {
3599    LONG ui=(int)mpz_get_si(&r->z);
3600    if ((((ui<<3)>>3)==ui)
3601    && (mpz_cmp_si(x->z,(long)ui)==0))
3602    {
3603      mpz_clear(&r->z);
3604      FREE_RNUMBER(r); //  omFreeBin((void *)r, rnumber_bin);
3605      r=INT_TO_SR(ui);
3606    }
3607  }
3608  return r;
3609}
3610#endif
3611#endif // not P_NUMBERS_H
3612#endif // LONGRAT_CC
Note: See TracBrowser for help on using the repository browser.