# source:git/kernel/rmodulo2m.cc@821a22

spielwiese
Last change on this file since 821a22 was 821a22, checked in by Oliver Wienand <wienand@…>, 16 years ago
*** empty log message *** git-svn-id: file:///usr/local/Singular/svn/trunk@10569 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: rmodulo2m.cc,v 1.15 2008-02-07 13:30:37 wienand Exp \$ */
5/*
6* ABSTRACT: numbers modulo 2^m
7*/
8
9#include <string.h>
10#include "mod2.h"
11
12#ifdef HAVE_RING2TOM
13#include <mylimits.h>
14#include "structs.h"
15#include "febase.h"
16#include "omalloc.h"
17#include "numbers.h"
18#include "longrat.h"
19#include "mpr_complex.h"
20#include "ring.h"
21#include "rmodulo2m.h"
22
23int nr2mExp;
24NATNUMBER nr2mModul;
25
26/*
27 * Multiply two numbers
28 */
29number nr2mMult (number a,number b)
30{
31  if (((NATNUMBER)a == 0) || ((NATNUMBER)b == 0))
32    return (number)0;
33  else
34    return nr2mMultM(a,b);
35}
36
37/*
38 * Give the smallest non unit k, such that a * x = k = b * y has a solution
39 */
40number nr2mLcm (number a,number b,ring r)
41{
42  NATNUMBER res = 0;
43  if ((NATNUMBER) a == 0) a = (number) 1;
44  if ((NATNUMBER) b == 0) b = (number) 1;
45  while ((NATNUMBER) a % 2 == 0)
46  {
47    a = (number) ((NATNUMBER) a / 2);
48    if ((NATNUMBER) b % 2 == 0) b = (number) ((NATNUMBER) b / 2);
49    res++;
50  }
51  while ((NATNUMBER) b % 2 == 0)
52  {
53    b = (number) ((NATNUMBER) b / 2);
54    res++;
55  }
56  return (number) (1L << res);  // (2**res)
57}
58
59/*
60 * Give the largest non unit k, such that a = x * k, b = y * k has
61 * a solution.
62 */
63number nr2mGcd (number a,number b,ring r)
64{
65  NATNUMBER res = 0;
66  if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1;
67  while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
68  {
69    a = (number) ((NATNUMBER) a / 2);
70    b = (number) ((NATNUMBER) b / 2);
71    res++;
72  }
73//  if ((NATNUMBER) b % 2 == 0)
74//  {
75//    return (number) ((1L << res));// * (NATNUMBER) a);  // (2**res)*a    a ist Einheit
76//  }
77//  else
78//  {
79    return (number) ((1L << res));// * (NATNUMBER) b);  // (2**res)*b    b ist Einheit
80//  }
81}
82
83/*
84 * Give the largest non unit k, such that a = x * k, b = y * k has
85 * a solution.
86 */
87number nr2mExtGcd (number a, number b, number *s, number *t)
88{
89  NATNUMBER res = 0;
90  if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1;
91  while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
92  {
93    a = (number) ((NATNUMBER) a / 2);
94    b = (number) ((NATNUMBER) b / 2);
95    res++;
96  }
97  if ((NATNUMBER) b % 2 == 0)
98  {
99    *t = NULL;
100    *s = nr2mInvers(a);
101    return (number) ((1L << res));// * (NATNUMBER) a);  // (2**res)*a    a ist Einheit
102  }
103  else
104  {
105    *s = NULL;
106    *t = nr2mInvers(b);
107    return (number) ((1L << res));// * (NATNUMBER) b);  // (2**res)*b    b ist Einheit
108  }
109}
110
111void nr2mPower (number a, int i, number * result)
112{
113  if (i==0)
114  {
115    //npInit(1,result);
116    *(NATNUMBER *)result = 1;
117  }
118  else if (i==1)
119  {
120    *result = a;
121  }
122  else
123  {
124    nr2mPower(a,i-1,result);
125    *result = nr2mMultM(a,*result);
126  }
127}
128
129/*
130 * create a number from int
131 */
132number nr2mInit (int i)
133{
134  long ii = i;
135  while (ii < 0) ii += nr2mModul;
136  while ((ii>1) && (ii >= nr2mModul)) ii -= nr2mModul;
137  return (number) ii;
138}
139
140/*
141 * convert a number to int (-p/2 .. p/2)
142 */
143int nr2mInt(number &n)
144{
145  if ((NATNUMBER)n > (nr2mModul >>1)) return (int)((NATNUMBER)n - nr2mModul);
146  else return (int)((NATNUMBER)n);
147}
148
149number nr2mAdd (number a, number b)
150{
152}
153
154number nr2mSub (number a, number b)
155{
156  return nr2mSubM(a,b);
157}
158
159BOOLEAN nr2mIsUnit (number a)
160{
161  return ((NATNUMBER) a % 2 == 1);
162}
163
164number  nr2mGetUnit (number k)
165{
166  if (k == NULL)
167    return (number) 1;
168  NATNUMBER tmp = (NATNUMBER) k;
169  while (tmp % 2 == 0)
170    tmp = tmp / 2;
171  return (number) tmp;
172}
173
174BOOLEAN nr2mIsZero (number  a)
175{
176  return 0 == (NATNUMBER)a;
177}
178
179BOOLEAN nr2mIsOne (number a)
180{
181  return 1 == (NATNUMBER)a;
182}
183
184BOOLEAN nr2mIsMOne (number a)
185{
186  return (nr2mModul == (NATNUMBER)a + 1) && (nr2mModul != 2);
187}
188
189BOOLEAN nr2mEqual (number a,number b)
190{
191  return nr2mEqualM(a,b);
192}
193
194BOOLEAN nr2mGreater (number a,number b)
195{
196  return nr2mDivBy(a, b);
197}
198
199BOOLEAN nr2mDivBy (number a,number b)
200{
201  return ((NATNUMBER) a % (NATNUMBER) b) == 0;
202  /*
203  if ((NATNUMBER) a == 0) return TRUE;
204  if ((NATNUMBER) b == 0) return FALSE;
205  while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
206  {
207    a = (number) ((NATNUMBER) a / 2);
208    b = (number) ((NATNUMBER) b / 2);
209}
210  return ((NATNUMBER) b % 2 == 1);
211  */
212}
213
214int nr2mComp(number as, number bs)
215{
216  NATNUMBER a = (NATNUMBER) as;
217  NATNUMBER b = (NATNUMBER) bs;
218  assume(a != 0 && b != 0);
219  while (a % 2 == 0 && b % 2 == 0)
220  {
221    a = a / 2;
222    b = b / 2;
223  }
224  if (a % 2 == 0)
225  {
226    return -1;
227  }
228  else
229  {
230    if (b % 2 == 1)
231    {
232      return 0;
233    }
234    else
235    {
236      return 1;
237    }
238  }
239}
240
241BOOLEAN nr2mGreaterZero (number k)
242{
243  return ((NATNUMBER) k !=0) && ((NATNUMBER) k <= (nr2mModul>>1));
244}
245
246//#ifdef HAVE_DIV_MOD
247#if 1 //ifdef HAVE_NTL // in ntl.a
248//extern void XGCD(long& d, long& s, long& t, long a, long b);
249#include <NTL/ZZ.h>
250#ifdef NTL_CLIENT
251NTL_CLIENT
252#endif
253#else
254void XGCD(long& d, long& s, long& t, long a, long b)
255{
256   long  u, v, u0, v0, u1, v1, u2, v2, q, r;
257
258   long aneg = 0, bneg = 0;
259
260   if (a < 0) {
261      a = -a;
262      aneg = 1;
263   }
264
265   if (b < 0) {
266      b = -b;
267      bneg = 1;
268   }
269
270   u1=1; v1=0;
271   u2=0; v2=1;
272   u = a; v = b;
273
274   while (v != 0) {
275      q = u / v;
276      r = u % v;
277      u = v;
278      v = r;
279      u0 = u2;
280      v0 = v2;
281      u2 =  u1 - q*u2;
282      v2 = v1- q*v2;
283      u1 = u0;
284      v1 = v0;
285   }
286
287   if (aneg)
288      u1 = -u1;
289
290   if (bneg)
291      v1 = -v1;
292
293   d = u;
294   s = u1;
295   t = v1;
296}
297#endif
298
299NATNUMBER InvMod(NATNUMBER a)
300{
301   long d, s, t;
302
303   XGCD(d, s, t, a, nr2mModul);
304   assume (d == 1);
305   if (s < 0)
306      return s + nr2mModul;
307   else
308      return s;
309}
310//#endif
311
312inline number nr2mInversM (number c)
313{
314  // Table !!!
315  NATNUMBER inv;
316  inv = InvMod((NATNUMBER)c);
317  return (number) inv;
318}
319
320number nr2mDiv (number a,number b)
321{
322  if ((NATNUMBER)a==0)
323    return (number)0;
324  else if ((NATNUMBER)b%2==0)
325  {
326    if ((NATNUMBER)b != 0)
327    {
328      while ((NATNUMBER) b%2 == 0 && (NATNUMBER) a%2 == 0)
329      {
330        a = (number) ((NATNUMBER) a / 2);
331        b = (number) ((NATNUMBER) b / 2);
332      }
333    }
334    if ((NATNUMBER) b%2 == 0)
335    {
336      WerrorS("div by zero divisor");
337      return (number)0;
338    }
339  }
340  return (number) nr2mMult(a, nr2mInversM(b));
341}
342
343number nr2mIntDiv (number a,number b)
344{
345  if ((NATNUMBER)a==0)
346  {
347    if ((NATNUMBER)b==0)
348      return (number) 1;
349    if ((NATNUMBER)b==1)
350      return (number) 0;
351    return (number) (nr2mModul / (NATNUMBER) b);
352  }
353  else
354  {
355    if ((NATNUMBER)b==0)
356      return (number) 0;
357    return (number) ((NATNUMBER) a / (NATNUMBER) b);
358  }
359}
360
361number  nr2mInvers (number c)
362{
363  if ((NATNUMBER)c%2==0)
364  {
365    WerrorS("division by zero divisor");
366    return (number)0;
367  }
368  return nr2mInversM(c);
369}
370
371number nr2mNeg (number c)
372{
373  if ((NATNUMBER)c==0) return c;
374  return nr2mNegM(c);
375}
376
377nMapFunc nr2mSetMap(ring src, ring dst)
378{
379  return NULL;      /* default */
380}
381
382
383/*
384 * set the exponent (allocate and init tables) (TODO)
385 */
386
387void nr2mSetExp(int m, ring r)
388{
389  if (m>1)
390  {
391    nr2mExp = m;
392    nr2mModul = 2;
393    for (int i = 1; i < m; i++) {
394      nr2mModul = nr2mModul * 2;
395    }
396  }
397  else
398  {
399    nr2mExp=0;
400    nr2mModul=0;
401  }
402//  PrintS("Modul: ");
403//  Print("%d\n", nr2mModul);
404}
405
406void nr2mInitExp(int m, ring r)
407{
408  int i, w;
409
410  if (m>1)
411  {
412    nr2mExp = m;
413    nr2mModul = 2;
414    for (int i = 1; i < m; i++) {
415      nr2mModul = nr2mModul * 2;
416
417    }
418  }
419  else
420  {
421    WarnS("nInitExp failed");
422  }
423}
424
425#ifdef LDEBUG
426BOOLEAN nr2mDBTest (number a, char *f, int l)
427{
428  if (((NATNUMBER)a<0) || ((NATNUMBER)a>nr2mModul))
429  {
430    return FALSE;
431  }
432  return TRUE;
433}
434#endif
435
436void nr2mWrite (number &a)
437{
438  if ((NATNUMBER)a > (nr2mModul >>1)) StringAppend("-%d",(int)(nr2mModul-((NATNUMBER)a)));
439  else                          StringAppend("%d",(int)((NATNUMBER)a));
440}
441
442char* nr2mEati(char *s, int *i)
443{
444
445  if (((*s) >= '0') && ((*s) <= '9'))
446  {
447    (*i) = 0;
448    do
449    {
450      (*i) *= 10;
451      (*i) += *s++ - '0';
452      if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) % nr2mModul;
453    }
454    while (((*s) >= '0') && ((*s) <= '9'));
455    if ((*i) >= nr2mModul) (*i) = (*i) % nr2mModul;
456  }
457  else (*i) = 1;
458  return s;
459}
460
461char * nr2mRead (char *s, number *a)
462{
463  int z;
464  int n=1;
465
466  s = nr2mEati(s, &z);
467  if ((*s) == '/')
468  {
469    s++;
470    s = nr2mEati(s, &n);
471  }
472  if (n == 1)
473    *a = (number)z;
474  else
475      *a = nr2mDiv((number)z,(number)n);
476  return s;
477}
478#endif
Note: See TracBrowser for help on using the repository browser.