source: git/libpolys/coeffs/rmodulon.cc @ 6cce775

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