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