rmodulo2m.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT: numbers modulo 2^m
6 */
7 #include <misc/auxiliary.h>
8 
9 #include <omalloc/omalloc.h>
10 
11 #include <misc/mylimits.h>
12 #include <reporter/reporter.h>
13 
14 #include "si_gmp.h"
15 #include "coeffs.h"
16 #include "numbers.h"
17 #include "longrat.h"
18 #include "mpr_complex.h"
19 
20 #include "rmodulo2m.h"
21 #include "rmodulon.h"
22 
23 #include <string.h>
24 
25 #ifdef HAVE_RINGS
26 
27 /// Our Type!
28 static const n_coeffType ID = n_Z2m;
29 
30 number nr2mCopy (number a, const coeffs r);
31 BOOLEAN nr2mGreaterZero (number k, const coeffs r);
32 number nr2mMult (number a, number b, const coeffs r);
33 number nr2mInit (long i, const coeffs r);
34 long nr2mInt (number &n, const coeffs r);
35 number nr2mAdd (number a, number b, const coeffs r);
36 number nr2mSub (number a, number b, const coeffs r);
37 void nr2mPower (number a, int i, number * result, const coeffs r);
38 BOOLEAN nr2mIsZero (number a, const coeffs r);
39 BOOLEAN nr2mIsOne (number a, const coeffs r);
40 BOOLEAN nr2mIsMOne (number a, const coeffs r);
41 BOOLEAN nr2mIsUnit (number a, const coeffs r);
42 number nr2mGetUnit (number a, const coeffs r);
43 number nr2mDiv (number a, number b, const coeffs r);
44 number nr2mIntDiv (number a,number b, const coeffs r);
45 number nr2mMod (number a,number b, const coeffs r);
46 number nr2mNeg (number c, const coeffs r);
47 number nr2mInvers (number c, const coeffs r);
48 BOOLEAN nr2mGreater (number a, number b, const coeffs r);
49 BOOLEAN nr2mDivBy (number a, number b, const coeffs r);
50 int nr2mDivComp (number a, number b, const coeffs r);
51 BOOLEAN nr2mEqual (number a, number b, const coeffs r);
52 number nr2mLcm (number a, number b, const coeffs r);
53 number nr2mGcd (number a, number b, const coeffs r);
54 number nr2mExtGcd (number a, number b, number *s, number *t, const coeffs r);
55 nMapFunc nr2mSetMap (const coeffs src, const coeffs dst);
56 void nr2mWrite (number a, const coeffs r);
57 const char * nr2mRead (const char *s, number *a, const coeffs r);
58 char * nr2mName (number n, const coeffs r);
59 void nr2mCoeffWrite (const coeffs r, BOOLEAN details);
60 coeffs nr2mQuot1(number c, const coeffs r);
61 #ifdef LDEBUG
62 BOOLEAN nr2mDBTest (number a, const char *f, const int l, const coeffs r);
63 #endif
64 void nr2mSetExp(int c, const coeffs r);
65 void nr2mInitExp(int c, const coeffs r);
66 
67 number nr2mMapQ(number from, const coeffs src, const coeffs dst);
68 
69 static inline number nr2mMultM(number a, number b, const coeffs r)
70 {
71  return (number)
72  ((((unsigned long) a) * ((unsigned long) b)) & ((unsigned long)r->mod2mMask));
73 }
74 
75 static inline number nr2mAddM(number a, number b, const coeffs r)
76 {
77  return (number)
78  ((((unsigned long) a) + ((unsigned long) b)) & ((unsigned long)r->mod2mMask));
79 }
80 
81 static inline number nr2mSubM(number a, number b, const coeffs r)
82 {
83  return (number)((unsigned long)a < (unsigned long)b ?
84  r->mod2mMask - (unsigned long)b + (unsigned long)a + 1:
85  (unsigned long)a - (unsigned long)b);
86 }
87 
88 #define nr2mNegM(A,r) (number)((r->mod2mMask - (unsigned long)(A) + 1) & r->mod2mMask)
89 #define nr2mEqualM(A,B) ((A)==(B))
90 
91 extern omBin gmp_nrz_bin; /* init in rintegers*/
92 
93 void nr2mCoeffWrite (const coeffs r, BOOLEAN /*details*/)
94 {
95  PrintS("// coeff. ring is : ");
96  Print("Z/2^%lu\n", r->modExponent);
97 }
98 
100 {
101  if (n==n_Z2m)
102  {
103  int m=(int)(long)(p);
104  unsigned long mm=r->mod2mMask;
105  if (((mm+1)>>m)==1L) return TRUE;
106  }
107  return FALSE;
108 }
109 
110 static char* nr2mCoeffString(const coeffs r)
111 {
112  char* s = (char*) omAlloc(11+11);
113  sprintf(s,"integer,2,%lu",r->modExponent);
114  return s;
115 }
116 
117 coeffs nr2mQuot1(number c, const coeffs r)
118 {
119  coeffs rr;
120  long ch = r->cfInt(c, r);
121  mpz_t a,b;
122  mpz_init_set(a, r->modNumber);
123  mpz_init_set_ui(b, ch);
124  mpz_ptr gcd;
125  gcd = (mpz_ptr) omAlloc(sizeof(mpz_t));
126  mpz_init(gcd);
127  mpz_gcd(gcd, a,b);
128  if(mpz_cmp_ui(gcd, 1) == 0)
129  {
130  WerrorS("constant in q-ideal is coprime to modulus in ground ring");
131  WerrorS("Unable to create qring!");
132  return NULL;
133  }
134  if(mpz_cmp_ui(gcd, 2) == 0)
135  {
136  rr = nInitChar(n_Zp, (void*)2);
137  }
138  else
139  {
140  ZnmInfo info;
141  info.base = r->modBase;
142  int kNew = 1;
143  mpz_t baseTokNew;
144  mpz_init(baseTokNew);
145  mpz_set(baseTokNew, r->modBase);
146  while(mpz_cmp(gcd, baseTokNew) > 0)
147  {
148  kNew++;
149  mpz_mul(baseTokNew, baseTokNew, r->modBase);
150  }
151  info.exp = kNew;
152  mpz_clear(baseTokNew);
153  rr = nInitChar(n_Z2m, (void*)(long)kNew);
154  }
155  return(rr);
156 }
157 
158 static number nr2mAnn(number b, const coeffs r);
159 /* for initializing function pointers */
161 {
162  assume( getCoeffType(r) == ID );
163  nr2mInitExp((int)(long)(p), r);
164 
165  r->is_field=FALSE;
166  r->is_domain=FALSE;
167  r->rep=n_rep_int;
168 
169  //r->cfKillChar = ndKillChar; /* dummy*/
170  r->nCoeffIsEqual = nr2mCoeffIsEqual;
171  r->cfCoeffString = nr2mCoeffString;
172 
173  r->modBase = (mpz_ptr) omAllocBin (gmp_nrz_bin);
174  mpz_init_set_si (r->modBase, 2L);
175  r->modNumber= (mpz_ptr) omAllocBin (gmp_nrz_bin);
176  mpz_init (r->modNumber);
177  mpz_pow_ui (r->modNumber, r->modBase, r->modExponent);
178 
179  /* next cast may yield an overflow as mod2mMask is an unsigned long */
180  r->ch = (int)r->mod2mMask + 1;
181 
182  r->cfInit = nr2mInit;
183  //r->cfCopy = ndCopy;
184  r->cfInt = nr2mInt;
185  r->cfAdd = nr2mAdd;
186  r->cfSub = nr2mSub;
187  r->cfMult = nr2mMult;
188  r->cfDiv = nr2mDiv;
189  r->cfAnn = nr2mAnn;
190  r->cfIntMod = nr2mMod;
191  r->cfExactDiv = nr2mDiv;
192  r->cfInpNeg = nr2mNeg;
193  r->cfInvers = nr2mInvers;
194  r->cfDivBy = nr2mDivBy;
195  r->cfDivComp = nr2mDivComp;
196  r->cfGreater = nr2mGreater;
197  r->cfEqual = nr2mEqual;
198  r->cfIsZero = nr2mIsZero;
199  r->cfIsOne = nr2mIsOne;
200  r->cfIsMOne = nr2mIsMOne;
201  r->cfGreaterZero = nr2mGreaterZero;
202  r->cfWriteLong = nr2mWrite;
203  r->cfRead = nr2mRead;
204  r->cfPower = nr2mPower;
205  r->cfSetMap = nr2mSetMap;
206 // r->cfNormalize = ndNormalize; // default
207  r->cfLcm = nr2mLcm;
208  r->cfGcd = nr2mGcd;
209  r->cfIsUnit = nr2mIsUnit;
210  r->cfGetUnit = nr2mGetUnit;
211  r->cfExtGcd = nr2mExtGcd;
212  r->cfCoeffWrite = nr2mCoeffWrite;
213  r->cfQuot1 = nr2mQuot1;
214 #ifdef LDEBUG
215  r->cfDBTest = nr2mDBTest;
216 #endif
217  r->has_simple_Alloc=TRUE;
218  return FALSE;
219 }
220 
221 /*
222  * Multiply two numbers
223  */
224 number nr2mMult(number a, number b, const coeffs r)
225 {
226  if (((unsigned long)a == 0) || ((unsigned long)b == 0))
227  return (number)0;
228  else
229  return nr2mMultM(a, b, r);
230 }
231 
232 /*
233  * Give the smallest k, such that a * x = k = b * y has a solution
234  */
235 number nr2mLcm(number a, number b, const coeffs)
236 {
237  unsigned long res = 0;
238  if ((unsigned long)a == 0) a = (number) 1;
239  if ((unsigned long)b == 0) b = (number) 1;
240  while ((unsigned long)a % 2 == 0)
241  {
242  a = (number)((unsigned long)a / 2);
243  if ((unsigned long)b % 2 == 0) b = (number)((unsigned long)b / 2);
244  res++;
245  }
246  while ((unsigned long)b % 2 == 0)
247  {
248  b = (number)((unsigned long)b / 2);
249  res++;
250  }
251  return (number)(1L << res); // (2**res)
252 }
253 
254 /*
255  * Give the largest k, such that a = x * k, b = y * k has
256  * a solution.
257  */
258 number nr2mGcd(number a, number b, const coeffs)
259 {
260  unsigned long res = 0;
261  if ((unsigned long)a == 0 && (unsigned long)b == 0) return (number)1;
262  while ((unsigned long)a % 2 == 0 && (unsigned long)b % 2 == 0)
263  {
264  a = (number)((unsigned long)a / 2);
265  b = (number)((unsigned long)b / 2);
266  res++;
267  }
268 // if ((unsigned long)b % 2 == 0)
269 // {
270 // return (number)((1L << res)); // * (unsigned long) a); // (2**res)*a a is a unit
271 // }
272 // else
273 // {
274  return (number)((1L << res)); // * (unsigned long) b); // (2**res)*b b is a unit
275 // }
276 }
277 
278 /*
279  * Give the largest k, such that a = x * k, b = y * k has
280  * a solution.
281  */
282 number nr2mExtGcd(number a, number b, number *s, number *t, const coeffs r)
283 {
284  unsigned long res = 0;
285  if ((unsigned long)a == 0 && (unsigned long)b == 0) return (number)1;
286  while ((unsigned long)a % 2 == 0 && (unsigned long)b % 2 == 0)
287  {
288  a = (number)((unsigned long)a / 2);
289  b = (number)((unsigned long)b / 2);
290  res++;
291  }
292  if ((unsigned long)b % 2 == 0)
293  {
294  *t = NULL;
295  *s = nr2mInvers(a,r);
296  return (number)((1L << res)); // * (unsigned long) a); // (2**res)*a a is a unit
297  }
298  else
299  {
300  *s = NULL;
301  *t = nr2mInvers(b,r);
302  return (number)((1L << res)); // * (unsigned long) b); // (2**res)*b b is a unit
303  }
304 }
305 
306 void nr2mPower(number a, int i, number * result, const coeffs r)
307 {
308  if (i == 0)
309  {
310  *(unsigned long *)result = 1;
311  }
312  else if (i == 1)
313  {
314  *result = a;
315  }
316  else
317  {
318  nr2mPower(a, i-1, result, r);
319  *result = nr2mMultM(a, *result, r);
320  }
321 }
322 
323 /*
324  * create a number from int
325  */
326 number nr2mInit(long i, const coeffs r)
327 {
328  if (i == 0) return (number)(unsigned long)i;
329 
330  long ii = i;
331  unsigned long j = (unsigned long)1;
332  if (ii < 0) { j = r->mod2mMask; ii = -ii; }
333  unsigned long k = (unsigned long)ii;
334  k = k & r->mod2mMask;
335  /* now we have: i = j * k mod 2^m */
336  return (number)nr2mMult((number)j, (number)k, r);
337 }
338 
339 /*
340  * convert a number to an int in ]-k/2 .. k/2],
341  * where k = 2^m; i.e., an int in ]-2^(m-1) .. 2^(m-1)];
342  */
343 long nr2mInt(number &n, const coeffs r)
344 {
345  unsigned long nn = (unsigned long)(unsigned long)n & r->mod2mMask;
346  unsigned long l = r->mod2mMask >> 1; l++; /* now: l = 2^(m-1) */
347  if ((unsigned long)nn > l)
348  return (long)((unsigned long)nn - r->mod2mMask - 1);
349  else
350  return (long)((unsigned long)nn);
351 }
352 
353 number nr2mAdd(number a, number b, const coeffs r)
354 {
355  return nr2mAddM(a, b, r);
356 }
357 
358 number nr2mSub(number a, number b, const coeffs r)
359 {
360  return nr2mSubM(a, b, r);
361 }
362 
363 BOOLEAN nr2mIsUnit(number a, const coeffs)
364 {
365  return ((unsigned long)a % 2 == 1);
366 }
367 
368 number nr2mGetUnit(number k, const coeffs)
369 {
370  if (k == NULL) return (number)1;
371  unsigned long erg = (unsigned long)k;
372  while (erg % 2 == 0) erg = erg / 2;
373  return (number)erg;
374 }
375 
376 BOOLEAN nr2mIsZero(number a, const coeffs)
377 {
378  return 0 == (unsigned long)a;
379 }
380 
381 BOOLEAN nr2mIsOne(number a, const coeffs)
382 {
383  return 1 == (unsigned long)a;
384 }
385 
386 BOOLEAN nr2mIsMOne(number a, const coeffs r)
387 {
388  return (r->mod2mMask == (unsigned long)a);
389 }
390 
391 BOOLEAN nr2mEqual(number a, number b, const coeffs)
392 {
393  return (a == b);
394 }
395 
396 BOOLEAN nr2mGreater(number a, number b, const coeffs r)
397 {
398  return nr2mDivBy(a, b,r);
399 }
400 
401 /* Is 'a' divisible by 'b'? There are two cases:
402  1) a = 0 mod 2^m; then TRUE iff b = 0 or b is a power of 2
403  2) a, b <> 0; then TRUE iff b/gcd(a, b) is a unit mod 2^m */
404 BOOLEAN nr2mDivBy (number a, number b, const coeffs r)
405 {
406  if (a == NULL)
407  {
408  unsigned long c = r->mod2mMask + 1;
409  if (c != 0) /* i.e., if no overflow */
410  return (c % (unsigned long)b) == 0;
411  else
412  {
413  /* overflow: we need to check whether b
414  is zero or a power of 2: */
415  c = (unsigned long)b;
416  while (c != 0)
417  {
418  if ((c % 2) != 0) return FALSE;
419  c = c >> 1;
420  }
421  return TRUE;
422  }
423  }
424  else
425  {
426  number n = nr2mGcd(a, b, r);
427  n = nr2mDiv(b, n, r);
428  return nr2mIsUnit(n, r);
429  }
430 }
431 
432 int nr2mDivComp(number as, number bs, const coeffs)
433 {
434  unsigned long a = (unsigned long)as;
435  unsigned long b = (unsigned long)bs;
436  assume(a != 0 && b != 0);
437  while (a % 2 == 0 && b % 2 == 0)
438  {
439  a = a / 2;
440  b = b / 2;
441  }
442  if (a % 2 == 0)
443  {
444  return -1;
445  }
446  else
447  {
448  if (b % 2 == 1)
449  {
450  return 2;
451  }
452  else
453  {
454  return 1;
455  }
456  }
457 }
458 
459 /* TRUE iff 0 < k <= 2^m / 2 */
461 {
462  if ((unsigned long)k == 0) return FALSE;
463  if ((unsigned long)k > ((r->mod2mMask >> 1) + 1)) return FALSE;
464  return TRUE;
465 }
466 
467 /* assumes that 'a' is odd, i.e., a unit in Z/2^m, and computes
468  the extended gcd of 'a' and 2^m, in order to find some 's'
469  and 't' such that a * s + 2^m * t = gcd(a, 2^m) = 1;
470  this code will always find a positive 's' */
471 void specialXGCD(unsigned long& s, unsigned long a, const coeffs r)
472 {
473  mpz_ptr u = (mpz_ptr)omAlloc(sizeof(mpz_t));
474  mpz_init_set_ui(u, a);
475  mpz_ptr u0 = (mpz_ptr)omAlloc(sizeof(mpz_t));
476  mpz_init(u0);
477  mpz_ptr u1 = (mpz_ptr)omAlloc(sizeof(mpz_t));
478  mpz_init_set_ui(u1, 1);
479  mpz_ptr u2 = (mpz_ptr)omAlloc(sizeof(mpz_t));
480  mpz_init(u2);
481  mpz_ptr v = (mpz_ptr)omAlloc(sizeof(mpz_t));
482  mpz_init_set_ui(v, r->mod2mMask);
483  mpz_add_ui(v, v, 1); /* now: v = 2^m */
484  mpz_ptr v0 = (mpz_ptr)omAlloc(sizeof(mpz_t));
485  mpz_init(v0);
486  mpz_ptr v1 = (mpz_ptr)omAlloc(sizeof(mpz_t));
487  mpz_init(v1);
488  mpz_ptr v2 = (mpz_ptr)omAlloc(sizeof(mpz_t));
489  mpz_init_set_ui(v2, 1);
490  mpz_ptr q = (mpz_ptr)omAlloc(sizeof(mpz_t));
491  mpz_init(q);
492  mpz_ptr rr = (mpz_ptr)omAlloc(sizeof(mpz_t));
493  mpz_init(rr);
494 
495  while (mpz_cmp_ui(v, 0) != 0) /* i.e., while v != 0 */
496  {
497  mpz_div(q, u, v);
498  mpz_mod(rr, u, v);
499  mpz_set(u, v);
500  mpz_set(v, rr);
501  mpz_set(u0, u2);
502  mpz_set(v0, v2);
503  mpz_mul(u2, u2, q); mpz_sub(u2, u1, u2); /* u2 = u1 - q * u2 */
504  mpz_mul(v2, v2, q); mpz_sub(v2, v1, v2); /* v2 = v1 - q * v2 */
505  mpz_set(u1, u0);
506  mpz_set(v1, v0);
507  }
508 
509  while (mpz_cmp_ui(u1, 0) < 0) /* i.e., while u1 < 0 */
510  {
511  /* we add 2^m = (2^m - 1) + 1 to u1: */
512  mpz_add_ui(u1, u1, r->mod2mMask);
513  mpz_add_ui(u1, u1, 1);
514  }
515  s = mpz_get_ui(u1); /* now: 0 <= s <= 2^m - 1 */
516 
517  mpz_clear(u); omFree((ADDRESS)u);
518  mpz_clear(u0); omFree((ADDRESS)u0);
519  mpz_clear(u1); omFree((ADDRESS)u1);
520  mpz_clear(u2); omFree((ADDRESS)u2);
521  mpz_clear(v); omFree((ADDRESS)v);
522  mpz_clear(v0); omFree((ADDRESS)v0);
523  mpz_clear(v1); omFree((ADDRESS)v1);
524  mpz_clear(v2); omFree((ADDRESS)v2);
525  mpz_clear(q); omFree((ADDRESS)q);
526  mpz_clear(rr); omFree((ADDRESS)rr);
527 }
528 
529 unsigned long InvMod(unsigned long a, const coeffs r)
530 {
531  assume((unsigned long)a % 2 != 0);
532  unsigned long s;
533  specialXGCD(s, a, r);
534  return s;
535 }
536 //#endif
537 
538 inline number nr2mInversM(number c, const coeffs r)
539 {
540  assume((unsigned long)c % 2 != 0);
541  // Table !!!
542  unsigned long inv;
543  inv = InvMod((unsigned long)c,r);
544  return (number)inv;
545 }
546 
547 number nr2mDiv(number a, number b, const coeffs r)
548 {
549  if ((unsigned long)a == 0) return (number)0;
550  else if ((unsigned long)b % 2 == 0)
551  {
552  if ((unsigned long)b != 0)
553  {
554  while (((unsigned long)b % 2 == 0) && ((unsigned long)a % 2 == 0))
555  {
556  a = (number)((unsigned long)a / 2);
557  b = (number)((unsigned long)b / 2);
558  }
559  }
560  if ((unsigned long)b % 2 == 0)
561  {
562  WerrorS("Division not possible, even by cancelling zero divisors.");
563  WerrorS("Result is integer division without remainder.");
564  return (number) ((unsigned long) a / (unsigned long) b);
565  }
566  }
567  return (number)nr2mMult(a, nr2mInversM(b,r),r);
568 }
569 
570 number nr2mMod(number a, number b, const coeffs r)
571 {
572  /*
573  We need to return the number rr which is uniquely determined by the
574  following two properties:
575  (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
576  (2) There exists some k in the integers Z such that a = k * b + rr.
577  Consider g := gcd(2^m, |b|). Note that then |b|/g is a unit in Z/2^m.
578  Now, there are three cases:
579  (a) g = 1
580  Then |b| is a unit in Z/2^m, i.e. |b| (and also b) divides a.
581  Thus rr = 0.
582  (b) g <> 1 and g divides a
583  Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
584  (c) g <> 1 and g does not divide a
585  Let's denote the division with remainder of a by g as follows:
586  a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
587  fulfills (1) and (2), i.e. rr := t is the correct result. Hence
588  in this third case, rr is the remainder of division of a by g in Z.
589  This algorithm is the same as for the case Z/n, except that we may
590  compute the gcd of |b| and 2^m "by hand": We just extract the highest
591  power of 2 (<= 2^m) that is contained in b.
592  */
593  assume((unsigned long) b != 0);
594  unsigned long g = 1;
595  unsigned long b_div = (unsigned long) b;
596 
597  /*
598  * b_div is unsigned, so that (b_div < 0) evaluates false at compile-time
599  *
600  if (b_div < 0) b_div = -b_div; // b_div now represents |b|, BUT b_div is unsigned!
601  */
602 
603  unsigned long rr = 0;
604  while ((g < r->mod2mMask ) && (b_div > 0) && (b_div % 2 == 0))
605  {
606  b_div = b_div >> 1;
607  g = g << 1;
608  } // g is now the gcd of 2^m and |b|
609 
610  if (g != 1) rr = (unsigned long)a % g;
611  return (number)rr;
612 }
613 
614 number nr2mIntDiv(number a, number b, const coeffs r)
615 {
616  if ((unsigned long)a == 0)
617  {
618  if ((unsigned long)b == 0)
619  return (number)1;
620  if ((unsigned long)b == 1)
621  return (number)0;
622  unsigned long c = r->mod2mMask + 1;
623  if (c != 0) /* i.e., if no overflow */
624  return (number)(c / (unsigned long)b);
625  else
626  {
627  /* overflow: c = 2^32 resp. 2^64, depending on platform */
628  mpz_ptr cc = (mpz_ptr)omAlloc(sizeof(mpz_t));
629  mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
630  mpz_div_ui(cc, cc, (unsigned long)(unsigned long)b);
631  unsigned long s = mpz_get_ui(cc);
632  mpz_clear(cc); omFree((ADDRESS)cc);
633  return (number)(unsigned long)s;
634  }
635  }
636  else
637  {
638  if ((unsigned long)b == 0)
639  return (number)0;
640  return (number)((unsigned long) a / (unsigned long) b);
641  }
642 }
643 
644 static number nr2mAnn(number b, const coeffs r)
645 {
646  if ((unsigned long)b == 0)
647  return NULL;
648  if ((unsigned long)b == 1)
649  return NULL;
650  unsigned long c = r->mod2mMask + 1;
651  if (c != 0) /* i.e., if no overflow */
652  return (number)(c / (unsigned long)b);
653  else
654  {
655  /* overflow: c = 2^32 resp. 2^64, depending on platform */
656  mpz_ptr cc = (mpz_ptr)omAlloc(sizeof(mpz_t));
657  mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
658  mpz_div_ui(cc, cc, (unsigned long)(unsigned long)b);
659  unsigned long s = mpz_get_ui(cc);
660  mpz_clear(cc); omFree((ADDRESS)cc);
661  return (number)(unsigned long)s;
662  }
663 }
664 
665 number nr2mInvers(number c, const coeffs r)
666 {
667  if ((unsigned long)c % 2 == 0)
668  {
669  WerrorS("division by zero divisor");
670  return (number)0;
671  }
672  return nr2mInversM(c, r);
673 }
674 
675 number nr2mNeg(number c, const coeffs r)
676 {
677  if ((unsigned long)c == 0) return c;
678  return nr2mNegM(c, r);
679 }
680 
681 number nr2mMapMachineInt(number from, const coeffs /*src*/, const coeffs dst)
682 {
683  unsigned long i = ((unsigned long)from) % dst->mod2mMask ;
684  return (number)i;
685 }
686 
687 number nr2mMapProject(number from, const coeffs /*src*/, const coeffs dst)
688 {
689  unsigned long i = ((unsigned long)from) % (dst->mod2mMask + 1);
690  return (number)i;
691 }
692 
693 number nr2mMapZp(number from, const coeffs /*src*/, const coeffs dst)
694 {
695  unsigned long j = (unsigned long)1;
696  long ii = (long)from;
697  if (ii < 0) { j = dst->mod2mMask; ii = -ii; }
698  unsigned long i = (unsigned long)ii;
699  i = i & dst->mod2mMask;
700  /* now we have: from = j * i mod 2^m */
701  return (number)nr2mMult((number)i, (number)j, dst);
702 }
703 
704 number nr2mMapQ(number from, const coeffs src, const coeffs dst)
705 {
706  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
707  mpz_init(erg);
708  mpz_ptr k = (mpz_ptr)omAlloc(sizeof(mpz_t));
709  mpz_init_set_ui(k, dst->mod2mMask);
710 
711  nlGMP(from, (number)erg, src); // FIXME? TODO? // extern void nlGMP(number &i, number n, const coeffs r); // to be replaced with n_MPZ(erg, from, src); // ?
712  mpz_and(erg, erg, k);
713  number res = (number)mpz_get_ui(erg);
714 
715  mpz_clear(erg); omFree((ADDRESS)erg);
716  mpz_clear(k); omFree((ADDRESS)k);
717 
718  return (number)res;
719 }
720 
721 number nr2mMapGMP(number from, const coeffs /*src*/, const coeffs dst)
722 {
723  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
724  mpz_init(erg);
725  mpz_ptr k = (mpz_ptr)omAlloc(sizeof(mpz_t));
726  mpz_init_set_ui(k, dst->mod2mMask);
727 
728  mpz_and(erg, (mpz_ptr)from, k);
729  number res = (number) mpz_get_ui(erg);
730 
731  mpz_clear(erg); omFree((ADDRESS)erg);
732  mpz_clear(k); omFree((ADDRESS)k);
733 
734  return (number)res;
735 }
736 
737 number nr2mMapZ(number from, const coeffs src, const coeffs dst)
738 {
739  if (SR_HDL(from) & SR_INT)
740  {
741  long f_i=SR_TO_INT(from);
742  return nr2mInit(f_i,dst);
743  }
744  return nr2mMapGMP(from,src,dst);
745 }
746 
747 nMapFunc nr2mSetMap(const coeffs src, const coeffs dst)
748 {
749  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
750  && (src->mod2mMask == dst->mod2mMask))
751  {
752  return ndCopyMap;
753  }
754  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
755  && (src->mod2mMask < dst->mod2mMask))
756  { /* i.e. map an integer mod 2^s into Z mod 2^t, where t < s */
757  return nr2mMapMachineInt;
758  }
759  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
760  && (src->mod2mMask > dst->mod2mMask))
761  { /* i.e. map an integer mod 2^s into Z mod 2^t, where t > s */
762  // to be done
763  return nr2mMapProject;
764  }
765  if ((src->rep==n_rep_gmp) && nCoeff_is_Ring_Z(src))
766  {
767  return nr2mMapGMP;
768  }
769  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
770  {
771  return nr2mMapZ;
772  }
773  if ((src->rep==n_rep_gap_rat) && nCoeff_is_Q(src))
774  {
775  return nr2mMapQ;
776  }
777  if ((src->rep==n_rep_int) && nCoeff_is_Zp(src) && (src->ch == 2))
778  {
779  return nr2mMapZp;
780  }
781  if ((src->rep==n_rep_gmp) &&
783  {
784  if (mpz_divisible_2exp_p(src->modNumber,dst->modExponent))
785  return nr2mMapGMP;
786  }
787  return NULL; // default
788 }
789 
790 /*
791  * set the exponent
792  */
793 
794 void nr2mSetExp(int m, coeffs r)
795 {
796  if (m > 1)
797  {
798  /* we want mod2mMask to be the bit pattern
799  '111..1' consisting of m one's: */
800  r->modExponent= m;
801  r->mod2mMask = 1;
802  for (int i = 1; i < m; i++) r->mod2mMask = (r->mod2mMask << 1) + 1;
803  }
804  else
805  {
806  r->modExponent= 2;
807  /* code unexpectedly called with m = 1; we continue with m = 2: */
808  r->mod2mMask = 3; /* i.e., '11' in binary representation */
809  }
810 }
811 
812 void nr2mInitExp(int m, coeffs r)
813 {
814  nr2mSetExp(m, r);
815  if (m < 2)
816  WarnS("nr2mInitExp unexpectedly called with m = 1 (we continue with Z/2^2");
817 }
818 
819 #ifdef LDEBUG
820 BOOLEAN nr2mDBTest (number a, const char *, const int, const coeffs r)
821 {
822  //if ((unsigned long)a < 0) return FALSE; // is unsigned!
823  if (((unsigned long)a & r->mod2mMask) != (unsigned long)a) return FALSE;
824  return TRUE;
825 }
826 #endif
827 
828 void nr2mWrite (number a, const coeffs r)
829 {
830  long i = nr2mInt(a, r);
831  StringAppend("%ld", i);
832 }
833 
834 static const char* nr2mEati(const char *s, int *i, const coeffs r)
835 {
836 
837  if (((*s) >= '0') && ((*s) <= '9'))
838  {
839  (*i) = 0;
840  do
841  {
842  (*i) *= 10;
843  (*i) += *s++ - '0';
844  if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & r->mod2mMask;
845  }
846  while (((*s) >= '0') && ((*s) <= '9'));
847  (*i) = (*i) & r->mod2mMask;
848  }
849  else (*i) = 1;
850  return s;
851 }
852 
853 const char * nr2mRead (const char *s, number *a, const coeffs r)
854 {
855  int z;
856  int n=1;
857 
858  s = nr2mEati(s, &z,r);
859  if ((*s) == '/')
860  {
861  s++;
862  s = nr2mEati(s, &n,r);
863  }
864  if (n == 1)
865  *a = (number)(long)z;
866  else
867  *a = nr2mDiv((number)(long)z,(number)(long)n,r);
868  return s;
869 }
870 #endif
871 /* #ifdef HAVE_RINGS */
mpz_ptr base
Definition: rmodulon.h:18
void nr2mInitExp(int c, const coeffs r)
Definition: rmodulo2m.cc:812
mpz_t z
Definition: longrat.h:48
number nr2mGcd(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:258
long nr2mInt(number &n, const coeffs r)
Definition: rmodulo2m.cc:343
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
number nr2mCopy(number a, const coeffs r)
number nr2mLcm(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:235
static number nr2mMultM(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:69
const CanonicalForm int s
Definition: facAbsFact.cc:55
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_ModN(const coeffs r)
Definition: coeffs.h:744
nMapFunc nr2mSetMap(const coeffs src, const coeffs dst)
Definition: rmodulo2m.cc:747
const poly a
Definition: syzextra.cc:212
number nr2mMapZ(number from, const coeffs src, const coeffs dst)
Definition: rmodulo2m.cc:737
omBin_t * omBin
Definition: omStructs.h:12
#define Print
Definition: emacs.cc:83
number nr2mInit(long i, const coeffs r)
Definition: rmodulo2m.cc:326
number nr2mAdd(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:353
static FORCE_INLINE BOOLEAN nCoeff_is_Zp(const coeffs r)
Definition: coeffs.h:818
number nr2mSub(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:358
static number nr2mAnn(number b, const coeffs r)
Definition: rmodulo2m.cc:644
number nr2mMult(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:224
const char * nr2mRead(const char *s, number *a, const coeffs r)
Definition: rmodulo2m.cc:853
only used if HAVE_RINGS is defined: ?
Definition: coeffs.h:45
#define FALSE
Definition: auxiliary.h:140
return P p
Definition: myNF.cc:203
number ndCopyMap(number a, const coeffs aRing, const coeffs r)
Definition: numbers.cc:239
BOOLEAN nr2mGreaterZero(number k, const coeffs r)
Definition: rmodulo2m.cc:460
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_Z(const coeffs r)
Definition: coeffs.h:750
{p < 2^31}
Definition: coeffs.h:30
#define nr2mNegM(A, r)
Definition: rmodulo2m.cc:88
const CanonicalForm CFMap CFMap int &both_non_zero int n
Definition: cfEzgcd.cc:52
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_2toM(const coeffs r)
Definition: coeffs.h:741
(), see rinteger.h, new impl.
Definition: coeffs.h:111
number nr2mGetUnit(number a, const coeffs r)
Definition: rmodulo2m.cc:368
#define TRUE
Definition: auxiliary.h:144
number nr2mMod(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:570
number nr2mMapZp(number from, const coeffs, const coeffs dst)
Definition: rmodulo2m.cc:693
void * ADDRESS
Definition: auxiliary.h:161
g
Definition: cfModGcd.cc:4031
void WerrorS(const char *s)
Definition: feFopen.cc:23
int k
Definition: cfEzgcd.cc:93
BOOLEAN nr2mIsMOne(number a, const coeffs r)
Definition: rmodulo2m.cc:386
BOOLEAN nr2mInitChar(coeffs r, void *p)
Definition: rmodulo2m.cc:160
number nr2mDiv(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:547
void nlGMP(number &i, number n, const coeffs r)
Definition: longrat.cc:1410
BOOLEAN nr2mGreater(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:396
static const n_coeffType ID
Our Type!
Definition: rmodulo2m.cc:28
static FORCE_INLINE BOOLEAN nCoeff_is_Q(const coeffs r)
Definition: coeffs.h:824
#define WarnS
Definition: emacs.cc:81
number nr2mMapMachineInt(number from, const coeffs, const coeffs dst)
Definition: rmodulo2m.cc:681
void nr2mWrite(number a, const coeffs r)
Definition: rmodulo2m.cc:828
coeffs nr2mQuot1(number c, const coeffs r)
Definition: rmodulo2m.cc:117
#define omAlloc(size)
Definition: omAllocDecl.h:210
number nr2mInvers(number c, const coeffs r)
Definition: rmodulo2m.cc:665
poly res
Definition: myNF.cc:322
number nr2mNeg(number c, const coeffs r)
Definition: rmodulo2m.cc:675
number nr2mIntDiv(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:614
const ring r
Definition: syzextra.cc:208
#define LDEBUG
Definition: mod2.h:323
Coefficient rings, fields and other domains suitable for Singular polynomials.
int nr2mDivComp(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:432
int j
Definition: myNF.cc:70
#define omFree(addr)
Definition: omAllocDecl.h:261
#define assume(x)
Definition: mod2.h:405
The main handler for Singular numbers which are suitable for Singular polynomials.
BOOLEAN nr2mDBTest(number a, const char *f, const int l, const coeffs r)
Definition: rmodulo2m.cc:820
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:72
static char * nr2mCoeffString(const coeffs r)
Definition: rmodulo2m.cc:110
const int MAX_INT_VAL
Definition: mylimits.h:12
All the auxiliary stuff.
int m
Definition: cfEzgcd.cc:119
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_PtoM(const coeffs r)
Definition: coeffs.h:747
unsigned long exp
Definition: rmodulon.h:18
#define info
Definition: libparse.cc:1254
#define StringAppend
Definition: emacs.cc:82
FILE * f
Definition: checklibs.c:7
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
number nr2mMapProject(number from, const coeffs, const coeffs dst)
Definition: rmodulo2m.cc:687
(mpz_ptr), see rmodulon,h
Definition: coeffs.h:114
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:422
static const char * nr2mEati(const char *s, int *i, const coeffs r)
Definition: rmodulo2m.cc:834
#define SR_TO_INT(SR)
Definition: longrat.h:67
(number), see longrat.h
Definition: coeffs.h:110
char * nr2mName(number n, const coeffs r)
void nr2mSetExp(int c, const coeffs r)
Definition: rmodulo2m.cc:794
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
number nr2mInversM(number c, const coeffs r)
Definition: rmodulo2m.cc:538
n_coeffType
Definition: coeffs.h:27
BOOLEAN nr2mDivBy(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:404
unsigned long InvMod(unsigned long a, const coeffs r)
Definition: rmodulo2m.cc:529
#define NULL
Definition: omList.c:10
BOOLEAN nr2mIsUnit(number a, const coeffs r)
Definition: rmodulo2m.cc:363
BOOLEAN nr2mIsZero(number a, const coeffs r)
Definition: rmodulo2m.cc:376
static number nr2mAddM(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:75
int gcd(int a, int b)
Definition: walkSupport.cc:839
#define SR_INT
Definition: longrat.h:65
number nr2mExtGcd(number a, number b, number *s, number *t, const coeffs r)
Definition: rmodulo2m.cc:282
number nr2mMapGMP(number from, const coeffs, const coeffs dst)
Definition: rmodulo2m.cc:721
number nr2mMapQ(number from, const coeffs src, const coeffs dst)
Definition: rmodulo2m.cc:704
static number nr2mSubM(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:81
void nr2mCoeffWrite(const coeffs r, BOOLEAN details)
Definition: rmodulo2m.cc:93
omBin gmp_nrz_bin
Definition: rintegers.cc:80
(int), see modulop.h
Definition: coeffs.h:109
#define SR_HDL(A)
Definition: tgb.cc:35
void specialXGCD(unsigned long &s, unsigned long a, const coeffs r)
Definition: rmodulo2m.cc:471
BOOLEAN nr2mEqual(number a, number b, const coeffs r)
Definition: rmodulo2m.cc:391
int BOOLEAN
Definition: auxiliary.h:131
const poly b
Definition: syzextra.cc:213
BOOLEAN nr2mIsOne(number a, const coeffs r)
Definition: rmodulo2m.cc:381
void nr2mPower(number a, int i, number *result, const coeffs r)
Definition: rmodulo2m.cc:306
int l
Definition: cfEzgcd.cc:94
return result
Definition: facAbsBiFact.cc:76
BOOLEAN nr2mCoeffIsEqual(const coeffs r, n_coeffType n, void *p)
Definition: rmodulo2m.cc:99
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:327