source: git/libpolys/coeffs/rintegers.cc @ e6cdad

spielwiese
Last change on this file since e6cdad was e6cdad, checked in by Hans Schoenemann <hannes@…>, 6 years ago
use xalloc if omalloc is disabled
  • Property mode set to 100644
File size: 40.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: numbers (integers)
6*/
7
8#include "misc/auxiliary.h"
9
10#include "factory/factory.h"
11
12#include "misc/mylimits.h"
13#include "reporter/reporter.h"
14
15#include "coeffs/coeffs.h"
16#include "coeffs/longrat.h"
17#include "coeffs/numbers.h"
18
19#include "coeffs/si_gmp.h"
20
21#include "coeffs/mpr_complex.h"
22#include "coeffs/rintegers.h"
23#include "coeffs/rmodulon.h"
24#include "coeffs/longrat.h"
25
26#include <string.h>
27
28#ifdef HAVE_RINGS
29
30omBin gmp_nrz_bin = omGetSpecBin(sizeof(mpz_t));
31
32#if SI_INTEGER_VARIANT == 2
33/*
34 * Multiply two numbers
35 */
36static number nrzMult (number a, number b, const coeffs)
37{
38  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
39  mpz_init(erg);
40  mpz_mul(erg, (mpz_ptr) a, (mpz_ptr) b);
41  return (number) erg;
42}
43
44/*
45 * Give the smallest non unit k, such that a * x = k = b * y has a solution
46 */
47static number nrzLcm (number a,number b,const coeffs)
48{
49  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
50  mpz_init(erg);
51  mpz_lcm(erg, (mpz_ptr) a, (mpz_ptr) b);
52  return (number) erg;
53}
54
55/*
56 * Give the largest non unit k, such that a = x * k, b = y * k has
57 * a solution.
58 */
59static number nrzGcd (number a,number b,const coeffs)
60{
61  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
62  mpz_init(erg);
63  mpz_gcd(erg, (mpz_ptr) a, (mpz_ptr) b);
64  return (number) erg;
65}
66
67/*
68 * Give the largest non unit k, such that a = x * k, b = y * k has
69 * a solution and r, s, s.t. k = s*a + t*b
70 */
71static number nrzExtGcd (number a, number b, number *s, number *t, const coeffs)
72{
73  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
74  mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
75  mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
76  mpz_init(erg);
77  mpz_init(bs);
78  mpz_init(bt);
79  mpz_gcdext(erg, bs, bt, (mpz_ptr) a, (mpz_ptr) b);
80  *s = (number) bs;
81  *t = (number) bt;
82  return (number) erg;
83}
84
85static number nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )
86{
87  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
88  mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
89  mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
90  mpz_init(erg);
91  mpz_init(bs);
92  mpz_init(bt);
93
94  mpz_gcdext(erg, bs, bt, (mpz_ptr)a, (mpz_ptr)b);
95
96  mpz_ptr bu = (mpz_ptr) omAllocBin(gmp_nrz_bin);
97  mpz_ptr bv = (mpz_ptr) omAllocBin(gmp_nrz_bin);
98
99  mpz_init_set(bu, (mpz_ptr) b);
100  mpz_init_set(bv, (mpz_ptr) a);
101
102  assume(mpz_cmp_si(erg, 0));
103
104  mpz_div(bu, bu, erg);
105  mpz_div(bv, bv, erg);
106
107  mpz_mul_si(bu, bu, -1);
108  *u = (number) bu;
109  *v = (number) bv;
110
111  *s = (number) bs;
112  *t = (number) bt;
113  return (number) erg;
114}
115
116static void nrzPower (number a, int i, number * result, const coeffs)
117{
118  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
119  mpz_init(erg);
120  mpz_pow_ui(erg, (mpz_ptr) a, i);
121  *result = (number) erg;
122}
123
124/*
125 * create a number from int
126 */
127number nrzInit (long i, const coeffs)
128{
129  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
130  mpz_init_set_si(erg, i);
131  return (number) erg;
132}
133
134void nrzDelete(number *a, const coeffs)
135{
136  if (*a != NULL)
137  {
138    mpz_clear((mpz_ptr) *a);
139    omFreeBin((ADDRESS) *a, gmp_nrz_bin);
140    *a = NULL;
141  }
142}
143
144static number nrzCopy(number a, const coeffs)
145{
146  if (a==NULL) return NULL;
147  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
148  mpz_init_set(erg, (mpz_ptr) a);
149  return (number) erg;
150}
151
152#if 0
153number nrzCopyMap(number a, const coeffs /*src*/, const coeffs dst)
154{
155  return nrzCopy(a,dst);
156}
157#endif
158
159int nrzSize(number a, const coeffs)
160{
161  mpz_ptr p=(mpz_ptr)a;
162  int s=p->_mp_alloc;
163  if (s==1) s=(mpz_cmp_ui(p,0)!=0);
164  return s;
165
166}
167
168/*
169 * convert a number to int
170 */
171static long nrzInt(number &n, const coeffs)
172{
173  return mpz_get_si( (mpz_ptr)n);
174}
175
176static number nrzAdd (number a, number b, const coeffs)
177{
178  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
179  mpz_init(erg);
180  mpz_add(erg, (mpz_ptr) a, (mpz_ptr) b);
181  return (number) erg;
182}
183
184static number nrzSub (number a, number b, const coeffs)
185{
186  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
187  mpz_init(erg);
188  mpz_sub(erg, (mpz_ptr) a, (mpz_ptr) b);
189  return (number) erg;
190}
191
192static number nrzGetUnit (number, const coeffs r)
193{
194  return nrzInit(1, r);
195}
196
197static BOOLEAN nrzIsUnit (number a, const coeffs)
198{
199  return 0 == mpz_cmpabs_ui((mpz_ptr) a, 1);
200}
201
202static BOOLEAN nrzIsZero (number  a, const coeffs)
203{
204  return 0 == mpz_cmpabs_ui((mpz_ptr) a, 0);
205}
206
207static BOOLEAN nrzIsOne (number a, const coeffs)
208{
209  return (0 == mpz_cmp_ui((mpz_ptr) a, 1));
210}
211
212static BOOLEAN nrzIsMOne (number a, const coeffs)
213{
214  return (0 == mpz_cmp_si((mpz_ptr) a, -1));
215}
216
217static BOOLEAN nrzEqual (number a,number b, const coeffs)
218{
219  return 0 == mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
220}
221
222static BOOLEAN nrzGreater (number a,number b, const coeffs)
223{
224  return 0 < mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
225}
226
227static BOOLEAN nrzGreaterZero (number k, const coeffs)
228{
229  return 0 < mpz_sgn1((mpz_ptr) k);
230}
231
232static BOOLEAN nrzDivBy (number a,number b, const coeffs)
233{
234  return mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b) != 0;
235}
236
237static int nrzDivComp(number a, number b, const coeffs r)
238{
239  if (nrzDivBy(a, b, r))
240  {
241    if (nrzDivBy(b, a, r)) return 2;
242    return -1;
243  }
244  if (nrzDivBy(b, a, r)) return 1;
245  return 0;
246}
247
248static number nrzDiv (number a,number b, const coeffs)
249{
250  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
251  mpz_init(erg);
252  mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);
253  mpz_init(r);
254  mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);
255  //if (!nrzIsZero((number) r, R))
256  //{
257  //  WerrorS("Division by non divisible element.");
258  //  WerrorS("Result is without remainder.");
259  //}
260  mpz_clear(r);
261  omFreeBin(r, gmp_nrz_bin);
262  return (number) erg;
263}
264
265static number nrzExactDiv (number a,number b, const coeffs)
266{
267  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
268  mpz_init(erg);
269  mpz_tdiv_q(erg, (mpz_ptr) a, (mpz_ptr) b);
270  return (number) erg;
271}
272
273static number nrzQuotRem (number a, number b, number * r, const coeffs )
274{
275  mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);
276  mpz_init(qq);
277  mpz_ptr rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);
278  mpz_init(rr);
279  mpz_tdiv_qr(qq, rr, (mpz_ptr) a, (mpz_ptr) b);
280  if (r==NULL)
281  {
282    mpz_clear(rr);
283    omFreeBin(rr,gmp_nrz_bin);
284  }
285  else
286  {
287    *r=(number)rr;
288  }
289  return (number) qq;
290}
291
292static number nrzIntMod (number a,number b, const coeffs)
293{
294  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
295  mpz_init(erg);
296  mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);
297  mpz_init(r);
298  mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);
299  mpz_clear(erg);
300  omFreeBin(erg, gmp_nrz_bin);
301  return (number) r;
302}
303
304static number  nrzInvers (number c, const coeffs r)
305{
306  if (!nrzIsUnit((number) c, r))
307  {
308    WerrorS("Non invertible element.");
309    return (number)NULL;
310  }
311  return nrzCopy(c,r);
312}
313
314static number nrzNeg (number c, const coeffs)
315{
316// nNeg inplace !!!
317  mpz_mul_si((mpz_ptr) c, (mpz_ptr) c, -1);
318  return c;
319}
320
321static number nrzMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
322{
323  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
324  mpz_init_set_ui(erg, (unsigned long) from);
325  return (number) erg;
326}
327
328static number nrzMapZp(number from, const coeffs /*src*/, const coeffs /*dst*/)
329{
330  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
331  mpz_init_set_si(erg, (long) from);
332  return (number) erg;
333}
334
335static number nrzMapQ(number from, const coeffs src, const coeffs /*dst*/)
336{
337  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
338  mpz_init(erg);
339  nlGMP(from, erg, src); // FIXME? TODO? // extern void   nlGMP(number &i, number n, const coeffs r); // to be replaced with n_MPZ(erg, from, src); // ?
340  return (number) erg;
341}
342
343static nMapFunc nrzSetMap(const coeffs src, const coeffs /*dst*/)
344{
345  /* dst = currRing */
346  /* dst = nrn */
347  if ((src->rep==n_rep_gmp)
348  && (nCoeff_is_Z(src) || nCoeff_is_Zn(src) || nCoeff_is_Ring_PtoM(src)))
349  {
350    return ndCopyMap; //nrzCopyMap;
351  }
352  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Z(src)*/)
353  {
354    return ndCopyMap; //nrzCopyMap;
355  }
356  if (nCoeff_is_Ring_2toM(src))
357  {
358    return nrzMapMachineInt;
359  }
360  if (nCoeff_is_Zp(src))
361  {
362    return nrzMapZp;
363  }
364  if (getCoeffType(src)==n_Q /*nCoeff_is_Q(src) or coeffs_BIGINT*/)
365  {
366    return nrzMapQ;
367  }
368  return NULL;      // default
369}
370
371/*
372 * set the exponent (allocate and init tables) (TODO)
373 */
374
375void nrzSetExp(int, coeffs)
376{
377}
378
379void nrzInitExp(int, coeffs)
380{
381}
382
383#ifdef LDEBUG
384static BOOLEAN nrzDBTest (number, const char *, const int, const coeffs)
385{
386  return TRUE;//TODO
387}
388#endif
389
390void nrzWrite (number a, const coeffs)
391{
392  char *s,*z;
393  if (a==NULL)
394  {
395    StringAppendS("o");
396  }
397  else
398  {
399    int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2;
400    s=(char*)omAlloc(l);
401    z=mpz_get_str(s,10,(mpz_ptr) a);
402    StringAppendS(z);
403    omFreeSize((ADDRESS)s,l);
404  }
405}
406
407/*2
408* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
409*/
410static const char * nlEatLongC(char *s, mpz_ptr i)
411{
412  const char * start=s;
413
414  if (*s<'0' || *s>'9')
415  {
416    mpz_set_ui(i,1);
417    return s;
418  }
419  while (*s >= '0' && *s <= '9') s++;
420  if (*s=='\0')
421  {
422    mpz_set_str(i,start,10);
423  }
424  else
425  {
426    char c=*s;
427    *s='\0';
428    mpz_set_str(i,start,10);
429    *s=c;
430  }
431  return s;
432}
433
434
435static CanonicalForm nrzConvSingNFactoryN(number n, BOOLEAN setChar, const coeffs /*r*/)
436{
437  if (setChar) setCharacteristic( 0 );
438
439  CanonicalForm term;
440  mpz_t num;
441  mpz_init_set(num, *((mpz_t*)n));
442  term = make_cf(num);
443  return term;
444}
445
446static number nrzConvFactoryNSingN(const CanonicalForm n, const coeffs r)
447{
448  if (n.isImm())
449    return nrzInit(n.intval(),r);
450  else
451  {
452    mpz_ptr m = (mpz_ptr) omAllocBin(gmp_nrz_bin);
453    gmp_numerator(n,m);
454    if (!n.den().isOne())
455    {
456      WarnS("denominator is not 1 in factory");
457    }
458    return (number) m;
459  }
460}
461
462static const char * nrzRead (const char *s, number *a, const coeffs)
463{
464  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
465  {
466    mpz_init(z);
467    s = nlEatLongC((char *) s, z);
468  }
469  *a = (number) z;
470  return s;
471}
472
473static void nrzCoeffWrite  (const coeffs, BOOLEAN /*details*/)
474{
475  PrintS("ZZ");
476}
477
478static char* nrzCoeffName(const coeffs)
479{
480  return (char*)"ZZ";
481}
482
483static char* nrzCoeffString(const coeffs cf)
484{
485  return omStrDup(nrzCoeffName(cf));
486}
487
488static coeffs nrzQuot1(number c, const coeffs r)
489{
490    long ch = r->cfInt(c, r);
491    mpz_t dummy;
492    mpz_init_set_ui(dummy, ch);
493    ZnmInfo info;
494    info.base = dummy;
495    info.exp = (unsigned long) 1;
496    coeffs rr = nInitChar(n_Zn, (void*)&info);
497    mpz_clear(dummy);
498    return(rr);
499}
500
501static number nrzInitMPZ(mpz_t m, const coeffs)
502{
503  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
504  mpz_init_set(z, m);
505  return (number)z;
506}
507
508static number nrzFarey(number r, number N, const coeffs R)
509{
510  number a0 = nrzCopy(N, R);
511  number b0 = nrzInit(0, R);
512  number a1 = nrzCopy(r, R);
513  number b1 = nrzInit(1, R);
514  number two = nrzInit(2, R);
515#if 0
516  PrintS("Farey start with ");
517  n_Print(r, R);
518  PrintS(" mod ");
519  n_Print(N, R);
520  PrintLn();
521#endif
522  while (1)
523  {
524    number as = nrzMult(a1, a1, R);
525    n_InpMult(as, two, R);
526    if (nrzGreater(N, as, R))
527    {
528      nrzDelete(&as, R);
529      break;
530    }
531    nrzDelete(&as, R);
532    number q = nrzDiv(a0, a1, R);
533    number t = nrzMult(a1, q, R),
534           s = nrzSub(a0, t, R);
535    nrzDelete(&a0, R);
536    a0 = a1;
537    a1 = s;
538    nrzDelete(&t, R);
539
540    t = nrzMult(b1, q, R);
541    s = nrzSub(b0, t, R);
542    nrzDelete(&b0, R);
543    b0 = b1;
544    b1 = s;
545    nrzDelete(&t, R);
546    nrzDelete(&q, R);
547  }
548  number as = nrzMult(b1, b1, R);
549  n_InpMult(as, two, R);
550  nrzDelete(&two, R);
551  if (nrzGreater(as, N, R))
552  {
553    nrzDelete(&a0, R);
554    nrzDelete(&a1, R);
555    nrzDelete(&b0, R);
556    nrzDelete(&b1, R);
557    nrzDelete(&as, R);
558    return NULL;
559  }
560  nrzDelete(&as, R);
561  nrzDelete(&a0, R);
562  nrzDelete(&b0, R);
563
564  number a, b, ab;
565  coeffs Q = nInitChar(n_Q, 0);
566  nMapFunc f = n_SetMap(R, Q);
567  a = f(a1, R, Q);
568  b = f(b1, R, Q);
569  ab = n_Div(a, b, Q);
570  n_Delete(&a, Q);
571  n_Delete(&b, Q);
572  nKillChar(Q);
573
574  nrzDelete(&a1, R);
575  nrzDelete(&b1, R);
576  return ab;
577}
578
579void nrzWriteFd(number n, const ssiInfo* d, const coeffs)
580{
581  mpz_out_str (d->f_write,SSI_BASE, (mpz_ptr)n);
582  fputc(' ',d->f_write);
583}
584
585number nrzReadFd(const ssiInfo *d, const coeffs)
586{
587  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
588  mpz_init(erg);
589  s_readmpz_base(d->f_read,erg,SSI_BASE);
590  return (number)erg;
591}
592
593BOOLEAN nrzInitChar(coeffs r,  void *)
594{
595  assume( getCoeffType(r) == n_Z );
596
597  r->is_field=FALSE;
598  r->is_domain=TRUE;
599  r->rep=n_rep_gmp;
600
601  //r->nCoeffIsEqual = ndCoeffIsEqual;
602  r->cfCoeffString = nrzCoeffString;
603  r->cfCoeffName = nrzCoeffName;
604  r->cfCoeffWrite = nrzCoeffWrite;
605  //r->cfKillChar = ndKillChar;
606  r->cfMult  = nrzMult;
607  r->cfSub   = nrzSub;
608  r->cfAdd   = nrzAdd;
609  r->cfDiv   = nrzDiv;
610  r->cfIntMod= nrzIntMod;
611  r->cfExactDiv= nrzExactDiv;
612  r->cfInit = nrzInit;
613  r->cfInitMPZ = nrzInitMPZ;
614  r->cfSize  = nrzSize;
615  r->cfInt  = nrzInt;
616  r->cfDivComp = nrzDivComp;
617  r->cfIsUnit = nrzIsUnit;
618  r->cfGetUnit = nrzGetUnit;
619  r->cfExtGcd = nrzExtGcd;
620  r->cfXExtGcd = nrzXExtGcd;
621  r->cfDivBy = nrzDivBy;
622  r->cfQuotRem = nrzQuotRem;
623  r->cfInpNeg   = nrzNeg;
624  r->cfInvers= nrzInvers;
625  r->cfCopy  = nrzCopy;
626  r->cfWriteLong = nrzWrite;
627  r->cfRead = nrzRead;
628  r->cfGreater = nrzGreater;
629  r->cfEqual = nrzEqual;
630  r->cfIsZero = nrzIsZero;
631  r->cfIsOne = nrzIsOne;
632  r->cfIsMOne = nrzIsMOne;
633  r->cfGreaterZero = nrzGreaterZero;
634  r->cfPower = nrzPower;
635  r->cfGcd  = nrzGcd;
636  r->cfLcm  = nrzLcm;
637  r->cfDelete= nrzDelete;
638  r->cfSetMap = nrzSetMap;
639  r->cfQuot1 = nrzQuot1;
640  r->convSingNFactoryN=nrzConvSingNFactoryN;
641  r->convFactoryNSingN=nrzConvFactoryNSingN;
642  r->cfChineseRemainder=nlChineseRemainderSym;
643  r->cfFarey=nrzFarey;
644  r->cfWriteFd=nrzWriteFd;
645  r->cfReadFd=nrzReadFd;
646  // debug stuff
647
648#ifdef LDEBUG
649  r->cfDBTest=nrzDBTest;
650#endif
651
652  r->ch = 0;
653  r->has_simple_Alloc=FALSE;
654  r->has_simple_Inverse=FALSE;
655  return FALSE;
656}
657
658#elif SI_INTEGER_VARIANT == 3
659
660//make sure that a small number is an immediate integer
661//bascially coped from longrat.cc nlShort3
662//TODO: is there any point in checking 0 first???
663//TODO: it is not clear that this works in 32/64 bit everywhere.
664//      too many hacks.
665#ifdef LDEBUG
666#define nrzTest(A) nrzDBTest(A,__FILE__,__LINE__,NULL)
667BOOLEAN nrzDBTest (number x, const char *f, const int l, const coeffs);
668#else
669#define nrzTest(A)
670#endif
671
672#undef CF_DEBUG
673static inline number nrz_short(number x)
674{
675#if CF_DEBUG
676  StringAppendS("short(");
677  nrzWrite(x, NULL);
678#endif
679  if (mpz_sgn1((mpz_ptr) x)==0)
680  {
681    mpz_clear((mpz_ptr)x);
682    omFreeBin(x, gmp_nrz_bin);
683#if CF_DEBUG
684    StringAppendS(")=0");
685#endif
686    return INT_TO_SR(0);
687  }
688  if (mpz_size1((mpz_ptr)x)<=MP_SMALL)
689  {
690    long ui=mpz_get_si((mpz_ptr)x);
691    if ((((ui<<3)>>3)==ui)
692    && (mpz_cmp_si((mpz_ptr)x,ui)==0))
693    {
694      mpz_clear((mpz_ptr)x);
695      omFreeBin(x, gmp_nrz_bin);
696#if CF_DEBUG
697    StringAppendS(")=imm");
698#endif
699      return INT_TO_SR(ui);
700    }
701  }
702#if CF_DEBUG
703  StringAppendS(")");
704#endif
705  return x;
706}
707
708
709static int nrzSize(number a, const coeffs)
710{
711  if (a==INT_TO_SR(0)) return 0;
712  if (n_Z_IS_SMALL(a)) return 1;
713  return ((mpz_ptr)a)->_mp_alloc;
714}
715
716
717/*
718 * Multiply two numbers
719 * check for 0, 1, -1 maybe
720 */
721#if CF_DEBUG
722number _nrzMult(number, number, const coeffs);
723number nrzMult(number a, number b, const coeffs R)
724{
725  StringSetS("Mult: ");
726  nrzWrite(a, R);
727  StringAppendS(" by ");
728  nrzWrite(b, R);
729  number c = _nrzMult(a, b, R);
730  StringAppendS(" is ");
731  nrzWrite(c, R);
732  char * s = StringEndS();
733  Print("%s\n", s);
734  omFree(s);
735  return c;
736}
737number _nrzMult (number a, number b, const coeffs R)
738#else
739number nrzMult (number a, number b, const coeffs R)
740#endif
741{
742  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
743  {
744  //from longrat.cc
745    if (SR_TO_INT(a)==0)
746      return a;
747    if (SR_TO_INT(b)==0)
748      return b;
749    long r=(long)((unsigned long)(SR_HDL(a)-1L))*((unsigned long)(SR_HDL(b)>>1));
750    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
751    {
752      number u=((number) ((r>>1)+SR_INT));
753    //  if (((((long)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
754      return nrzInit(SR_HDL(u)>>2, R);
755    }
756    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
757    mpz_init(erg);
758    mpz_set_si(erg, SR_TO_INT(a));
759    mpz_mul_si(erg, erg, SR_TO_INT(b));
760    nrzTest((number)erg);
761    return (number) erg;
762  }
763  else if (n_Z_IS_SMALL(a))
764  {
765    if (SR_TO_INT(a)==0)
766      return a;
767    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
768    mpz_init_set(erg, (mpz_ptr) b);
769    mpz_mul_si(erg, erg, SR_TO_INT(a));
770    nrzTest((number)erg);
771    return (number) erg;
772  }
773  else if (n_Z_IS_SMALL(b))
774  {
775    if (SR_TO_INT(b)==0)
776      return b;
777    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
778    mpz_init_set(erg, (mpz_ptr) a);
779    mpz_mul_si(erg, erg, SR_TO_INT(b));
780    nrzTest((number)erg);
781    return (number) erg;
782  }
783  else
784  {
785    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
786    mpz_init(erg);
787    mpz_mul(erg, (mpz_ptr) a, (mpz_ptr) b);
788    nrzTest((number)erg);
789    return (number) erg;
790  }
791}
792
793
794static long int_gcd(long a, long b)
795{
796  long r;
797  a = ABS(a);
798  b = ABS(b);
799  if (!a) return b;
800  if (!b) return a;
801  do
802  {
803    r = a % b;
804    a = b;
805    b = r;
806  } while (b);
807  return ABS(a); // % in c doeas not imply a signn
808                 // it would be unlikely to see a negative here
809                 // but who knows
810}
811
812/*
813 * Give the smallest non unit k, such that a * x = k = b * y has a solution
814 */
815static number nrzLcm (number a, number b, const coeffs R)
816{
817  #ifdef CF_DEBUG
818  PrintS("nrzLcm\n");
819  #endif
820  mpz_ptr erg;
821  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
822  {
823    long g = int_gcd(SR_TO_INT(a), SR_TO_INT(b));
824    return nrzMult(a, INT_TO_SR(SR_TO_INT(b)/g), R);
825  }
826  else
827  {
828    erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
829    if (n_Z_IS_SMALL(a))
830    {
831      mpz_init_set(erg, (mpz_ptr) b);
832      unsigned long g = mpz_gcd_ui(NULL, erg, (unsigned long) ABS(SR_TO_INT(a)));
833      mpz_mul_si(erg, erg, SR_TO_INT(a)/g);
834    }
835    else if (n_Z_IS_SMALL(b))
836    {
837      mpz_init_set(erg, (mpz_ptr) a);
838      unsigned long g = mpz_gcd_ui(NULL, erg, (unsigned long) ABS(SR_TO_INT(b)));
839      mpz_mul_si(erg, erg, SR_TO_INT(b)/g);
840    }
841    else
842    {
843      mpz_init(erg);
844      mpz_lcm(erg, (mpz_ptr) a, (mpz_ptr) b);
845    }
846  }
847  return (number) erg;
848}
849
850static number nrzCopy(number a, const coeffs)
851{
852  if (n_Z_IS_SMALL(a)) return a;
853  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
854  mpz_init_set(erg, (mpz_ptr) a);
855  return (number) erg;
856}
857
858/*
859 * Give the largest non unit k, such that a = x * k, b = y * k has
860 * a solution.
861 */
862static number nrzGcd (number a,number b,const coeffs R)
863{
864  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
865  {
866    long g = int_gcd(SR_TO_INT(a), SR_TO_INT(b));
867    return INT_TO_SR(g);
868  }
869  else if (n_Z_IS_SMALL(a))
870  {
871    if (a==INT_TO_SR(0))
872      return nrzCopy(b, R);
873    unsigned long g = mpz_gcd_ui(NULL, (mpz_ptr)b, (unsigned long) ABS(SR_TO_INT(a)));
874    return INT_TO_SR( g);
875  }
876  else if (n_Z_IS_SMALL(b))
877  {
878    if (b==INT_TO_SR(0))
879      return nrzCopy(a, R);
880    unsigned long g = mpz_gcd_ui(NULL, (mpz_ptr)a, (unsigned long) ABS(SR_TO_INT(b)));
881    return INT_TO_SR(g);
882  }
883  else
884  {
885    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
886    mpz_init(erg);
887    mpz_gcd(erg, (mpz_ptr) a, (mpz_ptr) b);
888    return (number) erg;
889  }
890}
891
892/*
893 * Give the largest non unit k, such that a = x * k, b = y * k has
894 * a solution and r, s, s.t. k = s*a + t*b
895 */
896static long int_extgcd(long a, long b, long * u, long* x, long * v, long* y)
897{
898  long q, r;
899  if (!a)
900  {
901    *u = 0;
902    *v = 1;
903    *x = -1;
904    *y = 0;
905    return b;
906  }
907  if (!b)
908  {
909    *u = 1;
910    *v = 0;
911    *x = 0;
912    *y = 1;
913    return a;
914  }
915  *u=1;
916  *v=0;
917  *x=0;
918  *y=1;
919  do
920  {
921    q = a/b;
922    r = a%b;
923    assume (q*b+r == a);
924    a = b;
925    b = r;
926
927    r = -(*v)*q+(*u);
928    (*u) =(*v);
929    (*v) = r;
930
931    r = -(*y)*q+(*x);
932    (*x) = (*y);
933    (*y) = r;
934  } while (b);
935
936  return a;
937}
938
939static number  nrzExtGcd (number a, number b, number *s, number *t, const coeffs)
940{
941  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
942  {
943    long u, v, x, y;
944    long g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &u, &v, &x, &y);
945    *s = INT_TO_SR(u);
946    *t = INT_TO_SR(v);
947    return INT_TO_SR(g);
948  }
949  else
950  {
951    mpz_t aa, bb;
952    if (n_Z_IS_SMALL(a))
953    {
954      mpz_init_set_si(aa, SR_TO_INT(a));
955    }
956    else
957    {
958      mpz_init_set(aa, (mpz_ptr) a);
959    }
960    if (n_Z_IS_SMALL(b))
961    {
962      mpz_init_set_si(bb, SR_TO_INT(b));
963    }
964    else
965    {
966      mpz_init_set(bb, (mpz_ptr) b);
967    }
968    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
969    mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
970    mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
971    mpz_init(erg);
972    mpz_init(bs);
973    mpz_init(bt);
974    mpz_gcdext(erg, bs, bt, aa, bb);
975    *s = nrz_short((number) bs);
976    *t = nrz_short((number) bt);
977    mpz_clear(aa);
978    mpz_clear(bb);
979    return nrz_short((number) erg);
980  }
981}
982#if CF_DEBUG
983static number _nrzXExtGcd(number, number, number *, number *, number *, number *, const coeffs);
984static number nrzXExtGcd(number a, number b, number *x, number * y, number * u, number * v, const coeffs R)
985{
986  char * s;
987  StringSetS("XExtGcd: ");
988  nrzWrite(a, R);
989  StringAppendS(" by ");
990  nrzWrite(b, R);
991  number c = _nrzXExtGcd(a, b, x, y, u, v, R);
992  StringAppendS(" is ");
993  nrzWrite(c, R);
994  StringAppendS("[[");
995  nrzWrite(*x, R);
996  StringAppendS(", ");
997  nrzWrite(*y, R);
998  StringAppendS("], ");
999  nrzWrite(*u, R);
1000  StringAppendS(", ");
1001  nrzWrite(*v, R);
1002  s=StringEndS();
1003  Print("%s]]\n", s);
1004  omFree(s);
1005  return c;
1006}
1007static number  _nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )
1008#else
1009static number  nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )
1010#endif
1011{
1012  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1013  {
1014    long uu, vv, x, y;
1015    long g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &uu, &vv, &x, &y);
1016    *s = INT_TO_SR(uu);
1017    *t = INT_TO_SR(vv);
1018    *u = INT_TO_SR(x);
1019    *v = INT_TO_SR(y);
1020    return INT_TO_SR(g);
1021  }
1022  else
1023  {
1024    mpz_t aa, bb;
1025    if (n_Z_IS_SMALL(a))
1026    {
1027      mpz_init_set_si(aa, SR_TO_INT(a));
1028    }
1029    else
1030    {
1031      mpz_init_set(aa, (mpz_ptr) a);
1032    }
1033    if (n_Z_IS_SMALL(b))
1034    {
1035      mpz_init_set_si(bb, SR_TO_INT(b));
1036    }
1037    else
1038    {
1039      mpz_init_set(bb, (mpz_ptr) b);
1040    }
1041    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1042    mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1043    mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1044    mpz_init(erg);
1045    mpz_init(bs);
1046    mpz_init(bt);
1047
1048    mpz_gcdext(erg, bs, bt, aa, bb);
1049
1050    mpz_ptr bu = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1051    mpz_ptr bv = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1052
1053    mpz_init_set(bu, (mpz_ptr) bb);
1054    mpz_init_set(bv, (mpz_ptr) aa);
1055
1056    mpz_clear(aa);
1057    mpz_clear(bb);
1058    assume(mpz_cmp_si(erg, 0));
1059
1060    mpz_div(bu, bu, erg);
1061    mpz_div(bv, bv, erg);
1062
1063    mpz_mul_si(bu, bu, -1);
1064    *u = nrz_short((number) bu);
1065    *v = nrz_short((number) bv);
1066
1067    *s = nrz_short((number) bs);
1068    *t = nrz_short((number) bt);
1069    return nrz_short((number) erg);
1070  }
1071}
1072#if CF_DEBUG
1073static number _nrzQuotRem(number, number, number *, const coeffs);
1074static number nrzQuotRem(number a, number b, number * r, const coeffs R)
1075{
1076  StringSetS("QuotRem: ");
1077  nrzWrite(a, R);
1078  StringAppendS(" by ");
1079  nrzWrite(b, R);
1080  number c = _nrzQuotRem(a, b, r, R);
1081  StringAppendS(" is ");
1082  nrzWrite(c, R);
1083  if (r) {
1084    StringAppendS("+R(");
1085    nrzWrite(*r, R);
1086    StringAppendS(")");
1087  }
1088  char * s = StringEndS();
1089  Print("%s\n", s);
1090  omFree(s);
1091  return c;
1092}
1093static number _nrzQuotRem (number a, number b, number * r, const coeffs )
1094#else
1095static number nrzQuotRem (number a, number b, number * r, const coeffs )
1096#endif
1097{
1098  assume(SR_TO_INT(b));
1099  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1100  {
1101    if (r)
1102      *r = INT_TO_SR(SR_TO_INT(a) % SR_TO_INT(b));
1103    return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
1104  }
1105  else if (n_Z_IS_SMALL(a))
1106  {
1107    //a is small, b is not, so q=0, r=a
1108    if (r)
1109      *r = a;
1110    return INT_TO_SR(0);
1111  }
1112  else if (n_Z_IS_SMALL(b))
1113  {
1114    unsigned long rr;
1115    mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1116    mpz_init(qq);
1117    mpz_t rrr;
1118    mpz_init(rrr);
1119    rr = mpz_divmod_ui(qq, rrr, (mpz_ptr) a, (unsigned long)ABS(SR_TO_INT(b)));
1120    mpz_clear(rrr);
1121
1122    if (r)
1123      *r = INT_TO_SR(rr);
1124    if (SR_TO_INT(b)<0)
1125    {
1126      mpz_mul_si(qq, qq, -1);
1127    }
1128    return nrz_short((number)qq);
1129  }
1130  mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin),
1131             rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1132  mpz_init(qq);
1133  mpz_init(rr);
1134  mpz_divmod(qq, rr, (mpz_ptr)a, (mpz_ptr)b);
1135  if (r)
1136    *r = (number) rr;
1137  else
1138  {
1139    mpz_clear(rr);
1140  }
1141  nrzTest((number)qq);
1142  return (number) qq;
1143}
1144
1145static void nrzPower (number a, int i, number * result, const coeffs)
1146{
1147  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1148  mpz_init(erg);
1149  mpz_t aa;
1150  if (n_Z_IS_SMALL(a))
1151    mpz_init_set_si(aa, SR_TO_INT(a));
1152  else
1153    mpz_init_set(aa, (mpz_ptr) a);
1154  mpz_pow_ui(erg, aa, i);
1155  *result = nrz_short((number) erg);
1156}
1157
1158/*
1159 * create a number from int
1160 * TODO: do not create an mpz if not necessary
1161 */
1162number nrzInit (long i, const coeffs)
1163{
1164  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1165  mpz_init_set_si(erg, i);
1166  return nrz_short((number) erg);
1167}
1168
1169static number nrzInitMPZ(mpz_t m, const coeffs)
1170{
1171  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1172  mpz_init_set(erg, m);
1173  return nrz_short((number) erg);
1174}
1175
1176
1177static void nrzDelete(number *a, const coeffs)
1178{
1179  if (*a == NULL) return;
1180  if (n_Z_IS_SMALL(*a)==0)
1181  {
1182    mpz_clear((mpz_ptr) *a);
1183    omFreeBin((ADDRESS) *a, gmp_nrz_bin);
1184  }
1185  *a = NULL;
1186}
1187
1188/*
1189 * convert a number to int
1190 */
1191static long nrzInt(number &n, const coeffs)
1192{
1193  if (n_Z_IS_SMALL(n)) return SR_TO_INT(n);
1194  return mpz_get_si( (mpz_ptr)n);
1195}
1196#if CF_DEBUG
1197static number _nrzAdd(number, number, const coeffs);
1198static number nrzAdd(number a, number b, const coeffs R)
1199{
1200  StringSetS("Add: ");
1201  nrzWrite(a, R);
1202  StringAppendS(" to ");
1203  nrzWrite(b, R);
1204  number c = _nrzAdd(a, b, R);
1205  StringAppendS(" is ");
1206  nrzWrite(c, R);
1207  char * s = StringEndS();
1208  Print("%s\n", s);
1209  omFree(s);
1210  return c;
1211}
1212static number _nrzAdd (number a, number b, const coeffs )
1213#else
1214static number nrzAdd (number a, number b, const coeffs )
1215#endif
1216{
1217  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1218  {
1219    long c = SR_TO_INT(a) + SR_TO_INT(b);
1220    if (INT_IS_SMALL(c))
1221      return INT_TO_SR(c);
1222    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1223    mpz_init_set_si(erg, c);
1224
1225    nrzTest((number)erg);
1226    return (number) erg;
1227  }
1228  else if (n_Z_IS_SMALL(a))
1229  {
1230    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1231    mpz_init(erg);
1232    if (SR_TO_INT(a)>0)
1233      mpz_add_ui(erg, (mpz_ptr) b, (unsigned long)SR_TO_INT(a));
1234    else
1235      mpz_sub_ui(erg, (mpz_ptr) b, (unsigned long)-(SR_TO_INT(a)));
1236    return nrz_short((number) erg);
1237  }
1238  else if (n_Z_IS_SMALL(b))
1239  {
1240    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1241    mpz_init(erg);
1242    if (SR_TO_INT(b)>0)
1243      mpz_add_ui(erg, (mpz_ptr) a, (unsigned long)SR_TO_INT(b));
1244    else
1245      mpz_sub_ui(erg, (mpz_ptr) a, (unsigned long)-(SR_TO_INT(b)));
1246    return nrz_short((number) erg);
1247  }
1248  else
1249  {
1250    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1251    mpz_init(erg);
1252    mpz_add(erg, (mpz_ptr) a, (mpz_ptr) b);
1253    return nrz_short((number) erg);
1254  }
1255}
1256
1257static number nrzSub (number a, number b,  const coeffs )
1258{
1259  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1260  {
1261    long c = SR_TO_INT(a) - SR_TO_INT(b);
1262    if (INT_IS_SMALL(c))
1263      return INT_TO_SR(c);
1264    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1265    mpz_init_set_si(erg, c);
1266    nrzTest((number)erg);
1267    return (number) erg;
1268  }
1269  else if (n_Z_IS_SMALL(a))
1270  {
1271    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1272    mpz_init(erg);
1273
1274    if (SR_TO_INT(a)>0)
1275      mpz_ui_sub(erg, (unsigned long)SR_TO_INT(a), (mpz_ptr) b);
1276    else
1277    {
1278      mpz_add_ui(erg, (mpz_ptr) b, (unsigned long)-SR_TO_INT(a));
1279      mpz_neg(erg, erg);
1280    }
1281    return nrz_short((number) erg);
1282  }
1283  else if (n_Z_IS_SMALL(b))
1284  {
1285    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1286    mpz_init(erg);
1287    if (SR_TO_INT(b)>0)
1288      mpz_sub_ui(erg, (mpz_ptr) a, (unsigned long)SR_TO_INT(b));
1289    else
1290      mpz_add_ui(erg, (mpz_ptr) a, (unsigned long)-SR_TO_INT(b));
1291    return nrz_short((number) erg);
1292  }
1293  else
1294  {
1295    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1296    mpz_init(erg);
1297    mpz_sub(erg, (mpz_ptr) a, (mpz_ptr) b);
1298    return nrz_short((number) erg);
1299  }
1300}
1301
1302static BOOLEAN nrzGreater (number a,number b, const coeffs)
1303{
1304  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1305    return ((long)a)>((long)b);
1306  else if (n_Z_IS_SMALL(a))
1307    return 0 > mpz_cmp_si((mpz_ptr)b,SR_TO_INT(a));
1308  else if (n_Z_IS_SMALL(b))
1309    return 0 < mpz_cmp_si((mpz_ptr)a,SR_TO_INT(b));
1310  return 0 < mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
1311}
1312
1313static BOOLEAN nrzGreaterZero (number k, const coeffs C)
1314{
1315  return nrzGreater(k, INT_TO_SR(0), C);
1316}
1317
1318static number  nrzGetUnit (number n, const coeffs r)
1319{
1320  if (nrzGreaterZero(n, r))
1321    return INT_TO_SR(1);
1322  /*else*/
1323    return INT_TO_SR(-1);
1324}
1325
1326static number nrzAnn(number n, const coeffs)
1327{
1328  if (SR_TO_INT(n))  // in Z: the annihilator of !=0 is 0
1329    return INT_TO_SR(0);
1330  else
1331    return INT_TO_SR(1);
1332}
1333
1334static BOOLEAN nrzIsUnit (number a, const coeffs)
1335{
1336  return ABS(SR_TO_INT(a))==1;
1337}
1338
1339static BOOLEAN nrzIsZero (number  a, const coeffs)
1340{
1341  return (a==INT_TO_SR(0));
1342}
1343
1344static BOOLEAN nrzIsOne (number a, const coeffs)
1345{
1346  return a==INT_TO_SR(1);
1347}
1348
1349static BOOLEAN nrzIsMOne (number a, const coeffs)
1350{
1351  return a==INT_TO_SR(-1);
1352}
1353
1354static BOOLEAN nrzEqual (number a,number b, const coeffs)
1355{
1356  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1357    return a==b;
1358  else if (n_Z_IS_SMALL(a) || n_Z_IS_SMALL(b))
1359    return FALSE;
1360  else
1361    return 0 == mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
1362}
1363
1364static BOOLEAN nrzDivBy (number a,number b, const coeffs)
1365{
1366  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1367  {
1368    return SR_TO_INT(a) %SR_TO_INT(b) ==0;
1369  }
1370  else if (n_Z_IS_SMALL(a))
1371  {
1372    return a==INT_TO_SR(0);
1373  }
1374  else if (n_Z_IS_SMALL(b))
1375  {
1376    return mpz_divisible_ui_p((mpz_ptr)a, (unsigned long)ABS(SR_TO_INT(b))) != 0;
1377  }
1378  else
1379    return mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b) != 0;
1380}
1381
1382static int nrzDivComp(number a, number b, const coeffs r)
1383{
1384  if (nrzDivBy(a, b, r))
1385  {
1386    if (nrzDivBy(b, a, r)) return 2;
1387    return -1;
1388  }
1389  if (nrzDivBy(b, a, r)) return 1;
1390  return 0;
1391}
1392
1393static number nrzDiv (number a,number b, const coeffs)
1394{
1395  assume(SR_TO_INT(b));
1396  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1397  {
1398    //if (SR_TO_INT(a) % SR_TO_INT(b))
1399    //{
1400    //  WerrorS("1:Division by non divisible element.");
1401    //  WerrorS("Result is without remainder.");
1402    //}
1403    return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
1404  }
1405  else if (n_Z_IS_SMALL(a))
1406  {
1407    //if (SR_TO_INT(a))
1408    //{
1409    //  WerrorS("2:Division by non divisible element.");
1410    //  WerrorS("Result is without remainder.");
1411    //}
1412    return INT_TO_SR(0);
1413  }
1414  else if (n_Z_IS_SMALL(b))
1415  {
1416    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1417    mpz_t r;
1418    mpz_init(r);
1419    mpz_init(erg);
1420    if (mpz_divmod_ui(erg, r, (mpz_ptr) a, (unsigned long)ABS(SR_TO_INT(b)))) {
1421    //  WerrorS("3:Division by non divisible element.");
1422    //  WerrorS("Result is without remainder.");
1423    }
1424    mpz_clear(r);
1425    if (SR_TO_INT(b)<0)
1426      mpz_neg(erg, erg);
1427    return nrz_short((number) erg);
1428  }
1429  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1430  mpz_init(erg);
1431  mpz_t r;
1432  mpz_init(r);
1433  mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);
1434#if CF_DEBUG
1435  StringSetS("division of");
1436  nrzWrite(a, R);
1437  StringAppendS(" by ");
1438  nrzWrite(b, R);
1439  StringAppendS(" is ");
1440  number du;
1441  nrzWrite(du = (number)erg, R);
1442  StringAppendS(" rest ");
1443  nrzWrite(du = (number)r, R);
1444  char * s = StringEndS();
1445  Print("%s\n", s);
1446  omFree(s);
1447#endif
1448
1449  if (mpz_sgn1(r)!=0)
1450  {
1451    //WerrorS("4:Division by non divisible element.");
1452    //WerrorS("Result is without remainder.");
1453  }
1454  mpz_clear(r);
1455  return nrz_short((number) erg);
1456}
1457
1458static number nrzExactDiv (number a,number b, const coeffs)
1459{
1460  assume(SR_TO_INT(b));
1461  mpz_t aa, bb;
1462  if (n_Z_IS_SMALL(a))
1463    mpz_init_set_si(aa, SR_TO_INT(a));
1464  else
1465    mpz_init_set(aa, (mpz_ptr) a);
1466  if (n_Z_IS_SMALL(b))
1467    mpz_init_set_si(bb, SR_TO_INT(b));
1468  else
1469    mpz_init_set(bb, (mpz_ptr) b);
1470  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1471  mpz_init(erg);
1472  mpz_tdiv_q(erg, (mpz_ptr) aa, (mpz_ptr) bb);
1473  mpz_clear(aa);
1474  mpz_clear(bb);
1475  nrzTest((number)erg);
1476  return (number) erg;
1477}
1478
1479static number nrzIntMod (number a,number b, const coeffs)
1480{
1481  mpz_t aa, bb;
1482  assume(SR_TO_INT(b));
1483  if (n_Z_IS_SMALL(a))
1484    mpz_init_set_si(aa, SR_TO_INT(a));
1485  else
1486    mpz_init_set(aa, (mpz_ptr) a);
1487  if (n_Z_IS_SMALL(b))
1488    mpz_init_set_si(bb, SR_TO_INT(b));
1489  else
1490    mpz_init_set(bb, (mpz_ptr) b);
1491
1492  mpz_t erg;
1493  mpz_init(erg);
1494  mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1495  mpz_init(r);
1496  mpz_tdiv_qr(erg, r, (mpz_ptr) aa, (mpz_ptr) bb);
1497  mpz_clear(erg);
1498  mpz_clear(aa);
1499  mpz_clear(bb);
1500
1501  return nrz_short((number) r);
1502}
1503
1504static number  nrzInvers (number c, const coeffs r)
1505{
1506  if (!nrzIsUnit((number) c, r))
1507  {
1508    WerrorS("Non invertible element.");
1509    return (number)NULL;
1510  }
1511  return c; // has to be 1 or -1....
1512}
1513
1514static number nrzNeg (number c, const coeffs)
1515{
1516// nNeg inplace !!!
1517  if (n_Z_IS_SMALL(c))
1518    return INT_TO_SR(-SR_TO_INT(c));
1519  mpz_mul_si((mpz_ptr) c, (mpz_ptr) c, -1);
1520  return c;
1521}
1522
1523static number nrzFarey(number r, number N, const coeffs R)
1524{
1525  number a0 = nrzCopy(N, R);
1526  number b0 = nrzInit(0, R);
1527  number a1 = nrzCopy(r, R);
1528  number b1 = nrzInit(1, R);
1529  number two = nrzInit(2, R);
1530#if 0
1531  PrintS("Farey start with ");
1532  n_Print(r, R);
1533  PrintS(" mod ");
1534  n_Print(N, R);
1535  PrintLn();
1536#endif
1537  while (1)
1538  {
1539    number as = nrzMult(a1, a1, R);
1540    n_InpMult(as, two, R);
1541    if (nrzGreater(N, as, R))
1542    {
1543      nrzDelete(&as, R);
1544      break;
1545    }
1546    nrzDelete(&as, R);
1547    number q = nrzDiv(a0, a1, R);
1548    number t = nrzMult(a1, q, R),
1549           s = nrzSub(a0, t, R);
1550    nrzDelete(&a0, R);
1551    a0 = a1;
1552    a1 = s;
1553    nrzDelete(&t, R);
1554
1555    t = nrzMult(b1, q, R);
1556    s = nrzSub(b0, t, R);
1557    nrzDelete(&b0, R);
1558    b0 = b1;
1559    b1 = s;
1560    nrzDelete(&t, R);
1561    nrzDelete(&q, R);
1562  }
1563  number as = nrzMult(b1, b1, R);
1564  n_InpMult(as, two, R);
1565  nrzDelete(&two, R);
1566  if (nrzGreater(as, N, R))
1567  {
1568    nrzDelete(&a0, R);
1569    nrzDelete(&a1, R);
1570    nrzDelete(&b0, R);
1571    nrzDelete(&b1, R);
1572    nrzDelete(&as, R);
1573    return NULL;
1574  }
1575  nrzDelete(&as, R);
1576  nrzDelete(&a0, R);
1577  nrzDelete(&b0, R);
1578
1579  number a, b, ab;
1580  coeffs Q = nInitChar(n_Q, 0);
1581  nMapFunc f = n_SetMap(R, Q);
1582  a = f(a1, R, Q);
1583  b = f(b1, R, Q);
1584  ab = n_Div(a, b, Q);
1585  n_Delete(&a, Q);
1586  n_Delete(&b, Q);
1587  nKillChar(Q);
1588
1589  nrzDelete(&a1, R);
1590  nrzDelete(&b1, R);
1591  return ab;
1592}
1593
1594static number nrzMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
1595{
1596  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1597  mpz_init_set_ui(erg, (unsigned long) from);
1598  return nrz_short((number) erg);
1599}
1600
1601static number nrzMapZp(number from, const coeffs /*src*/, const coeffs /*dst*/)
1602{
1603  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1604  mpz_init_set_si(erg, (long) from);
1605  return nrz_short((number) erg);
1606}
1607
1608static number nrzModNMap(number from, const coeffs /* src */, const coeffs /*dst*/)
1609{
1610  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1611  mpz_init_set(erg, (mpz_ptr) from);
1612  return nrz_short((number) erg);
1613}
1614
1615static number nrzMapQ(number from, const coeffs /* src */, const coeffs dst)
1616{
1617  if (SR_HDL(from) & SR_INT)
1618    return nrzInit(SR_TO_INT(from),dst);
1619  if (from->s!=3)
1620  {
1621    WerrorS("rational in map to integer");
1622    return NULL;
1623  }
1624  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1625  mpz_init_set(erg, from->z);
1626  return nrz_short((number) erg);
1627}
1628
1629static nMapFunc nrzSetMap(const coeffs src, const coeffs /*dst*/)
1630{
1631  /* dst = rintegers */
1632  if (src->rep==n_rep_gmp) //nCoeff_is_Zn(src) || nCoeff_is_Ring_PtoM(src))
1633    return nrzModNMap;
1634
1635  if ((src->rep==n_rep_gap_gmp) && nCoeff_is_Z(src))
1636  {
1637    return ndCopyMap; //nrzCopyMap;
1638  }
1639  if (src->rep==n_rep_gap_rat) /*&& nCoeff_is_Z(src)) Q, bigint*/
1640  {
1641    return nrzMapQ;
1642  }
1643  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src))
1644  {
1645    return nrzMapMachineInt;
1646  }
1647  if ((src->rep==n_rep_int) && nCoeff_is_Zp(src))
1648  {
1649    return nrzMapZp;
1650  }
1651  return NULL;      // default
1652}
1653
1654
1655/*
1656 * set the exponent (allocate and init tables) (TODO)
1657 */
1658
1659void nrzSetExp(int, coeffs)
1660{
1661}
1662
1663void nrzInitExp(int, coeffs)
1664{
1665}
1666
1667#ifdef LDEBUG
1668BOOLEAN nrzDBTest (number x, const char *f, const int l, const coeffs)
1669{
1670  if (SR_HDL(x) & SR_INT) return TRUE;
1671  if (mpz_sgn1((mpz_ptr) x)==0)
1672  {
1673    Print("gmp-0 %s:%d\n",f,l);
1674    return FALSE;
1675  }
1676  if (mpz_size1((mpz_ptr)x)<=MP_SMALL)
1677  {
1678    long ui=mpz_get_si((mpz_ptr)x);
1679    if ((((ui<<3)>>3)==ui)
1680    && (mpz_cmp_si((mpz_ptr)x,ui)==0))
1681    {
1682      Print("gmp-small %s:%d\n",f,l);
1683      return FALSE;
1684    }
1685  }
1686  return TRUE;
1687}
1688#endif
1689
1690void nrzWrite (number a, const coeffs)
1691{
1692  char *s,*z;
1693  if (a==NULL)
1694  {
1695    StringAppendS("o");
1696  }
1697  else
1698  {
1699    if (n_Z_IS_SMALL(a))
1700    {
1701      StringAppend("%d", SR_TO_INT(a));
1702    }
1703    else
1704    {
1705      int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2;
1706      s=(char*)omAlloc(l);
1707      z=mpz_get_str(s,10,(mpz_ptr) a);
1708      StringAppendS(z);
1709      omFreeSize((ADDRESS)s,l);
1710    }
1711  }
1712}
1713
1714/*2
1715* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
1716*/
1717static const char * nlEatLongC(char *s, mpz_ptr i)
1718{
1719  const char * start=s;
1720
1721  if (*s<'0' || *s>'9')
1722  {
1723    mpz_set_ui(i,1);
1724    return s;
1725  }
1726  while (*s >= '0' && *s <= '9') s++;
1727  if (*s=='\0')
1728  {
1729    mpz_set_str(i,start,10);
1730  }
1731  else
1732  {
1733    char c=*s;
1734    *s='\0';
1735    mpz_set_str(i,start,10);
1736    *s=c;
1737  }
1738  return s;
1739}
1740
1741static const char * nrzRead (const char *s, number *a, const coeffs)
1742{
1743  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1744  {
1745    mpz_init(z);
1746    s = nlEatLongC((char *) s, z);
1747  }
1748  *a = nrz_short((number) z);
1749  return s;
1750}
1751
1752static void nrzCoeffWrite  (const coeffs, BOOLEAN /*details*/)
1753{
1754  //PrintS("// ZZ\n");
1755  PrintS("//   coeff. ring is : Integers\n");
1756}
1757
1758static char* nrzCoeffString(const coeffs)
1759{
1760  return omStrDup("integer");
1761}
1762
1763static CanonicalForm nrzConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs /*r*/ )
1764{
1765  if (setChar) setCharacteristic( 0 );
1766
1767  CanonicalForm term;
1768  if ( n_Z_IS_SMALL(n))
1769  {
1770    term = SR_TO_INT(n);
1771  }
1772  else
1773  {
1774    mpz_t dummy;
1775    mpz_init_set( dummy,n->z );
1776    term = make_cf( dummy );
1777  }
1778  return term;
1779}
1780
1781static number nrzConvFactoryNSingN( const CanonicalForm n, const coeffs r)
1782{
1783  if (n.isImm())
1784  {
1785    return nrzInit(n.intval(),r);
1786  }
1787  else
1788  {
1789    if ( !n.den().isOne() )
1790    {
1791     WerrorS("rational in conversion to integer");
1792     return NULL;
1793    }
1794    mpz_ptr z = (mpz_ptr) omAlloc0Bin(gmp_nrz_bin);
1795    gmp_numerator( n,z);
1796    return nrz_short((number)z);
1797  }
1798}
1799
1800static void nrzMPZ(mpz_t res, number &a, const coeffs)
1801{
1802  if (n_Z_IS_SMALL(a))
1803    mpz_init_set_si(res, SR_TO_INT(a));
1804  else
1805    mpz_init_set(res, (mpz_ptr) a);
1806}
1807
1808static coeffs nrzQuot1(number c, const coeffs r)
1809{
1810    mpz_t dummy;
1811    if(n_Z_IS_SMALL(c))
1812    {
1813      long ch = r->cfInt(c, r);
1814      mpz_init_set_ui(dummy, ch);
1815    }
1816    else
1817    {
1818      mpz_init_set(dummy, (mpz_ptr)c);
1819    }
1820    ZnmInfo info;
1821    info.base = dummy;
1822    info.exp = (unsigned long) 1;
1823    coeffs rr = nInitChar(n_Zn, (void*)&info);
1824    mpz_clear(dummy);
1825    return(rr);
1826}
1827
1828number nlReadFd(const ssiInfo *d, const coeffs);
1829void nlWriteFd(number n, const ssiInfo* d, const coeffs);
1830
1831BOOLEAN nrzInitChar(coeffs r,  void *)
1832{
1833  assume( getCoeffType(r) == n_Z );
1834
1835  r->is_field=FALSE;
1836  r->is_domain=TRUE;
1837  r->rep=n_rep_gap_gmp;
1838
1839  //r->nCoeffIsEqual = ndCoeffIsEqual;
1840  r->cfCoeffString = nrzCoeffString;
1841  //r->cfKillChar = ndKillChar;
1842  r->cfMult  = nrzMult;
1843  r->cfSub   = nrzSub;
1844  r->cfAdd   = nrzAdd;
1845  r->cfDiv   = nrzDiv;
1846  r->cfIntMod= nrzIntMod;
1847  r->cfExactDiv= nrzExactDiv;
1848  r->cfInit = nrzInit;
1849  r->cfInitMPZ = nrzInitMPZ;
1850  r->cfSize  = nrzSize;
1851  r->cfInt  = nrzInt;
1852  //#ifdef HAVE_RINGS
1853  r->cfDivComp = nrzDivComp; // only for ring stuff
1854  r->cfIsUnit = nrzIsUnit; // only for ring stuff
1855  r->cfGetUnit = nrzGetUnit; // only for ring stuff
1856  r->cfAnn = nrzAnn;
1857  r->cfExtGcd = nrzExtGcd; // only for ring stuff
1858  r->cfXExtGcd = nrzXExtGcd; // only for ring stuff
1859  r->cfQuotRem = nrzQuotRem;
1860  r->cfDivBy = nrzDivBy; // only for ring stuff
1861  //#endif
1862  r->cfInpNeg   = nrzNeg;
1863  r->cfInvers= nrzInvers;
1864  r->cfCopy  = nrzCopy;
1865  r->cfWriteLong = nrzWrite;
1866  r->cfRead = nrzRead;
1867  r->cfGreater = nrzGreater;
1868  r->cfEqual = nrzEqual;
1869  r->cfIsZero = nrzIsZero;
1870  r->cfIsOne = nrzIsOne;
1871  r->cfIsMOne = nrzIsMOne;
1872  r->cfGreaterZero = nrzGreaterZero;
1873  r->cfPower = nrzPower;
1874  r->cfGcd  = nrzGcd;
1875  r->cfLcm  = nrzLcm;
1876  r->cfDelete= nrzDelete;
1877  r->cfSetMap = nrzSetMap;
1878  r->cfCoeffWrite = nrzCoeffWrite;
1879  r->convSingNFactoryN = nrzConvSingNFactoryN;
1880  r->convFactoryNSingN = nrzConvFactoryNSingN;
1881  r->cfMPZ = nrzMPZ;
1882  r->cfFarey = nrzFarey;
1883  r->cfWriteFd=nlWriteFd;
1884  r->cfReadFd=nlReadFd;
1885
1886  r->cfQuot1 = nrzQuot1;
1887  // requires conversion to factory:
1888  r->cfChineseRemainder=nlChineseRemainderSym;
1889  // debug stuff
1890
1891#ifdef LDEBUG
1892  r->cfDBTest=nrzDBTest;
1893#endif
1894
1895  r->ch = 0;
1896  r->has_simple_Alloc=FALSE;
1897  r->has_simple_Inverse=FALSE;
1898  return FALSE;
1899}
1900
1901#elif SI_INTEGER_VARIANT == 1
1902BOOLEAN nrzInitChar(coeffs r,  void *)
1903{
1904  return nlInitChar(r,(void*)1);
1905}
1906#else
1907#error set SI_INTEGER_VARIANT
1908#endif
1909#endif
Note: See TracBrowser for help on using the repository browser.