# source:git/kernel/rmodulo2m.cc@c90b43

spielwiese
Last change on this file since c90b43 was c90b43, checked in by Hans Schönemann <hannes@…>, 14 years ago
*hannes: HAVE_RINGS only git-svn-id: file:///usr/local/Singular/svn/trunk@11781 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 10.0 KB
Line
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* \$Id: rmodulo2m.cc,v 1.24 2009-05-06 12:53:49 Singular Exp \$ */
5/*
6* ABSTRACT: numbers modulo 2^m
7*/
8
9#include <string.h>
10#include "mod2.h"
11
12#ifdef HAVE_RINGS
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#include "si_gmp.h"
23
24int nr2mExp;
25NATNUMBER nr2mModul;
26
27/*
28 * Multiply two numbers
29 */
30number nr2mMult (number a,number b)
31{
32  if (((NATNUMBER)a == 0) || ((NATNUMBER)b == 0))
33    return (number)0;
34  else
35    return nr2mMultM(a,b);
36}
37
38/*
39 * Give the smallest non unit k, such that a * x = k = b * y has a solution
40 */
41number nr2mLcm (number a,number b,ring r)
42{
43  NATNUMBER res = 0;
44  if ((NATNUMBER) a == 0) a = (number) 1;
45  if ((NATNUMBER) b == 0) b = (number) 1;
46  while ((NATNUMBER) a % 2 == 0)
47  {
48    a = (number) ((NATNUMBER) a / 2);
49    if ((NATNUMBER) b % 2 == 0) b = (number) ((NATNUMBER) b / 2);
50    res++;
51  }
52  while ((NATNUMBER) b % 2 == 0)
53  {
54    b = (number) ((NATNUMBER) b / 2);
55    res++;
56  }
57  return (number) (1L << res);  // (2**res)
58}
59
60/*
61 * Give the largest non unit k, such that a = x * k, b = y * k has
62 * a solution.
63 */
64number nr2mGcd (number a,number b,ring r)
65{
66  NATNUMBER res = 0;
67  if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1;
68  while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
69  {
70    a = (number) ((NATNUMBER) a / 2);
71    b = (number) ((NATNUMBER) b / 2);
72    res++;
73  }
74//  if ((NATNUMBER) b % 2 == 0)
75//  {
76//    return (number) ((1L << res));// * (NATNUMBER) a);  // (2**res)*a    a ist Einheit
77//  }
78//  else
79//  {
80    return (number) ((1L << res));// * (NATNUMBER) b);  // (2**res)*b    b ist Einheit
81//  }
82}
83
84/*
85 * Give the largest non unit k, such that a = x * k, b = y * k has
86 * a solution.
87 */
88number nr2mExtGcd (number a, number b, number *s, number *t)
89{
90  NATNUMBER res = 0;
91  if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1;
92  while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
93  {
94    a = (number) ((NATNUMBER) a / 2);
95    b = (number) ((NATNUMBER) b / 2);
96    res++;
97  }
98  if ((NATNUMBER) b % 2 == 0)
99  {
100    *t = NULL;
101    *s = nr2mInvers(a);
102    return (number) ((1L << res));// * (NATNUMBER) a);  // (2**res)*a    a ist Einheit
103  }
104  else
105  {
106    *s = NULL;
107    *t = nr2mInvers(b);
108    return (number) ((1L << res));// * (NATNUMBER) b);  // (2**res)*b    b ist Einheit
109  }
110}
111
112void nr2mPower (number a, int i, number * result)
113{
114  if (i==0)
115  {
116    //npInit(1,result);
117    *(NATNUMBER *)result = 1;
118  }
119  else if (i==1)
120  {
121    *result = a;
122  }
123  else
124  {
125    nr2mPower(a,i-1,result);
126    *result = nr2mMultM(a,*result);
127  }
128}
129
130/*
131 * create a number from int
132 */
133number nr2mInit (int i)
134{
135  long ii = i;
136  while (ii < 0) ii += nr2mModul;
137  while ((ii>1) && (ii >= nr2mModul)) ii -= nr2mModul;
138  return (number) ii;
139}
140
141/*
142 * convert a number to int (-p/2 .. p/2)
143 */
144int nr2mInt(number &n)
145{
146  if ((NATNUMBER)n > (nr2mModul >>1)) return (int)((NATNUMBER)n - nr2mModul);
147  else return (int)((NATNUMBER)n);
148}
149
150number nr2mAdd (number a, number b)
151{
153}
154
155number nr2mSub (number a, number b)
156{
157  return nr2mSubM(a,b);
158}
159
160BOOLEAN nr2mIsUnit (number a)
161{
162  return ((NATNUMBER) a % 2 == 1);
163}
164
165number  nr2mGetUnit (number k)
166{
167  if (k == NULL)
168    return (number) 1;
169  NATNUMBER tmp = (NATNUMBER) k;
170  while (tmp % 2 == 0)
171    tmp = tmp / 2;
172  return (number) tmp;
173}
174
175BOOLEAN nr2mIsZero (number  a)
176{
177  return 0 == (NATNUMBER)a;
178}
179
180BOOLEAN nr2mIsOne (number a)
181{
182  return 1 == (NATNUMBER)a;
183}
184
185BOOLEAN nr2mIsMOne (number a)
186{
187  return (nr2mModul == (NATNUMBER)a + 1) && (nr2mModul != 2);
188}
189
190BOOLEAN nr2mEqual (number a,number b)
191{
192  return nr2mEqualM(a,b);
193}
194
195BOOLEAN nr2mGreater (number a,number b)
196{
197  return nr2mDivBy(a, b);
198}
199
200BOOLEAN nr2mDivBy (number a,number b)
201{
202  if (a == NULL)
203    return (nr2mModul % (NATNUMBER) b) == 0;
204  else
205    return ((NATNUMBER) a % (NATNUMBER) b) == 0;
206  /*
207  if ((NATNUMBER) a == 0) return TRUE;
208  if ((NATNUMBER) b == 0) return FALSE;
209  while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
210  {
211    a = (number) ((NATNUMBER) a / 2);
212    b = (number) ((NATNUMBER) b / 2);
213}
214  return ((NATNUMBER) b % 2 == 1);
215  */
216}
217
218int nr2mDivComp(number as, number bs)
219{
220  NATNUMBER a = (NATNUMBER) as;
221  NATNUMBER b = (NATNUMBER) bs;
222  assume(a != 0 && b != 0);
223  while (a % 2 == 0 && b % 2 == 0)
224  {
225    a = a / 2;
226    b = b / 2;
227  }
228  if (a % 2 == 0)
229  {
230    return -1;
231  }
232  else
233  {
234    if (b % 2 == 1)
235    {
236      return 0;
237    }
238    else
239    {
240      return 1;
241    }
242  }
243}
244
245BOOLEAN nr2mGreaterZero (number k)
246{
247  return ((NATNUMBER) k !=0) && ((NATNUMBER) k <= (nr2mModul>>1));
248}
249
250//#ifdef HAVE_DIV_MOD
251#if 1 //ifdef HAVE_NTL // in ntl.a
252//extern void XGCD(long& d, long& s, long& t, long a, long b);
253#include <NTL/ZZ.h>
254#ifdef NTL_CLIENT
255NTL_CLIENT
256#endif
257#else
258void XGCD(long& d, long& s, long& t, long a, long b)
259{
260   long  u, v, u0, v0, u1, v1, u2, v2, q, r;
261
262   long aneg = 0, bneg = 0;
263
264   if (a < 0) {
265      a = -a;
266      aneg = 1;
267   }
268
269   if (b < 0) {
270      b = -b;
271      bneg = 1;
272   }
273
274   u1=1; v1=0;
275   u2=0; v2=1;
276   u = a; v = b;
277
278   while (v != 0) {
279      q = u / v;
280      r = u % v;
281      u = v;
282      v = r;
283      u0 = u2;
284      v0 = v2;
285      u2 =  u1 - q*u2;
286      v2 = v1- q*v2;
287      u1 = u0;
288      v1 = v0;
289   }
290
291   if (aneg)
292      u1 = -u1;
293
294   if (bneg)
295      v1 = -v1;
296
297   d = u;
298   s = u1;
299   t = v1;
300}
301#endif
302
303NATNUMBER InvMod(NATNUMBER a)
304{
305   long d, s, t;
306
307   XGCD(d, s, t, a, nr2mModul);
308   assume (d == 1);
309   if (s < 0)
310      return s + nr2mModul;
311   else
312      return s;
313}
314//#endif
315
316inline number nr2mInversM (number c)
317{
318  // Table !!!
319  NATNUMBER inv;
320  inv = InvMod((NATNUMBER)c);
321  return (number) inv;
322}
323
324number nr2mDiv (number a,number b)
325{
326  if ((NATNUMBER)a==0)
327    return (number)0;
328  else if ((NATNUMBER)b%2==0)
329  {
330    if ((NATNUMBER)b != 0)
331    {
332      while ((NATNUMBER) b%2 == 0 && (NATNUMBER) a%2 == 0)
333      {
334        a = (number) ((NATNUMBER) a / 2);
335        b = (number) ((NATNUMBER) b / 2);
336      }
337    }
338    if ((NATNUMBER) b%2 == 0)
339    {
340      WarnS("Division not possible, even by cancelling zero divisors.");
341      WarnS("Result is integer division without remainder.");
342      return (number) ((NATNUMBER) a / (NATNUMBER) b);
343    }
344  }
345  return (number) nr2mMult(a, nr2mInversM(b));
346}
347
348number nr2mIntDiv (number a,number b)
349{
350  if ((NATNUMBER)a==0)
351  {
352    if ((NATNUMBER)b==0)
353      return (number) 1;
354    if ((NATNUMBER)b==1)
355      return (number) 0;
356    return (number) (nr2mModul / (NATNUMBER) b);
357  }
358  else
359  {
360    if ((NATNUMBER)b==0)
361      return (number) 0;
362    return (number) ((NATNUMBER) a / (NATNUMBER) b);
363  }
364}
365
366number  nr2mInvers (number c)
367{
368  if ((NATNUMBER)c%2==0)
369  {
370    WerrorS("division by zero divisor");
371    return (number)0;
372  }
373  return nr2mInversM(c);
374}
375
376number nr2mNeg (number c)
377{
378  if ((NATNUMBER)c==0) return c;
379  return nr2mNegM(c);
380}
381
382number nr2mMapMachineInt(number from)
383{
384  NATNUMBER i = ((NATNUMBER) from) % nr2mModul;
385  return (number) i;
386}
387
388number nr2mMapZp(number from)
389{
390  long ii = (long) from;
391  while (ii < 0) ii += nr2mModul;
392  while ((ii>1) && (ii >= nr2mModul)) ii -= nr2mModul;
393  return (number) ii;
394}
395
396number nr2mMapQ(number from)
397{
398  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
399  mpz_init(erg);
400
401  nlGMP(from, (number) erg);
402  mpz_mod_ui(erg, erg, nr2mModul);
403  number r = (number) mpz_get_ui(erg);
404
405  mpz_clear(erg);
407  return (number) r;
408}
409
410number nr2mMapGMP(number from)
411{
412  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
413  mpz_init(erg);
414
415  mpz_mod_ui(erg, (int_number) from, nr2mModul);
416  number r = (number) mpz_get_ui(erg);
417
418  mpz_clear(erg);
420  return (number) r;
421}
422
423nMapFunc nr2mSetMap(ring src, ring dst)
424{
425  if (rField_is_Ring_2toM(src)
426     && (src->ringflagb >= dst->ringflagb))
427  {
428    return nr2mMapMachineInt;
429  }
430  if (rField_is_Ring_Z(src))
431  {
432    return nr2mMapGMP;
433  }
434  if (rField_is_Q(src))
435  {
436    return nr2mMapQ;
437  }
438  if (rField_is_Zp(src)
439     && (src->ch == 2)
440     && (dst->ringflagb == 1))
441  {
442    return nr2mMapZp;
443  }
444  if (rField_is_Ring_PtoM(src) || rField_is_Ring_ModN(src))
445  {
446    // Computing the n of Z/n
447    int_number modul = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
448    mpz_init(modul);
449    mpz_set(modul, src->ringflaga);
450    mpz_pow_ui(modul, modul, src->ringflagb);
451    if (mpz_divisible_2exp_p(modul, dst->ringflagb))
452    {
453      mpz_clear(modul);
455      return nr2mMapGMP;
456    }
457    mpz_clear(modul);
459  }
460  return NULL;      // default
461}
462
463/*
464 * set the exponent (allocate and init tables) (TODO)
465 */
466
467void nr2mSetExp(int m, ring r)
468{
469  if (m>1)
470  {
471    nr2mExp = m;
472    nr2mModul = 2;
473    for (int i = 1; i < m; i++) {
474      nr2mModul = nr2mModul * 2;
475    }
476  }
477  else
478  {
479    nr2mExp=2;
480    nr2mModul=4;
481  }
482}
483
484void nr2mInitExp(int m, ring r)
485{
486  nr2mSetExp(m, r);
487  if (m<2) WarnS("nInitExp failed: using Z/2^2");
488}
489
490#ifdef LDEBUG
491BOOLEAN nr2mDBTest (number a, const char *f, const int l)
492{
493  if (((NATNUMBER)a<0) || ((NATNUMBER)a>nr2mModul))
494  {
495    return FALSE;
496  }
497  return TRUE;
498}
499#endif
500
501void nr2mWrite (number &a)
502{
503  if ((NATNUMBER)a > (nr2mModul >>1)) StringAppend("-%d",(int)(nr2mModul-((NATNUMBER)a)));
504  else                          StringAppend("%d",(int)((NATNUMBER)a));
505}
506
507static const char* nr2mEati(const char *s, int *i)
508{
509
510  if (((*s) >= '0') && ((*s) <= '9'))
511  {
512    (*i) = 0;
513    do
514    {
515      (*i) *= 10;
516      (*i) += *s++ - '0';
517      if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) % nr2mModul;
518    }
519    while (((*s) >= '0') && ((*s) <= '9'));
520    if ((*i) >= nr2mModul) (*i) = (*i) % nr2mModul;
521  }
522  else (*i) = 1;
523  return s;
524}
525
526const char * nr2mRead (const char *s, number *a)
527{
528  int z;
529  int n=1;
530
531  s = nr2mEati(s, &z);
532  if ((*s) == '/')
533  {
534    s++;
535    s = nr2mEati(s, &n);
536  }
537  if (n == 1)
538    *a = (number)z;
539  else
540      *a = nr2mDiv((number)z,(number)n);
541  return s;
542}
543#endif
Note: See TracBrowser for help on using the repository browser.