rmodulon.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT: numbers modulo n
6 */
7 #include <misc/auxiliary.h>
8 #include <omalloc/omalloc.h>
9 
10 #include <misc/mylimits.h>
11 #include <reporter/reporter.h>
12 
13 #include "si_gmp.h"
14 #include "coeffs.h"
15 #include "numbers.h"
16 
17 #include "mpr_complex.h"
18 
19 #include "longrat.h"
20 #include "rmodulon.h"
21 
22 #include <string.h>
23 
24 #ifdef HAVE_RINGS
25 
26 /// Our Type!
27 static const n_coeffType ID = n_Zn;
28 static const n_coeffType ID2 = n_Znm;
29 
30 number nrnCopy (number a, const coeffs r);
31 int nrnSize (number a, const coeffs r);
32 void nrnDelete (number *a, const coeffs r);
33 BOOLEAN nrnGreaterZero (number k, const coeffs r);
34 number nrnMult (number a, number b, const coeffs r);
35 number nrnInit (long i, const coeffs r);
36 long nrnInt (number &n, const coeffs r);
37 number nrnAdd (number a, number b, const coeffs r);
38 number nrnSub (number a, number b, const coeffs r);
39 void nrnPower (number a, int i, number * result, const coeffs r);
40 BOOLEAN nrnIsZero (number a, const coeffs r);
41 BOOLEAN nrnIsOne (number a, const coeffs r);
42 BOOLEAN nrnIsMOne (number a, const coeffs r);
43 BOOLEAN nrnIsUnit (number a, const coeffs r);
44 number nrnGetUnit (number a, const coeffs r);
45 number nrnAnn (number a, const coeffs r);
46 number nrnDiv (number a, number b, const coeffs r);
47 number nrnMod (number a,number b, const coeffs r);
48 number nrnIntDiv (number a,number b, const coeffs r);
49 number nrnNeg (number c, const coeffs r);
50 number nrnInvers (number c, const coeffs r);
51 BOOLEAN nrnGreater (number a, number b, const coeffs r);
52 BOOLEAN nrnDivBy (number a, number b, const coeffs r);
53 int nrnDivComp (number a, number b, const coeffs r);
54 BOOLEAN nrnEqual (number a, number b, const coeffs r);
55 number nrnLcm (number a,number b, const coeffs r);
56 number nrnGcd (number a,number b, const coeffs r);
57 number nrnExtGcd (number a, number b, number *s, number *t, const coeffs r);
58 number nrnXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs r);
59 number nrnQuotRem (number a, number b, number *s, const coeffs r);
60 nMapFunc nrnSetMap (const coeffs src, const coeffs dst);
61 #if SI_INTEGER_VARIANT==2
62 // FIXME? TODO? // extern void nrzWrite (number &a, const coeffs r); // FIXME
63 # define nrnWrite nrzWrite
64 #else
65 void nrnWrite (number &a, const coeffs);
66 #endif
67 const char * nrnRead (const char *s, number *a, const coeffs r);
68 void nrnCoeffWrite (const coeffs r, BOOLEAN details);
69 #ifdef LDEBUG
70 BOOLEAN nrnDBTest (number a, const char *f, const int l, const coeffs r);
71 #endif
72 void nrnSetExp(unsigned long c, const coeffs r);
73 void nrnInitExp(unsigned long c, const coeffs r);
74 coeffs nrnQuot1(number c, const coeffs r);
75 
76 number nrnMapQ(number from, const coeffs src, const coeffs dst);
77 
78 
79 extern omBin gmp_nrz_bin;
80 
81 void nrnCoeffWrite (const coeffs r, BOOLEAN /*details*/)
82 {
83  long l = (long)mpz_sizeinbase(r->modBase, 10) + 2;
84  char* s = (char*) omAlloc(l);
85  s= mpz_get_str (s, 10, r->modBase);
86  if (nCoeff_is_Ring_ModN(r)) Print("// coeff. ring is : Z/%s\n", s);
87  else if (nCoeff_is_Ring_PtoM(r)) Print("// coeff. ring is : Z/%s^%lu\n", s, r->modExponent);
88  omFreeSize((ADDRESS)s, l);
89 }
90 
91 static BOOLEAN nrnCoeffsEqual(const coeffs r, n_coeffType n, void * parameter)
92 {
93  /* test, if r is an instance of nInitCoeffs(n,parameter) */
94  return (n==n_Zn) && (mpz_cmp_ui(r->modNumber,(long)parameter)==0);
95 }
96 
97 static char* nrnCoeffString(const coeffs r)
98 {
99  long l = (long)mpz_sizeinbase(r->modBase, 10) +2;
100  char* b = (char*) omAlloc(l);
101  b= mpz_get_str (b, 10, r->modBase);
102  char* s = (char*) omAlloc(7+2+10+l);
103  if (nCoeff_is_Ring_ModN(r)) sprintf(s,"integer,%s",b);
104  else /*if (nCoeff_is_Ring_PtoM(r))*/ sprintf(s,"integer,%s^%lu",b,r->modExponent);
105  omFreeSize(b,l);
106  return s;
107 }
108 
109 static void nrnKillChar(coeffs r)
110 {
111  mpz_clear(r->modNumber);
112  mpz_clear(r->modBase);
113  omFreeBin((void *) r->modBase, gmp_nrz_bin);
114  omFreeBin((void *) r->modNumber, gmp_nrz_bin);
115 }
116 
117 coeffs nrnQuot1(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(r->modExponent == 1)
135  {
136  ZnmInfo info;
137  info.base = gcd;
138  info.exp = (unsigned long) 1;
139  rr = nInitChar(n_Zn, (void*)&info);
140  }
141  else
142  {
143  ZnmInfo info;
144  info.base = r->modBase;
145  int kNew = 1;
146  mpz_t baseTokNew;
147  mpz_init(baseTokNew);
148  mpz_set(baseTokNew, r->modBase);
149  while(mpz_cmp(gcd, baseTokNew) > 0)
150  {
151  kNew++;
152  mpz_mul(baseTokNew, baseTokNew, r->modBase);
153  }
154  //printf("\nkNew = %i\n",kNew);
155  info.exp = kNew;
156  mpz_clear(baseTokNew);
157  rr = nInitChar(n_Znm, (void*)&info);
158  }
159  return(rr);
160 }
161 
162 /* for initializing function pointers */
164 {
165  assume( (getCoeffType(r) == ID) || (getCoeffType (r) == ID2) );
166  ZnmInfo * info= (ZnmInfo *) p;
167  r->modBase= (mpz_ptr)nrnCopy((number)info->base, r); //this circumvents the problem
168  //in bigintmat.cc where we cannot create a "legal" nrn that can be freed.
169  //If we take a copy, we can do whatever we want.
170 
171  nrnInitExp (info->exp, r);
172 
173  /* next computation may yield wrong characteristic as r->modNumber
174  is a GMP number */
175  r->ch = mpz_get_ui(r->modNumber);
176 
177  r->is_field=FALSE;
178  r->is_domain=FALSE;
179  r->rep=n_rep_gmp;
180 
181 
182  r->cfCoeffString = nrnCoeffString;
183 
184  r->cfInit = nrnInit;
185  r->cfDelete = nrnDelete;
186  r->cfCopy = nrnCopy;
187  r->cfSize = nrnSize;
188  r->cfInt = nrnInt;
189  r->cfAdd = nrnAdd;
190  r->cfSub = nrnSub;
191  r->cfMult = nrnMult;
192  r->cfDiv = nrnDiv;
193  r->cfAnn = nrnAnn;
194  r->cfIntMod = nrnMod;
195  r->cfExactDiv = nrnDiv;
196  r->cfInpNeg = nrnNeg;
197  r->cfInvers = nrnInvers;
198  r->cfDivBy = nrnDivBy;
199  r->cfDivComp = nrnDivComp;
200  r->cfGreater = nrnGreater;
201  r->cfEqual = nrnEqual;
202  r->cfIsZero = nrnIsZero;
203  r->cfIsOne = nrnIsOne;
204  r->cfIsMOne = nrnIsMOne;
205  r->cfGreaterZero = nrnGreaterZero;
206  r->cfWriteLong = nrnWrite;
207  r->cfRead = nrnRead;
208  r->cfPower = nrnPower;
209  r->cfSetMap = nrnSetMap;
210  //r->cfNormalize = ndNormalize;
211  r->cfLcm = nrnLcm;
212  r->cfGcd = nrnGcd;
213  r->cfIsUnit = nrnIsUnit;
214  r->cfGetUnit = nrnGetUnit;
215  r->cfExtGcd = nrnExtGcd;
216  r->cfXExtGcd = nrnXExtGcd;
217  r->cfQuotRem = nrnQuotRem;
218  r->cfCoeffWrite = nrnCoeffWrite;
219  r->nCoeffIsEqual = nrnCoeffsEqual;
220  r->cfKillChar = nrnKillChar;
221  r->cfQuot1 = nrnQuot1;
222 #ifdef LDEBUG
223  r->cfDBTest = nrnDBTest;
224 #endif
225  return FALSE;
226 }
227 
228 /*
229  * create a number from int
230  */
231 number nrnInit(long i, const coeffs r)
232 {
233  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
234  mpz_init_set_si(erg, i);
235  mpz_mod(erg, erg, r->modNumber);
236  return (number) erg;
237 }
238 
239 void nrnDelete(number *a, const coeffs)
240 {
241  if (*a == NULL) return;
242  mpz_clear((mpz_ptr) *a);
243  omFreeBin((void *) *a, gmp_nrz_bin);
244  *a = NULL;
245 }
246 
247 number nrnCopy(number a, const coeffs)
248 {
249  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
250  mpz_init_set(erg, (mpz_ptr) a);
251  return (number) erg;
252 }
253 
254 int nrnSize(number a, const coeffs)
255 {
256  if (a == NULL) return 0;
257  return sizeof(mpz_t);
258 }
259 
260 /*
261  * convert a number to int
262  */
263 long nrnInt(number &n, const coeffs)
264 {
265  return mpz_get_si((mpz_ptr) n);
266 }
267 
268 /*
269  * Multiply two numbers
270  */
271 number nrnMult(number a, number b, const coeffs r)
272 {
273  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
274  mpz_init(erg);
275  mpz_mul(erg, (mpz_ptr)a, (mpz_ptr) b);
276  mpz_mod(erg, erg, r->modNumber);
277  return (number) erg;
278 }
279 
280 void nrnPower(number a, int i, number * result, const coeffs r)
281 {
282  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
283  mpz_init(erg);
284  mpz_powm_ui(erg, (mpz_ptr)a, i, r->modNumber);
285  *result = (number) erg;
286 }
287 
288 number nrnAdd(number a, number b, const coeffs r)
289 {
290  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
291  mpz_init(erg);
292  mpz_add(erg, (mpz_ptr)a, (mpz_ptr) b);
293  mpz_mod(erg, erg, r->modNumber);
294  return (number) erg;
295 }
296 
297 number nrnSub(number a, number b, const coeffs r)
298 {
299  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
300  mpz_init(erg);
301  mpz_sub(erg, (mpz_ptr)a, (mpz_ptr) b);
302  mpz_mod(erg, erg, r->modNumber);
303  return (number) erg;
304 }
305 
306 number nrnNeg(number c, const coeffs r)
307 {
308  if( !nrnIsZero(c, r) )
309  // Attention: This method operates in-place.
310  mpz_sub((mpz_ptr)c, r->modNumber, (mpz_ptr)c);
311  return c;
312 }
313 
314 number nrnInvers(number c, const coeffs r)
315 {
316  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
317  mpz_init(erg);
318  mpz_invert(erg, (mpz_ptr)c, r->modNumber);
319  return (number) erg;
320 }
321 
322 /*
323  * Give the smallest k, such that a * x = k = b * y has a solution
324  * TODO: lcm(gcd,gcd) better than gcd(lcm) ?
325  */
326 number nrnLcm(number a, number b, const coeffs r)
327 {
328  number erg = nrnGcd(NULL, a, r);
329  number tmp = nrnGcd(NULL, b, r);
330  mpz_lcm((mpz_ptr)erg, (mpz_ptr)erg, (mpz_ptr)tmp);
331  nrnDelete(&tmp, r);
332  return (number)erg;
333 }
334 
335 /*
336  * Give the largest k, such that a = x * k, b = y * k has
337  * a solution.
338  */
339 number nrnGcd(number a, number b, const coeffs r)
340 {
341  if ((a == NULL) && (b == NULL)) return nrnInit(0,r);
342  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
343  mpz_init_set(erg, r->modNumber);
344  if (a != NULL) mpz_gcd(erg, erg, (mpz_ptr)a);
345  if (b != NULL) mpz_gcd(erg, erg, (mpz_ptr)b);
346  return (number)erg;
347 }
348 
349 /* Not needed any more, but may have room for improvement
350  number nrnGcd3(number a,number b, number c,ring r)
351 {
352  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
353  mpz_init(erg);
354  if (a == NULL) a = (number)r->modNumber;
355  if (b == NULL) b = (number)r->modNumber;
356  if (c == NULL) c = (number)r->modNumber;
357  mpz_gcd(erg, (mpz_ptr)a, (mpz_ptr)b);
358  mpz_gcd(erg, erg, (mpz_ptr)c);
359  mpz_gcd(erg, erg, r->modNumber);
360  return (number)erg;
361 }
362 */
363 
364 /*
365  * Give the largest k, such that a = x * k, b = y * k has
366  * a solution and r, s, s.t. k = s*a + t*b
367  * CF: careful: ExtGcd is wrong as implemented (or at least may not
368  * give you what you want:
369  * ExtGcd(5, 10 modulo 12):
370  * the gcdext will return 5 = 1*5 + 0*10
371  * however, mod 12, the gcd should be 1
372  */
373 number nrnExtGcd(number a, number b, number *s, number *t, const coeffs r)
374 {
375  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
376  mpz_ptr bs = (mpz_ptr)omAllocBin(gmp_nrz_bin);
377  mpz_ptr bt = (mpz_ptr)omAllocBin(gmp_nrz_bin);
378  mpz_init(erg);
379  mpz_init(bs);
380  mpz_init(bt);
381  mpz_gcdext(erg, bs, bt, (mpz_ptr)a, (mpz_ptr)b);
382  mpz_mod(bs, bs, r->modNumber);
383  mpz_mod(bt, bt, r->modNumber);
384  *s = (number)bs;
385  *t = (number)bt;
386  return (number)erg;
387 }
388 /* XExtGcd returns a unimodular matrix ((s,t)(u,v)) sth.
389  * (a,b)^t ((st)(uv)) = (g,0)^t
390  * Beware, the ExtGcd will not necessaairly do this.
391  * Problem: if g = as+bt then (in Z/nZ) it follows NOT that
392  * 1 = (a/g)s + (b/g) t
393  * due to the zero divisors.
394  */
395 
396 //#define CF_DEB;
397 number nrnXExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
398 {
399  number xx;
400 #ifdef CF_DEB
401  StringSetS("XExtGcd of ");
402  nrnWrite(a, r);
403  StringAppendS("\t");
404  nrnWrite(b, r);
405  StringAppendS(" modulo ");
406  nrnWrite(xx = (number)r->modNumber, r);
407  Print("%s\n", StringEndS());
408 #endif
409 
410  mpz_ptr one = (mpz_ptr)omAllocBin(gmp_nrz_bin);
411  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
412  mpz_ptr bs = (mpz_ptr)omAllocBin(gmp_nrz_bin);
413  mpz_ptr bt = (mpz_ptr)omAllocBin(gmp_nrz_bin);
414  mpz_ptr bu = (mpz_ptr)omAllocBin(gmp_nrz_bin);
415  mpz_ptr bv = (mpz_ptr)omAllocBin(gmp_nrz_bin);
416  mpz_init(erg);
417  mpz_init(one);
418  mpz_init_set(bs, (mpz_ptr) a);
419  mpz_init_set(bt, (mpz_ptr) b);
420  mpz_init(bu);
421  mpz_init(bv);
422  mpz_gcd(erg, bs, bt);
423 
424 #ifdef CF_DEB
425  StringSetS("1st gcd:");
426  nrnWrite(xx= (number)erg, r);
427 #endif
428 
429  mpz_gcd(erg, erg, r->modNumber);
430 
431  mpz_div(bs, bs, erg);
432  mpz_div(bt, bt, erg);
433 
434 #ifdef CF_DEB
435  Print("%s\n", StringEndS());
436  StringSetS("xgcd: ");
437 #endif
438 
439  mpz_gcdext(one, bu, bv, bs, bt);
440  number ui = nrnGetUnit(xx = (number) one, r);
441 #ifdef CF_DEB
442  n_Write(xx, r);
443  StringAppendS("\t");
444  n_Write(ui, r);
445  Print("%s\n", StringEndS());
446 #endif
447  nrnDelete(&xx, r);
448  if (!nrnIsOne(ui, r)) {
449 #ifdef CF_DEB
450  Print("Scaling\n");
451 #endif
452  number uii = nrnInvers(ui, r);
453  nrnDelete(&ui, r);
454  ui = uii;
455  mpz_ptr uu = (mpz_ptr)omAllocBin(gmp_nrz_bin);
456  mpz_init_set(uu, (mpz_ptr)ui);
457  mpz_mul(bu, bu, uu);
458  mpz_mul(bv, bv, uu);
459  mpz_clear(uu);
460  omFreeBin(uu, gmp_nrz_bin);
461  }
462  nrnDelete(&ui, r);
463 #ifdef CF_DEB
464  StringSetS("xgcd");
465  nrnWrite(xx= (number)bs, r);
466  StringAppendS("*");
467  nrnWrite(xx= (number)bu, r);
468  StringAppendS(" + ");
469  nrnWrite(xx= (number)bt, r);
470  StringAppendS("*");
471  nrnWrite(xx= (number)bv, r);
472  Print("%s\n", StringEndS());
473 #endif
474 
475  mpz_mod(bs, bs, r->modNumber);
476  mpz_mod(bt, bt, r->modNumber);
477  mpz_mod(bu, bu, r->modNumber);
478  mpz_mod(bv, bv, r->modNumber);
479  *s = (number)bu;
480  *t = (number)bv;
481  *u = (number)bt;
482  *u = nrnNeg(*u, r);
483  *v = (number)bs;
484  return (number)erg;
485 }
486 
487 
488 BOOLEAN nrnIsZero(number a, const coeffs)
489 {
490 #ifdef LDEBUG
491  if (a == NULL) return FALSE;
492 #endif
493  return 0 == mpz_cmpabs_ui((mpz_ptr)a, 0);
494 }
495 
496 BOOLEAN nrnIsOne(number a, const coeffs)
497 {
498 #ifdef LDEBUG
499  if (a == NULL) return FALSE;
500 #endif
501  return 0 == mpz_cmp_si((mpz_ptr)a, 1);
502 }
503 
504 BOOLEAN nrnIsMOne(number a, const coeffs r)
505 {
506 #ifdef LDEBUG
507  if (a == NULL) return FALSE;
508 #endif
509  mpz_t t; mpz_init_set(t, (mpz_ptr)a);
510  mpz_add_ui(t, t, 1);
511  bool erg = (0 == mpz_cmp(t, r->modNumber));
512  mpz_clear(t);
513  return erg;
514 }
515 
516 BOOLEAN nrnEqual(number a, number b, const coeffs)
517 {
518  return 0 == mpz_cmp((mpz_ptr)a, (mpz_ptr)b);
519 }
520 
521 BOOLEAN nrnGreater(number a, number b, const coeffs)
522 {
523  return 0 < mpz_cmp((mpz_ptr)a, (mpz_ptr)b);
524 }
525 
527 {
528  return 0 < mpz_cmp_si((mpz_ptr)k, 0);
529 }
530 
531 BOOLEAN nrnIsUnit(number a, const coeffs r)
532 {
533  number tmp = nrnGcd(a, (number)r->modNumber, r);
534  bool res = nrnIsOne(tmp, r);
535  nrnDelete(&tmp, NULL);
536  return res;
537 }
538 
539 number nrnGetUnit(number k, const coeffs r)
540 {
541  if (mpz_divisible_p(r->modNumber, (mpz_ptr)k)) return nrnInit(1,r);
542 
543  mpz_ptr unit = (mpz_ptr)nrnGcd(k, 0, r);
544  mpz_tdiv_q(unit, (mpz_ptr)k, unit);
545  mpz_ptr gcd = (mpz_ptr)nrnGcd((number)unit, 0, r);
546  if (!nrnIsOne((number)gcd,r))
547  {
548  mpz_ptr ctmp;
549  // tmp := unit^2
550  mpz_ptr tmp = (mpz_ptr) nrnMult((number) unit,(number) unit,r);
551  // gcd_new := gcd(tmp, 0)
552  mpz_ptr gcd_new = (mpz_ptr) nrnGcd((number) tmp, 0, r);
553  while (!nrnEqual((number) gcd_new,(number) gcd,r))
554  {
555  // gcd := gcd_new
556  ctmp = gcd;
557  gcd = gcd_new;
558  gcd_new = ctmp;
559  // tmp := tmp * unit
560  mpz_mul(tmp, tmp, unit);
561  mpz_mod(tmp, tmp, r->modNumber);
562  // gcd_new := gcd(tmp, 0)
563  mpz_gcd(gcd_new, tmp, r->modNumber);
564  }
565  // unit := unit + modNumber / gcd_new
566  mpz_tdiv_q(tmp, r->modNumber, gcd_new);
567  mpz_add(unit, unit, tmp);
568  mpz_mod(unit, unit, r->modNumber);
569  nrnDelete((number*) &gcd_new, NULL);
570  nrnDelete((number*) &tmp, NULL);
571  }
572  nrnDelete((number*) &gcd, NULL);
573  return (number)unit;
574 }
575 
576 number nrnAnn(number k, const coeffs r)
577 {
578  mpz_ptr tmp = (mpz_ptr) omAllocBin(gmp_nrz_bin);
579  mpz_init(tmp);
580  mpz_gcd(tmp, (mpz_ptr) k, r->modNumber);
581  if (mpz_cmp_si(tmp, 1)==0) {
582  mpz_set_si(tmp, 0);
583  return (number) tmp;
584  }
585  mpz_divexact(tmp, r->modNumber, tmp);
586  return (number) tmp;
587 }
588 
589 BOOLEAN nrnDivBy(number a, number b, const coeffs r)
590 {
591  if (a == NULL)
592  return mpz_divisible_p(r->modNumber, (mpz_ptr)b);
593  else
594  { /* b divides a iff b/gcd(a, b) is a unit in the given ring: */
595  number n = nrnGcd(a, b, r);
596  mpz_tdiv_q((mpz_ptr)n, (mpz_ptr)b, (mpz_ptr)n);
597  bool result = nrnIsUnit(n, r);
598  nrnDelete(&n, NULL);
599  return result;
600  }
601 }
602 
603 int nrnDivComp(number a, number b, const coeffs r)
604 {
605  if (nrnEqual(a, b,r)) return 2;
606  if (mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b)) return -1;
607  if (mpz_divisible_p((mpz_ptr) b, (mpz_ptr) a)) return 1;
608  return 0;
609 }
610 
611 number nrnDiv(number a, number b, const coeffs r)
612 {
613  if (a == NULL) a = (number)r->modNumber;
614  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
615  mpz_init(erg);
616  if (mpz_divisible_p((mpz_ptr)a, (mpz_ptr)b))
617  {
618  mpz_divexact(erg, (mpz_ptr)a, (mpz_ptr)b);
619  return (number)erg;
620  }
621  else
622  {
623  mpz_ptr gcd = (mpz_ptr)nrnGcd(a, b, r);
624  mpz_divexact(erg, (mpz_ptr)b, gcd);
625  if (!nrnIsUnit((number)erg, r))
626  {
627  WerrorS("Division not possible, even by cancelling zero divisors.");
628  WerrorS("Result is integer division without remainder.");
629  mpz_tdiv_q(erg, (mpz_ptr) a, (mpz_ptr) b);
630  nrnDelete((number*) &gcd, NULL);
631  return (number)erg;
632  }
633  // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
634  mpz_ptr tmp = (mpz_ptr)nrnInvers((number) erg,r);
635  mpz_divexact(erg, (mpz_ptr)a, gcd);
636  mpz_mul(erg, erg, tmp);
637  nrnDelete((number*) &gcd, NULL);
638  nrnDelete((number*) &tmp, NULL);
639  mpz_mod(erg, erg, r->modNumber);
640  return (number)erg;
641  }
642 }
643 
644 number nrnMod(number a, number b, const coeffs r)
645 {
646  /*
647  We need to return the number rr which is uniquely determined by the
648  following two properties:
649  (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
650  (2) There exists some k in the integers Z such that a = k * b + rr.
651  Consider g := gcd(n, |b|). Note that then |b|/g is a unit in Z/n.
652  Now, there are three cases:
653  (a) g = 1
654  Then |b| is a unit in Z/n, i.e. |b| (and also b) divides a.
655  Thus rr = 0.
656  (b) g <> 1 and g divides a
657  Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
658  (c) g <> 1 and g does not divide a
659  Then denote the division with remainder of a by g as this:
660  a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
661  fulfills (1) and (2), i.e. rr := t is the correct result. Hence
662  in this third case, rr is the remainder of division of a by g in Z.
663  Remark: according to mpz_mod: a,b are always non-negative
664  */
665  mpz_ptr g = (mpz_ptr)omAllocBin(gmp_nrz_bin);
666  mpz_ptr rr = (mpz_ptr)omAllocBin(gmp_nrz_bin);
667  mpz_init(g);
668  mpz_init_set_si(rr, 0);
669  mpz_gcd(g, (mpz_ptr)r->modNumber, (mpz_ptr)b); // g is now as above
670  if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(rr, (mpz_ptr)a, g); // the case g <> 1
671  mpz_clear(g);
673  return (number)rr;
674 }
675 
676 number nrnIntDiv(number a, number b, const coeffs r)
677 {
678  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
679  mpz_init(erg);
680  if (a == NULL) a = (number)r->modNumber;
681  mpz_tdiv_q(erg, (mpz_ptr)a, (mpz_ptr)b);
682  return (number)erg;
683 }
684 
685 /* CF: note that Z/nZ has (at least) two distinct euclidean structures
686  * 1st phi(a) := (a mod n) which is just the structure directly
687  * inherited from Z
688  * 2nd phi(a) := gcd(a, n)
689  * The 1st version is probably faster as everything just comes from Z,
690  * but the 2nd version behaves nicely wrt. to quotient operations
691  * and HNF and such. In agreement with nrnMod we imlement the 2nd here
692  *
693  * For quotrem note that if b exactly divides a, then
694  * min(v_p(a), v_p(n)) >= min(v_p(b), v_p(n))
695  * so if we divide a and b by g:= gcd(a,b,n), then b becomes a
696  * unit mod n/g.
697  * Thus we 1st compute the remainder (similar to nrnMod) and then
698  * the exact quotient.
699  */
700 number nrnQuotRem(number a, number b, number * rem, const coeffs r)
701 {
702  mpz_t g, aa, bb;
703  mpz_ptr qq = (mpz_ptr)omAllocBin(gmp_nrz_bin);
704  mpz_ptr rr = (mpz_ptr)omAllocBin(gmp_nrz_bin);
705  mpz_init(qq);
706  mpz_init(rr);
707  mpz_init(g);
708  mpz_init_set(aa, (mpz_ptr)a);
709  mpz_init_set(bb, (mpz_ptr)b);
710 
711  mpz_gcd(g, bb, r->modNumber);
712  mpz_mod(rr, aa, g);
713  mpz_sub(aa, aa, rr);
714  mpz_gcd(g, aa, g);
715  mpz_div(aa, aa, g);
716  mpz_div(bb, bb, g);
717  mpz_div(g, r->modNumber, g);
718  mpz_invert(g, bb, g);
719  mpz_mul(qq, aa, g);
720  if (rem)
721  *rem = (number)rr;
722  else {
723  mpz_clear(rr);
724  omFreeBin(rr, gmp_nrz_bin);
725  }
726  mpz_clear(g);
727  mpz_clear(aa);
728  mpz_clear(bb);
729  return (number) qq;
730 }
731 
732 /*
733  * Helper function for computing the module
734  */
735 
736 mpz_ptr nrnMapCoef = NULL;
737 
738 number nrnMapModN(number from, const coeffs /*src*/, const coeffs dst)
739 {
740  return nrnMult(from, (number) nrnMapCoef, dst);
741 }
742 
743 number nrnMap2toM(number from, const coeffs /*src*/, const coeffs dst)
744 {
745  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
746  mpz_init(erg);
747  mpz_mul_ui(erg, nrnMapCoef, (unsigned long)from);
748  mpz_mod(erg, erg, dst->modNumber);
749  return (number)erg;
750 }
751 
752 number nrnMapZp(number from, const coeffs /*src*/, const coeffs dst)
753 {
754  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
755  mpz_init(erg);
756  // TODO: use npInt(...)
757  mpz_mul_si(erg, nrnMapCoef, (unsigned long)from);
758  mpz_mod(erg, erg, dst->modNumber);
759  return (number)erg;
760 }
761 
762 number nrnMapGMP(number from, const coeffs /*src*/, const coeffs dst)
763 {
764  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
765  mpz_init(erg);
766  mpz_mod(erg, (mpz_ptr)from, dst->modNumber);
767  return (number)erg;
768 }
769 
770 #if SI_INTEGER_VARIANT==3
771 number nrnMapZ(number from, const coeffs /*src*/, const coeffs dst)
772 {
773  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
774  if (n_Z_IS_SMALL(from))
775  mpz_init_set_si(erg, SR_TO_INT(from));
776  else
777  mpz_init_set(erg, (mpz_ptr) from);
778  mpz_mod(erg, erg, dst->modNumber);
779  return (number)erg;
780 }
781 #elif SI_INTEGER_VARIANT==2
782 
783 number nrnMapZ(number from, const coeffs src, const coeffs dst)
784 {
785  if (SR_HDL(from) & SR_INT)
786  {
787  long f_i=SR_TO_INT(from);
788  return nrnInit(f_i,dst);
789  }
790  return nrnMapGMP(from,src,dst);
791 }
792 #elif SI_INTEGER_VARIANT==1
793 number nrnMapZ(number from, const coeffs src, const coeffs dst)
794 {
795  return nrnMapQ(from,src,dst);
796 }
797 #endif
798 #if SI_INTEGER_VARIANT!=2
799 void nrnWrite (number &a, const coeffs)
800 {
801  char *s,*z;
802  if (a==NULL)
803  {
804  StringAppendS("o");
805  }
806  else
807  {
808  int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2;
809  s=(char*)omAlloc(l);
810  z=mpz_get_str(s,10,(mpz_ptr) a);
811  StringAppendS(z);
812  omFreeSize((ADDRESS)s,l);
813  }
814 }
815 #endif
816 
817 number nrnMapQ(number from, const coeffs src, const coeffs dst)
818 {
819  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
820  mpz_init(erg);
821  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); // ?
822  mpz_mod(erg, erg, dst->modNumber);
823  return (number)erg;
824 }
825 
826 nMapFunc nrnSetMap(const coeffs src, const coeffs dst)
827 {
828  /* dst = nrn */
829  if ((src->rep==n_rep_gmp) && nCoeff_is_Ring_Z(src))
830  {
831  return nrnMapZ;
832  }
833  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
834  {
835  return nrnMapZ;
836  }
837  if (src->rep==n_rep_gap_rat) /*&& nCoeff_is_Q(src)) or Z*/
838  {
839  return nrnMapQ;
840  }
841  // Some type of Z/n ring / field
842  if (nCoeff_is_Ring_ModN(src) || nCoeff_is_Ring_PtoM(src) ||
843  nCoeff_is_Ring_2toM(src) || nCoeff_is_Zp(src))
844  {
845  if ( (!nCoeff_is_Zp(src))
846  && (mpz_cmp(src->modBase, dst->modBase) == 0)
847  && (src->modExponent == dst->modExponent)) return nrnMapGMP;
848  else
849  {
850  mpz_ptr nrnMapModul = (mpz_ptr) omAllocBin(gmp_nrz_bin);
851  // Computing the n of Z/n
852  if (nCoeff_is_Zp(src))
853  {
854  mpz_init_set_si(nrnMapModul, src->ch);
855  }
856  else
857  {
858  mpz_init(nrnMapModul);
859  mpz_set(nrnMapModul, src->modNumber);
860  }
861  // nrnMapCoef = 1 in dst if dst is a subring of src
862  // nrnMapCoef = 0 in dst / src if src is a subring of dst
863  if (nrnMapCoef == NULL)
864  {
865  nrnMapCoef = (mpz_ptr) omAllocBin(gmp_nrz_bin);
866  mpz_init(nrnMapCoef);
867  }
868  if (mpz_divisible_p(nrnMapModul, dst->modNumber))
869  {
870  mpz_set_si(nrnMapCoef, 1);
871  }
872  else
873  if (nrnDivBy(NULL, (number) nrnMapModul,dst))
874  {
875  mpz_divexact(nrnMapCoef, dst->modNumber, nrnMapModul);
876  mpz_ptr tmp = dst->modNumber;
877  dst->modNumber = nrnMapModul;
878  if (!nrnIsUnit((number) nrnMapCoef,dst))
879  {
880  dst->modNumber = tmp;
881  nrnDelete((number*) &nrnMapModul, dst);
882  return NULL;
883  }
884  mpz_ptr inv = (mpz_ptr) nrnInvers((number) nrnMapCoef,dst);
885  dst->modNumber = tmp;
886  mpz_mul(nrnMapCoef, nrnMapCoef, inv);
887  mpz_mod(nrnMapCoef, nrnMapCoef, dst->modNumber);
888  nrnDelete((number*) &inv, dst);
889  }
890  else
891  {
892  nrnDelete((number*) &nrnMapModul, dst);
893  return NULL;
894  }
895  nrnDelete((number*) &nrnMapModul, dst);
896  if (nCoeff_is_Ring_2toM(src))
897  return nrnMap2toM;
898  else if (nCoeff_is_Zp(src))
899  return nrnMapZp;
900  else
901  return nrnMapModN;
902  }
903  }
904  return NULL; // default
905 }
906 
907 /*
908  * set the exponent (allocate and init tables) (TODO)
909  */
910 
911 void nrnSetExp(unsigned long m, coeffs r)
912 {
913  /* clean up former stuff */
914  if (r->modNumber != NULL) mpz_clear(r->modNumber);
915 
916  r->modExponent= m;
917  r->modNumber = (mpz_ptr)omAllocBin(gmp_nrz_bin);
918  mpz_init_set (r->modNumber, r->modBase);
919  mpz_pow_ui (r->modNumber, r->modNumber, m);
920 }
921 
922 /* We expect this ring to be Z/n^m for some m > 0 and for some n > 2 which is not a prime. */
923 void nrnInitExp(unsigned long m, coeffs r)
924 {
925  nrnSetExp(m, r);
926  assume (r->modNumber != NULL);
927 //CF: in geenral, the modulus is computed somewhere. I don't want to
928 // check it's size before I construct the best ring.
929 // if (mpz_cmp_ui(r->modNumber,2) <= 0)
930 // WarnS("nrnInitExp failed (m in Z/m too small)");
931 }
932 
933 #ifdef LDEBUG
934 BOOLEAN nrnDBTest (number a, const char *, const int, const coeffs r)
935 {
936  if (a==NULL) return TRUE;
937  if ( (mpz_cmp_si((mpz_ptr) a, 0) < 0) || (mpz_cmp((mpz_ptr) a, r->modNumber) > 0) )
938  {
939  return FALSE;
940  }
941  return TRUE;
942 }
943 #endif
944 
945 /*2
946 * extracts a long integer from s, returns the rest (COPY FROM longrat0.cc)
947 */
948 static const char * nlCPEatLongC(char *s, mpz_ptr i)
949 {
950  const char * start=s;
951  if (!(*s >= '0' && *s <= '9'))
952  {
953  mpz_init_set_si(i, 1);
954  return s;
955  }
956  mpz_init(i);
957  while (*s >= '0' && *s <= '9') s++;
958  if (*s=='\0')
959  {
960  mpz_set_str(i,start,10);
961  }
962  else
963  {
964  char c=*s;
965  *s='\0';
966  mpz_set_str(i,start,10);
967  *s=c;
968  }
969  return s;
970 }
971 
972 const char * nrnRead (const char *s, number *a, const coeffs r)
973 {
974  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
975  {
976  s = nlCPEatLongC((char *)s, z);
977  }
978  mpz_mod(z, z, r->modNumber);
979  *a = (number) z;
980  return s;
981 }
982 #endif
983 /* #ifdef HAVE_RINGS */
mpz_ptr base
Definition: rmodulon.h:18
void nrnSetExp(unsigned long c, const coeffs r)
Definition: rmodulon.cc:911
mpz_t z
Definition: longrat.h:48
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
const CanonicalForm int s
Definition: facAbsFact.cc:55
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_ModN(const coeffs r)
Definition: coeffs.h:744
void mpz_mul_si(mpz_ptr r, mpz_srcptr s, long int si)
Definition: longrat.cc:181
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:574
const poly a
Definition: syzextra.cc:212
omBin_t * omBin
Definition: omStructs.h:12
#define Print
Definition: emacs.cc:83
only used if HAVE_RINGS is defined: ?
Definition: coeffs.h:43
static FORCE_INLINE BOOLEAN nCoeff_is_Zp(const coeffs r)
Definition: coeffs.h:818
number nrnMult(number a, number b, const coeffs r)
Definition: rmodulon.cc:271
number nrnQuotRem(number a, number b, number *s, const coeffs r)
Definition: rmodulon.cc:700
#define FALSE
Definition: auxiliary.h:140
return P p
Definition: myNF.cc:203
f
Definition: cfModGcd.cc:4022
static const char * nlCPEatLongC(char *s, mpz_ptr i)
Definition: rmodulon.cc:948
number nrnExtGcd(number a, number b, number *s, number *t, const coeffs r)
Definition: rmodulon.cc:373
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_Z(const coeffs r)
Definition: coeffs.h:750
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
number nrnLcm(number a, number b, const coeffs r)
Definition: rmodulon.cc:326
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_2toM(const coeffs r)
Definition: coeffs.h:741
(), see rinteger.h, new impl.
Definition: coeffs.h:111
BOOLEAN nrnIsMOne(number a, const coeffs r)
Definition: rmodulon.cc:504
BOOLEAN nrnDivBy(number a, number b, const coeffs r)
Definition: rmodulon.cc:589
#define TRUE
Definition: auxiliary.h:144
mpz_ptr nrnMapCoef
Definition: rmodulon.cc:736
void * ADDRESS
Definition: auxiliary.h:161
g
Definition: cfModGcd.cc:4031
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
char * StringEndS()
Definition: reporter.cc:151
void nlGMP(number &i, number n, const coeffs r)
Definition: longrat.cc:1410
static BOOLEAN nrnCoeffsEqual(const coeffs r, n_coeffType n, void *parameter)
Definition: rmodulon.cc:91
#define omAlloc(size)
Definition: omAllocDecl.h:210
number nrnMod(number a, number b, const coeffs r)
Definition: rmodulon.cc:644
number nrnMapZ(number from, const coeffs src, const coeffs dst)
Definition: rmodulon.cc:783
BOOLEAN nrnGreaterZero(number k, const coeffs r)
Definition: rmodulon.cc:526
static const n_coeffType ID2
Definition: rmodulon.cc:28
void nrnInitExp(unsigned long c, const coeffs r)
Definition: rmodulon.cc:923
number nrnDiv(number a, number b, const coeffs r)
Definition: rmodulon.cc:611
poly res
Definition: myNF.cc:322
void nrnDelete(number *a, const coeffs r)
Definition: rmodulon.cc:239
number nrnAnn(number a, const coeffs r)
Definition: rmodulon.cc:576
mpz_t n
Definition: longrat.h:49
const ring r
Definition: syzextra.cc:208
static char * nrnCoeffString(const coeffs r)
Definition: rmodulon.cc:97
Coefficient rings, fields and other domains suitable for Singular polynomials.
number nrnMap2toM(number from, const coeffs, const coeffs dst)
Definition: rmodulon.cc:743
BOOLEAN nrnIsUnit(number a, const coeffs r)
Definition: rmodulon.cc:531
only used if HAVE_RINGS is defined: ?
Definition: coeffs.h:44
BOOLEAN nrnInitChar(coeffs r, void *p)
Definition: rmodulon.cc:163
const char * nrnRead(const char *s, number *a, const coeffs r)
Definition: rmodulon.cc:972
number nrnMapQ(number from, const coeffs src, const coeffs dst)
Definition: rmodulon.cc:817
#define assume(x)
Definition: mod2.h:405
The main handler for Singular numbers which are suitable for Singular polynomials.
void StringSetS(const char *st)
Definition: reporter.cc:128
void StringAppendS(const char *st)
Definition: reporter.cc:107
const ExtensionInfo & info
< [in] sqrfree poly
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:72
BOOLEAN nrnDBTest(number a, const char *f, const int l, const coeffs r)
Definition: rmodulon.cc:934
static FORCE_INLINE void n_Write(number &n, const coeffs r, const BOOLEAN bShortOut=TRUE)
Definition: coeffs.h:592
int nrnSize(number a, const coeffs r)
Definition: rmodulon.cc:254
omBin gmp_nrz_bin
Definition: rintegers.cc:80
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
number nrnMapGMP(number from, const coeffs, const coeffs dst)
Definition: rmodulon.cc:762
int i
Definition: cfEzgcd.cc:123
(mpz_ptr), see rmodulon,h
Definition: coeffs.h:114
static const n_coeffType ID
Our Type!
Definition: rmodulon.cc:27
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:422
void nrnPower(number a, int i, number *result, const coeffs r)
Definition: rmodulon.cc:280
number nrnAdd(number a, number b, const coeffs r)
Definition: rmodulon.cc:288
#define SR_TO_INT(SR)
Definition: longrat.h:67
(number), see longrat.h
Definition: coeffs.h:110
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
BOOLEAN nrnGreater(number a, number b, const coeffs r)
Definition: rmodulon.cc:521
static void nrnKillChar(coeffs r)
Definition: rmodulon.cc:109
n_coeffType
Definition: coeffs.h:27
#define NULL
Definition: omList.c:10
int gcd(int a, int b)
Definition: walkSupport.cc:839
#define SR_INT
Definition: longrat.h:65
number nrnXExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
Definition: rmodulon.cc:397
number nrnCopy(number a, const coeffs r)
Definition: rmodulon.cc:247
number nrnGetUnit(number a, const coeffs r)
Definition: rmodulon.cc:539
number nrnMapZp(number from, const coeffs, const coeffs dst)
Definition: rmodulon.cc:752
number nrnSub(number a, number b, const coeffs r)
Definition: rmodulon.cc:297
BOOLEAN nrnIsOne(number a, const coeffs r)
Definition: rmodulon.cc:496
BOOLEAN nrnIsZero(number a, const coeffs r)
Definition: rmodulon.cc:488
int nrnDivComp(number a, number b, const coeffs r)
Definition: rmodulon.cc:603
nMapFunc nrnSetMap(const coeffs src, const coeffs dst)
Definition: rmodulon.cc:826
#define SR_HDL(A)
Definition: tgb.cc:35
BOOLEAN nrnEqual(number a, number b, const coeffs r)
Definition: rmodulon.cc:516
#define nrnWrite
Definition: rmodulon.cc:63
number nrnInit(long i, const coeffs r)
Definition: rmodulon.cc:231
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
number nrnInvers(number c, const coeffs r)
Definition: rmodulon.cc:314
int BOOLEAN
Definition: auxiliary.h:131
const poly b
Definition: syzextra.cc:213
void nrnCoeffWrite(const coeffs r, BOOLEAN details)
Definition: rmodulon.cc:81
int l
Definition: cfEzgcd.cc:94
return result
Definition: facAbsBiFact.cc:76
long nrnInt(number &n, const coeffs r)
Definition: rmodulon.cc:263
number nrnIntDiv(number a, number b, const coeffs r)
Definition: rmodulon.cc:676
number nrnNeg(number c, const coeffs r)
Definition: rmodulon.cc:306
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:327
number nrnGcd(number a, number b, const coeffs r)
Definition: rmodulon.cc:339
coeffs nrnQuot1(number c, const coeffs r)
Definition: rmodulon.cc:117
number nrnMapModN(number from, const coeffs, const coeffs dst)
Definition: rmodulon.cc:738