# source:git/kernel/rmodulo2m.cc@6ea941

spielwiese
Last change on this file since 6ea941 was 6ea941, checked in by Frank Seelisch <seelisch@…>, 14 years ago
corrected a modulo b in Z/n and Z/2^m git-svn-id: file:///usr/local/Singular/svn/trunk@11940 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 11.4 KB
Line
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* \$Id: rmodulo2m.cc,v 1.25 2009-07-03 13:14:10 seelisch 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 nr2mMod (number a, number b)
349{
350  /*
351    We need to return the number r which is uniquely determined by the
352    following two properties:
353      (1) 0 <= r < |b| (with respect to '<' and '<=' performed in Z x Z)
354      (2) There exists some k in the integers Z such that a = k * b + r.
355    Consider g := gcd(2^m, |b|). Note that then |b|/g is a unit in Z/2^m.
356    Now, there are three cases:
357      (a) g = 1
358          Then |b| is a unit in Z/2^m, i.e. |b| (and also b) divides a.
359          Thus r = 0.
360      (b) g <> 1 and g divides a
361          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again r = 0.
362      (c) g <> 1 and g does not divide a
363          Let's denote the division with remainder of a by g as follows:
364          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
365          fulfills (1) and (2), i.e. r := t is the correct result. Hence
366          in this third case, r is the remainder of division of a by g in Z.
367    This algorithm is the same as for the case Z/n, except that we may
368    compute the gcd of |b| and 2^m "by hand": We just extract the highest
369    power of 2 (<= 2^m) that is contained in b.
370  */
371  NATNUMBER g = 1;
372  NATNUMBER b_div = b;
373  if (b_div < 0) b_div = - b_div; // b_div now represents |b|
374  NATNUMBER r = 0;
375  while ((g < nr2mModul) && (b_div > 0) && (b_div % 2 == 0))
376  {
377    b_div = b_div >> 1;
378    g = g << 1;
379  } // g is now the gcd of 2^m and |b|
380
381  if (g != 1) r = (NATNUMBER)a % g;
382  return (number)r;
383}
384
385number nr2mIntDiv (number a,number b)
386{
387  if ((NATNUMBER)a==0)
388  {
389    if ((NATNUMBER)b==0)
390      return (number) 1;
391    if ((NATNUMBER)b==1)
392      return (number) 0;
393    return (number) (nr2mModul / (NATNUMBER) b);
394  }
395  else
396  {
397    if ((NATNUMBER)b==0)
398      return (number) 0;
399    return (number) ((NATNUMBER) a / (NATNUMBER) b);
400  }
401}
402
403number  nr2mInvers (number c)
404{
405  if ((NATNUMBER)c%2==0)
406  {
407    WerrorS("division by zero divisor");
408    return (number)0;
409  }
410  return nr2mInversM(c);
411}
412
413number nr2mNeg (number c)
414{
415  if ((NATNUMBER)c==0) return c;
416  return nr2mNegM(c);
417}
418
419number nr2mMapMachineInt(number from)
420{
421  NATNUMBER i = ((NATNUMBER) from) % nr2mModul;
422  return (number) i;
423}
424
425number nr2mMapZp(number from)
426{
427  long ii = (long) from;
428  while (ii < 0) ii += nr2mModul;
429  while ((ii>1) && (ii >= nr2mModul)) ii -= nr2mModul;
430  return (number) ii;
431}
432
433number nr2mMapQ(number from)
434{
435  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
436  mpz_init(erg);
437
438  nlGMP(from, (number) erg);
439  mpz_mod_ui(erg, erg, nr2mModul);
440  number r = (number) mpz_get_ui(erg);
441
442  mpz_clear(erg);
444  return (number) r;
445}
446
447number nr2mMapGMP(number from)
448{
449  int_number erg = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
450  mpz_init(erg);
451
452  mpz_mod_ui(erg, (int_number) from, nr2mModul);
453  number r = (number) mpz_get_ui(erg);
454
455  mpz_clear(erg);
457  return (number) r;
458}
459
460nMapFunc nr2mSetMap(ring src, ring dst)
461{
462  if (rField_is_Ring_2toM(src)
463     && (src->ringflagb >= dst->ringflagb))
464  {
465    return nr2mMapMachineInt;
466  }
467  if (rField_is_Ring_Z(src))
468  {
469    return nr2mMapGMP;
470  }
471  if (rField_is_Q(src))
472  {
473    return nr2mMapQ;
474  }
475  if (rField_is_Zp(src)
476     && (src->ch == 2)
477     && (dst->ringflagb == 1))
478  {
479    return nr2mMapZp;
480  }
481  if (rField_is_Ring_PtoM(src) || rField_is_Ring_ModN(src))
482  {
483    // Computing the n of Z/n
484    int_number modul = (int_number) omAlloc(sizeof(MP_INT)); // evtl. spaeter mit bin
485    mpz_init(modul);
486    mpz_set(modul, src->ringflaga);
487    mpz_pow_ui(modul, modul, src->ringflagb);
488    if (mpz_divisible_2exp_p(modul, dst->ringflagb))
489    {
490      mpz_clear(modul);
492      return nr2mMapGMP;
493    }
494    mpz_clear(modul);
496  }
497  return NULL;      // default
498}
499
500/*
501 * set the exponent (allocate and init tables) (TODO)
502 */
503
504void nr2mSetExp(int m, ring r)
505{
506  if (m>1)
507  {
508    nr2mExp = m;
509    nr2mModul = 2;
510    for (int i = 1; i < m; i++) {
511      nr2mModul = nr2mModul * 2;
512    }
513  }
514  else
515  {
516    nr2mExp=2;
517    nr2mModul=4;
518  }
519}
520
521void nr2mInitExp(int m, ring r)
522{
523  nr2mSetExp(m, r);
524  if (m<2) WarnS("nInitExp failed: using Z/2^2");
525}
526
527#ifdef LDEBUG
528BOOLEAN nr2mDBTest (number a, const char *f, const int l)
529{
530  if (((NATNUMBER)a<0) || ((NATNUMBER)a>nr2mModul))
531  {
532    return FALSE;
533  }
534  return TRUE;
535}
536#endif
537
538void nr2mWrite (number &a)
539{
540  if ((NATNUMBER)a > (nr2mModul >>1)) StringAppend("-%d",(int)(nr2mModul-((NATNUMBER)a)));
541  else                          StringAppend("%d",(int)((NATNUMBER)a));
542}
543
544static const char* nr2mEati(const char *s, int *i)
545{
546
547  if (((*s) >= '0') && ((*s) <= '9'))
548  {
549    (*i) = 0;
550    do
551    {
552      (*i) *= 10;
553      (*i) += *s++ - '0';
554      if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) % nr2mModul;
555    }
556    while (((*s) >= '0') && ((*s) <= '9'));
557    if ((*i) >= nr2mModul) (*i) = (*i) % nr2mModul;
558  }
559  else (*i) = 1;
560  return s;
561}
562
563const char * nr2mRead (const char *s, number *a)
564{
565  int z;
566  int n=1;
567
568  s = nr2mEati(s, &z);
569  if ((*s) == '/')
570  {
571    s++;
572    s = nr2mEati(s, &n);
573  }
574  if (n == 1)
575    *a = (number)z;
576  else
577      *a = nr2mDiv((number)z,(number)n);
578  return s;
579}
580#endif
Note: See TracBrowser for help on using the repository browser.