source: git/libpolys/coeffs/rmodulon.cc @ 0acf3e

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