source: git/libpolys/coeffs/longrat.cc @ 778b36

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