source: git/kernel/rmodulon.cc @ b1c0a9

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