source: git/libpolys/coeffs/rmodulon.cc @ 6a1aa7

spielwiese
Last change on this file since 6a1aa7 was 6a1aa7, checked in by Hans Schoenemann <hannes@…>, 9 years ago
preparation for more integer representations n_Init_bigint: substituted by n_setMap/nMap introduced cf->rep for different data representation introduced several (experimental) integer variants via #define in libpolys/misc/auxiliary.h.in
  • Property mode set to 100644
File size: 18.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: numbers modulo n
6*/
7
8
9
10
11#include <misc/auxiliary.h>
12
13#ifdef HAVE_RINGS
14
15#include <misc/mylimits.h>
16#include <coeffs/coeffs.h>
17#include <reporter/reporter.h>
18#include <omalloc/omalloc.h>
19#include <coeffs/numbers.h>
20#include <coeffs/longrat.h>
21#include <coeffs/mpr_complex.h>
22#include <coeffs/rmodulon.h>
23#include "si_gmp.h"
24
25#include <string.h>
26
27/// Our Type!
28static const n_coeffType ID = n_Zn;
29static const n_coeffType ID2 = n_Znm;
30
31extern omBin gmp_nrz_bin;
32
33void    nrnCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
34{
35  long l = (long)mpz_sizeinbase(r->modBase, 10) + 2;
36  char* s = (char*) omAlloc(l);
37  s= mpz_get_str (s, 10, r->modBase);
38  if (nCoeff_is_Ring_ModN(r)) Print("//   coeff. ring is : Z/%s\n", s);
39  else if (nCoeff_is_Ring_PtoM(r)) Print("//   coeff. ring is : Z/%s^%lu\n", s, r->modExponent);
40  omFreeSize((ADDRESS)s, l);
41}
42
43static BOOLEAN nrnCoeffsEqual(const coeffs r, n_coeffType n, void * parameter)
44{
45  /* test, if r is an instance of nInitCoeffs(n,parameter) */
46  return (n==n_Zn) && (mpz_cmp_ui(r->modNumber,(long)parameter)==0);
47}
48
49static char* nrnCoeffString(const coeffs r)
50{
51  long l = (long)mpz_sizeinbase(r->modBase, 10) +2;
52  char* b = (char*) omAlloc(l);
53  b= mpz_get_str (b, 10, r->modBase);
54  char* s = (char*) omAlloc(7+2+10+l);
55  if (nCoeff_is_Ring_ModN(r)) sprintf(s,"integer,%s",b);
56  else /*if (nCoeff_is_Ring_PtoM(r))*/ sprintf(s,"integer,%s^%lu",b,r->modExponent);
57  omFreeSize(b,l);
58  return s;
59}
60
61coeffs nrnQuot1(number c, const coeffs r)
62{
63    coeffs rr;
64    int ch = r->cfInt(c, r);
65    mpz_t a,b;
66    mpz_init_set(a, r->modNumber);
67    mpz_init_set_ui(b, ch);
68    int_number gcd;
69    gcd = (int_number) omAlloc(sizeof(mpz_t));
70    mpz_init(gcd);
71    mpz_gcd(gcd, a,b);
72    if(mpz_cmp_ui(gcd, 1) == 0)
73        {
74            WerrorS("constant in q-ideal is coprime to modulus in ground ring");
75            WerrorS("Unable to create qring!");
76            return NULL;
77        }
78    if(r->modExponent == 1)
79    {
80        ZnmInfo info;
81        info.base = gcd;
82        info.exp = (unsigned long) 1;
83        rr = nInitChar(n_Zn, (void*)&info);
84    }
85    else
86    {
87        ZnmInfo info;
88        info.base = r->modBase;
89        int kNew = 1;
90        mpz_t baseTokNew;
91        mpz_init(baseTokNew);
92        mpz_set(baseTokNew, r->modBase);
93        while(mpz_cmp(gcd, baseTokNew) > 0)
94        {
95          kNew++;
96          mpz_mul(baseTokNew, baseTokNew, r->modBase);
97        }
98        //printf("\nkNew = %i\n",kNew);
99        info.exp = kNew;
100        mpz_clear(baseTokNew);
101        rr = nInitChar(n_Znm, (void*)&info);
102    }
103    return(rr);
104}
105
106static number nrnAnn(number b, const coeffs r);
107/* for initializing function pointers */
108BOOLEAN nrnInitChar (coeffs r, void* p)
109{
110  assume( (getCoeffType(r) == ID) || (getCoeffType (r) == ID2) );
111  ZnmInfo * info= (ZnmInfo *) p;
112  r->modBase= info->base;
113
114  nrnInitExp (info->exp, r);
115
116  /* next computation may yield wrong characteristic as r->modNumber
117     is a GMP number */
118  r->ch = mpz_get_ui(r->modNumber);
119
120  r->is_field=FALSE;
121  r->is_domain=FALSE;
122  r->rep=n_rep_gmp;
123
124
125  r->cfCoeffString = nrnCoeffString;
126
127  r->cfInit        = nrnInit;
128  r->cfDelete      = nrnDelete;
129  r->cfCopy        = nrnCopy;
130  r->cfSize        = nrnSize;
131  r->cfInt         = nrnInt;
132  r->cfAdd         = nrnAdd;
133  r->cfSub         = nrnSub;
134  r->cfMult        = nrnMult;
135  r->cfDiv         = nrnDiv;
136  r->cfAnn         = nrnAnn;
137  r->cfIntMod      = nrnMod;
138  r->cfExactDiv    = nrnDiv;
139  r->cfInpNeg         = nrnNeg;
140  r->cfInvers      = nrnInvers;
141  r->cfDivBy       = nrnDivBy;
142  r->cfDivComp     = nrnDivComp;
143  r->cfGreater     = nrnGreater;
144  r->cfEqual       = nrnEqual;
145  r->cfIsZero      = nrnIsZero;
146  r->cfIsOne       = nrnIsOne;
147  r->cfIsMOne      = nrnIsMOne;
148  r->cfGreaterZero = nrnGreaterZero;
149  r->cfWriteLong   = nrnWrite;
150  r->cfRead        = nrnRead;
151  r->cfPower       = nrnPower;
152  r->cfSetMap      = nrnSetMap;
153  r->cfNormalize   = ndNormalize;
154  r->cfLcm         = nrnLcm;
155  r->cfGcd         = nrnGcd;
156  r->cfIsUnit      = nrnIsUnit;
157  r->cfGetUnit     = nrnGetUnit;
158  r->cfExtGcd      = nrnExtGcd;
159  r->cfName        = ndName;
160  r->cfCoeffWrite  = nrnCoeffWrite;
161  r->nCoeffIsEqual = nrnCoeffsEqual;
162  r->cfKillChar    = ndKillChar;
163  r->cfQuot1       = nrnQuot1;
164#ifdef LDEBUG
165  r->cfDBTest      = nrnDBTest;
166#endif
167  return FALSE;
168}
169
170/*
171 * create a number from int
172 */
173number nrnInit(long i, const coeffs r)
174{
175  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
176  mpz_init_set_si(erg, i);
177  mpz_mod(erg, erg, r->modNumber);
178  return (number) erg;
179}
180
181void nrnDelete(number *a, const coeffs)
182{
183  if (*a == NULL) return;
184  mpz_clear((int_number) *a);
185  omFreeBin((void *) *a, gmp_nrz_bin);
186  *a = NULL;
187}
188
189number nrnCopy(number a, const coeffs)
190{
191  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
192  mpz_init_set(erg, (int_number) a);
193  return (number) erg;
194}
195
196int nrnSize(number a, const coeffs)
197{
198  if (a == NULL) return 0;
199  return sizeof(mpz_t);
200}
201
202/*
203 * convert a number to int
204 */
205int nrnInt(number &n, const coeffs)
206{
207  return (int)mpz_get_si((int_number) n);
208}
209
210/*
211 * Multiply two numbers
212 */
213number nrnMult(number a, number b, const coeffs r)
214{
215  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
216  mpz_init(erg);
217  mpz_mul(erg, (int_number)a, (int_number) b);
218  mpz_mod(erg, erg, r->modNumber);
219  return (number) erg;
220}
221
222void nrnPower(number a, int i, number * result, const coeffs r)
223{
224  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
225  mpz_init(erg);
226  mpz_powm_ui(erg, (int_number)a, i, r->modNumber);
227  *result = (number) erg;
228}
229
230number nrnAdd(number a, number b, const coeffs r)
231{
232  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
233  mpz_init(erg);
234  mpz_add(erg, (int_number)a, (int_number) b);
235  mpz_mod(erg, erg, r->modNumber);
236  return (number) erg;
237}
238
239number nrnSub(number a, number b, const coeffs r)
240{
241  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
242  mpz_init(erg);
243  mpz_sub(erg, (int_number)a, (int_number) b);
244  mpz_mod(erg, erg, r->modNumber);
245  return (number) erg;
246}
247
248number nrnNeg(number c, const coeffs r)
249{
250  if( !nrnIsZero(c, r) )
251    // Attention: This method operates in-place.
252    mpz_sub((int_number)c, r->modNumber, (int_number)c);
253  return c;
254}
255
256number nrnInvers(number c, const coeffs r)
257{
258  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
259  mpz_init(erg);
260  mpz_invert(erg, (int_number)c, r->modNumber);
261  return (number) erg;
262}
263
264/*
265 * Give the smallest k, such that a * x = k = b * y has a solution
266 * TODO: lcm(gcd,gcd) better than gcd(lcm) ?
267 */
268number nrnLcm(number a, number b, const coeffs r)
269{
270  number erg = nrnGcd(NULL, a, r);
271  number tmp = nrnGcd(NULL, b, r);
272  mpz_lcm((int_number)erg, (int_number)erg, (int_number)tmp);
273  nrnDelete(&tmp, NULL);
274  return (number)erg;
275}
276
277/*
278 * Give the largest k, such that a = x * k, b = y * k has
279 * a solution.
280 */
281number nrnGcd(number a, number b, const coeffs r)
282{
283  if ((a == NULL) && (b == NULL)) return nrnInit(0,r);
284  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
285  mpz_init_set(erg, r->modNumber);
286  if (a != NULL) mpz_gcd(erg, erg, (int_number)a);
287  if (b != NULL) mpz_gcd(erg, erg, (int_number)b);
288  return (number)erg;
289}
290
291/* Not needed any more, but may have room for improvement
292   number nrnGcd3(number a,number b, number c,ring r)
293{
294  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
295  mpz_init(erg);
296  if (a == NULL) a = (number)r->modNumber;
297  if (b == NULL) b = (number)r->modNumber;
298  if (c == NULL) c = (number)r->modNumber;
299  mpz_gcd(erg, (int_number)a, (int_number)b);
300  mpz_gcd(erg, erg, (int_number)c);
301  mpz_gcd(erg, erg, r->modNumber);
302  return (number)erg;
303}
304*/
305
306/*
307 * Give the largest k, such that a = x * k, b = y * k has
308 * a solution and r, s, s.t. k = s*a + t*b
309 */
310number nrnExtGcd(number a, number b, number *s, number *t, const coeffs r)
311{
312  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
313  int_number bs  = (int_number)omAllocBin(gmp_nrz_bin);
314  int_number bt  = (int_number)omAllocBin(gmp_nrz_bin);
315  mpz_init(erg);
316  mpz_init(bs);
317  mpz_init(bt);
318  mpz_gcdext(erg, bs, bt, (int_number)a, (int_number)b);
319  mpz_mod(bs, bs, r->modNumber);
320  mpz_mod(bt, bt, r->modNumber);
321  *s = (number)bs;
322  *t = (number)bt;
323  return (number)erg;
324}
325
326BOOLEAN nrnIsZero(number a, const coeffs)
327{
328#ifdef LDEBUG
329  if (a == NULL) return FALSE;
330#endif
331  return 0 == mpz_cmpabs_ui((int_number)a, 0);
332}
333
334BOOLEAN nrnIsOne(number a, const coeffs)
335{
336#ifdef LDEBUG
337  if (a == NULL) return FALSE;
338#endif
339  return 0 == mpz_cmp_si((int_number)a, 1);
340}
341
342BOOLEAN nrnIsMOne(number a, const coeffs r)
343{
344#ifdef LDEBUG
345  if (a == NULL) return FALSE;
346#endif
347  mpz_t t; mpz_init_set(t, (int_number)a);
348  mpz_add_ui(t, t, 1);
349  bool erg = (0 == mpz_cmp(t, r->modNumber));
350  mpz_clear(t);
351  return erg;
352}
353
354BOOLEAN nrnEqual(number a, number b, const coeffs)
355{
356  return 0 == mpz_cmp((int_number)a, (int_number)b);
357}
358
359BOOLEAN nrnGreater(number a, number b, const coeffs)
360{
361  return 0 < mpz_cmp((int_number)a, (int_number)b);
362}
363
364BOOLEAN nrnGreaterZero(number k, const coeffs)
365{
366  return 0 < mpz_cmp_si((int_number)k, 0);
367}
368
369BOOLEAN nrnIsUnit(number a, const coeffs r)
370{
371  number tmp = nrnGcd(a, (number)r->modNumber, r);
372  bool res = nrnIsOne(tmp, r);
373  nrnDelete(&tmp, NULL);
374  return res;
375}
376
377number nrnGetUnit(number k, const coeffs r)
378{
379  if (mpz_divisible_p(r->modNumber, (int_number)k)) return nrnInit(1,r);
380
381  int_number unit = (int_number)nrnGcd(k, 0, r);
382  mpz_tdiv_q(unit, (int_number)k, unit);
383  int_number gcd = (int_number)nrnGcd((number)unit, 0, r);
384  if (!nrnIsOne((number)gcd,r))
385  {
386    int_number ctmp;
387    // tmp := unit^2
388    int_number tmp = (int_number) nrnMult((number) unit,(number) unit,r);
389    // gcd_new := gcd(tmp, 0)
390    int_number gcd_new = (int_number) nrnGcd((number) tmp, 0, r);
391    while (!nrnEqual((number) gcd_new,(number) gcd,r))
392    {
393      // gcd := gcd_new
394      ctmp = gcd;
395      gcd = gcd_new;
396      gcd_new = ctmp;
397      // tmp := tmp * unit
398      mpz_mul(tmp, tmp, unit);
399      mpz_mod(tmp, tmp, r->modNumber);
400      // gcd_new := gcd(tmp, 0)
401      mpz_gcd(gcd_new, tmp, r->modNumber);
402    }
403    // unit := unit + modNumber / gcd_new
404    mpz_tdiv_q(tmp, r->modNumber, gcd_new);
405    mpz_add(unit, unit, tmp);
406    mpz_mod(unit, unit, r->modNumber);
407    nrnDelete((number*) &gcd_new, NULL);
408    nrnDelete((number*) &tmp, NULL);
409  }
410  nrnDelete((number*) &gcd, NULL);
411  return (number)unit;
412}
413
414BOOLEAN nrnDivBy(number a, number b, const coeffs r)
415{
416  if (a == NULL)
417    return mpz_divisible_p(r->modNumber, (int_number)b);
418  else
419  { /* b divides a iff b/gcd(a, b) is a unit in the given ring: */
420    number n = nrnGcd(a, b, r);
421    mpz_tdiv_q((int_number)n, (int_number)b, (int_number)n);
422    bool result = nrnIsUnit(n, r);
423    nrnDelete(&n, NULL);
424    return result;
425  }
426}
427
428int nrnDivComp(number a, number b, const coeffs r)
429{
430  if (nrnEqual(a, b,r)) return 2;
431  if (mpz_divisible_p((int_number) a, (int_number) b)) return -1;
432  if (mpz_divisible_p((int_number) b, (int_number) a)) return 1;
433  return 0;
434}
435
436number nrnDiv(number a, number b, const coeffs r)
437{
438  if (a == NULL) a = (number)r->modNumber;
439  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
440  mpz_init(erg);
441  if (mpz_divisible_p((int_number)a, (int_number)b))
442  {
443    mpz_divexact(erg, (int_number)a, (int_number)b);
444    return (number)erg;
445  }
446  else
447  {
448    int_number gcd = (int_number)nrnGcd(a, b, r);
449    mpz_divexact(erg, (int_number)b, gcd);
450    if (!nrnIsUnit((number)erg, r))
451    {
452      WerrorS("Division not possible, even by cancelling zero divisors.");
453      WerrorS("Result is integer division without remainder.");
454      mpz_tdiv_q(erg, (int_number) a, (int_number) b);
455      nrnDelete((number*) &gcd, NULL);
456      return (number)erg;
457    }
458    // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
459    int_number tmp = (int_number)nrnInvers((number) erg,r);
460    mpz_divexact(erg, (int_number)a, gcd);
461    mpz_mul(erg, erg, tmp);
462    nrnDelete((number*) &gcd, NULL);
463    nrnDelete((number*) &tmp, NULL);
464    mpz_mod(erg, erg, r->modNumber);
465    return (number)erg;
466  }
467}
468
469number nrnMod(number a, number b, const coeffs r)
470{
471  /*
472    We need to return the number rr which is uniquely determined by the
473    following two properties:
474      (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
475      (2) There exists some k in the integers Z such that a = k * b + rr.
476    Consider g := gcd(n, |b|). Note that then |b|/g is a unit in Z/n.
477    Now, there are three cases:
478      (a) g = 1
479          Then |b| is a unit in Z/n, i.e. |b| (and also b) divides a.
480          Thus rr = 0.
481      (b) g <> 1 and g divides a
482          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
483      (c) g <> 1 and g does not divide a
484          Then denote the division with remainder of a by g as this:
485          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
486          fulfills (1) and (2), i.e. rr := t is the correct result. Hence
487          in this third case, rr is the remainder of division of a by g in Z.
488     Remark: according to mpz_mod: a,b are always non-negative
489  */
490  int_number g = (int_number)omAllocBin(gmp_nrz_bin);
491  int_number rr = (int_number)omAllocBin(gmp_nrz_bin);
492  mpz_init(g);
493  mpz_init_set_si(rr, 0);
494  mpz_gcd(g, (int_number)r->modNumber, (int_number)b); // g is now as above
495  if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(rr, (int_number)a, g); // the case g <> 1
496  mpz_clear(g);
497  omFreeBin(g, gmp_nrz_bin);
498  return (number)rr;
499}
500
501number nrnIntDiv(number a, number b, const coeffs r)
502{
503  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
504  mpz_init(erg);
505  if (a == NULL) a = (number)r->modNumber;
506  mpz_tdiv_q(erg, (int_number)a, (int_number)b);
507  return (number)erg;
508}
509
510static number nrnAnn(number b, const coeffs r)
511{
512  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
513  mpz_init(erg);
514  mpz_tdiv_q(erg, (int_number)r->modNumber, (int_number)b);
515  return (number)erg;
516}
517
518/*
519 * Helper function for computing the module
520 */
521
522int_number nrnMapCoef = NULL;
523
524number nrnMapModN(number from, const coeffs /*src*/, const coeffs dst)
525{
526  return nrnMult(from, (number) nrnMapCoef, dst);
527}
528
529number nrnMap2toM(number from, const coeffs /*src*/, const coeffs dst)
530{
531  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
532  mpz_init(erg);
533  mpz_mul_ui(erg, nrnMapCoef, (NATNUMBER)from);
534  mpz_mod(erg, erg, dst->modNumber);
535  return (number)erg;
536}
537
538number nrnMapZp(number from, const coeffs /*src*/, const coeffs dst)
539{
540  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
541  mpz_init(erg);
542  // TODO: use npInt(...)
543  mpz_mul_si(erg, nrnMapCoef, (NATNUMBER)from);
544  mpz_mod(erg, erg, dst->modNumber);
545  return (number)erg;
546}
547
548number nrnMapGMP(number from, const coeffs /*src*/, const coeffs dst)
549{
550  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
551  mpz_init(erg);
552  mpz_mod(erg, (int_number)from, dst->modNumber);
553  return (number)erg;
554}
555
556number nrnMapZ(number from, const coeffs src, const coeffs dst)
557{
558  if (SR_HDL(from) & SR_INT)
559  {
560    long f_i=SR_TO_INT(from);
561    return nrnInit(f_i,dst);
562  }
563  return nrnMapGMP(from,src,dst);
564}
565
566number nrnMapQ(number from, const coeffs src, const coeffs dst)
567{
568  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
569  mpz_init(erg);
570  nlGMP(from, (number)erg, src);
571  mpz_mod(erg, erg, dst->modNumber);
572  return (number)erg;
573}
574
575nMapFunc nrnSetMap(const coeffs src, const coeffs dst)
576{
577  /* dst = nrn */
578  if ((src->rep==n_rep_gmp) && nCoeff_is_Ring_Z(src))
579  {
580    return nrnMapGMP;
581  }
582  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
583  {
584    return nrnMapZ;
585  }
586  if ((src->rep==n_rep_gap_rat) && nCoeff_is_Q(src))
587  {
588    return nrnMapQ;
589  }
590  // Some type of Z/n ring / field
591  if (nCoeff_is_Ring_ModN(src) || nCoeff_is_Ring_PtoM(src) ||
592      nCoeff_is_Ring_2toM(src) || nCoeff_is_Zp(src))
593  {
594    if (   (!nCoeff_is_Zp(src))
595        && (mpz_cmp(src->modBase, dst->modBase) == 0)
596        && (src->modExponent == dst->modExponent)) return nrnMapGMP;
597    else
598    {
599      int_number nrnMapModul = (int_number) omAllocBin(gmp_nrz_bin);
600      // Computing the n of Z/n
601      if (nCoeff_is_Zp(src))
602      {
603        mpz_init_set_si(nrnMapModul, src->ch);
604      }
605      else
606      {
607        mpz_init(nrnMapModul);
608        mpz_set(nrnMapModul, src->modNumber);
609      }
610      // nrnMapCoef = 1 in dst       if dst is a subring of src
611      // nrnMapCoef = 0 in dst / src if src is a subring of dst
612      if (nrnMapCoef == NULL)
613      {
614        nrnMapCoef = (int_number) omAllocBin(gmp_nrz_bin);
615        mpz_init(nrnMapCoef);
616      }
617      if (mpz_divisible_p(nrnMapModul, dst->modNumber))
618      {
619        mpz_set_si(nrnMapCoef, 1);
620      }
621      else
622      if (nrnDivBy(NULL, (number) nrnMapModul,dst))
623      {
624        mpz_divexact(nrnMapCoef, dst->modNumber, nrnMapModul);
625        int_number tmp = dst->modNumber;
626        dst->modNumber = nrnMapModul;
627        if (!nrnIsUnit((number) nrnMapCoef,dst))
628        {
629          dst->modNumber = tmp;
630          nrnDelete((number*) &nrnMapModul, dst);
631          return NULL;
632        }
633        int_number inv = (int_number) nrnInvers((number) nrnMapCoef,dst);
634        dst->modNumber = tmp;
635        mpz_mul(nrnMapCoef, nrnMapCoef, inv);
636        mpz_mod(nrnMapCoef, nrnMapCoef, dst->modNumber);
637        nrnDelete((number*) &inv, dst);
638      }
639      else
640      {
641        nrnDelete((number*) &nrnMapModul, dst);
642        return NULL;
643      }
644      nrnDelete((number*) &nrnMapModul, dst);
645      if (nCoeff_is_Ring_2toM(src))
646        return nrnMap2toM;
647      else if (nCoeff_is_Zp(src))
648        return nrnMapZp;
649      else
650        return nrnMapModN;
651    }
652  }
653  return NULL;      // default
654}
655
656/*
657 * set the exponent (allocate and init tables) (TODO)
658 */
659
660void nrnSetExp(unsigned long m, coeffs r)
661{
662  /* clean up former stuff */
663  if (r->modNumber != NULL) mpz_clear(r->modNumber);
664
665  r->modExponent= m;
666  r->modNumber = (int_number)omAllocBin(gmp_nrz_bin);
667  mpz_init_set (r->modNumber, r->modBase);
668  mpz_pow_ui (r->modNumber, r->modNumber, m);
669}
670
671/* We expect this ring to be Z/n^m for some m > 0 and for some n > 2 which is not a prime. */
672void nrnInitExp(unsigned long m, coeffs r)
673{
674  nrnSetExp(m, r);
675  assume (r->modNumber != NULL);
676  if (mpz_cmp_ui(r->modNumber,2) <= 0)
677    WarnS("nrnInitExp failed (m in Z/m too small)");
678}
679
680#ifdef LDEBUG
681BOOLEAN nrnDBTest (number a, const char *, const int, const coeffs r)
682{
683  if (a==NULL) return TRUE;
684  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, r->modNumber) > 0) )
685  {
686    return FALSE;
687  }
688  return TRUE;
689}
690#endif
691
692/*2
693* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
694*/
695static const char * nlCPEatLongC(char *s, mpz_ptr i)
696{
697  const char * start=s;
698  if (!(*s >= '0' && *s <= '9'))
699  {
700    mpz_init_set_si(i, 1);
701    return s;
702  }
703  mpz_init(i);
704  while (*s >= '0' && *s <= '9') s++;
705  if (*s=='\0')
706  {
707    mpz_set_str(i,start,10);
708  }
709  else
710  {
711    char c=*s;
712    *s='\0';
713    mpz_set_str(i,start,10);
714    *s=c;
715  }
716  return s;
717}
718
719const char * nrnRead (const char *s, number *a, const coeffs r)
720{
721  int_number z = (int_number) omAllocBin(gmp_nrz_bin);
722  {
723    s = nlCPEatLongC((char *)s, z);
724  }
725  mpz_mod(z, z, r->modNumber);
726  *a = (number) z;
727  return s;
728}
729#endif
730/* #ifdef HAVE_RINGS */
Note: See TracBrowser for help on using the repository browser.