source: git/libpolys/coeffs/rmodulon.cc

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