source: git/kernel/rmodulon.cc @ 01d4d3

spielwiese
Last change on this file since 01d4d3 was 01d4d3, checked in by Hans Schönemann <hannes@…>, 14 years ago
debug stuff git-svn-id: file:///usr/local/Singular/svn/trunk@12362 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: numbers modulo n
7*/
8
9#include <string.h>
10#include "mod2.h"
11#include <mylimits.h>
12#include "structs.h"
13#include "febase.h"
14#include "omalloc.h"
15#include "numbers.h"
16#include "longrat.h"
17#include "mpr_complex.h"
18#include "ring.h"
19#include "rmodulon.h"
20#include "si_gmp.h"
21
22#ifdef HAVE_RINGS
23  extern omBin gmp_nrz_bin;
24
25int_number nrnMinusOne = NULL;
26unsigned long nrnExponent = 0;
27
28/*
29 * create a number from int
30 */
31number nrnInit (int i, const ring r)
32{
33  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
34  mpz_init_set_si(erg, i);
35  mpz_mod(erg, erg, r->nrnModul);
36  return (number) erg;
37}
38
39void nrnDelete(number *a, const ring r)
40{
41  if (*a == NULL) return;
42  mpz_clear((int_number) *a);
43  omFreeBin((ADDRESS) *a, gmp_nrz_bin);
44  *a = NULL;
45}
46
47number nrnCopy(number a)
48{
49  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
50  mpz_init_set(erg, (int_number) a);
51  return (number) erg;
52}
53
54number cfrnCopy(number a, const ring r)
55{
56  return nrnCopy(a);
57}
58
59int nrnSize(number a)
60{
61  if (a == NULL) return 0;
62  return sizeof(MP_INT);
63}
64
65/*
66 * convert a number to int (-p/2 .. p/2)
67 */
68int nrnInt(number &n, const ring r)
69{
70  return (int) mpz_get_si( (int_number) n);
71}
72
73/*
74 * Multiply two numbers
75 */
76number nrnMult (number a, number b)
77{
78  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
79  mpz_init(erg);
80  mpz_mul(erg, (int_number) a, (int_number) b);
81  mpz_mod(erg, erg, currRing->nrnModul);
82  return (number) erg;
83}
84
85void nrnPower (number a, int i, number * result)
86{
87  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
88  mpz_init(erg);
89  mpz_powm_ui(erg, (int_number) a, i, currRing->nrnModul);
90  *result = (number) erg;
91}
92
93number nrnAdd (number a, number b)
94{
95  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
96  mpz_init(erg);
97  mpz_add(erg, (int_number) a, (int_number) b);
98  mpz_mod(erg, erg, currRing->nrnModul);
99  return (number) erg;
100}
101
102number nrnSub (number a, number b)
103{
104  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
105  mpz_init(erg);
106  mpz_sub(erg, (int_number) a, (int_number) b);
107  mpz_mod(erg, erg, currRing->nrnModul);
108  return (number) erg;
109}
110
111number nrnNeg (number c)
112{
113// nNeg inplace !!!
114  mpz_sub((int_number) c, currRing->nrnModul, (int_number) c);
115  return c;
116}
117
118number  nrnInvers (number c)
119{
120  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
121  mpz_init(erg);
122  mpz_invert(erg, (int_number) c, currRing->nrnModul);
123  return (number) erg;
124}
125
126/*
127 * Give the smallest non unit k, such that a * x = k = b * y has a solution
128 * TODO: lcm(gcd,gcd) besser als gcd(lcm) ?
129 */
130number nrnLcm (number a,number b,ring r)
131{
132  number erg = nrnGcd(NULL, a, r);
133  number tmp = nrnGcd(NULL, b, r);
134  mpz_lcm((int_number) erg, (int_number) erg, (int_number) tmp);
135  nrnDelete(&tmp, NULL);
136  return (number) erg;
137}
138
139/*
140 * Give the largest non unit k, such that a = x * k, b = y * k has
141 * a solution.
142 */
143number nrnGcd (number a,number b,ring r)
144{
145  if ((a == NULL) && (b == NULL)) return nrnInit(0,r);
146  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
147  mpz_init_set(erg, r->nrnModul);
148  if (a != NULL) mpz_gcd(erg, erg, (int_number) a);
149  if (b != NULL) mpz_gcd(erg, erg, (int_number) b);
150  return (number) erg;
151}
152
153/* Not needed any more, but may have room for improvement
154number nrnGcd3 (number a,number b, number c,ring r)
155{
156  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
157  mpz_init(erg);
158  if (a == NULL) a = (number) r->nrnModul;
159  if (b == NULL) b = (number) r->nrnModul;
160  if (c == NULL) c = (number) r->nrnModul;
161  mpz_gcd(erg, (int_number) a, (int_number) b);
162  mpz_gcd(erg, erg, (int_number) c);
163  mpz_gcd(erg, erg, r->nrnModul);
164  return (number) erg;
165}
166*/
167
168/*
169 * Give the largest non unit k, such that a = x * k, b = y * k has
170 * a solution and r, s, s.t. k = s*a + t*b
171 */
172number  nrnExtGcd (number a, number b, number *s, number *t)
173{
174  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
175  int_number bs = (int_number) omAllocBin(gmp_nrz_bin);
176  int_number bt = (int_number) omAllocBin(gmp_nrz_bin);
177  mpz_init(erg);
178  mpz_init(bs);
179  mpz_init(bt);
180  mpz_gcdext(erg, bs, bt, (int_number) a, (int_number) b);
181  mpz_mod(bs, bs, currRing->nrnModul);
182  mpz_mod(bt, bt, currRing->nrnModul);
183  *s = (number) bs;
184  *t = (number) bt;
185  return (number) erg;
186}
187
188BOOLEAN nrnIsZero (number  a)
189{
190#ifdef LDEBUG
191  if (a==NULL) return FALSE;
192#endif
193  return 0 == mpz_cmpabs_ui((int_number) a, 0);
194}
195
196BOOLEAN nrnIsOne (number a)
197{
198  return 0 == mpz_cmp_si((int_number) a, 1);
199}
200
201BOOLEAN nrnIsMOne (number a)
202{
203  return 0 == mpz_cmp((int_number) a, nrnMinusOne);
204}
205
206BOOLEAN nrnEqual (number a,number b)
207{
208  return 0 == mpz_cmp((int_number) a, (int_number) b);
209}
210
211BOOLEAN nrnGreater (number a,number b)
212{
213  return 0 < mpz_cmp((int_number) a, (int_number) b);
214}
215
216BOOLEAN nrnGreaterZero (number k)
217{
218  return 0 < mpz_cmp_si((int_number) k, 0);
219}
220
221BOOLEAN nrnIsUnit (number a)
222{
223  number tmp = nrnGcd(a, (number) currRing->nrnModul, currRing);
224  bool res = nrnIsOne(tmp);
225  nrnDelete(&tmp, NULL);
226  return res;
227}
228
229number  nrnGetUnit (number k)
230{
231  if (mpz_divisible_p(currRing->nrnModul, (int_number) k)) return nrnInit(1,currRing);
232
233  int_number unit = (int_number) nrnGcd(k, 0, currRing);
234  mpz_tdiv_q(unit, (int_number) k, unit);
235  int_number gcd = (int_number) nrnGcd((number) unit, 0, currRing);
236  if (!nrnIsOne((number) gcd))
237  {
238    int_number ctmp;
239    // tmp := unit^2
240    int_number tmp = (int_number) nrnMult((number) unit,(number) unit);
241    // gcd_new := gcd(tmp, 0)
242    int_number gcd_new = (int_number) nrnGcd((number) tmp, 0, currRing);
243    while (!nrnEqual((number) gcd_new,(number) gcd))
244    {
245      // gcd := gcd_new
246      ctmp = gcd;
247      gcd = gcd_new;
248      gcd_new = ctmp;
249      // tmp := tmp * unit
250      mpz_mul(tmp, tmp, unit);
251      mpz_mod(tmp, tmp, currRing->nrnModul);
252      // gcd_new := gcd(tmp, 0)
253      mpz_gcd(gcd_new, tmp, currRing->nrnModul);
254    }
255    // unit := unit + nrnModul / gcd_new
256    mpz_tdiv_q(tmp, currRing->nrnModul, gcd_new);
257    mpz_add(unit, unit, tmp);
258    mpz_mod(unit, unit, currRing->nrnModul);
259    nrnDelete((number*) &gcd_new, NULL);
260    nrnDelete((number*) &tmp, NULL);
261  }
262  nrnDelete((number*) &gcd, NULL);
263  return (number) unit;
264}
265
266BOOLEAN nrnDivBy (number a,number b)
267{
268  if (a == NULL)
269    return mpz_divisible_p(currRing->nrnModul, (int_number) b);
270  else
271    return mpz_divisible_p((int_number) a, (int_number) b);
272  /*
273  number bs = nrnGcd(a, b, currRing);
274  mpz_tdiv_q((int_number) bs, (int_number) b, (int_number) bs);
275  bool res = nrnIsUnit(bs);
276  nrnDelete(&bs, NULL);
277  return res;
278  */
279}
280
281int nrnDivComp(number a, number b)
282{
283  if (nrnEqual(a, b)) return 0;
284  if (mpz_divisible_p((int_number) a, (int_number) b)) return -1;
285  if (mpz_divisible_p((int_number) b, (int_number) a)) return 1;
286  return 2;
287}
288
289number nrnDiv (number a,number b)
290{
291  if (a == NULL) a = (number) currRing->nrnModul;
292  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
293  mpz_init(erg);
294  if (mpz_divisible_p((int_number) a, (int_number) b))
295  {
296    mpz_divexact(erg, (int_number) a, (int_number) b);
297    return (number) erg;
298  }
299  else
300  {
301    int_number gcd = (int_number) nrnGcd(a, b, currRing);
302    mpz_divexact(erg, (int_number) b, gcd);
303    if (!nrnIsUnit((number) erg))
304    {
305      WerrorS("Division not possible, even by cancelling zero divisors.");
306      WerrorS("Result is integer division without remainder.");
307      mpz_tdiv_q(erg, (int_number) a, (int_number) b);
308      nrnDelete((number*) &gcd, NULL);
309      return (number) erg;
310    }
311    // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
312    int_number tmp = (int_number) nrnInvers((number) erg);
313    mpz_divexact(erg, (int_number) a, gcd);
314    mpz_mul(erg, erg, tmp);
315    nrnDelete((number*) &gcd, NULL);
316    nrnDelete((number*) &tmp, NULL);
317    mpz_mod(erg, erg, currRing->nrnModul);
318    return (number) erg;
319  }
320}
321
322number nrnMod (number a, number b)
323{
324  /*
325    We need to return the number r which is uniquely determined by the
326    following two properties:
327      (1) 0 <= r < |b| (with respect to '<' and '<=' performed in Z x Z)
328      (2) There exists some k in the integers Z such that a = k * b + r.
329    Consider g := gcd(n, |b|). Note that then |b|/g is a unit in Z/n.
330    Now, there are three cases:
331      (a) g = 1
332          Then |b| is a unit in Z/n, i.e. |b| (and also b) divides a.
333          Thus r = 0.
334      (b) g <> 1 and g divides a
335          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again r = 0.
336      (c) g <> 1 and g does not divide a
337          Then denote the division with remainder of a by g as this:
338          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
339          fulfills (1) and (2), i.e. r := t is the correct result. Hence
340          in this third case, r is the remainder of division of a by g in Z.
341     Remark: according to mpz_mod: a,b are always non-negative
342  */
343  int_number g = (int_number) omAllocBin(gmp_nrz_bin);
344  int_number r = (int_number) omAllocBin(gmp_nrz_bin);
345  mpz_init(g);
346  mpz_init_set_si(r,(long)0);
347  mpz_gcd(g, (int_number) currRing->nrnModul, (int_number)b); // g is now as above
348  if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(r, (int_number)a, g); // the case g <> 1
349  mpz_clear(g);
350  omFreeBin(g, gmp_nrz_bin);
351  return (number)r;
352}
353
354number nrnIntDiv (number a,number b)
355{
356  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
357  mpz_init(erg);
358  if (a == NULL) a = (number) currRing->nrnModul;
359  mpz_tdiv_q(erg, (int_number) a, (int_number) b);
360  return (number) erg;
361}
362
363/*
364 * Helper function for computing the module
365 */
366
367int_number nrnMapCoef = NULL;
368
369number nrnMapModN(number from)
370{
371  return nrnMult(from, (number) nrnMapCoef);
372}
373
374number nrnMap2toM(number from)
375{
376  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
377  mpz_init(erg);
378  mpz_mul_ui(erg, nrnMapCoef, (NATNUMBER) from);
379  mpz_mod(erg, erg, currRing->nrnModul);
380  return (number) erg;
381}
382
383number nrnMapZp(number from)
384{
385  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
386  mpz_init(erg);
387  mpz_mul_si(erg, nrnMapCoef, (NATNUMBER) from);
388  mpz_mod(erg, erg, currRing->nrnModul);
389  return (number) erg;
390}
391
392number nrnMapGMP(number from)
393{
394  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
395  mpz_init(erg);
396  mpz_mod(erg, (int_number) from, currRing->nrnModul);
397  return (number) erg;
398}
399
400number nrnMapQ(number from)
401{
402  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
403  mpz_init(erg);
404  nlGMP(from, (number) erg);
405  mpz_mod(erg, erg, currRing->nrnModul);
406  return (number) erg;
407}
408
409nMapFunc nrnSetMap(ring src, ring dst)
410{
411  /* dst = currRing */
412  if (rField_is_Ring_Z(src))
413  {
414    return nrnMapGMP;
415  }
416  if (rField_is_Q(src))
417  {
418    return nrnMapQ;
419  }
420  // Some type of Z/n ring / field
421  if (rField_is_Ring_ModN(src) || rField_is_Ring_PtoM(src) || rField_is_Ring_2toM(src) || rField_is_Zp(src))
422  {
423    if (   (src->ringtype > 0)
424        && (mpz_cmp(src->ringflaga, dst->ringflaga) == 0)
425        && (src->ringflagb == dst->ringflagb)) return nrnMapGMP;
426    else
427    {
428      int_number nrnMapModul = (int_number) omAllocBin(gmp_nrz_bin);
429      // Computing the n of Z/n
430      if (rField_is_Zp(src))
431      {
432        mpz_init_set_si(nrnMapModul, src->ch);
433      }
434      else
435      {
436        mpz_init(nrnMapModul);
437        mpz_set(nrnMapModul, src->ringflaga);
438        mpz_pow_ui(nrnMapModul, nrnMapModul, src->ringflagb);
439      }
440      // nrnMapCoef = 1 in dst       if dst is a subring of src
441      // nrnMapCoef = 0 in dst / src if src is a subring of dst
442      if (nrnMapCoef == NULL)
443      {
444        nrnMapCoef = (int_number) omAllocBin(gmp_nrz_bin);
445        mpz_init(nrnMapCoef);
446      }
447      if (mpz_divisible_p(nrnMapModul, currRing->nrnModul))
448      {
449        mpz_set_si(nrnMapCoef, 1);
450      }
451      else
452      if (nrnDivBy(NULL, (number) nrnMapModul))
453      {
454        mpz_divexact(nrnMapCoef, currRing->nrnModul, nrnMapModul);
455        int_number tmp = currRing->nrnModul;
456        currRing->nrnModul = nrnMapModul;
457        if (!nrnIsUnit((number) nrnMapCoef))
458        {
459          currRing->nrnModul = tmp;
460          nrnDelete((number*) &nrnMapModul, currRing);
461          return NULL;
462        }
463        int_number inv = (int_number) nrnInvers((number) nrnMapCoef);
464        currRing->nrnModul = tmp;
465        mpz_mul(nrnMapCoef, nrnMapCoef, inv);
466        mpz_mod(nrnMapCoef, nrnMapCoef, currRing->nrnModul);
467        nrnDelete((number*) &inv, currRing);
468      }
469      else
470      {
471        nrnDelete((number*) &nrnMapModul, currRing);
472        return NULL;
473      }
474      nrnDelete((number*) &nrnMapModul, currRing);
475      if (rField_is_Ring_2toM(src))
476        return nrnMap2toM;
477      else if (rField_is_Zp(src))
478        return nrnMapZp;
479      else
480        return nrnMapModN;
481    }
482  }
483  return NULL;      // default
484}
485
486/*
487 * set the exponent (allocate and init tables) (TODO)
488 */
489
490void nrnSetExp(int m, ring r)
491{
492  if ((r->nrnModul != NULL) && (mpz_cmp(r->nrnModul, r->ringflaga) == 0) && (nrnExponent == r->ringflagb)) return;
493
494  nrnExponent = r->ringflagb;
495  if (r->nrnModul == NULL)
496  {
497    r->nrnModul = (int_number) omAllocBin(gmp_nrz_bin);
498    mpz_init(r->nrnModul);
499    nrnMinusOne = (int_number) omAllocBin(gmp_nrz_bin);
500    mpz_init(nrnMinusOne);
501  }
502  mpz_set(r->nrnModul, r->ringflaga);
503  mpz_pow_ui(r->nrnModul, r->nrnModul, nrnExponent);
504  mpz_sub_ui(nrnMinusOne, r->nrnModul, 1);
505}
506
507void nrnInitExp(int m, ring r)
508{
509  nrnSetExp(m, r);
510
511  if (mpz_cmp_ui(r->nrnModul,2) <= 0)
512  {
513    WarnS("nrnInitExp failed");
514  }
515}
516
517#ifdef LDEBUG
518BOOLEAN nrnDBTest (number a, const char *f, const int l)
519{
520  if (a==NULL) return TRUE;
521  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, currRing->nrnModul) > 0) )
522  {
523    return FALSE;
524  }
525  return TRUE;
526}
527#endif
528
529void nrnWrite (number &a)
530{
531  char *s,*z;
532  if (a==NULL)
533  {
534    StringAppendS("o");
535  }
536  else
537  {
538    int l=mpz_sizeinbase((int_number) a, 10);
539    if (a->s<2) l=si_max(l,mpz_sizeinbase((int_number) a,10));
540    l+=2;
541    s=(char*)omAlloc(l);
542    z=mpz_get_str(s,10,(int_number) a);
543    StringAppendS(z);
544    omFreeSize((ADDRESS)s,l);
545  }
546}
547
548/*2
549* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
550*/
551static const char * nlCPEatLongC(char *s, MP_INT *i)
552{
553  const char * start=s;
554  if (!(*s >= '0' && *s <= '9'))
555  {
556    mpz_init_set_si(i, 1);
557    return s;
558  }
559  mpz_init(i);
560  while (*s >= '0' && *s <= '9') s++;
561  if (*s=='\0')
562  {
563    mpz_set_str(i,start,10);
564  }
565  else
566  {
567    char c=*s;
568    *s='\0';
569    mpz_set_str(i,start,10);
570    *s=c;
571  }
572  return s;
573}
574
575const char * nrnRead (const char *s, number *a)
576{
577  int_number z = (int_number) omAllocBin(gmp_nrz_bin);
578  {
579    s = nlCPEatLongC((char *)s, z);
580  }
581  mpz_mod(z, z, currRing->nrnModul);
582  *a = (number) z;
583  return s;
584}
585#endif
Note: See TracBrowser for help on using the repository browser.