source: git/kernel/rmodulon.cc @ c3e83f

spielwiese
Last change on this file since c3e83f was c3e83f, checked in by Oliver Wienand <wienand@…>, 17 years ago
nrnMapSet liften git-svn-id: file:///usr/local/Singular/svn/trunk@10150 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: rmodulon.cc,v 1.7 2007-06-26 18:34:16 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
21#ifdef HAVE_RINGMODN
22
23NATNUMBER nrnModul;
24
25NATNUMBER XSGCD3(NATNUMBER a, NATNUMBER b, NATNUMBER c)
26{
27  while ((a != 0 && b != 0) || (a != 0 && c != 0) || (b != 0 && c != 0))
28  {
29    if (a > b)
30    {
31      if (b > c) a = a % b;       // a > b > c
32      else
33      {
34        if (a > c) a = a % c;     // a > b, c > b, a > c ==> a > c > b
35        else c = c % a;           // c > a > b
36      }
37    }
38    else
39    {
40      if (a > c) b = b % a;       // a > b > c
41      else
42      {
43        if (b > c) b = b % c;     // a > b, c > b, a > c ==> a > c > b
44        else c = c % b;           // c > a > b
45      }
46    }
47  }
48  return a + b + c;
49}
50
51NATNUMBER XSGCD2(NATNUMBER a, NATNUMBER b)
52{
53  NATNUMBER TMP;
54  while (a != 0 && b != 0)
55  {
56     TMP = a % b;
57     a = b;
58     b = TMP;
59  }
60  return a;
61}
62
63/*
64 * Multiply two numbers
65 */
66number nrnMult (number a,number b)
67{
68  if (((NATNUMBER)a == 0) || ((NATNUMBER)b == 0))
69    return (number)0;
70  else
71    return nrnMultM(a,b);
72}
73
74/*
75 * Give the smallest non unit k, such that a * x = k = b * y has a solution
76 */
77number nrnLcm (number a,number b,ring r)
78{
79  NATNUMBER erg = XSGCD2(nrnModul, (NATNUMBER) a) * XSGCD2(nrnModul, (NATNUMBER) b) / XSGCD3(nrnModul, (NATNUMBER) a, (NATNUMBER) b);
80  if (erg == nrnModul) return NULL;   // Schneller als return erg % nrnModul ?
81  return (number) erg;
82}
83
84/*
85 * Give the largest non unit k, such that a = x * k, b = y * k has
86 * a solution.
87 */
88number nrnGcd (number a,number b,ring r)
89{
90  return (number) XSGCD3(nrnModul, (NATNUMBER) a, (NATNUMBER) b);
91}
92
93void nrnPower (number a, int i, number * result)
94{
95  if (i==0)
96  {
97    //npInit(1,result);
98    *(NATNUMBER *)result = 1;
99  }
100  else if (i==1)
101  {
102    *result = a;
103  }
104  else
105  {
106    nrnPower(a,i-1,result);
107    *result = nrnMultM(a,*result);
108  }
109}
110
111/*
112 * create a number from int
113 */
114number nrnInit (int i)
115{
116  long ii = i;
117  while (ii < 0) ii += nrnModul;
118  while ((ii>1) && (ii >= nrnModul)) ii -= nrnModul;
119  return (number) ii;
120}
121
122/*
123 * convert a number to int (-p/2 .. p/2)
124 */
125int nrnInt(number &n)
126{
127  if ((NATNUMBER)n > (nrnModul >>1)) return (int)((NATNUMBER)n - nrnModul);
128  else return (int)((NATNUMBER)n);
129}
130
131number nrnAdd (number a, number b)
132{
133  return nrnAddM(a,b);
134}
135
136number nrnSub (number a, number b)
137{
138  return nrnSubM(a,b);
139}
140
141number  nrnGetUnit (number k)
142{
143  number unit = nrnIntDiv(k, nrnGcd(k, 0, currRing));
144  number gcd = nrnGcd(unit, 0, currRing);
145  if (!nrnIsOne(gcd))
146  {
147    number tmp = nrnMult(unit, unit);
148    number gcd_new = nrnGcd(tmp, 0, currRing);
149    while (gcd_new != gcd)
150    {
151      gcd = gcd_new;
152      tmp = nrnMult(tmp, unit);
153      gcd_new = nrnGcd(tmp, 0, currRing);
154    }
155    unit = nrnAdd(unit, nrnIntDiv(0, gcd_new));
156  }
157//  Print("k = %d ; unit = %d ; gcd = %d", k, unit, gcd);
158  return unit;
159}
160
161BOOLEAN nrnIsZero (number  a)
162{
163  return 0 == (NATNUMBER)a;
164}
165
166BOOLEAN nrnIsOne (number a)
167{
168  return 1 == (NATNUMBER)a;
169}
170
171BOOLEAN nrnIsUnit (number a)
172{
173  return nrnIsOne(nrnGcd(0, a, currRing));
174}
175
176BOOLEAN nrnIsMOne (number a)
177{
178  return nrnModul == (NATNUMBER)a + 1;
179}
180
181BOOLEAN nrnEqual (number a,number b)
182{
183  return nrnEqualM(a,b);
184}
185
186BOOLEAN nrnGreater (number a,number b)
187{
188  nrnDivBy(a, b);
189}
190
191int nrnComp(number a, number b)
192{
193   NATNUMBER bs = XSGCD2((NATNUMBER) b, nrnModul);
194   NATNUMBER as = XSGCD2((NATNUMBER) a, nrnModul);
195   if (bs == as) return 0;
196   if (as % bs == 0) return -1;
197   if (bs % as == 0) return 1;
198   return 2;
199}
200
201BOOLEAN nrnDivBy (number a,number b)
202{
203  return (XSGCD2((NATNUMBER) b / XSGCD2((NATNUMBER) a, (NATNUMBER) b), nrnModul) == 1);
204}
205
206BOOLEAN nrnGreaterZero (number k)
207{
208  return (((NATNUMBER) k) != 0) && ((NATNUMBER) k <= (nrnModul>>1));
209}
210
211//#ifdef HAVE_DIV_MOD
212#if 1 //ifdef HAVE_NTL // in ntl.a
213//extern void XGCD(long& d, long& s, long& t, long a, long b);
214#include <NTL/ZZ.h>
215#ifdef NTL_CLIENT
216NTL_CLIENT
217#endif
218#else
219void XGCD(long& d, long& s, long& t, long a, long b)
220{
221   long  u, v, u0, v0, u1, v1, u2, v2, q, r;
222
223   long aneg = 0, bneg = 0;
224
225   if (a < 0) {
226      a = -a;
227      aneg = 1;
228   }
229
230   if (b < 0) {
231      b = -b;
232      bneg = 1;
233   }
234
235   u1=1; v1=0;
236   u2=0; v2=1;
237   u = a; v = b;
238
239   while (v != 0) {
240      q = u / v;
241      r = u % v;
242      u = v;
243      v = r;
244      u0 = u2;
245      v0 = v2;
246      u2 =  u1 - q*u2;
247      v2 = v1- q*v2;
248      u1 = u0;
249      v1 = v0;
250   }
251
252   if (aneg)
253      u1 = -u1;
254
255   if (bneg)
256      v1 = -v1;
257
258   d = u;
259   s = u1;
260   t = v1;
261}
262#endif
263//#endif
264
265/*
266 * Give the largest non unit k, such that a = x * k, b = y * k has
267 * a solution and r, s, s.t. k = s*a + t*b
268 */
269number  nrnExtGcd (number a, number b, number *s, number *t)
270{
271  long bs;
272  long bt;
273  long d;
274  XGCD(d, bs, bt, (long) a, (long) b);
275  *s = nrnInit(bs);
276  *t = nrnInit(bt);
277  return (number) d;
278}
279
280NATNUMBER InvModN(NATNUMBER a)
281{
282   long d, s, t;
283
284   // TODO in chain wird XSGCD2 aufgerufen
285   XGCD(d, s, t, a, nrnModul);
286   assume (d == 1);
287   if (s < 0)
288      return s + nrnModul;
289   else
290      return s;
291}
292
293inline number nrnInversM (number c)
294{
295  // Table !!!
296  NATNUMBER inv;
297  inv = InvModN((NATNUMBER)c);
298  return (number) inv;
299}
300
301number nrnDiv (number a,number b)
302{
303  NATNUMBER tmp = XSGCD3(nrnModul, (NATNUMBER) a, (NATNUMBER) b);
304  a = (number) ((NATNUMBER) a / tmp);
305  b = (number) ((NATNUMBER) b / tmp);
306  if (XSGCD2(nrnModul, (NATNUMBER) b) == 1)
307  {
308    return (number) nrnMult(a, nrnInversM(b));
309  }
310  else
311  {
312    WerrorS("div by zero divisor");
313    return (number)0;
314  }
315}
316
317number nrnIntDiv (number a,number b)
318{
319  if ((NATNUMBER)a==0)
320  {
321    if ((NATNUMBER)b==0)
322      return (number) 1;
323    if ((NATNUMBER)b==1)
324      return (number) 0;
325    return (number) ( nrnModul / (NATNUMBER) b);
326  }
327  else
328  {
329    if ((NATNUMBER)b==0)
330      return (number) 0;
331    return (number) ((NATNUMBER) a / (NATNUMBER) b);
332  }
333}
334
335number  nrnInvers (number c)
336{
337  if (XSGCD2(nrnModul, (NATNUMBER) c) != 1)
338  {
339    WerrorS("division by zero divisor");
340    return (number)0;
341  }
342  return nrnInversM(c);
343}
344
345number nrnNeg (number c)
346{
347  if ((NATNUMBER)c==0) return c;
348  return nrnNegM(c);
349}
350
351NATNUMBER nrnMapModul;
352NATNUMBER nrnMapCoef;
353
354number nrnMapModN(number from)
355{
356  NATNUMBER i = (nrnMapCoef * (NATNUMBER) from) % nrnModul;
357  return (number) i;
358}
359
360nMapFunc nrnSetMap(ring src, ring dst)
361{
362  if (rField_is_Ring_ModN(src))
363  {
364    if (src->ringflaga == dst->ringflaga) return ndCopy;
365    else
366    {
367      nrnMapModul = (NATNUMBER) src->ringflaga;
368      if (nrnModul % nrnMapModul == 0)
369      {
370        nrnMapCoef = (nrnModul / nrnMapModul);
371        NATNUMBER tmp = nrnModul;
372        nrnModul = nrnMapModul;
373        nrnMapCoef *= (NATNUMBER) nrnInvers((number) (nrnMapCoef % nrnMapModul));
374        nrnModul = tmp;
375      }
376      else
377        nrnMapCoef = 1;
378      return nrnMapModN;
379    }
380  }
381  return NULL;      /* default */
382}
383
384
385/*
386 * set the exponent (allocate and init tables) (TODO)
387 */
388
389void nrnSetExp(int m, ring r)
390{
391  nrnModul = m;
392//  PrintS("Modul: ");
393//  Print("%d\n", nrnModul);
394}
395
396void nrnInitExp(int m, ring r)
397{
398  nrnModul = m;
399
400  if (m < 2)
401  {
402    WarnS("nInitChar failed");
403  }
404}
405
406#ifdef LDEBUG
407BOOLEAN nrnDBTest (number a, char *f, int l)
408{
409  if (((NATNUMBER)a<0) || ((NATNUMBER)a>nrnModul))
410  {
411    return FALSE;
412  }
413  return TRUE;
414}
415#endif
416
417void nrnWrite (number &a)
418{
419  if ((NATNUMBER)a > (nrnModul >>1)) StringAppend("-%d",(int)(nrnModul-((NATNUMBER)a)));
420  else                          StringAppend("%d",(int)((NATNUMBER)a));
421}
422
423char* nrnEati(char *s, int *i)
424{
425
426  if (((*s) >= '0') && ((*s) <= '9'))
427  {
428    (*i) = 0;
429    do
430    {
431      (*i) *= 10;
432      (*i) += *s++ - '0';
433      if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) % nrnModul;
434    }
435    while (((*s) >= '0') && ((*s) <= '9'));
436    if ((*i) >= nrnModul) (*i) = (*i) % nrnModul;
437  }
438  else (*i) = 1;
439  return s;
440}
441
442char * nrnRead (char *s, number *a)
443{
444  int z;
445  int n=1;
446
447  s = nrnEati(s, &z);
448  if ((*s) == '/')
449  {
450    s++;
451    s = nrnEati(s, &n);
452  }
453  if (n == 1)
454    *a = (number)z;
455  else
456      *a = nrnDiv((number)z,(number)n);
457  return s;
458}
459#endif
Note: See TracBrowser for help on using the repository browser.