source: git/kernel/rmodulon.cc @ 528f5b7

spielwiese
Last change on this file since 528f5b7 was 91d286, checked in by Oliver Wienand <wienand@…>, 13 years ago
* Error in component check while creating gcd polys (fixed #300) * Further changed return value of nDivComp to match return values of pDivComp git-svn-id: file:///usr/local/Singular/svn/trunk@13717 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.3 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 <kernel/mod2.h>
11#include <omalloc/mylimits.h>
12#include <kernel/structs.h>
13#include <kernel/febase.h>
14#include <omalloc/omalloc.h>
15#include <kernel/numbers.h>
16#include <kernel/longrat.h>
17#include <kernel/mpr_complex.h>
18#include <kernel/ring.h>
19#include <kernel/rmodulon.h>
20#include <kernel/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(mpz_t);
63}
64
65/*
66 * convert a number to int
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#ifdef LDEBUG
199  if (a == NULL) return FALSE;
200#endif
201  return 0 == mpz_cmp_si((int_number) a, 1);
202}
203
204BOOLEAN nrnIsMOne (number a)
205{
206#ifdef LDEBUG
207  if (a == NULL) return FALSE;
208#endif
209  return 0 == mpz_cmp((int_number) a, nrnMinusOne);
210}
211
212BOOLEAN nrnEqual (number a,number b)
213{
214  return 0 == mpz_cmp((int_number) a, (int_number) b);
215}
216
217BOOLEAN nrnGreater (number a,number b)
218{
219  return 0 < mpz_cmp((int_number) a, (int_number) b);
220}
221
222BOOLEAN nrnGreaterZero (number k)
223{
224  return 0 < mpz_cmp_si((int_number) k, 0);
225}
226
227BOOLEAN nrnIsUnit (number a)
228{
229  number tmp = nrnGcd(a, (number) currRing->nrnModul, currRing);
230  bool res = nrnIsOne(tmp);
231  nrnDelete(&tmp, NULL);
232  return res;
233}
234
235number  nrnGetUnit (number k)
236{
237  if (mpz_divisible_p(currRing->nrnModul, (int_number) k)) return nrnInit(1,currRing);
238
239  int_number unit = (int_number) nrnGcd(k, 0, currRing);
240  mpz_tdiv_q(unit, (int_number) k, unit);
241  int_number gcd = (int_number) nrnGcd((number) unit, 0, currRing);
242  if (!nrnIsOne((number) gcd))
243  {
244    int_number ctmp;
245    // tmp := unit^2
246    int_number tmp = (int_number) nrnMult((number) unit,(number) unit);
247    // gcd_new := gcd(tmp, 0)
248    int_number gcd_new = (int_number) nrnGcd((number) tmp, 0, currRing);
249    while (!nrnEqual((number) gcd_new,(number) gcd))
250    {
251      // gcd := gcd_new
252      ctmp = gcd;
253      gcd = gcd_new;
254      gcd_new = ctmp;
255      // tmp := tmp * unit
256      mpz_mul(tmp, tmp, unit);
257      mpz_mod(tmp, tmp, currRing->nrnModul);
258      // gcd_new := gcd(tmp, 0)
259      mpz_gcd(gcd_new, tmp, currRing->nrnModul);
260    }
261    // unit := unit + nrnModul / gcd_new
262    mpz_tdiv_q(tmp, currRing->nrnModul, gcd_new);
263    mpz_add(unit, unit, tmp);
264    mpz_mod(unit, unit, currRing->nrnModul);
265    nrnDelete((number*) &gcd_new, NULL);
266    nrnDelete((number*) &tmp, NULL);
267  }
268  nrnDelete((number*) &gcd, NULL);
269  return (number) unit;
270}
271
272BOOLEAN nrnDivBy (number a,number b)
273{
274  if (a == NULL)
275    return mpz_divisible_p(currRing->nrnModul, (int_number) b);
276  else
277    return mpz_divisible_p((int_number) a, (int_number) b);
278  /*
279  number bs = nrnGcd(a, b, currRing);
280  mpz_tdiv_q((int_number) bs, (int_number) b, (int_number) bs);
281  bool res = nrnIsUnit(bs);
282  nrnDelete(&bs, NULL);
283  return res;
284  */
285}
286
287int nrnDivComp(number a, number b)
288{
289  if (nrnEqual(a, b)) return 2;
290  if (mpz_divisible_p((int_number) a, (int_number) b)) return -1;
291  if (mpz_divisible_p((int_number) b, (int_number) a)) return 1;
292  return 0;
293}
294
295number nrnDiv (number a,number b)
296{
297  if (a == NULL) a = (number) currRing->nrnModul;
298  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
299  mpz_init(erg);
300  if (mpz_divisible_p((int_number) a, (int_number) b))
301  {
302    mpz_divexact(erg, (int_number) a, (int_number) b);
303    return (number) erg;
304  }
305  else
306  {
307    int_number gcd = (int_number) nrnGcd(a, b, currRing);
308    mpz_divexact(erg, (int_number) b, gcd);
309    if (!nrnIsUnit((number) erg))
310    {
311      WerrorS("Division not possible, even by cancelling zero divisors.");
312      WerrorS("Result is integer division without remainder.");
313      mpz_tdiv_q(erg, (int_number) a, (int_number) b);
314      nrnDelete((number*) &gcd, NULL);
315      return (number) erg;
316    }
317    // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
318    int_number tmp = (int_number) nrnInvers((number) erg);
319    mpz_divexact(erg, (int_number) a, gcd);
320    mpz_mul(erg, erg, tmp);
321    nrnDelete((number*) &gcd, NULL);
322    nrnDelete((number*) &tmp, NULL);
323    mpz_mod(erg, erg, currRing->nrnModul);
324    return (number) erg;
325  }
326}
327
328number nrnMod (number a, number b)
329{
330  /*
331    We need to return the number r which is uniquely determined by the
332    following two properties:
333      (1) 0 <= r < |b| (with respect to '<' and '<=' performed in Z x Z)
334      (2) There exists some k in the integers Z such that a = k * b + r.
335    Consider g := gcd(n, |b|). Note that then |b|/g is a unit in Z/n.
336    Now, there are three cases:
337      (a) g = 1
338          Then |b| is a unit in Z/n, i.e. |b| (and also b) divides a.
339          Thus r = 0.
340      (b) g <> 1 and g divides a
341          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again r = 0.
342      (c) g <> 1 and g does not divide a
343          Then denote the division with remainder of a by g as this:
344          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
345          fulfills (1) and (2), i.e. r := t is the correct result. Hence
346          in this third case, r is the remainder of division of a by g in Z.
347     Remark: according to mpz_mod: a,b are always non-negative
348  */
349  int_number g = (int_number) omAllocBin(gmp_nrz_bin);
350  int_number r = (int_number) omAllocBin(gmp_nrz_bin);
351  mpz_init(g);
352  mpz_init_set_si(r,(long)0);
353  mpz_gcd(g, (int_number) currRing->nrnModul, (int_number)b); // g is now as above
354  if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(r, (int_number)a, g); // the case g <> 1
355  mpz_clear(g);
356  omFreeBin(g, gmp_nrz_bin);
357  return (number)r;
358}
359
360number nrnIntDiv (number a,number b)
361{
362  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
363  mpz_init(erg);
364  if (a == NULL) a = (number) currRing->nrnModul;
365  mpz_tdiv_q(erg, (int_number) a, (int_number) b);
366  return (number) erg;
367}
368
369/*
370 * Helper function for computing the module
371 */
372
373int_number nrnMapCoef = NULL;
374
375number nrnMapModN(number from)
376{
377  return nrnMult(from, (number) nrnMapCoef);
378}
379
380number nrnMap2toM(number from)
381{
382  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
383  mpz_init(erg);
384  mpz_mul_ui(erg, nrnMapCoef, (NATNUMBER) from);
385  mpz_mod(erg, erg, currRing->nrnModul);
386  return (number) erg;
387}
388
389number nrnMapZp(number from)
390{
391  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
392  mpz_init(erg);
393  mpz_mul_si(erg, nrnMapCoef, (NATNUMBER) from);
394  mpz_mod(erg, erg, currRing->nrnModul);
395  return (number) erg;
396}
397
398number nrnMapGMP(number from)
399{
400  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
401  mpz_init(erg);
402  mpz_mod(erg, (int_number) from, currRing->nrnModul);
403  return (number) erg;
404}
405
406number nrnMapQ(number from)
407{
408  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
409  mpz_init(erg);
410  nlGMP(from, (number) erg);
411  mpz_mod(erg, erg, currRing->nrnModul);
412  return (number) erg;
413}
414
415nMapFunc nrnSetMap(const ring src, const ring dst)
416{
417  /* dst = currRing */
418  if (rField_is_Ring_Z(src))
419  {
420    return nrnMapGMP;
421  }
422  if (rField_is_Q(src))
423  {
424    return nrnMapQ;
425  }
426  // Some type of Z/n ring / field
427  if (rField_is_Ring_ModN(src) || rField_is_Ring_PtoM(src) || rField_is_Ring_2toM(src) || rField_is_Zp(src))
428  {
429    if (   (src->ringtype > 0)
430        && (mpz_cmp(src->ringflaga, dst->ringflaga) == 0)
431        && (src->ringflagb == dst->ringflagb)) return nrnMapGMP;
432    else
433    {
434      int_number nrnMapModul = (int_number) omAllocBin(gmp_nrz_bin);
435      // Computing the n of Z/n
436      if (rField_is_Zp(src))
437      {
438        mpz_init_set_si(nrnMapModul, src->ch);
439      }
440      else
441      {
442        mpz_init(nrnMapModul);
443        mpz_set(nrnMapModul, src->ringflaga);
444        mpz_pow_ui(nrnMapModul, nrnMapModul, src->ringflagb);
445      }
446      // nrnMapCoef = 1 in dst       if dst is a subring of src
447      // nrnMapCoef = 0 in dst / src if src is a subring of dst
448      if (nrnMapCoef == NULL)
449      {
450        nrnMapCoef = (int_number) omAllocBin(gmp_nrz_bin);
451        mpz_init(nrnMapCoef);
452      }
453      if (mpz_divisible_p(nrnMapModul, currRing->nrnModul))
454      {
455        mpz_set_si(nrnMapCoef, 1);
456      }
457      else
458      if (nrnDivBy(NULL, (number) nrnMapModul))
459      {
460        mpz_divexact(nrnMapCoef, currRing->nrnModul, nrnMapModul);
461        int_number tmp = currRing->nrnModul;
462        currRing->nrnModul = nrnMapModul;
463        if (!nrnIsUnit((number) nrnMapCoef))
464        {
465          currRing->nrnModul = tmp;
466          nrnDelete((number*) &nrnMapModul, currRing);
467          return NULL;
468        }
469        int_number inv = (int_number) nrnInvers((number) nrnMapCoef);
470        currRing->nrnModul = tmp;
471        mpz_mul(nrnMapCoef, nrnMapCoef, inv);
472        mpz_mod(nrnMapCoef, nrnMapCoef, currRing->nrnModul);
473        nrnDelete((number*) &inv, currRing);
474      }
475      else
476      {
477        nrnDelete((number*) &nrnMapModul, currRing);
478        return NULL;
479      }
480      nrnDelete((number*) &nrnMapModul, currRing);
481      if (rField_is_Ring_2toM(src))
482        return nrnMap2toM;
483      else if (rField_is_Zp(src))
484        return nrnMapZp;
485      else
486        return nrnMapModN;
487    }
488  }
489  return NULL;      // default
490}
491
492/*
493 * set the exponent (allocate and init tables) (TODO)
494 */
495
496void nrnSetExp(int m, ring r)
497{
498  if ((r->nrnModul != NULL) && (mpz_cmp(r->nrnModul, r->ringflaga) == 0) && (nrnExponent == r->ringflagb)) return;
499
500  nrnExponent = r->ringflagb;
501  if (r->nrnModul == NULL)
502  {
503    r->nrnModul = (int_number) omAllocBin(gmp_nrz_bin);
504    mpz_init(r->nrnModul);
505    nrnMinusOne = (int_number) omAllocBin(gmp_nrz_bin);
506    mpz_init(nrnMinusOne);
507  }
508  mpz_set(r->nrnModul, r->ringflaga);
509  mpz_pow_ui(r->nrnModul, r->nrnModul, nrnExponent);
510  mpz_sub_ui(nrnMinusOne, r->nrnModul, 1);
511}
512
513void nrnInitExp(int m, ring r)
514{
515  nrnSetExp(m, r);
516
517  if (mpz_cmp_ui(r->nrnModul,2) <= 0)
518  {
519    WarnS("nrnInitExp failed");
520  }
521}
522
523#ifdef LDEBUG
524BOOLEAN nrnDBTest (number a, const char *f, const int l)
525{
526  if (a==NULL) return TRUE;
527  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, currRing->nrnModul) > 0) )
528  {
529    return FALSE;
530  }
531  return TRUE;
532}
533#endif
534
535/*2
536* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
537*/
538static const char * nlCPEatLongC(char *s, mpz_ptr i)
539{
540  const char * start=s;
541  if (!(*s >= '0' && *s <= '9'))
542  {
543    mpz_init_set_si(i, 1);
544    return s;
545  }
546  mpz_init(i);
547  while (*s >= '0' && *s <= '9') s++;
548  if (*s=='\0')
549  {
550    mpz_set_str(i,start,10);
551  }
552  else
553  {
554    char c=*s;
555    *s='\0';
556    mpz_set_str(i,start,10);
557    *s=c;
558  }
559  return s;
560}
561
562const char * nrnRead (const char *s, number *a)
563{
564  int_number z = (int_number) omAllocBin(gmp_nrz_bin);
565  {
566    s = nlCPEatLongC((char *)s, z);
567  }
568  mpz_mod(z, z, currRing->nrnModul);
569  *a = (number) z;
570  return s;
571}
572#endif
Note: See TracBrowser for help on using the repository browser.