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

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