source: git/kernel/rmodulon.cc @ 093f30e

spielwiese
Last change on this file since 093f30e was 093f30e, checked in by Oliver Wienand <wienand@…>, 16 years ago
kstd1.cc, kstd2.cc: redRing2toM --> redRing kutil.*: extendedspoly unit case p_Mult_q.cc: assume domain for bucket mult pp_Mult_mm__T.cc: free tmp memory rmodulo2m.cc rmodulon.cc git-svn-id: file:///usr/local/Singular/svn/trunk@10573 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: rmodulon.cc,v 1.22 2008-02-08 10:11:30 wienand 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 "gmp.h"
21
22#ifdef HAVE_RINGMODN
23
24typedef MP_INT *int_number;
25omBin gmp_nrn_bin = omGetSpecBin(sizeof(MP_INT));
26
27int_number nrnModul = NULL;
28int_number nrnMinusOne = NULL;
29unsigned long nrnExponent = 0;
30unsigned long long nrnBase = 0;
31
32/*
33 * create a number from int
34 */
35number nrnInit (int i)
36{
37  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
38  mpz_init_set_si(erg, i);
39  mpz_mod(erg, erg, nrnModul);
40  return (number) erg;
41}
42
43void nrnDelete(number *a, const ring r)
44{
45  if (*a == NULL) return;
46  mpz_clear((int_number) *a);
47  omFreeBin((ADDRESS) *a, gmp_nrn_bin);
48  *a = NULL;
49}
50
51number nrnCopy(number a)
52{
53  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
54  mpz_init_set(erg, (int_number) a);
55  return (number) erg;
56}
57
58number cfrnCopy(number a, const ring r)
59{
60  return nrnCopy(a);
61}
62
63int nrnSize(number a)
64{
65  if (a == NULL) return 0;
66  return sizeof(MP_INT);
67}
68
69/*
70 * convert a number to int (-p/2 .. p/2)
71 */
72int nrnInt(number &n)
73{
74  return (int) mpz_get_si( (int_number) &n);
75}
76
77/*
78 * Multiply two numbers
79 */
80number nrnMult (number a, number b)
81{
82  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
83  mpz_init(erg);
84  mpz_mul(erg, (int_number) a, (int_number) b);
85  mpz_mod(erg, erg, nrnModul);
86  return (number) erg;
87}
88
89void nrnPower (number a, int i, number * result)
90{
91  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
92  mpz_init(erg);
93  mpz_powm_ui(erg, (int_number) a, i, nrnModul);
94  *result = (number) erg;
95}
96
97number nrnAdd (number a, number b)
98{
99  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
100  mpz_init(erg);
101  mpz_add(erg, (int_number) a, (int_number) b);
102  mpz_mod(erg, erg, nrnModul);
103  return (number) erg;
104}
105
106number nrnSub (number a, number b)
107{
108  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
109  mpz_init(erg);
110  mpz_sub(erg, (int_number) a, (int_number) b);
111  mpz_mod(erg, erg, nrnModul);
112  return (number) erg;
113}
114
115number nrnNeg (number c)
116{
117// nNeg inplace !!! TODO
118//  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
119//  mpz_init(erg);
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); // evtl. spaeter mit 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); // evtl. spaeter mit 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); // evtl. spaeter mit 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); // evtl. spaeter mit bin
181  int_number bs = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
182  int_number bt = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit 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 nrnComp(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); // evtl. spaeter mit 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 zero.");
310      mpz_set_ui(erg, 0);
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); // evtl. spaeter mit 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
334nMapFunc nrnSetMap(ring src, ring dst)
335{
336  return NULL;      /* default */
337}
338
339/*
340 * set the exponent (allocate and init tables) (TODO)
341 */
342
343void mpz_set_ull(int_number res, unsigned long long xx)
344{
345  unsigned long h = xx >> 32;
346  mpz_set_ui (res, h);
347  mpz_mul_2exp (res, res, 32);
348  mpz_add_ui (res, res, (unsigned long) xx);
349}
350
351void nrnSetExp(int m, ring r)
352{
353  if ((nrnBase == r->ringflaga) && (nrnExponent == r->ringflagb)) return;
354  nrnBase = r->ringflaga;
355  nrnExponent = r->ringflagb;
356  if (nrnModul == NULL)
357  {
358    nrnModul = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
359    mpz_init(nrnModul);
360    nrnMinusOne = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
361    mpz_init(nrnMinusOne);
362  }
363  mpz_set_ull(nrnModul, nrnBase);
364  mpz_pow_ui(nrnModul, nrnModul, nrnExponent);
365  mpz_sub_ui(nrnMinusOne, nrnModul, 1);
366}
367
368void nrnInitExp(int m, ring r)
369{
370  nrnSetExp(m, r);
371
372  if (mpz_cmp_ui(nrnModul,2) <= 0)
373  {
374    WarnS("nrnInitExp failed");
375  }
376}
377
378#ifdef LDEBUG
379BOOLEAN nrnDBTest (number a, char *f, int l)
380{
381  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, nrnModul) > 0) )
382  {
383    return FALSE;
384  }
385  return TRUE;
386}
387#endif
388
389void nrnWrite (number &a)
390{
391  char *s,*z;
392  if (a==NULL)
393  {
394    StringAppendS("o");
395  }
396  else
397  {
398    int l=mpz_sizeinbase((int_number) a, 10);
399    if (a->s<2) l=si_max(l,mpz_sizeinbase((int_number) a,10));
400    l+=2;
401    s=(char*)omAlloc(l);
402    z=mpz_get_str(s,10,(int_number) a);
403    StringAppendS(z);
404    omFreeSize((ADDRESS)s,l);
405  }
406}
407
408/*2
409* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
410*/
411char * nlCPEatLongC(char *s, MP_INT *i)
412{
413  char * start=s;
414  if (!(*s >= '0' && *s <= '9'))
415  {
416    mpz_init_set_si(i, 1);
417    return s;
418  }
419  mpz_init(i);
420  while (*s >= '0' && *s <= '9') s++;
421  if (*s=='\0')
422  {
423    mpz_set_str(i,start,10);
424  }
425  else
426  {
427    char c=*s;
428    *s='\0';
429    mpz_set_str(i,start,10);
430    *s=c;
431  }
432  return s;
433}
434
435char * nrnRead (char *s, number *a)
436{
437  int_number z = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
438  {
439    s = nlCPEatLongC(s, z);
440  }
441  *a = (number) z;
442  return s;
443}
444#endif
Note: See TracBrowser for help on using the repository browser.