source: git/libpolys/coeffs/rintegers.cc @ 780c1b

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