source: git/libpolys/coeffs/rintegers.cc @ 4b5b36

spielwiese
Last change on this file since 4b5b36 was 4b5b36, checked in by Hans Schoenemann <hannes@…>, 5 years ago
add: ssi for poly over Z/n
  • Property mode set to 100644
File size: 40.7 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_Ring_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_Ring_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    return (number) m;
456  }
457}
458
459static const char * nrzRead (const char *s, number *a, const coeffs)
460{
461  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
462  {
463    mpz_init(z);
464    s = nlEatLongC((char *) s, z);
465  }
466  *a = (number) z;
467  return s;
468}
469
470static void nrzCoeffWrite  (const coeffs, BOOLEAN /*details*/)
471{
472  PrintS("ZZ");
473}
474
475static char* nrzCoeffName(const coeffs)
476{
477  return (char*)"ZZ";
478}
479
480static char* nrzCoeffString(const coeffs cf)
481{
482  return omStrDup(nrzCoeffName(cf));
483}
484
485static coeffs nrzQuot1(number c, const coeffs r)
486{
487    long ch = r->cfInt(c, r);
488    mpz_t dummy;
489    mpz_init_set_ui(dummy, ch);
490    ZnmInfo info;
491    info.base = dummy;
492    info.exp = (unsigned long) 1;
493    coeffs rr = nInitChar(n_Zn, (void*)&info);
494    mpz_clear(dummy);
495    return(rr);
496}
497
498static number nrzInitMPZ(mpz_t m, const coeffs)
499{
500  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
501  mpz_init_set(z, m);
502  return (number)z;
503}
504
505static number nrzFarey(number r, number N, const coeffs R)
506{
507  number a0 = nrzCopy(N, R);
508  number b0 = nrzInit(0, R);
509  number a1 = nrzCopy(r, R);
510  number b1 = nrzInit(1, R);
511  number two = nrzInit(2, R);
512#if 0
513  PrintS("Farey start with ");
514  n_Print(r, R);
515  PrintS(" mod ");
516  n_Print(N, R);
517  PrintLn();
518#endif
519  while (1)
520  {
521    number as = nrzMult(a1, a1, R);
522    n_InpMult(as, two, R);
523    if (nrzGreater(N, as, R))
524    {
525      nrzDelete(&as, R);
526      break;
527    }
528    nrzDelete(&as, R);
529    number q = nrzDiv(a0, a1, R);
530    number t = nrzMult(a1, q, R),
531           s = nrzSub(a0, t, R);
532    nrzDelete(&a0, R);
533    a0 = a1;
534    a1 = s;
535    nrzDelete(&t, R);
536
537    t = nrzMult(b1, q, R);
538    s = nrzSub(b0, t, R);
539    nrzDelete(&b0, R);
540    b0 = b1;
541    b1 = s;
542    nrzDelete(&t, R);
543    nrzDelete(&q, R);
544  }
545  number as = nrzMult(b1, b1, R);
546  n_InpMult(as, two, R);
547  nrzDelete(&two, R);
548  if (nrzGreater(as, N, R))
549  {
550    nrzDelete(&a0, R);
551    nrzDelete(&a1, R);
552    nrzDelete(&b0, R);
553    nrzDelete(&b1, R);
554    nrzDelete(&as, R);
555    return NULL;
556  }
557  nrzDelete(&as, R);
558  nrzDelete(&a0, R);
559  nrzDelete(&b0, R);
560
561  number a, b, ab;
562  coeffs Q = nInitChar(n_Q, 0);
563  nMapFunc f = n_SetMap(R, Q);
564  a = f(a1, R, Q);
565  b = f(b1, R, Q);
566  ab = n_Div(a, b, Q);
567  n_Delete(&a, Q);
568  n_Delete(&b, Q);
569  nKillChar(Q);
570
571  nrzDelete(&a1, R);
572  nrzDelete(&b1, R);
573  return ab;
574}
575
576void nrzWriteFd(number n, const ssiInfo* d, const coeffs)
577{
578  mpz_out_str (d->f_write,SSI_BASE, (mpz_ptr)n);
579  fputc(' ',d->f_write);
580}
581
582number nrzReadFd(const ssiInfo *d, const coeffs)
583{
584  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
585  mpz_init(erg);
586  s_readmpz_base(d->f_read,erg,SSI_BASE);
587  return (number)erg;
588}
589
590BOOLEAN nrzInitChar(coeffs r,  void *)
591{
592  assume( getCoeffType(r) == n_Z );
593
594  r->is_field=FALSE;
595  r->is_domain=TRUE;
596  r->rep=n_rep_gmp;
597
598  //r->nCoeffIsEqual = ndCoeffIsEqual;
599  r->cfCoeffString = nrzCoeffString;
600  r->cfCoeffName = nrzCoeffName;
601  r->cfCoeffWrite = nrzCoeffWrite;
602  //r->cfKillChar = ndKillChar;
603  r->cfMult  = nrzMult;
604  r->cfSub   = nrzSub;
605  r->cfAdd   = nrzAdd;
606  r->cfDiv   = nrzDiv;
607  r->cfIntMod= nrzIntMod;
608  r->cfExactDiv= nrzExactDiv;
609  r->cfInit = nrzInit;
610  r->cfInitMPZ = nrzInitMPZ;
611  r->cfSize  = nrzSize;
612  r->cfInt  = nrzInt;
613  r->cfDivComp = nrzDivComp;
614  r->cfIsUnit = nrzIsUnit;
615  r->cfGetUnit = nrzGetUnit;
616  r->cfExtGcd = nrzExtGcd;
617  r->cfXExtGcd = nrzXExtGcd;
618  r->cfDivBy = nrzDivBy;
619  r->cfQuotRem = nrzQuotRem;
620  r->cfInpNeg   = nrzNeg;
621  r->cfInvers= nrzInvers;
622  r->cfCopy  = nrzCopy;
623  r->cfWriteLong = nrzWrite;
624  r->cfRead = nrzRead;
625  r->cfGreater = nrzGreater;
626  r->cfEqual = nrzEqual;
627  r->cfIsZero = nrzIsZero;
628  r->cfIsOne = nrzIsOne;
629  r->cfIsMOne = nrzIsMOne;
630  r->cfGreaterZero = nrzGreaterZero;
631  r->cfPower = nrzPower;
632  r->cfGcd  = nrzGcd;
633  r->cfLcm  = nrzLcm;
634  r->cfDelete= nrzDelete;
635  r->cfSetMap = nrzSetMap;
636  r->cfQuot1 = nrzQuot1;
637  r->convSingNFactoryN=nrzConvSingNFactoryN;
638  r->convFactoryNSingN=nrzConvFactoryNSingN;
639  r->cfChineseRemainder=nlChineseRemainderSym;
640  r->cfFarey=nrzFarey;
641  r->cfWriteFd=nrzWriteFd;
642  r->cfReadFd=nrzReadFd;
643  // debug stuff
644
645#ifdef LDEBUG
646  r->cfDBTest=nrzDBTest;
647#endif
648
649  r->ch = 0;
650  r->has_simple_Alloc=FALSE;
651  r->has_simple_Inverse=FALSE;
652  return FALSE;
653}
654
655#elif SI_INTEGER_VARIANT == 3
656
657//make sure that a small number is an immediate integer
658//bascially coped from longrat.cc nlShort3
659//TODO: is there any point in checking 0 first???
660//TODO: it is not clear that this works in 32/64 bit everywhere.
661//      too many hacks.
662#ifdef LDEBUG
663#define nrzTest(A) nrzDBTest(A,__FILE__,__LINE__,NULL)
664BOOLEAN nrzDBTest (number x, const char *f, const int l, const coeffs);
665#else
666#define nrzTest(A)
667#endif
668
669#undef CF_DEBUG
670static inline number nrz_short(number x)
671{
672#if CF_DEBUG
673  StringAppendS("short(");
674  nrzWrite(x, NULL);
675#endif
676  if (mpz_sgn1((mpz_ptr) x)==0)
677  {
678    mpz_clear((mpz_ptr)x);
679    omFreeBin(x, gmp_nrz_bin);
680#if CF_DEBUG
681    StringAppendS(")=0");
682#endif
683    return INT_TO_SR(0);
684  }
685  if (mpz_size1((mpz_ptr)x)<=MP_SMALL)
686  {
687    long ui=mpz_get_si((mpz_ptr)x);
688    if ((((ui<<3)>>3)==ui)
689    && (mpz_cmp_si((mpz_ptr)x,ui)==0))
690    {
691      mpz_clear((mpz_ptr)x);
692      omFreeBin(x, gmp_nrz_bin);
693#if CF_DEBUG
694    StringAppendS(")=imm");
695#endif
696      return INT_TO_SR(ui);
697    }
698  }
699#if CF_DEBUG
700  StringAppendS(")");
701#endif
702  return x;
703}
704
705
706static int nrzSize(number a, const coeffs)
707{
708  if (a==INT_TO_SR(0)) return 0;
709  if (n_Z_IS_SMALL(a)) return 1;
710  return ((mpz_ptr)a)->_mp_alloc;
711}
712
713
714/*
715 * Multiply two numbers
716 * check for 0, 1, -1 maybe
717 */
718#if CF_DEBUG
719number _nrzMult(number, number, const coeffs);
720number nrzMult(number a, number b, const coeffs R)
721{
722  StringSetS("Mult: ");
723  nrzWrite(a, R);
724  StringAppendS(" by ");
725  nrzWrite(b, R);
726  number c = _nrzMult(a, b, R);
727  StringAppendS(" is ");
728  nrzWrite(c, R);
729  char * s = StringEndS();
730  Print("%s\n", s);
731  omFree(s);
732  return c;
733}
734number _nrzMult (number a, number b, const coeffs R)
735#else
736number nrzMult (number a, number b, const coeffs R)
737#endif
738{
739  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
740  {
741  //from longrat.cc
742    if (SR_TO_INT(a)==0)
743      return a;
744    if (SR_TO_INT(b)==0)
745      return b;
746    long r=(long)((unsigned long)(SR_HDL(a)-1L))*((unsigned long)(SR_HDL(b)>>1));
747    if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
748    {
749      number u=((number) ((r>>1)+SR_INT));
750    //  if (((((long)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
751      return nrzInit(SR_HDL(u)>>2, R);
752    }
753    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
754    mpz_init(erg);
755    mpz_set_si(erg, SR_TO_INT(a));
756    mpz_mul_si(erg, erg, SR_TO_INT(b));
757    nrzTest((number)erg);
758    return (number) erg;
759  }
760  else if (n_Z_IS_SMALL(a))
761  {
762    if (SR_TO_INT(a)==0)
763      return a;
764    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
765    mpz_init_set(erg, (mpz_ptr) b);
766    mpz_mul_si(erg, erg, SR_TO_INT(a));
767    nrzTest((number)erg);
768    return (number) erg;
769  }
770  else if (n_Z_IS_SMALL(b))
771  {
772    if (SR_TO_INT(b)==0)
773      return b;
774    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
775    mpz_init_set(erg, (mpz_ptr) a);
776    mpz_mul_si(erg, erg, SR_TO_INT(b));
777    nrzTest((number)erg);
778    return (number) erg;
779  }
780  else
781  {
782    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
783    mpz_init(erg);
784    mpz_mul(erg, (mpz_ptr) a, (mpz_ptr) b);
785    nrzTest((number)erg);
786    return (number) erg;
787  }
788}
789
790
791static long int_gcd(long a, long b)
792{
793  long r;
794  a = ABS(a);
795  b = ABS(b);
796  if (!a) return b;
797  if (!b) return a;
798  do
799  {
800    r = a % b;
801    a = b;
802    b = r;
803  } while (b);
804  return ABS(a); // % in c doeas not imply a signn
805                 // it would be unlikely to see a negative here
806                 // but who knows
807}
808
809/*
810 * Give the smallest non unit k, such that a * x = k = b * y has a solution
811 */
812static number nrzLcm (number a, number b, const coeffs R)
813{
814  #ifdef CF_DEBUG
815  PrintS("nrzLcm\n");
816  #endif
817  mpz_ptr erg;
818  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
819  {
820    long g = int_gcd(SR_TO_INT(a), SR_TO_INT(b));
821    return nrzMult(a, INT_TO_SR(SR_TO_INT(b)/g), R);
822  }
823  else
824  {
825    erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
826    if (n_Z_IS_SMALL(a))
827    {
828      mpz_init_set(erg, (mpz_ptr) b);
829      unsigned long g = mpz_gcd_ui(NULL, erg, (unsigned long) ABS(SR_TO_INT(a)));
830      mpz_mul_si(erg, erg, SR_TO_INT(a)/g);
831    }
832    else if (n_Z_IS_SMALL(b))
833    {
834      mpz_init_set(erg, (mpz_ptr) a);
835      unsigned long g = mpz_gcd_ui(NULL, erg, (unsigned long) ABS(SR_TO_INT(b)));
836      mpz_mul_si(erg, erg, SR_TO_INT(b)/g);
837    }
838    else
839    {
840      mpz_init(erg);
841      mpz_lcm(erg, (mpz_ptr) a, (mpz_ptr) b);
842    }
843  }
844  return (number) erg;
845}
846
847static number nrzCopy(number a, const coeffs)
848{
849  if (n_Z_IS_SMALL(a)) return a;
850  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
851  mpz_init_set(erg, (mpz_ptr) a);
852  return (number) erg;
853}
854
855/*
856 * Give the largest non unit k, such that a = x * k, b = y * k has
857 * a solution.
858 */
859static number nrzGcd (number a,number b,const coeffs R)
860{
861  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
862  {
863    long g = int_gcd(SR_TO_INT(a), SR_TO_INT(b));
864    return INT_TO_SR(g);
865  }
866  else if (n_Z_IS_SMALL(a))
867  {
868    if (a==INT_TO_SR(0))
869      return nrzCopy(b, R);
870    unsigned long g = mpz_gcd_ui(NULL, (mpz_ptr)b, (unsigned long) ABS(SR_TO_INT(a)));
871    return INT_TO_SR( g);
872  }
873  else if (n_Z_IS_SMALL(b))
874  {
875    if (b==INT_TO_SR(0))
876      return nrzCopy(a, R);
877    unsigned long g = mpz_gcd_ui(NULL, (mpz_ptr)a, (unsigned long) ABS(SR_TO_INT(b)));
878    return INT_TO_SR(g);
879  }
880  else
881  {
882    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
883    mpz_init(erg);
884    mpz_gcd(erg, (mpz_ptr) a, (mpz_ptr) b);
885    return (number) erg;
886  }
887}
888
889/*
890 * Give the largest non unit k, such that a = x * k, b = y * k has
891 * a solution and r, s, s.t. k = s*a + t*b
892 */
893static long int_extgcd(long a, long b, long * u, long* x, long * v, long* y)
894{
895  long q, r;
896  if (!a)
897  {
898    *u = 0;
899    *v = 1;
900    *x = -1;
901    *y = 0;
902    return b;
903  }
904  if (!b)
905  {
906    *u = 1;
907    *v = 0;
908    *x = 0;
909    *y = 1;
910    return a;
911  }
912  *u=1;
913  *v=0;
914  *x=0;
915  *y=1;
916  do
917  {
918    q = a/b;
919    r = a%b;
920    assume (q*b+r == a);
921    a = b;
922    b = r;
923
924    r = -(*v)*q+(*u);
925    (*u) =(*v);
926    (*v) = r;
927
928    r = -(*y)*q+(*x);
929    (*x) = (*y);
930    (*y) = r;
931  } while (b);
932
933  return a;
934}
935
936static number  nrzExtGcd (number a, number b, number *s, number *t, const coeffs)
937{
938  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
939  {
940    long u, v, x, y;
941    long g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &u, &v, &x, &y);
942    *s = INT_TO_SR(u);
943    *t = INT_TO_SR(v);
944    return INT_TO_SR(g);
945  }
946  else
947  {
948    mpz_t aa, bb;
949    if (n_Z_IS_SMALL(a))
950    {
951      mpz_init_set_si(aa, SR_TO_INT(a));
952    }
953    else
954    {
955      mpz_init_set(aa, (mpz_ptr) a);
956    }
957    if (n_Z_IS_SMALL(b))
958    {
959      mpz_init_set_si(bb, SR_TO_INT(b));
960    }
961    else
962    {
963      mpz_init_set(bb, (mpz_ptr) b);
964    }
965    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
966    mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
967    mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
968    mpz_init(erg);
969    mpz_init(bs);
970    mpz_init(bt);
971    mpz_gcdext(erg, bs, bt, aa, bb);
972    *s = nrz_short((number) bs);
973    *t = nrz_short((number) bt);
974    mpz_clear(aa);
975    mpz_clear(bb);
976    return nrz_short((number) erg);
977  }
978}
979#if CF_DEBUG
980static number _nrzXExtGcd(number, number, number *, number *, number *, number *, const coeffs);
981static number nrzXExtGcd(number a, number b, number *x, number * y, number * u, number * v, const coeffs R)
982{
983  char * s;
984  StringSetS("XExtGcd: ");
985  nrzWrite(a, R);
986  StringAppendS(" by ");
987  nrzWrite(b, R);
988  number c = _nrzXExtGcd(a, b, x, y, u, v, R);
989  StringAppendS(" is ");
990  nrzWrite(c, R);
991  StringAppendS("[[");
992  nrzWrite(*x, R);
993  StringAppendS(", ");
994  nrzWrite(*y, R);
995  StringAppendS("], ");
996  nrzWrite(*u, R);
997  StringAppendS(", ");
998  nrzWrite(*v, R);
999  s=StringEndS();
1000  Print("%s]]\n", s);
1001  omFree(s);
1002  return c;
1003}
1004static number  _nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )
1005#else
1006static number  nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )
1007#endif
1008{
1009  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1010  {
1011    long uu, vv, x, y;
1012    long g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &uu, &vv, &x, &y);
1013    *s = INT_TO_SR(uu);
1014    *t = INT_TO_SR(vv);
1015    *u = INT_TO_SR(x);
1016    *v = INT_TO_SR(y);
1017    return INT_TO_SR(g);
1018  }
1019  else
1020  {
1021    mpz_t aa, bb;
1022    if (n_Z_IS_SMALL(a))
1023    {
1024      mpz_init_set_si(aa, SR_TO_INT(a));
1025    }
1026    else
1027    {
1028      mpz_init_set(aa, (mpz_ptr) a);
1029    }
1030    if (n_Z_IS_SMALL(b))
1031    {
1032      mpz_init_set_si(bb, SR_TO_INT(b));
1033    }
1034    else
1035    {
1036      mpz_init_set(bb, (mpz_ptr) b);
1037    }
1038    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1039    mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1040    mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1041    mpz_init(erg);
1042    mpz_init(bs);
1043    mpz_init(bt);
1044
1045    mpz_gcdext(erg, bs, bt, aa, bb);
1046
1047    mpz_ptr bu = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1048    mpz_ptr bv = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1049
1050    mpz_init_set(bu, (mpz_ptr) bb);
1051    mpz_init_set(bv, (mpz_ptr) aa);
1052
1053    mpz_clear(aa);
1054    mpz_clear(bb);
1055    assume(mpz_cmp_si(erg, 0));
1056
1057    mpz_div(bu, bu, erg);
1058    mpz_div(bv, bv, erg);
1059
1060    mpz_mul_si(bu, bu, -1);
1061    *u = nrz_short((number) bu);
1062    *v = nrz_short((number) bv);
1063
1064    *s = nrz_short((number) bs);
1065    *t = nrz_short((number) bt);
1066    return nrz_short((number) erg);
1067  }
1068}
1069#if CF_DEBUG
1070static number _nrzQuotRem(number, number, number *, const coeffs);
1071static number nrzQuotRem(number a, number b, number * r, const coeffs R)
1072{
1073  StringSetS("QuotRem: ");
1074  nrzWrite(a, R);
1075  StringAppendS(" by ");
1076  nrzWrite(b, R);
1077  number c = _nrzQuotRem(a, b, r, R);
1078  StringAppendS(" is ");
1079  nrzWrite(c, R);
1080  if (r) {
1081    StringAppendS("+R(");
1082    nrzWrite(*r, R);
1083    StringAppendS(")");
1084  }
1085  char * s = StringEndS();
1086  Print("%s\n", s);
1087  omFree(s);
1088  return c;
1089}
1090static number _nrzQuotRem (number a, number b, number * r, const coeffs )
1091#else
1092static number nrzQuotRem (number a, number b, number * r, const coeffs )
1093#endif
1094{
1095  assume(SR_TO_INT(b));
1096  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1097  {
1098    if (r)
1099      *r = INT_TO_SR(SR_TO_INT(a) % SR_TO_INT(b));
1100    return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
1101  }
1102  else if (n_Z_IS_SMALL(a))
1103  {
1104    //a is small, b is not, so q=0, r=a
1105    if (r)
1106      *r = a;
1107    return INT_TO_SR(0);
1108  }
1109  else if (n_Z_IS_SMALL(b))
1110  {
1111    unsigned long rr;
1112    mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1113    mpz_init(qq);
1114    mpz_t rrr;
1115    mpz_init(rrr);
1116    rr = mpz_divmod_ui(qq, rrr, (mpz_ptr) a, (unsigned long)ABS(SR_TO_INT(b)));
1117    mpz_clear(rrr);
1118
1119    if (r)
1120      *r = INT_TO_SR(rr);
1121    if (SR_TO_INT(b)<0)
1122    {
1123      mpz_mul_si(qq, qq, -1);
1124    }
1125    return nrz_short((number)qq);
1126  }
1127  mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin),
1128             rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1129  mpz_init(qq);
1130  mpz_init(rr);
1131  mpz_divmod(qq, rr, (mpz_ptr)a, (mpz_ptr)b);
1132  if (r)
1133    *r = (number) rr;
1134  else
1135  {
1136    mpz_clear(rr);
1137  }
1138  nrzTest((number)qq);
1139  return (number) qq;
1140}
1141
1142static void nrzPower (number a, int i, number * result, const coeffs)
1143{
1144  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1145  mpz_init(erg);
1146  mpz_t aa;
1147  if (n_Z_IS_SMALL(a))
1148    mpz_init_set_si(aa, SR_TO_INT(a));
1149  else
1150    mpz_init_set(aa, (mpz_ptr) a);
1151  mpz_pow_ui(erg, aa, i);
1152  *result = nrz_short((number) erg);
1153}
1154
1155/*
1156 * create a number from int
1157 * TODO: do not create an mpz if not necessary
1158 */
1159number nrzInit (long i, const coeffs)
1160{
1161  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1162  mpz_init_set_si(erg, i);
1163  return nrz_short((number) erg);
1164}
1165
1166static number nrzInitMPZ(mpz_t m, const coeffs)
1167{
1168  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1169  mpz_init_set(erg, m);
1170  return nrz_short((number) erg);
1171}
1172
1173
1174static void nrzDelete(number *a, const coeffs)
1175{
1176  if (*a == NULL) return;
1177  if (n_Z_IS_SMALL(*a)==0)
1178  {
1179    mpz_clear((mpz_ptr) *a);
1180    omFreeBin((ADDRESS) *a, gmp_nrz_bin);
1181  }
1182  *a = NULL;
1183}
1184
1185/*
1186 * convert a number to int
1187 */
1188static long nrzInt(number &n, const coeffs)
1189{
1190  if (n_Z_IS_SMALL(n)) return SR_TO_INT(n);
1191  return mpz_get_si( (mpz_ptr)n);
1192}
1193#if CF_DEBUG
1194static number _nrzAdd(number, number, const coeffs);
1195static number nrzAdd(number a, number b, const coeffs R)
1196{
1197  StringSetS("Add: ");
1198  nrzWrite(a, R);
1199  StringAppendS(" to ");
1200  nrzWrite(b, R);
1201  number c = _nrzAdd(a, b, R);
1202  StringAppendS(" is ");
1203  nrzWrite(c, R);
1204  char * s = StringEndS();
1205  Print("%s\n", s);
1206  omFree(s);
1207  return c;
1208}
1209static number _nrzAdd (number a, number b, const coeffs )
1210#else
1211static number nrzAdd (number a, number b, const coeffs )
1212#endif
1213{
1214  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1215  {
1216    long c = SR_TO_INT(a) + SR_TO_INT(b);
1217    if (INT_IS_SMALL(c))
1218      return INT_TO_SR(c);
1219    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1220    mpz_init_set_si(erg, c);
1221
1222    nrzTest((number)erg);
1223    return (number) erg;
1224  }
1225  else if (n_Z_IS_SMALL(a))
1226  {
1227    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1228    mpz_init(erg);
1229    if (SR_TO_INT(a)>0)
1230      mpz_add_ui(erg, (mpz_ptr) b, (unsigned long)SR_TO_INT(a));
1231    else
1232      mpz_sub_ui(erg, (mpz_ptr) b, (unsigned long)-(SR_TO_INT(a)));
1233    return nrz_short((number) erg);
1234  }
1235  else if (n_Z_IS_SMALL(b))
1236  {
1237    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1238    mpz_init(erg);
1239    if (SR_TO_INT(b)>0)
1240      mpz_add_ui(erg, (mpz_ptr) a, (unsigned long)SR_TO_INT(b));
1241    else
1242      mpz_sub_ui(erg, (mpz_ptr) a, (unsigned long)-(SR_TO_INT(b)));
1243    return nrz_short((number) erg);
1244  }
1245  else
1246  {
1247    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1248    mpz_init(erg);
1249    mpz_add(erg, (mpz_ptr) a, (mpz_ptr) b);
1250    return nrz_short((number) erg);
1251  }
1252}
1253
1254static number nrzSub (number a, number b,  const coeffs )
1255{
1256  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1257  {
1258    long c = SR_TO_INT(a) - SR_TO_INT(b);
1259    if (INT_IS_SMALL(c))
1260      return INT_TO_SR(c);
1261    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1262    mpz_init_set_si(erg, c);
1263    nrzTest((number)erg);
1264    return (number) erg;
1265  }
1266  else if (n_Z_IS_SMALL(a))
1267  {
1268    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1269    mpz_init(erg);
1270
1271    if (SR_TO_INT(a)>0)
1272      mpz_ui_sub(erg, (unsigned long)SR_TO_INT(a), (mpz_ptr) b);
1273    else
1274    {
1275      mpz_add_ui(erg, (mpz_ptr) b, (unsigned long)-SR_TO_INT(a));
1276      mpz_neg(erg, erg);
1277    }
1278    return nrz_short((number) erg);
1279  }
1280  else if (n_Z_IS_SMALL(b))
1281  {
1282    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1283    mpz_init(erg);
1284    if (SR_TO_INT(b)>0)
1285      mpz_sub_ui(erg, (mpz_ptr) a, (unsigned long)SR_TO_INT(b));
1286    else
1287      mpz_add_ui(erg, (mpz_ptr) a, (unsigned long)-SR_TO_INT(b));
1288    return nrz_short((number) erg);
1289  }
1290  else
1291  {
1292    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1293    mpz_init(erg);
1294    mpz_sub(erg, (mpz_ptr) a, (mpz_ptr) b);
1295    return nrz_short((number) erg);
1296  }
1297}
1298
1299static BOOLEAN nrzGreater (number a,number b, const coeffs)
1300{
1301  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1302    return ((long)a)>((long)b);
1303  else if (n_Z_IS_SMALL(a))
1304    return 0 > mpz_cmp_si((mpz_ptr)b,SR_TO_INT(a));
1305  else if (n_Z_IS_SMALL(b))
1306    return 0 < mpz_cmp_si((mpz_ptr)a,SR_TO_INT(b));
1307  return 0 < mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
1308}
1309
1310static BOOLEAN nrzGreaterZero (number k, const coeffs C)
1311{
1312  return nrzGreater(k, INT_TO_SR(0), C);
1313}
1314
1315static number  nrzGetUnit (number n, const coeffs r)
1316{
1317  if (nrzGreaterZero(n, r))
1318    return INT_TO_SR(1);
1319  /*else*/
1320    return INT_TO_SR(-1);
1321}
1322
1323static number nrzAnn(number n, const coeffs)
1324{
1325  if (SR_TO_INT(n))  // in Z: the annihilator of !=0 is 0
1326    return INT_TO_SR(0);
1327  else
1328    return INT_TO_SR(1);
1329}
1330
1331static BOOLEAN nrzIsUnit (number a, const coeffs)
1332{
1333  return ABS(SR_TO_INT(a))==1;
1334}
1335
1336static BOOLEAN nrzIsZero (number  a, const coeffs)
1337{
1338  return (a==INT_TO_SR(0));
1339}
1340
1341static BOOLEAN nrzIsOne (number a, const coeffs)
1342{
1343  return a==INT_TO_SR(1);
1344}
1345
1346static BOOLEAN nrzIsMOne (number a, const coeffs)
1347{
1348  return a==INT_TO_SR(-1);
1349}
1350
1351static BOOLEAN nrzEqual (number a,number b, const coeffs)
1352{
1353  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1354    return a==b;
1355  else if (n_Z_IS_SMALL(a) || n_Z_IS_SMALL(b))
1356    return FALSE;
1357  else
1358    return 0 == mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
1359}
1360
1361static BOOLEAN nrzDivBy (number a,number b, const coeffs)
1362{
1363  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1364  {
1365    return SR_TO_INT(a) %SR_TO_INT(b) ==0;
1366  }
1367  else if (n_Z_IS_SMALL(a))
1368  {
1369    return a==INT_TO_SR(0);
1370  }
1371  else if (n_Z_IS_SMALL(b))
1372  {
1373    return mpz_divisible_ui_p((mpz_ptr)a, (unsigned long)ABS(SR_TO_INT(b))) != 0;
1374  }
1375  else
1376    return mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b) != 0;
1377}
1378
1379static int nrzDivComp(number a, number b, const coeffs r)
1380{
1381  if (nrzDivBy(a, b, r))
1382  {
1383    if (nrzDivBy(b, a, r)) return 2;
1384    return -1;
1385  }
1386  if (nrzDivBy(b, a, r)) return 1;
1387  return 0;
1388}
1389
1390static number nrzDiv (number a,number b, const coeffs)
1391{
1392  assume(SR_TO_INT(b));
1393  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
1394  {
1395    //if (SR_TO_INT(a) % SR_TO_INT(b))
1396    //{
1397    //  WerrorS("1:Division by non divisible element.");
1398    //  WerrorS("Result is without remainder.");
1399    //}
1400    return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
1401  }
1402  else if (n_Z_IS_SMALL(a))
1403  {
1404    //if (SR_TO_INT(a))
1405    //{
1406    //  WerrorS("2:Division by non divisible element.");
1407    //  WerrorS("Result is without remainder.");
1408    //}
1409    return INT_TO_SR(0);
1410  }
1411  else if (n_Z_IS_SMALL(b))
1412  {
1413    mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1414    mpz_t r;
1415    mpz_init(r);
1416    mpz_init(erg);
1417    if (mpz_divmod_ui(erg, r, (mpz_ptr) a, (unsigned long)ABS(SR_TO_INT(b)))) {
1418    //  WerrorS("3:Division by non divisible element.");
1419    //  WerrorS("Result is without remainder.");
1420    }
1421    mpz_clear(r);
1422    if (SR_TO_INT(b)<0)
1423      mpz_neg(erg, erg);
1424    return nrz_short((number) erg);
1425  }
1426  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1427  mpz_init(erg);
1428  mpz_t r;
1429  mpz_init(r);
1430  mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);
1431#if CF_DEBUG
1432  StringSetS("division of");
1433  nrzWrite(a, R);
1434  StringAppendS(" by ");
1435  nrzWrite(b, R);
1436  StringAppendS(" is ");
1437  number du;
1438  nrzWrite(du = (number)erg, R);
1439  StringAppendS(" rest ");
1440  nrzWrite(du = (number)r, R);
1441  char * s = StringEndS();
1442  Print("%s\n", s);
1443  omFree(s);
1444#endif
1445
1446  if (mpz_sgn1(r)!=0)
1447  {
1448    //WerrorS("4:Division by non divisible element.");
1449    //WerrorS("Result is without remainder.");
1450  }
1451  mpz_clear(r);
1452  return nrz_short((number) erg);
1453}
1454
1455static number nrzExactDiv (number a,number b, const coeffs)
1456{
1457  assume(SR_TO_INT(b));
1458  mpz_t aa, bb;
1459  if (n_Z_IS_SMALL(a))
1460    mpz_init_set_si(aa, SR_TO_INT(a));
1461  else
1462    mpz_init_set(aa, (mpz_ptr) a);
1463  if (n_Z_IS_SMALL(b))
1464    mpz_init_set_si(bb, SR_TO_INT(b));
1465  else
1466    mpz_init_set(bb, (mpz_ptr) b);
1467  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1468  mpz_init(erg);
1469  mpz_tdiv_q(erg, (mpz_ptr) aa, (mpz_ptr) bb);
1470  mpz_clear(aa);
1471  mpz_clear(bb);
1472  nrzTest((number)erg);
1473  return (number) erg;
1474}
1475
1476static number nrzIntMod (number a,number b, const coeffs)
1477{
1478  mpz_t aa, bb;
1479  assume(SR_TO_INT(b));
1480  if (n_Z_IS_SMALL(a))
1481    mpz_init_set_si(aa, SR_TO_INT(a));
1482  else
1483    mpz_init_set(aa, (mpz_ptr) a);
1484  if (n_Z_IS_SMALL(b))
1485    mpz_init_set_si(bb, SR_TO_INT(b));
1486  else
1487    mpz_init_set(bb, (mpz_ptr) b);
1488
1489  mpz_t erg;
1490  mpz_init(erg);
1491  mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1492  mpz_init(r);
1493  mpz_tdiv_qr(erg, r, (mpz_ptr) aa, (mpz_ptr) bb);
1494  mpz_clear(erg);
1495  mpz_clear(aa);
1496  mpz_clear(bb);
1497
1498  return nrz_short((number) r);
1499}
1500
1501static number  nrzInvers (number c, const coeffs r)
1502{
1503  if (!nrzIsUnit((number) c, r))
1504  {
1505    WerrorS("Non invertible element.");
1506    return (number)NULL;
1507  }
1508  return c; // has to be 1 or -1....
1509}
1510
1511static number nrzNeg (number c, const coeffs)
1512{
1513// nNeg inplace !!!
1514  if (n_Z_IS_SMALL(c))
1515    return INT_TO_SR(-SR_TO_INT(c));
1516  mpz_mul_si((mpz_ptr) c, (mpz_ptr) c, -1);
1517  return c;
1518}
1519
1520static number nrzFarey(number r, number N, const coeffs R)
1521{
1522  number a0 = nrzCopy(N, R);
1523  number b0 = nrzInit(0, R);
1524  number a1 = nrzCopy(r, R);
1525  number b1 = nrzInit(1, R);
1526  number two = nrzInit(2, R);
1527#if 0
1528  PrintS("Farey start with ");
1529  n_Print(r, R);
1530  PrintS(" mod ");
1531  n_Print(N, R);
1532  PrintLn();
1533#endif
1534  while (1)
1535  {
1536    number as = nrzMult(a1, a1, R);
1537    n_InpMult(as, two, R);
1538    if (nrzGreater(N, as, R))
1539    {
1540      nrzDelete(&as, R);
1541      break;
1542    }
1543    nrzDelete(&as, R);
1544    number q = nrzDiv(a0, a1, R);
1545    number t = nrzMult(a1, q, R),
1546           s = nrzSub(a0, t, R);
1547    nrzDelete(&a0, R);
1548    a0 = a1;
1549    a1 = s;
1550    nrzDelete(&t, R);
1551
1552    t = nrzMult(b1, q, R);
1553    s = nrzSub(b0, t, R);
1554    nrzDelete(&b0, R);
1555    b0 = b1;
1556    b1 = s;
1557    nrzDelete(&t, R);
1558    nrzDelete(&q, R);
1559  }
1560  number as = nrzMult(b1, b1, R);
1561  n_InpMult(as, two, R);
1562  nrzDelete(&two, R);
1563  if (nrzGreater(as, N, R))
1564  {
1565    nrzDelete(&a0, R);
1566    nrzDelete(&a1, R);
1567    nrzDelete(&b0, R);
1568    nrzDelete(&b1, R);
1569    nrzDelete(&as, R);
1570    return NULL;
1571  }
1572  nrzDelete(&as, R);
1573  nrzDelete(&a0, R);
1574  nrzDelete(&b0, R);
1575
1576  number a, b, ab;
1577  coeffs Q = nInitChar(n_Q, 0);
1578  nMapFunc f = n_SetMap(R, Q);
1579  a = f(a1, R, Q);
1580  b = f(b1, R, Q);
1581  ab = n_Div(a, b, Q);
1582  n_Delete(&a, Q);
1583  n_Delete(&b, Q);
1584  nKillChar(Q);
1585
1586  nrzDelete(&a1, R);
1587  nrzDelete(&b1, R);
1588  return ab;
1589}
1590
1591static number nrzMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
1592{
1593  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1594  mpz_init_set_ui(erg, (unsigned long) from);
1595  return nrz_short((number) erg);
1596}
1597
1598static number nrzMapZp(number from, const coeffs /*src*/, const coeffs /*dst*/)
1599{
1600  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1601  mpz_init_set_si(erg, (long) from);
1602  return nrz_short((number) erg);
1603}
1604
1605static number nrzModNMap(number from, const coeffs /* src */, const coeffs /*dst*/)
1606{
1607  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1608  mpz_init_set(erg, (mpz_ptr) from);
1609  return nrz_short((number) erg);
1610}
1611
1612static number nrzMapQ(number from, const coeffs /* src */, const coeffs dst)
1613{
1614  if (SR_HDL(from) & SR_INT)
1615    return nrzInit(SR_TO_INT(from),dst);
1616  if (from->s!=3)
1617  {
1618    WerrorS("rational in map to integer");
1619    return NULL;
1620  }
1621  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1622  mpz_init_set(erg, from->z);
1623  return nrz_short((number) erg);
1624}
1625
1626static nMapFunc nrzSetMap(const coeffs src, const coeffs /*dst*/)
1627{
1628  /* dst = rintegers */
1629  if (src->rep==n_rep_gmp) //nCoeff_is_Ring_ModN(src) || nCoeff_is_Ring_PtoM(src))
1630    return nrzModNMap;
1631
1632  if ((src->rep==n_rep_gap_gmp) && nCoeff_is_Ring_Z(src))
1633  {
1634    return ndCopyMap; //nrzCopyMap;
1635  }
1636  if (src->rep==n_rep_gap_rat) /*&& nCoeff_is_Ring_Z(src)) Q, bigint*/
1637  {
1638    return nrzMapQ;
1639  }
1640  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src))
1641  {
1642    return nrzMapMachineInt;
1643  }
1644  if ((src->rep==n_rep_int) && nCoeff_is_Zp(src))
1645  {
1646    return nrzMapZp;
1647  }
1648  return NULL;      // default
1649}
1650
1651
1652/*
1653 * set the exponent (allocate and init tables) (TODO)
1654 */
1655
1656void nrzSetExp(int, coeffs)
1657{
1658}
1659
1660void nrzInitExp(int, coeffs)
1661{
1662}
1663
1664#ifdef LDEBUG
1665BOOLEAN nrzDBTest (number x, const char *f, const int l, const coeffs)
1666{
1667  if (SR_HDL(x) & SR_INT) return TRUE;
1668  if (mpz_sgn1((mpz_ptr) x)==0)
1669  {
1670    Print("gmp-0 %s:%d\n",f,l);
1671    return FALSE;
1672  }
1673  if (mpz_size1((mpz_ptr)x)<=MP_SMALL)
1674  {
1675    long ui=mpz_get_si((mpz_ptr)x);
1676    if ((((ui<<3)>>3)==ui)
1677    && (mpz_cmp_si((mpz_ptr)x,ui)==0))
1678    {
1679      Print("gmp-small %s:%d\n",f,l);
1680      return FALSE;
1681    }
1682  }
1683  return TRUE;
1684}
1685#endif
1686
1687void nrzWrite (number a, const coeffs)
1688{
1689  char *s,*z;
1690  if (a==NULL)
1691  {
1692    StringAppendS("o");
1693  }
1694  else
1695  {
1696    if (n_Z_IS_SMALL(a))
1697    {
1698      StringAppend("%d", SR_TO_INT(a));
1699    }
1700    else
1701    {
1702      int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2;
1703      s=(char*)omAlloc(l);
1704      z=mpz_get_str(s,10,(mpz_ptr) a);
1705      StringAppendS(z);
1706      omFreeSize((ADDRESS)s,l);
1707    }
1708  }
1709}
1710
1711/*2
1712* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
1713*/
1714static const char * nlEatLongC(char *s, mpz_ptr i)
1715{
1716  const char * start=s;
1717
1718  if (*s<'0' || *s>'9')
1719  {
1720    mpz_set_ui(i,1);
1721    return s;
1722  }
1723  while (*s >= '0' && *s <= '9') s++;
1724  if (*s=='\0')
1725  {
1726    mpz_set_str(i,start,10);
1727  }
1728  else
1729  {
1730    char c=*s;
1731    *s='\0';
1732    mpz_set_str(i,start,10);
1733    *s=c;
1734  }
1735  return s;
1736}
1737
1738static const char * nrzRead (const char *s, number *a, const coeffs)
1739{
1740  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
1741  {
1742    mpz_init(z);
1743    s = nlEatLongC((char *) s, z);
1744  }
1745  *a = nrz_short((number) z);
1746  return s;
1747}
1748
1749static void nrzCoeffWrite  (const coeffs, BOOLEAN /*details*/)
1750{
1751  //PrintS("// ZZ\n");
1752  PrintS("//   coeff. ring is : Integers\n");
1753}
1754
1755static char* nrzCoeffString(const coeffs)
1756{
1757  return omStrDup("integer");
1758}
1759
1760static CanonicalForm nrzConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs /*r*/ )
1761{
1762  if (setChar) setCharacteristic( 0 );
1763
1764  CanonicalForm term;
1765  if ( n_Z_IS_SMALL(n))
1766  {
1767    term = SR_TO_INT(n);
1768  }
1769  else
1770  {
1771    mpz_t dummy;
1772    mpz_init_set( dummy,n->z );
1773    term = make_cf( dummy );
1774  }
1775  return term;
1776}
1777
1778static number nrzConvFactoryNSingN( const CanonicalForm n, const coeffs r)
1779{
1780  if (n.isImm())
1781  {
1782    return nrzInit(n.intval(),r);
1783  }
1784  else
1785  {
1786    if ( !n.den().isOne() )
1787    {
1788     WerrorS("rational in conversion to integer");
1789     return NULL;
1790    }
1791    mpz_ptr z = (mpz_ptr) omAlloc0Bin(gmp_nrz_bin);
1792    gmp_numerator( n,z);
1793    return nrz_short((number)z);
1794  }
1795}
1796
1797static void nrzMPZ(mpz_t res, number &a, const coeffs)
1798{
1799  if (n_Z_IS_SMALL(a))
1800    mpz_init_set_si(res, SR_TO_INT(a));
1801  else
1802    mpz_init_set(res, (mpz_ptr) a);
1803}
1804
1805static coeffs nrzQuot1(number c, const coeffs r)
1806{
1807    mpz_t dummy;
1808    if(n_Z_IS_SMALL(c))
1809    {
1810      long ch = r->cfInt(c, r);
1811      mpz_init_set_ui(dummy, ch);
1812    }
1813    else
1814    {
1815      mpz_init_set(dummy, (mpz_ptr)c);
1816    }
1817    ZnmInfo info;
1818    info.base = dummy;
1819    info.exp = (unsigned long) 1;
1820    coeffs rr = nInitChar(n_Zn, (void*)&info);
1821    mpz_clear(dummy);
1822    return(rr);
1823}
1824
1825number nlReadFd(const ssiInfo *d, const coeffs);
1826void nlWriteFd(number n, const ssiInfo* d, const coeffs);
1827
1828BOOLEAN nrzInitChar(coeffs r,  void *)
1829{
1830  assume( getCoeffType(r) == n_Z );
1831
1832  r->is_field=FALSE;
1833  r->is_domain=TRUE;
1834  r->rep=n_rep_gap_gmp;
1835
1836  //r->nCoeffIsEqual = ndCoeffIsEqual;
1837  r->cfCoeffString = nrzCoeffString;
1838  //r->cfKillChar = ndKillChar;
1839  r->cfMult  = nrzMult;
1840  r->cfSub   = nrzSub;
1841  r->cfAdd   = nrzAdd;
1842  r->cfDiv   = nrzDiv;
1843  r->cfIntMod= nrzIntMod;
1844  r->cfExactDiv= nrzExactDiv;
1845  r->cfInit = nrzInit;
1846  r->cfInitMPZ = nrzInitMPZ;
1847  r->cfSize  = nrzSize;
1848  r->cfInt  = nrzInt;
1849  //#ifdef HAVE_RINGS
1850  r->cfDivComp = nrzDivComp; // only for ring stuff
1851  r->cfIsUnit = nrzIsUnit; // only for ring stuff
1852  r->cfGetUnit = nrzGetUnit; // only for ring stuff
1853  r->cfAnn = nrzAnn;
1854  r->cfExtGcd = nrzExtGcd; // only for ring stuff
1855  r->cfXExtGcd = nrzXExtGcd; // only for ring stuff
1856  r->cfQuotRem = nrzQuotRem;
1857  r->cfDivBy = nrzDivBy; // only for ring stuff
1858  //#endif
1859  r->cfInpNeg   = nrzNeg;
1860  r->cfInvers= nrzInvers;
1861  r->cfCopy  = nrzCopy;
1862  r->cfWriteLong = nrzWrite;
1863  r->cfRead = nrzRead;
1864  r->cfGreater = nrzGreater;
1865  r->cfEqual = nrzEqual;
1866  r->cfIsZero = nrzIsZero;
1867  r->cfIsOne = nrzIsOne;
1868  r->cfIsMOne = nrzIsMOne;
1869  r->cfGreaterZero = nrzGreaterZero;
1870  r->cfPower = nrzPower;
1871  r->cfGcd  = nrzGcd;
1872  r->cfLcm  = nrzLcm;
1873  r->cfDelete= nrzDelete;
1874  r->cfSetMap = nrzSetMap;
1875  r->cfCoeffWrite = nrzCoeffWrite;
1876  r->convSingNFactoryN = nrzConvSingNFactoryN;
1877  r->convFactoryNSingN = nrzConvFactoryNSingN;
1878  r->cfMPZ = nrzMPZ;
1879  r->cfFarey = nrzFarey;
1880  r->cfWriteFd=nlWriteFd;
1881  r->cfReadFd=nlReadFd;
1882
1883  r->cfQuot1 = nrzQuot1;
1884  // requires conversion to factory:
1885  r->cfChineseRemainder=nlChineseRemainderSym;
1886  // debug stuff
1887
1888#ifdef LDEBUG
1889  r->cfDBTest=nrzDBTest;
1890#endif
1891
1892  r->ch = 0;
1893  r->has_simple_Alloc=FALSE;
1894  r->has_simple_Inverse=FALSE;
1895  return FALSE;
1896}
1897
1898#elif SI_INTEGER_VARIANT == 1
1899BOOLEAN nrzInitChar(coeffs r,  void *)
1900{
1901  return nlInitChar(r,(void*)1);
1902}
1903#else
1904#error set SI_INTEGER_VARIANT
1905#endif
1906#endif
Note: See TracBrowser for help on using the repository browser.