source: git/libpolys/coeffs/rintegers.cc @ 353a42

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