source: git/libpolys/coeffs/rintegers.cc @ 5a0d2ae

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