source: git/libpolys/coeffs/rintegers.cc @ 12e275

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