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

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