source: git/kernel/rmodulon.cc @ bac8611

spielwiese
Last change on this file since bac8611 was bac8611, checked in by Oliver Wienand <wienand@…>, 16 years ago
ring.h: no simple alloc for modn kutil.cc: delete generated coeffs in pair creation number.cc: copy, delete, size for rmodulon and rintegers kspoly.cc: pSetCoeff durch pSetCoeff0 ersetzt rint*, rmod*: delete, etc. git-svn-id: file:///usr/local/Singular/svn/trunk@10550 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: rmodulon.cc,v 1.15 2008-02-01 15:11:34 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  mpz_clear((int_number) *a);
46  omFreeBin((ADDRESS) *a, gmp_nrn_bin);
47  *a = NULL;
48}
49
50number nrnCopy(number a)
51{
52  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
53  mpz_init_set(erg, (int_number) a);
54  return (number) erg;
55}
56
57number cfrnCopy(number a, const ring r)
58{
59  return nrnCopy(a);
60}
61
62int nrnSize(number a)
63{
64  if (a == NULL) return 0;
65  return sizeof(MP_INT);
66}
67
68/*
69 * convert a number to int (-p/2 .. p/2)
70 */
71int nrnInt(number &n)
72{
73  return (int) mpz_get_si( (int_number) &n);
74}
75
76/*
77 * Multiply two numbers
78 */
79number nrnMult (number a, number b)
80{
81  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
82  mpz_init(erg);
83  mpz_mul(erg, (int_number) a, (int_number) b);
84  mpz_mod(erg, erg, nrnModul);
85  return (number) erg;
86}
87
88void nrnPower (number a, int i, number * result)
89{
90  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
91  mpz_init(erg);
92  mpz_powm_ui(erg, (int_number) a, i, nrnModul);
93  *result = (number) erg;
94}
95
96number nrnAdd (number a, number b)
97{
98  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
99  mpz_init(erg);
100  mpz_add(erg, (int_number) a, (int_number) b);
101  mpz_mod(erg, erg, nrnModul);
102  return (number) erg;
103}
104
105number nrnSub (number a, number b)
106{
107  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
108  mpz_init(erg);
109  mpz_sub(erg, (int_number) a, (int_number) b);
110  mpz_mod(erg, erg, nrnModul);
111  return (number) erg;
112}
113
114number nrnNeg (number c)
115{
116  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
117  mpz_init(erg);
118  mpz_sub(erg, nrnModul, (int_number) c);
119  return (number) erg;
120}
121
122number  nrnInvers (number c)
123{
124  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
125  mpz_init(erg);
126  mpz_invert(erg, (int_number) c, nrnModul);
127  return (number) erg;
128}
129
130/*
131 * Give the smallest non unit k, such that a * x = k = b * y has a solution
132 * TODO: lcm(gcd,gcd) besser als gcd(lcm) ?
133 */
134number nrnLcm (number a,number b,ring r)
135{
136  number erg = nrnGcd(NULL, a, NULL);
137  number tmp = nrnGcd(NULL, b, NULL);
138  mpz_lcm((int_number) erg, (int_number) erg, (int_number) tmp);
139  nrnDelete(&tmp, NULL);
140  return (number) erg;
141}
142
143/*
144 * Give the largest non unit k, such that a = x * k, b = y * k has
145 * a solution.
146 */
147number nrnGcd (number a,number b,ring r)
148{
149  if ((a == NULL) && (b == NULL)) return nrnInit(0);
150  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
151  mpz_init_set(erg, nrnModul);
152  if (a != NULL) mpz_gcd(erg, erg, (int_number) a);
153  if (b != NULL) mpz_gcd(erg, erg, (int_number) b);
154  return (number) erg;
155}
156
157/* Not needed any more, but may have room for improvement
158number nrnGcd3 (number a,number b, number c,ring r)
159{
160  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
161  mpz_init(erg);
162  if (a == NULL) a = (number) nrnModul;
163  if (b == NULL) b = (number) nrnModul;
164  if (c == NULL) c = (number) nrnModul;
165  mpz_gcd(erg, (int_number) a, (int_number) b);
166  mpz_gcd(erg, erg, (int_number) c);
167  mpz_gcd(erg, erg, nrnModul);
168  return (number) erg;
169}
170*/
171
172/*
173 * Give the largest non unit k, such that a = x * k, b = y * k has
174 * a solution and r, s, s.t. k = s*a + t*b
175 */
176number  nrnExtGcd (number a, number b, number *s, number *t)
177{
178  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
179  int_number bs = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
180  int_number bt = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
181  mpz_init(erg);
182  mpz_init(bs);
183  mpz_init(bt);
184  mpz_gcdext(erg, bs, bt, (int_number) a, (int_number) b);
185  mpz_mod(bs, bs, nrnModul);
186  mpz_mod(bt, bt, nrnModul);
187  *s = (number) bs;
188  *t = (number) bt;
189  return (number) erg;
190}
191
192BOOLEAN nrnIsZero (number  a)
193{
194  return 0 == mpz_cmpabs_ui((int_number) a, 0);
195}
196
197BOOLEAN nrnIsOne (number a)
198{
199  return 0 == mpz_cmp_si((int_number) a, 1);
200}
201
202BOOLEAN nrnIsMOne (number a)
203{
204  return 0 == mpz_cmp((int_number) a, nrnMinusOne);
205}
206
207BOOLEAN nrnEqual (number a,number b)
208{
209  return 0 == mpz_cmp((int_number) a, (int_number) b);
210}
211
212BOOLEAN nrnGreater (number a,number b)
213{
214  return 0 < mpz_cmp((int_number) a, (int_number) b);
215}
216
217BOOLEAN nrnGreaterZero (number k)
218{
219  return 0 <= mpz_cmp_si((int_number) k, 0);
220}
221
222BOOLEAN nrnIsUnit (number a)
223{
224  number tmp = nrnGcd(a, (number) nrnModul, NULL);
225  bool res = nrnIsOne(tmp);
226  nrnDelete(&tmp, NULL);
227  return res;
228}
229
230number  nrnGetUnit (number k)
231{
232  int_number unit = (int_number) nrnGcd(k, 0, currRing);
233  mpz_tdiv_q(unit, (int_number) k, unit);
234  int_number gcd = (int_number) nrnGcd((number) unit, 0, currRing);
235  if (!nrnIsOne((number) gcd))
236  {
237    int_number ctmp;
238    // tmp := unit^2
239    int_number tmp = (int_number) nrnMult((number) unit,(number) unit);
240    // gcd_new := gcd(tmp, 0)
241    int_number gcd_new = (int_number) nrnGcd((number) tmp, 0, currRing);
242    while (!nrnEqual((number) gcd_new,(number) gcd))
243    {
244      // gcd := gcd_new
245      ctmp = gcd;
246      gcd = gcd_new;
247      gcd_new = ctmp;
248      // tmp := tmp * unit
249      mpz_mul(tmp, tmp, unit);
250      mpz_mod(tmp, tmp, nrnModul);
251      // gcd_new := gcd(tmp, 0)
252      mpz_gcd(gcd_new, tmp, nrnModul);
253    }
254    // unit := unit + nrnModul / gcd_new
255    mpz_tdiv_q(tmp, nrnModul, gcd_new);
256    mpz_add(unit, unit, tmp);
257    mpz_mod(unit, unit, nrnModul);
258    nrnDelete((number*) &gcd_new, NULL);
259    nrnDelete((number*) &tmp, NULL);
260  }
261  nrnDelete((number*) &gcd, NULL);
262  return (number) unit;
263}
264
265BOOLEAN nrnDivBy (number a,number b)
266{
267  number bs = nrnGcd(a, b, NULL);
268  mpz_tdiv_q((int_number) bs, (int_number) b, (int_number) bs);
269  bool res = nrnIsUnit(bs);
270  nrnDelete(&bs, NULL);
271  return res;
272}
273
274int nrnComp(number a, number b)
275{
276  if (nrnEqual(a, b)) return 0;
277  if (mpz_divisible_p((int_number) a, (int_number) b)) return -1;
278  if (mpz_divisible_p((int_number) b, (int_number) a)) return 1;
279  return 2;
280}
281
282number nrnDiv (number a,number b)
283{
284  if (a == NULL) a = (number) nrnModul;
285  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
286  mpz_init(erg);
287  if (mpz_divisible_p((int_number) a, (int_number) b))
288  {
289    mpz_divexact(erg, (int_number) a, (int_number) b);
290    return (number) erg;
291  }
292  else
293  {
294    int_number gcd = (int_number) nrnGcd(a, b, NULL);
295    mpz_divexact(erg, (int_number) b, gcd);
296    if (!nrnIsUnit((number) erg))
297    {
298      WarnS("Division not possible, even by cancelling zero divisors.");
299      WarnS("Result is zero.");
300      mpz_set_ui(erg, 0);
301      nrnDelete((number*) &gcd, NULL);
302      return (number) erg;
303    }
304    // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
305    int_number tmp = (int_number) nrnInvers((number) erg);
306    mpz_divexact(erg, (int_number) a, gcd);
307    mpz_mul(erg, erg, tmp);
308    nrnDelete((number*) &gcd, NULL);
309    nrnDelete((number*) &tmp, NULL);
310    mpz_mod(erg, erg, nrnModul);
311    return (number) erg;
312  }
313}
314
315number nrnIntDiv (number a,number b)
316{
317  int_number erg = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
318  mpz_init(erg);
319  if (a == NULL) a = (number) nrnModul;
320  mpz_tdiv_q(erg, (int_number) a, (int_number) b);
321  return (number) erg;
322}
323
324nMapFunc nrnSetMap(ring src, ring dst)
325{
326  return NULL;      /* default */
327}
328
329/*
330 * set the exponent (allocate and init tables) (TODO)
331 */
332
333void mpz_set_ull(int_number res, unsigned long long xx)
334{
335  unsigned long h = xx >> 32;
336  mpz_set_ui (res, h);
337  mpz_mul_2exp (res, res, 32);
338  mpz_add_ui (res, res, (unsigned long) xx);
339}
340
341void nrnSetExp(int m, ring r)
342{
343  if ((nrnBase == r->ringflaga) && (nrnExponent == r->ringflagb)) return;
344  nrnBase = r->ringflaga;
345  nrnExponent = r->ringflagb;
346  if (nrnModul != NULL)
347  {
348    nrnDelete((number*) &nrnModul, NULL);
349    nrnDelete((number*) &nrnMinusOne, NULL);
350  }
351  nrnModul = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
352  mpz_init(nrnModul);
353  mpz_set_ull(nrnModul, nrnBase);
354  mpz_pow_ui(nrnModul, nrnModul, nrnExponent);
355
356  nrnMinusOne = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
357  mpz_init(nrnMinusOne);
358  mpz_sub_ui(nrnMinusOne, nrnModul, 1);
359}
360
361void nrnInitExp(int m, ring r)
362{
363  nrnSetExp(m, r);
364  if (mpz_cmp_ui(nrnModul,2) <= 0)
365  {
366    WarnS("nInitChar failed");
367  }
368}
369
370#ifdef LDEBUG
371BOOLEAN nrnDBTest (number a, char *f, int l)
372{
373  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, nrnModul) > 0) )
374  {
375    return FALSE;
376  }
377  return TRUE;
378}
379#endif
380
381void nrnWrite (number &a)
382{
383  char *s,*z;
384  if (a==NULL)
385  {
386    StringAppendS("o");
387  }
388  else
389  {
390    int l=mpz_sizeinbase((int_number) a, 10);
391    if (a->s<2) l=si_max(l,mpz_sizeinbase((int_number) a,10));
392    l+=2;
393    s=(char*)omAlloc(l);
394    z=mpz_get_str(s,10,(int_number) a);
395    StringAppendS(z);
396    omFreeSize((ADDRESS)s,l);
397  }
398}
399
400/*2
401* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
402*/
403char * nlCPEatLongC(char *s, MP_INT *i)
404{
405  char * start=s;
406  if (!(*s >= '0' && *s <= '9'))
407  {
408    mpz_init_set_si(i, 1);
409    return s;
410  }
411  mpz_init(i);
412  while (*s >= '0' && *s <= '9') s++;
413  if (*s=='\0')
414  {
415    mpz_set_str(i,start,10);
416  }
417  else
418  {
419    char c=*s;
420    *s='\0';
421    mpz_set_str(i,start,10);
422    *s=c;
423  }
424  return s;
425}
426
427char * nrnRead (char *s, number *a)
428{
429  int_number z = (int_number) omAllocBin(gmp_nrn_bin); // evtl. spaeter mit bin
430  {
431    s = nlCPEatLongC(s, z);
432  }
433  *a = (number) z;
434  return s;
435}
436#endif
Note: See TracBrowser for help on using the repository browser.