source: git/libpolys/coeffs/rmodulon.cc @ 66ce6d

fieker-DuValspielwiese
Last change on this file since 66ce6d was 66ce6d, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: got rid of ringtype
  • Property mode set to 100644
File size: 16.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: numbers modulo n
6*/
7
8#include "config.h"
9#include <misc/auxiliary.h>
10
11#ifdef HAVE_RINGS
12
13#include <misc/mylimits.h>
14#include <coeffs/coeffs.h>
15#include <reporter/reporter.h>
16#include <omalloc/omalloc.h>
17#include <coeffs/numbers.h>
18#include <coeffs/longrat.h>
19#include <coeffs/mpr_complex.h>
20#include <coeffs/rmodulon.h>
21#include "si_gmp.h"
22
23#include <string.h>
24
25/// Our Type!
26static const n_coeffType ID = n_Zn;
27static const n_coeffType ID2 = n_Znm;
28
29extern omBin gmp_nrz_bin;
30
31void    nrnCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
32{
33  long l = (long)mpz_sizeinbase(r->modBase, 10) + 2;
34  char* s = (char*) omAlloc(l);
35  s= mpz_get_str (s, 10, r->modBase);
36  if (nCoeff_is_Ring_ModN(r)) Print("//   coeff. ring is : Z/%s\n", s);
37  else if (nCoeff_is_Ring_PtoM(r)) Print("//   coeff. ring is : Z/%s^%lu\n", s, r->modExponent);
38  omFreeSize((ADDRESS)s, l);
39}
40
41static BOOLEAN nrnCoeffsEqual(const coeffs r, n_coeffType n, void * parameter)
42{
43  /* test, if r is an instance of nInitCoeffs(n,parameter) */
44  return (n==n_Zn) && (mpz_cmp_ui(r->modNumber,(long)parameter)==0);
45}
46
47
48/* for initializing function pointers */
49BOOLEAN nrnInitChar (coeffs r, void* p)
50{
51  assume( (getCoeffType(r) == ID) || (getCoeffType (r) == ID2) );
52  ZnmInfo * info= (ZnmInfo *) p;
53  r->modBase= info->base;
54
55  nrnInitExp (info->exp, r);
56
57  /* next computation may yield wrong characteristic as r->modNumber
58     is a GMP number */
59  r->ch = mpz_get_ui(r->modNumber);
60
61  r->cfInit        = nrnInit;
62  r->cfDelete      = nrnDelete;
63  r->cfCopy        = nrnCopy;
64  r->cfSize        = nrnSize;
65  r->cfInt         = nrnInt;
66  r->cfAdd         = nrnAdd;
67  r->cfSub         = nrnSub;
68  r->cfMult        = nrnMult;
69  r->cfDiv         = nrnDiv;
70  r->cfIntDiv      = nrnIntDiv;
71  r->cfIntMod      = nrnMod;
72  r->cfExactDiv    = nrnDiv;
73  r->cfNeg         = nrnNeg;
74  r->cfInvers      = nrnInvers;
75  r->cfDivBy       = nrnDivBy;
76  r->cfDivComp     = nrnDivComp;
77  r->cfGreater     = nrnGreater;
78  r->cfEqual       = nrnEqual;
79  r->cfIsZero      = nrnIsZero;
80  r->cfIsOne       = nrnIsOne;
81  r->cfIsMOne      = nrnIsMOne;
82  r->cfGreaterZero = nrnGreaterZero;
83  r->cfWriteLong       = nrnWrite;
84  r->cfRead        = nrnRead;
85  r->cfPower       = nrnPower;
86  r->cfSetMap      = nrnSetMap;
87  r->cfNormalize   = ndNormalize;
88  r->cfLcm         = nrnLcm;
89  r->cfGcd         = nrnGcd;
90  r->cfIsUnit      = nrnIsUnit;
91  r->cfGetUnit     = nrnGetUnit;
92  r->cfExtGcd      = nrnExtGcd;
93  r->cfName        = ndName;
94  r->cfCoeffWrite  = nrnCoeffWrite;
95  r->nCoeffIsEqual = nrnCoeffsEqual;
96  r->cfInit_bigint = nrnMapQ;
97
98#ifdef LDEBUG
99  r->cfDBTest      = nrnDBTest;
100#endif
101  return FALSE;
102}
103
104/*
105 * create a number from int
106 */
107number nrnInit(long i, const coeffs r)
108{
109  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
110  mpz_init_set_si(erg, i);
111  mpz_mod(erg, erg, r->modNumber);
112  return (number) erg;
113}
114
115void nrnDelete(number *a, const coeffs)
116{
117  if (*a == NULL) return;
118  mpz_clear((int_number) *a);
119  omFreeBin((void *) *a, gmp_nrz_bin);
120  *a = NULL;
121}
122
123number nrnCopy(number a, const coeffs)
124{
125  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
126  mpz_init_set(erg, (int_number) a);
127  return (number) erg;
128}
129
130int nrnSize(number a, const coeffs)
131{
132  if (a == NULL) return 0;
133  return sizeof(mpz_t);
134}
135
136/*
137 * convert a number to int
138 */
139int nrnInt(number &n, const coeffs)
140{
141  return (int)mpz_get_si((int_number) n);
142}
143
144/*
145 * Multiply two numbers
146 */
147number nrnMult(number a, number b, const coeffs r)
148{
149  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
150  mpz_init(erg);
151  mpz_mul(erg, (int_number)a, (int_number) b);
152  mpz_mod(erg, erg, r->modNumber);
153  return (number) erg;
154}
155
156void nrnPower(number a, int i, number * result, const coeffs r)
157{
158  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
159  mpz_init(erg);
160  mpz_powm_ui(erg, (int_number)a, i, r->modNumber);
161  *result = (number) erg;
162}
163
164number nrnAdd(number a, number b, const coeffs r)
165{
166  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
167  mpz_init(erg);
168  mpz_add(erg, (int_number)a, (int_number) b);
169  mpz_mod(erg, erg, r->modNumber);
170  return (number) erg;
171}
172
173number nrnSub(number a, number b, const coeffs r)
174{
175  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
176  mpz_init(erg);
177  mpz_sub(erg, (int_number)a, (int_number) b);
178  mpz_mod(erg, erg, r->modNumber);
179  return (number) erg;
180}
181
182number nrnNeg(number c, const coeffs r)
183{
184  if( !nrnIsZero(c, r) )
185    // Attention: This method operates in-place.
186    mpz_sub((int_number)c, r->modNumber, (int_number)c);
187  return c;
188}
189
190number nrnInvers(number c, const coeffs r)
191{
192  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
193  mpz_init(erg);
194  mpz_invert(erg, (int_number)c, r->modNumber);
195  return (number) erg;
196}
197
198/*
199 * Give the smallest k, such that a * x = k = b * y has a solution
200 * TODO: lcm(gcd,gcd) better than gcd(lcm) ?
201 */
202number nrnLcm(number a, number b, const coeffs r)
203{
204  number erg = nrnGcd(NULL, a, r);
205  number tmp = nrnGcd(NULL, b, r);
206  mpz_lcm((int_number)erg, (int_number)erg, (int_number)tmp);
207  nrnDelete(&tmp, NULL);
208  return (number)erg;
209}
210
211/*
212 * Give the largest k, such that a = x * k, b = y * k has
213 * a solution.
214 */
215number nrnGcd(number a, number b, const coeffs r)
216{
217  if ((a == NULL) && (b == NULL)) return nrnInit(0,r);
218  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
219  mpz_init_set(erg, r->modNumber);
220  if (a != NULL) mpz_gcd(erg, erg, (int_number)a);
221  if (b != NULL) mpz_gcd(erg, erg, (int_number)b);
222  return (number)erg;
223}
224
225/* Not needed any more, but may have room for improvement
226   number nrnGcd3(number a,number b, number c,ring r)
227{
228  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
229  mpz_init(erg);
230  if (a == NULL) a = (number)r->modNumber;
231  if (b == NULL) b = (number)r->modNumber;
232  if (c == NULL) c = (number)r->modNumber;
233  mpz_gcd(erg, (int_number)a, (int_number)b);
234  mpz_gcd(erg, erg, (int_number)c);
235  mpz_gcd(erg, erg, r->modNumber);
236  return (number)erg;
237}
238*/
239
240/*
241 * Give the largest k, such that a = x * k, b = y * k has
242 * a solution and r, s, s.t. k = s*a + t*b
243 */
244number nrnExtGcd(number a, number b, number *s, number *t, const coeffs r)
245{
246  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
247  int_number bs  = (int_number)omAllocBin(gmp_nrz_bin);
248  int_number bt  = (int_number)omAllocBin(gmp_nrz_bin);
249  mpz_init(erg);
250  mpz_init(bs);
251  mpz_init(bt);
252  mpz_gcdext(erg, bs, bt, (int_number)a, (int_number)b);
253  mpz_mod(bs, bs, r->modNumber);
254  mpz_mod(bt, bt, r->modNumber);
255  *s = (number)bs;
256  *t = (number)bt;
257  return (number)erg;
258}
259
260BOOLEAN nrnIsZero(number a, const coeffs)
261{
262#ifdef LDEBUG
263  if (a == NULL) return FALSE;
264#endif
265  return 0 == mpz_cmpabs_ui((int_number)a, 0);
266}
267
268BOOLEAN nrnIsOne(number a, const coeffs)
269{
270#ifdef LDEBUG
271  if (a == NULL) return FALSE;
272#endif
273  return 0 == mpz_cmp_si((int_number)a, 1);
274}
275
276BOOLEAN nrnIsMOne(number a, const coeffs r)
277{
278#ifdef LDEBUG
279  if (a == NULL) return FALSE;
280#endif
281  mpz_t t; mpz_init_set(t, (int_number)a);
282  mpz_add_ui(t, t, 1);
283  bool erg = (0 == mpz_cmp(t, r->modNumber));
284  mpz_clear(t);
285  return erg;
286}
287
288BOOLEAN nrnEqual(number a, number b, const coeffs)
289{
290  return 0 == mpz_cmp((int_number)a, (int_number)b);
291}
292
293BOOLEAN nrnGreater(number a, number b, const coeffs)
294{
295  return 0 < mpz_cmp((int_number)a, (int_number)b);
296}
297
298BOOLEAN nrnGreaterZero(number k, const coeffs)
299{
300  return 0 < mpz_cmp_si((int_number)k, 0);
301}
302
303BOOLEAN nrnIsUnit(number a, const coeffs r)
304{
305  number tmp = nrnGcd(a, (number)r->modNumber, r);
306  bool res = nrnIsOne(tmp, r);
307  nrnDelete(&tmp, NULL);
308  return res;
309}
310
311number nrnGetUnit(number k, const coeffs r)
312{
313  if (mpz_divisible_p(r->modNumber, (int_number)k)) return nrnInit(1,r);
314
315  int_number unit = (int_number)nrnGcd(k, 0, r);
316  mpz_tdiv_q(unit, (int_number)k, unit);
317  int_number gcd = (int_number)nrnGcd((number)unit, 0, r);
318  if (!nrnIsOne((number)gcd,r))
319  {
320    int_number ctmp;
321    // tmp := unit^2
322    int_number tmp = (int_number) nrnMult((number) unit,(number) unit,r);
323    // gcd_new := gcd(tmp, 0)
324    int_number gcd_new = (int_number) nrnGcd((number) tmp, 0, r);
325    while (!nrnEqual((number) gcd_new,(number) gcd,r))
326    {
327      // gcd := gcd_new
328      ctmp = gcd;
329      gcd = gcd_new;
330      gcd_new = ctmp;
331      // tmp := tmp * unit
332      mpz_mul(tmp, tmp, unit);
333      mpz_mod(tmp, tmp, r->modNumber);
334      // gcd_new := gcd(tmp, 0)
335      mpz_gcd(gcd_new, tmp, r->modNumber);
336    }
337    // unit := unit + modNumber / gcd_new
338    mpz_tdiv_q(tmp, r->modNumber, gcd_new);
339    mpz_add(unit, unit, tmp);
340    mpz_mod(unit, unit, r->modNumber);
341    nrnDelete((number*) &gcd_new, NULL);
342    nrnDelete((number*) &tmp, NULL);
343  }
344  nrnDelete((number*) &gcd, NULL);
345  return (number)unit;
346}
347
348BOOLEAN nrnDivBy(number a, number b, const coeffs r)
349{
350  if (a == NULL)
351    return mpz_divisible_p(r->modNumber, (int_number)b);
352  else
353  { /* b divides a iff b/gcd(a, b) is a unit in the given ring: */
354    number n = nrnGcd(a, b, r);
355    mpz_tdiv_q((int_number)n, (int_number)b, (int_number)n);
356    bool result = nrnIsUnit(n, r);
357    nrnDelete(&n, NULL);
358    return result;
359  }
360}
361
362int nrnDivComp(number a, number b, const coeffs r)
363{
364  if (nrnEqual(a, b,r)) return 2;
365  if (mpz_divisible_p((int_number) a, (int_number) b)) return -1;
366  if (mpz_divisible_p((int_number) b, (int_number) a)) return 1;
367  return 0;
368}
369
370number nrnDiv(number a, number b, const coeffs r)
371{
372  if (a == NULL) a = (number)r->modNumber;
373  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
374  mpz_init(erg);
375  if (mpz_divisible_p((int_number)a, (int_number)b))
376  {
377    mpz_divexact(erg, (int_number)a, (int_number)b);
378    return (number)erg;
379  }
380  else
381  {
382    int_number gcd = (int_number)nrnGcd(a, b, r);
383    mpz_divexact(erg, (int_number)b, gcd);
384    if (!nrnIsUnit((number)erg, r))
385    {
386      WerrorS("Division not possible, even by cancelling zero divisors.");
387      WerrorS("Result is integer division without remainder.");
388      mpz_tdiv_q(erg, (int_number) a, (int_number) b);
389      nrnDelete((number*) &gcd, NULL);
390      return (number)erg;
391    }
392    // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
393    int_number tmp = (int_number)nrnInvers((number) erg,r);
394    mpz_divexact(erg, (int_number)a, gcd);
395    mpz_mul(erg, erg, tmp);
396    nrnDelete((number*) &gcd, NULL);
397    nrnDelete((number*) &tmp, NULL);
398    mpz_mod(erg, erg, r->modNumber);
399    return (number)erg;
400  }
401}
402
403number nrnMod(number a, number b, const coeffs r)
404{
405  /*
406    We need to return the number rr which is uniquely determined by the
407    following two properties:
408      (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
409      (2) There exists some k in the integers Z such that a = k * b + rr.
410    Consider g := gcd(n, |b|). Note that then |b|/g is a unit in Z/n.
411    Now, there are three cases:
412      (a) g = 1
413          Then |b| is a unit in Z/n, i.e. |b| (and also b) divides a.
414          Thus rr = 0.
415      (b) g <> 1 and g divides a
416          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
417      (c) g <> 1 and g does not divide a
418          Then denote the division with remainder of a by g as this:
419          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
420          fulfills (1) and (2), i.e. rr := t is the correct result. Hence
421          in this third case, rr is the remainder of division of a by g in Z.
422     Remark: according to mpz_mod: a,b are always non-negative
423  */
424  int_number g = (int_number)omAllocBin(gmp_nrz_bin);
425  int_number rr = (int_number)omAllocBin(gmp_nrz_bin);
426  mpz_init(g);
427  mpz_init_set_si(rr, 0);
428  mpz_gcd(g, (int_number)r->modNumber, (int_number)b); // g is now as above
429  if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(rr, (int_number)a, g); // the case g <> 1
430  mpz_clear(g);
431  omFreeBin(g, gmp_nrz_bin);
432  return (number)rr;
433}
434
435number nrnIntDiv(number a, number b, const coeffs r)
436{
437  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
438  mpz_init(erg);
439  if (a == NULL) a = (number)r->modNumber;
440  mpz_tdiv_q(erg, (int_number)a, (int_number)b);
441  return (number)erg;
442}
443
444/*
445 * Helper function for computing the module
446 */
447
448int_number nrnMapCoef = NULL;
449
450number nrnMapModN(number from, const coeffs /*src*/, const coeffs dst)
451{
452  return nrnMult(from, (number) nrnMapCoef, dst);
453}
454
455number nrnMap2toM(number from, const coeffs /*src*/, const coeffs dst)
456{
457  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
458  mpz_init(erg);
459  mpz_mul_ui(erg, nrnMapCoef, (NATNUMBER)from);
460  mpz_mod(erg, erg, dst->modNumber);
461  return (number)erg;
462}
463
464number nrnMapZp(number from, const coeffs /*src*/, const coeffs dst)
465{
466  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
467  mpz_init(erg);
468  // TODO: use npInt(...)
469  mpz_mul_si(erg, nrnMapCoef, (NATNUMBER)from);
470  mpz_mod(erg, erg, dst->modNumber);
471  return (number)erg;
472}
473
474number nrnMapGMP(number from, const coeffs /*src*/, const coeffs dst)
475{
476  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
477  mpz_init(erg);
478  mpz_mod(erg, (int_number)from, dst->modNumber);
479  return (number)erg;
480}
481
482number nrnMapQ(number from, const coeffs src, const coeffs dst)
483{
484  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
485  mpz_init(erg);
486  nlGMP(from, (number)erg, src);
487  mpz_mod(erg, erg, dst->modNumber);
488  return (number)erg;
489}
490
491nMapFunc nrnSetMap(const coeffs src, const coeffs dst)
492{
493  /* dst = currRing->cf */
494  if (nCoeff_is_Ring_Z(src))
495  {
496    return nrnMapGMP;
497  }
498  if (nCoeff_is_Q(src))
499  {
500    return nrnMapQ;
501  }
502  // Some type of Z/n ring / field
503  if (nCoeff_is_Ring_ModN(src) || nCoeff_is_Ring_PtoM(src) ||
504      nCoeff_is_Ring_2toM(src) || nCoeff_is_Zp(src))
505  {
506    if (   (!nCoeff_is_Zp(src))
507        && (mpz_cmp(src->modBase, dst->modBase) == 0)
508        && (src->modExponent == dst->modExponent)) return nrnMapGMP;
509    else
510    {
511      int_number nrnMapModul = (int_number) omAllocBin(gmp_nrz_bin);
512      // Computing the n of Z/n
513      if (nCoeff_is_Zp(src))
514      {
515        mpz_init_set_si(nrnMapModul, src->ch);
516      }
517      else
518      {
519        mpz_init(nrnMapModul);
520        mpz_set(nrnMapModul, src->modNumber);
521      }
522      // nrnMapCoef = 1 in dst       if dst is a subring of src
523      // nrnMapCoef = 0 in dst / src if src is a subring of dst
524      if (nrnMapCoef == NULL)
525      {
526        nrnMapCoef = (int_number) omAllocBin(gmp_nrz_bin);
527        mpz_init(nrnMapCoef);
528      }
529      if (mpz_divisible_p(nrnMapModul, dst->modNumber))
530      {
531        mpz_set_si(nrnMapCoef, 1);
532      }
533      else
534      if (nrnDivBy(NULL, (number) nrnMapModul,dst))
535      {
536        mpz_divexact(nrnMapCoef, dst->modNumber, nrnMapModul);
537        int_number tmp = dst->modNumber;
538        dst->modNumber = nrnMapModul;
539        if (!nrnIsUnit((number) nrnMapCoef,dst))
540        {
541          dst->modNumber = tmp;
542          nrnDelete((number*) &nrnMapModul, dst);
543          return NULL;
544        }
545        int_number inv = (int_number) nrnInvers((number) nrnMapCoef,dst);
546        dst->modNumber = tmp;
547        mpz_mul(nrnMapCoef, nrnMapCoef, inv);
548        mpz_mod(nrnMapCoef, nrnMapCoef, dst->modNumber);
549        nrnDelete((number*) &inv, dst);
550      }
551      else
552      {
553        nrnDelete((number*) &nrnMapModul, dst);
554        return NULL;
555      }
556      nrnDelete((number*) &nrnMapModul, dst);
557      if (nCoeff_is_Ring_2toM(src))
558        return nrnMap2toM;
559      else if (nCoeff_is_Zp(src))
560        return nrnMapZp;
561      else
562        return nrnMapModN;
563    }
564  }
565  return NULL;      // default
566}
567
568/*
569 * set the exponent (allocate and init tables) (TODO)
570 */
571
572void nrnSetExp(unsigned long m, coeffs r)
573{
574  /* clean up former stuff */
575  if (r->modNumber != NULL) mpz_clear(r->modNumber);
576
577  r->modExponent= m;
578  r->modNumber = (int_number)omAllocBin(gmp_nrz_bin);
579  mpz_init_set (r->modNumber, r->modBase);
580  mpz_pow_ui (r->modNumber, r->modNumber, m);
581}
582
583/* We expect this ring to be Z/n^m for some m > 0 and for some n > 2 which is not a prime. */
584void nrnInitExp(unsigned long m, coeffs r)
585{
586  nrnSetExp(m, r);
587  assume (r->modNumber != NULL);
588  if (mpz_cmp_ui(r->modNumber,2) <= 0)
589    WarnS("nrnInitExp failed (m in Z/m too small)");
590}
591
592#ifdef LDEBUG
593BOOLEAN nrnDBTest (number a, const char *, const int, const coeffs r)
594{
595  if (a==NULL) return TRUE;
596  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, r->modNumber) > 0) )
597  {
598    return FALSE;
599  }
600  return TRUE;
601}
602#endif
603
604/*2
605* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
606*/
607static const char * nlCPEatLongC(char *s, mpz_ptr i)
608{
609  const char * start=s;
610  if (!(*s >= '0' && *s <= '9'))
611  {
612    mpz_init_set_si(i, 1);
613    return s;
614  }
615  mpz_init(i);
616  while (*s >= '0' && *s <= '9') s++;
617  if (*s=='\0')
618  {
619    mpz_set_str(i,start,10);
620  }
621  else
622  {
623    char c=*s;
624    *s='\0';
625    mpz_set_str(i,start,10);
626    *s=c;
627  }
628  return s;
629}
630
631const char * nrnRead (const char *s, number *a, const coeffs r)
632{
633  int_number z = (int_number) omAllocBin(gmp_nrz_bin);
634  {
635    s = nlCPEatLongC((char *)s, z);
636  }
637  mpz_mod(z, z, r->modNumber);
638  *a = (number) z;
639  return s;
640}
641#endif
642/* #ifdef HAVE_RINGS */
Note: See TracBrowser for help on using the repository browser.