source: git/kernel/rmodulon.cc @ 8e56ad

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