source: git/kernel/rmodulon.cc @ 12ea9d

fieker-DuValspielwiese
Last change on this file since 12ea9d was 12ea9d, checked in by Oliver Wienand <wienand@…>, 16 years ago
numbers.cc: n^m ringtyp ring.cc: ringausgabe für long long okay ring.h, structs.h: ull für base rmodulon.cc: init, set für llu git-svn-id: file:///usr/local/Singular/svn/trunk@10537 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: rmodulon.cc,v 1.10 2008-01-30 16:14:25 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;
25
26int_number nrnModul = NULL;
27int_number nrnMinusOne = NULL;
28unsigned long nrnExponent = 0;
29unsigned long long nrnBase = 0;
30
31/*
32 * Multiply two numbers
33 */
34number nrnMult (number a, number b)
35{
36  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
37  mpz_init(erg);
38  mpz_mul(erg, (int_number) a, (int_number) b);
39  mpz_mod(erg, erg, nrnModul);
40  return (number) erg;
41}
42
43/*
44 * Give the smallest non unit k, such that a * x = k = b * y has a solution
45 */
46number nrnLcm (number a,number b,ring r)
47{
48  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
49  mpz_init(erg);
50  mpz_lcm(erg, (int_number) a, (int_number) b);
51  mpz_mod(erg, erg, nrnModul);
52  return (number) erg;
53}
54
55/*
56 * Give the largest non unit k, such that a = x * k, b = y * k has
57 * a solution.
58 */
59number nrnGcd (number a,number b,ring r)
60{
61  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
62  mpz_init(erg);
63  mpz_gcd(erg, (int_number) a, (int_number) b);
64  return (number) erg;
65}
66
67/*
68 * Give the largest non unit k, such that a = x * k, b = y * k has
69 * a solution and r, s, s.t. k = s*a + t*b
70 */
71number  nrnExtGcd (number a, number b, number *s, number *t)
72{
73  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
74  int_number bs = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
75  int_number bt = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
76  mpz_init(erg);
77  mpz_init(bs);
78  mpz_init(bt);
79  mpz_gcdext(erg, bs, bt, (int_number) a, (int_number) b);
80  *s = (number) bs;
81  *t = (number) bt;
82  return (number) erg;
83}
84
85void nrnPower (number a, int i, number * result)
86{
87  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
88  mpz_init(erg);
89  mpz_pow_ui(erg, (int_number) a, i);
90  mpz_mod(erg, erg, nrnModul);
91//  int_number tmp = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
92//  mpz_init(tmp);
93//  mpz_clear(erg);
94//  omFree(erg);
95//  erg = tmp;
96  *result = (number) erg;
97}
98
99/*
100 * create a number from int
101 */
102number nrnInit (int i)
103{
104  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
105  mpz_init_set_si(erg, i);
106  return (number) erg;
107}
108
109void nrnDelete(number *a, const ring r)
110{
111  mpz_clear((int_number) *a);
112  omFree((ADDRESS) *a);
113}
114
115/*
116 * convert a number to int (-p/2 .. p/2)
117 */
118int nrnInt(number &n)
119{
120  return (int) mpz_get_si( (int_number) &n);
121}
122
123number nrnAdd (number a, number b)
124{
125  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
126  mpz_init(erg);
127  mpz_add(erg, (int_number) a, (int_number) b);
128  mpz_mod(erg, erg, nrnModul);
129  return (number) erg;
130}
131
132number nrnSub (number a, number b)
133{
134  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
135  mpz_init(erg);
136  mpz_sub(erg, (int_number) a, (int_number) b);
137  mpz_mod(erg, erg, nrnModul);
138  return (number) erg;
139}
140
141number  nrnGetUnit (number a)
142{
143  return nrnDiv((number) nrnModul, nrnGcd(a, (number) nrnModul, NULL));
144}
145
146BOOLEAN nrnIsUnit (number a)
147{
148  return nrnIsOne(nrnGcd(a, (number) nrnModul, NULL));
149}
150
151BOOLEAN nrnIsZero (number  a)
152{
153  return 0 == mpz_cmpabs_ui((int_number) a, 0);
154}
155
156BOOLEAN nrnIsOne (number a)
157{
158  return 0 == mpz_cmp_si((int_number) a, 1);
159}
160
161BOOLEAN nrnIsMOne (number a)
162{
163  return 0 == mpz_cmp((int_number) a, nrnMinusOne);
164}
165
166BOOLEAN nrnEqual (number a,number b)
167{
168  return 0 == mpz_cmp((int_number) a, (int_number) b);
169}
170
171BOOLEAN nrnGreater (number a,number b)
172{
173  return 0 < mpz_cmp((int_number) a, (int_number) b);
174}
175
176BOOLEAN nrnGreaterZero (number k)
177{
178  return 0 <= mpz_cmp_si((int_number) k, 0);
179}
180
181int nrnComp(number a, number b)
182{
183  if (nrnEqual(a, b)) return 0;
184  if (nrnDivBy(a, b)) return -1;
185  if (nrnDivBy(b, a)) return 1;
186  return 2;
187}
188
189BOOLEAN nrnDivBy (number a,number b)
190{
191  return mpz_divisible_p((int_number) a, (int_number) b) != 0;
192}
193
194number nrnDiv (number a,number b)
195{
196  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
197  mpz_init(erg);
198  if (nrnDivBy(a, b))
199  {
200    mpz_divexact(erg, (int_number) a, (int_number) b);
201    return (number) erg;
202  }
203  else
204  {
205    int_number gcd = (int_number) nrnGcd(a, b, NULL);
206    int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
207    mpz_init(erg);
208    mpz_divexact(erg, (int_number) b, gcd);
209    if (!nrnIsUnit((number) erg))
210    {
211      WarnS("Division not possible, even by cancelling zero divisors.");
212      WarnS("Result is zero.");
213      mpz_set_ui(erg, 0);
214      mpz_clear(gcd);
215      omFree(gcd);
216      return (number) erg;
217    }
218    int_number tmp = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
219    mpz_init(tmp);
220    tmp = (int_number) nrnInvers((number) erg);
221    mpz_divexact(erg, (int_number) a, gcd);
222    mpz_mul(erg, erg, tmp);
223    mpz_clear(gcd);
224    omFree(gcd);
225    mpz_clear(tmp);
226    omFree(tmp);
227    mpz_mod(erg, erg, nrnModul);
228    return (number) erg;
229  }
230}
231
232number nrnIntDiv (number a,number b)
233{
234  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
235  mpz_init(erg);
236  mpz_tdiv_q(erg, (int_number) a, (int_number) b);
237  return (number) erg;
238}
239
240number  nrnInvers (number c)
241{
242  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
243  mpz_init(erg);
244  mpz_invert(erg, (int_number) c, nrnModul);
245  return (number) erg;
246}
247
248number nrnNeg (number c)
249{
250  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
251  mpz_init(erg);
252  mpz_sub(erg, nrnModul, (int_number) c);
253  return (number) erg;
254}
255
256nMapFunc nrnSetMap(ring src, ring dst)
257{
258  return NULL;      /* default */
259}
260
261
262/*
263 * set the exponent (allocate and init tables) (TODO)
264 */
265
266void mpz_set_ull(int_number res, unsigned long long xx)
267{
268  unsigned long h = xx >> 32;
269  mpz_set_ui (res, h);
270  mpz_mul_2exp (res, res, 32);
271  mpz_add_ui (res, res, (unsigned long) xx);
272}
273
274void nrnSetExp(int m, ring r)
275{
276  if ((nrnBase == r->ringflaga) && (nrnExponent == r->ringflagb)) return;
277  nrnBase = r->ringflaga;
278  nrnExponent = r->ringflagb;
279  if (nrnModul != NULL)
280  {
281    mpz_clear(nrnModul);
282    omFree(nrnModul);
283    mpz_clear(nrnMinusOne);
284    omFree(nrnMinusOne);
285  }
286  nrnModul = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
287  mpz_init(nrnModul);
288  mpz_set_ull(nrnModul, nrnBase);
289  mpz_pow_ui(nrnModul, nrnModul, nrnExponent);
290
291  nrnMinusOne = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
292  mpz_init(nrnMinusOne);
293  mpz_sub_ui(nrnMinusOne, nrnModul, 1);
294}
295
296void nrnInitExp(int m, ring r)
297{
298  nrnSetExp(m, r);
299  if (mpz_cmp_ui(nrnModul,2) <= 0)
300  {
301    WarnS("nInitChar failed");
302  }
303}
304
305#ifdef LDEBUG
306BOOLEAN nrnDBTest (number a, char *f, int l)
307{
308  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, nrnModul) > 0) )
309  {
310    return FALSE;
311  }
312  return TRUE;
313}
314#endif
315
316void nrnWrite (number &a)
317{
318  char *s,*z;
319  if (a==NULL)
320  {
321    StringAppendS("o");
322  }
323  else
324  {
325    int l=mpz_sizeinbase((int_number) a, 10);
326    if (a->s<2) l=si_max(l,mpz_sizeinbase((int_number) a,10));
327    l+=2;
328    s=(char*)omAlloc(l);
329    z=mpz_get_str(s,10,(int_number) a);
330    StringAppendS(z);
331    omFreeSize((ADDRESS)s,l);
332  }
333}
334
335/*2
336* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
337*/
338char * nlCPEatLongC(char *s, MP_INT *i)
339{
340  char * start=s;
341
342  while (*s >= '0' && *s <= '9') s++;
343  if (*s=='\0')
344  {
345    mpz_set_str(i,start,10);
346  }
347  else
348  {
349    char c=*s;
350    *s='\0';
351    mpz_set_str(i,start,10);
352    *s=c;
353  }
354  return s;
355}
356
357char * nrnRead (char *s, number *a)
358{
359  int_number z = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
360  {
361    mpz_init(z);
362    s = nlCPEatLongC(s, z);
363  }
364  *a = (number) z;
365  return s;
366}
367#endif
Note: See TracBrowser for help on using the repository browser.