source: git/libpolys/coeffs/rintegers.cc @ 8d1432e

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