source: git/libpolys/coeffs/rmodulon.cc @ 45cc512

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