source: git/libpolys/coeffs/rmodulon.cc

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