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

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