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

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