source: git/libpolys/coeffs/rintegers.cc @ 7f7808

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