Changeset bd8dc07 in git


Ignore:
Timestamp:
Nov 4, 2016, 12:14:16 PM (7 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
23ca7b7a1c132bcc62bab0c48743fe235bbe1764
Parents:
446a69d5ed9308c77490b1c0193e81aba162e1fb
Message:
fix: nr2mCoeffName, more static, no forward decl. in rmodulo2m.cc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • libpolys/coeffs/rmodulo2m.cc

    r446a69 rbd8dc07  
    2525#ifdef HAVE_RINGS
    2626
    27 number  nr2mCopy        (number a, const coeffs r);
    28 BOOLEAN nr2mGreaterZero (number k, const coeffs r);
    29 number  nr2mMult        (number a, number b, const coeffs r);
    30 number  nr2mInit        (long i, const coeffs r);
    31 long    nr2mInt         (number &n, const coeffs r);
    32 number  nr2mAdd         (number a, number b, const coeffs r);
    33 number  nr2mSub         (number a, number b, const coeffs r);
    34 void    nr2mPower       (number a, int i, number * result, const coeffs r);
    35 BOOLEAN nr2mIsZero      (number a, const coeffs r);
    36 BOOLEAN nr2mIsOne       (number a, const coeffs r);
    37 BOOLEAN nr2mIsMOne      (number a, const coeffs r);
    38 BOOLEAN nr2mIsUnit      (number a, const coeffs r);
    39 number  nr2mGetUnit     (number a, const coeffs r);
    40 number  nr2mDiv         (number a, number b, const coeffs r);
    41 number  nr2mIntDiv      (number a,number b, const coeffs r);
    42 number  nr2mMod         (number a,number b, const coeffs r);
    43 number  nr2mNeg         (number c, const coeffs r);
    44 number  nr2mInvers      (number c, const coeffs r);
    45 BOOLEAN nr2mGreater     (number a, number b, const coeffs r);
    46 BOOLEAN nr2mDivBy       (number a, number b, const coeffs r);
    47 int     nr2mDivComp     (number a, number b, const coeffs r);
    48 BOOLEAN nr2mEqual       (number a, number b, const coeffs r);
    49 number  nr2mLcm         (number a, number b, const coeffs r);
    50 number  nr2mGcd         (number a, number b, const coeffs r);
    51 number  nr2mExtGcd      (number a, number b, number *s, number *t, const coeffs r);
    52 nMapFunc nr2mSetMap     (const coeffs src, const coeffs dst);
    53 void    nr2mWrite       (number a, const coeffs r);
    54 const char *  nr2mRead  (const char *s, number *a, const coeffs r);
    55 char *  nr2mName        (number n, const coeffs r);
    56 void    nr2mCoeffWrite  (const coeffs r, BOOLEAN details);
    57 coeffs  nr2mQuot1(number c, const coeffs r);
    58 #ifdef LDEBUG
    59 BOOLEAN nr2mDBTest      (number a, const char *f, const int l, const coeffs r);
    60 #endif
    61 void    nr2mSetExp(int c, const coeffs r);
    62 void    nr2mInitExp(int c, const coeffs r);
    63 
    64 number nr2mMapQ(number from, const coeffs src, const coeffs dst);
    65 
    6627static inline number nr2mMultM(number a, number b, const coeffs r)
    6728{
    … …  
    8849extern omBin gmp_nrz_bin; /* init in rintegers*/
    8950
    90 void    nr2mCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
     51static char* nr2mCoeffName(const coeffs cf)
     52{
     53  static char n2mCoeffName_buf[22];
     54#ifdef SINGULAR_4_1
     55  snprintf(n2mCoeffName_buf,21,"ZZ/(2^%lu)",cf->modExponent);
     56#else
     57  snprintf(n2mCoeffName_buf,21,"integer,2,%lu",cf->modExponent);
     58#endif
     59  return n2mCoeffName_buf;
     60}
     61
     62static void    nr2mCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
    9163{
    9264  PrintS("//   coeff. ring is : ");
    … …  
    9466}
    9567
    96 BOOLEAN nr2mCoeffIsEqual(const coeffs r, n_coeffType n, void * p)
     68static BOOLEAN nr2mCoeffIsEqual(const coeffs r, n_coeffType n, void * p)
    9769{
    9870  if (n==n_Z2m)
    … …  
    11789}
    11890
    119 coeffs nr2mQuot1(number c, const coeffs r)
     91static coeffs nr2mQuot1(number c, const coeffs r)
    12092{
    12193    coeffs rr;
    … …  
    158130}
    159131
     132/* TRUE iff 0 < k <= 2^m / 2 */
     133static BOOLEAN nr2mGreaterZero(number k, const coeffs r)
     134{
     135  if ((unsigned long)k == 0) return FALSE;
     136  if ((unsigned long)k > ((r->mod2mMask >> 1) + 1)) return FALSE;
     137  return TRUE;
     138}
     139
     140/*
     141 * Multiply two numbers
     142 */
     143static number nr2mMult(number a, number b, const coeffs r)
     144{
     145  if (((unsigned long)a == 0) || ((unsigned long)b == 0))
     146    return (number)0;
     147  else
     148    return nr2mMultM(a, b, r);
     149}
     150
    160151static number nr2mAnn(number b, const coeffs r);
     152/*
     153 * Give the smallest k, such that a * x = k = b * y has a solution
     154 */
     155static number nr2mLcm(number a, number b, const coeffs)
     156{
     157  unsigned long res = 0;
     158  if ((unsigned long)a == 0) a = (number) 1;
     159  if ((unsigned long)b == 0) b = (number) 1;
     160  while ((unsigned long)a % 2 == 0)
     161  {
     162    a = (number)((unsigned long)a / 2);
     163    if ((unsigned long)b % 2 == 0) b = (number)((unsigned long)b / 2);
     164    res++;
     165  }
     166  while ((unsigned long)b % 2 == 0)
     167  {
     168    b = (number)((unsigned long)b / 2);
     169    res++;
     170  }
     171  return (number)(1L << res);  // (2**res)
     172}
     173
     174/*
     175 * Give the largest k, such that a = x * k, b = y * k has
     176 * a solution.
     177 */
     178static number nr2mGcd(number a, number b, const coeffs)
     179{
     180  unsigned long res = 0;
     181  if ((unsigned long)a == 0 && (unsigned long)b == 0) return (number)1;
     182  while ((unsigned long)a % 2 == 0 && (unsigned long)b % 2 == 0)
     183  {
     184    a = (number)((unsigned long)a / 2);
     185    b = (number)((unsigned long)b / 2);
     186    res++;
     187  }
     188//  if ((unsigned long)b % 2 == 0)
     189//  {
     190//    return (number)((1L << res)); // * (unsigned long) a);  // (2**res)*a    a is a unit
     191//  }
     192//  else
     193//  {
     194    return (number)((1L << res)); // * (unsigned long) b);  // (2**res)*b    b is a unit
     195//  }
     196}
     197
     198/* assumes that 'a' is odd, i.e., a unit in Z/2^m, and computes
     199   the extended gcd of 'a' and 2^m, in order to find some 's'
     200   and 't' such that a * s + 2^m * t = gcd(a, 2^m) = 1;
     201   this code will always find a positive 's' */
     202static void specialXGCD(unsigned long& s, unsigned long a, const coeffs r)
     203{
     204  mpz_ptr u = (mpz_ptr)omAlloc(sizeof(mpz_t));
     205  mpz_init_set_ui(u, a);
     206  mpz_ptr u0 = (mpz_ptr)omAlloc(sizeof(mpz_t));
     207  mpz_init(u0);
     208  mpz_ptr u1 = (mpz_ptr)omAlloc(sizeof(mpz_t));
     209  mpz_init_set_ui(u1, 1);
     210  mpz_ptr u2 = (mpz_ptr)omAlloc(sizeof(mpz_t));
     211  mpz_init(u2);
     212  mpz_ptr v = (mpz_ptr)omAlloc(sizeof(mpz_t));
     213  mpz_init_set_ui(v, r->mod2mMask);
     214  mpz_add_ui(v, v, 1); /* now: v = 2^m */
     215  mpz_ptr v0 = (mpz_ptr)omAlloc(sizeof(mpz_t));
     216  mpz_init(v0);
     217  mpz_ptr v1 = (mpz_ptr)omAlloc(sizeof(mpz_t));
     218  mpz_init(v1);
     219  mpz_ptr v2 = (mpz_ptr)omAlloc(sizeof(mpz_t));
     220  mpz_init_set_ui(v2, 1);
     221  mpz_ptr q = (mpz_ptr)omAlloc(sizeof(mpz_t));
     222  mpz_init(q);
     223  mpz_ptr rr = (mpz_ptr)omAlloc(sizeof(mpz_t));
     224  mpz_init(rr);
     225
     226  while (mpz_cmp_ui(v, 0) != 0) /* i.e., while v != 0 */
     227  {
     228    mpz_div(q, u, v);
     229    mpz_mod(rr, u, v);
     230    mpz_set(u, v);
     231    mpz_set(v, rr);
     232    mpz_set(u0, u2);
     233    mpz_set(v0, v2);
     234    mpz_mul(u2, u2, q); mpz_sub(u2, u1, u2); /* u2 = u1 - q * u2 */
     235    mpz_mul(v2, v2, q); mpz_sub(v2, v1, v2); /* v2 = v1 - q * v2 */
     236    mpz_set(u1, u0);
     237    mpz_set(v1, v0);
     238  }
     239
     240  while (mpz_cmp_ui(u1, 0) < 0) /* i.e., while u1 < 0 */
     241  {
     242    /* we add 2^m = (2^m - 1) + 1 to u1: */
     243    mpz_add_ui(u1, u1, r->mod2mMask);
     244    mpz_add_ui(u1, u1, 1);
     245  }
     246  s = mpz_get_ui(u1); /* now: 0 <= s <= 2^m - 1 */
     247
     248  mpz_clear(u);  omFree((ADDRESS)u);
     249  mpz_clear(u0); omFree((ADDRESS)u0);
     250  mpz_clear(u1); omFree((ADDRESS)u1);
     251  mpz_clear(u2); omFree((ADDRESS)u2);
     252  mpz_clear(v);  omFree((ADDRESS)v);
     253  mpz_clear(v0); omFree((ADDRESS)v0);
     254  mpz_clear(v1); omFree((ADDRESS)v1);
     255  mpz_clear(v2); omFree((ADDRESS)v2);
     256  mpz_clear(q); omFree((ADDRESS)q);
     257  mpz_clear(rr); omFree((ADDRESS)rr);
     258}
     259
     260static unsigned long InvMod(unsigned long a, const coeffs r)
     261{
     262  assume((unsigned long)a % 2 != 0);
     263  unsigned long s;
     264  specialXGCD(s, a, r);
     265  return s;
     266}
     267
     268static inline number nr2mInversM(number c, const coeffs r)
     269{
     270  assume((unsigned long)c % 2 != 0);
     271  // Table !!!
     272  unsigned long inv;
     273  inv = InvMod((unsigned long)c,r);
     274  return (number)inv;
     275}
     276
     277static number nr2mInvers(number c, const coeffs r)
     278{
     279  if ((unsigned long)c % 2 == 0)
     280  {
     281    WerrorS("division by zero divisor");
     282    return (number)0;
     283  }
     284  return nr2mInversM(c, r);
     285}
     286
     287/*
     288 * Give the largest k, such that a = x * k, b = y * k has
     289 * a solution.
     290 */
     291static number nr2mExtGcd(number a, number b, number *s, number *t, const coeffs r)
     292{
     293  unsigned long res = 0;
     294  if ((unsigned long)a == 0 && (unsigned long)b == 0) return (number)1;
     295  while ((unsigned long)a % 2 == 0 && (unsigned long)b % 2 == 0)
     296  {
     297    a = (number)((unsigned long)a / 2);
     298    b = (number)((unsigned long)b / 2);
     299    res++;
     300  }
     301  if ((unsigned long)b % 2 == 0)
     302  {
     303    *t = NULL;
     304    *s = nr2mInvers(a,r);
     305    return (number)((1L << res)); // * (unsigned long) a);  // (2**res)*a    a is a unit
     306  }
     307  else
     308  {
     309    *s = NULL;
     310    *t = nr2mInvers(b,r);
     311    return (number)((1L << res)); // * (unsigned long) b);  // (2**res)*b    b is a unit
     312  }
     313}
     314
     315static void nr2mPower(number a, int i, number * result, const coeffs r)
     316{
     317  if (i == 0)
     318  {
     319    *(unsigned long *)result = 1;
     320  }
     321  else if (i == 1)
     322  {
     323    *result = a;
     324  }
     325  else
     326  {
     327    nr2mPower(a, i-1, result, r);
     328    *result = nr2mMultM(a, *result, r);
     329  }
     330}
     331
     332/*
     333 * create a number from int
     334 */
     335static number nr2mInit(long i, const coeffs r)
     336{
     337  if (i == 0) return (number)(unsigned long)i;
     338
     339  long ii = i;
     340  unsigned long j = (unsigned long)1;
     341  if (ii < 0) { j = r->mod2mMask; ii = -ii; }
     342  unsigned long k = (unsigned long)ii;
     343  k = k & r->mod2mMask;
     344  /* now we have: i = j * k mod 2^m */
     345  return (number)nr2mMult((number)j, (number)k, r);
     346}
     347
     348/*
     349 * convert a number to an int in ]-k/2 .. k/2],
     350 * where k = 2^m; i.e., an int in ]-2^(m-1) .. 2^(m-1)];
     351 */
     352static long nr2mInt(number &n, const coeffs r)
     353{
     354  unsigned long nn = (unsigned long)(unsigned long)n & r->mod2mMask;
     355  unsigned long l = r->mod2mMask >> 1; l++; /* now: l = 2^(m-1) */
     356  if ((unsigned long)nn > l)
     357    return (long)((unsigned long)nn - r->mod2mMask - 1);
     358  else
     359    return (long)((unsigned long)nn);
     360}
     361
     362static number nr2mAdd(number a, number b, const coeffs r)
     363{
     364  return nr2mAddM(a, b, r);
     365}
     366
     367static number nr2mSub(number a, number b, const coeffs r)
     368{
     369  return nr2mSubM(a, b, r);
     370}
     371
     372static BOOLEAN nr2mIsUnit(number a, const coeffs)
     373{
     374  return ((unsigned long)a % 2 == 1);
     375}
     376
     377static number nr2mGetUnit(number k, const coeffs)
     378{
     379  if (k == NULL) return (number)1;
     380  unsigned long erg = (unsigned long)k;
     381  while (erg % 2 == 0) erg = erg / 2;
     382  return (number)erg;
     383}
     384
     385static BOOLEAN nr2mIsZero(number a, const coeffs)
     386{
     387  return 0 == (unsigned long)a;
     388}
     389
     390static BOOLEAN nr2mIsOne(number a, const coeffs)
     391{
     392  return 1 == (unsigned long)a;
     393}
     394
     395static BOOLEAN nr2mIsMOne(number a, const coeffs r)
     396{
     397  return ((r->mod2mMask  == (unsigned long)a) &&(1L!=(long)a))/*for char 2^1*/;
     398}
     399
     400static BOOLEAN nr2mEqual(number a, number b, const coeffs)
     401{
     402  return (a == b);
     403}
     404
     405static number nr2mDiv(number a, number b, const coeffs r)
     406{
     407  if ((unsigned long)a == 0) return (number)0;
     408  else if ((unsigned long)b % 2 == 0)
     409  {
     410    if ((unsigned long)b != 0)
     411    {
     412      while (((unsigned long)b % 2 == 0) && ((unsigned long)a % 2 == 0))
     413      {
     414        a = (number)((unsigned long)a / 2);
     415        b = (number)((unsigned long)b / 2);
     416      }
     417    }
     418    if ((unsigned long)b % 2 == 0)
     419    {
     420      WerrorS("Division not possible, even by cancelling zero divisors.");
     421      WerrorS("Result is integer division without remainder.");
     422      return (number) ((unsigned long) a / (unsigned long) b);
     423    }
     424  }
     425  return (number)nr2mMult(a, nr2mInversM(b,r),r);
     426}
     427
     428/* Is 'a' divisible by 'b'? There are two cases:
     429   1) a = 0 mod 2^m; then TRUE iff b = 0 or b is a power of 2
     430   2) a, b <> 0; then TRUE iff b/gcd(a, b) is a unit mod 2^m */
     431static BOOLEAN nr2mDivBy (number a, number b, const coeffs r)
     432{
     433  if (a == NULL)
     434  {
     435    unsigned long c = r->mod2mMask + 1;
     436    if (c != 0) /* i.e., if no overflow */
     437      return (c % (unsigned long)b) == 0;
     438    else
     439    {
     440      /* overflow: we need to check whether b
     441         is zero or a power of 2: */
     442      c = (unsigned long)b;
     443      while (c != 0)
     444      {
     445        if ((c % 2) != 0) return FALSE;
     446        c = c >> 1;
     447      }
     448      return TRUE;
     449    }
     450  }
     451  else
     452  {
     453    number n = nr2mGcd(a, b, r);
     454    n = nr2mDiv(b, n, r);
     455    return nr2mIsUnit(n, r);
     456  }
     457}
     458
     459static BOOLEAN nr2mGreater(number a, number b, const coeffs r)
     460{
     461  return nr2mDivBy(a, b,r);
     462}
     463
     464static int nr2mDivComp(number as, number bs, const coeffs)
     465{
     466  unsigned long a = (unsigned long)as;
     467  unsigned long b = (unsigned long)bs;
     468  assume(a != 0 && b != 0);
     469  while (a % 2 == 0 && b % 2 == 0)
     470  {
     471    a = a / 2;
     472    b = b / 2;
     473  }
     474  if (a % 2 == 0)
     475  {
     476    return -1;
     477  }
     478  else
     479  {
     480    if (b % 2 == 1)
     481    {
     482      return 2;
     483    }
     484    else
     485    {
     486      return 1;
     487    }
     488  }
     489}
     490
     491static number nr2mMod(number a, number b, const coeffs r)
     492{
     493  /*
     494    We need to return the number rr which is uniquely determined by the
     495    following two properties:
     496      (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
     497      (2) There exists some k in the integers Z such that a = k * b + rr.
     498    Consider g := gcd(2^m, |b|). Note that then |b|/g is a unit in Z/2^m.
     499    Now, there are three cases:
     500      (a) g = 1
     501          Then |b| is a unit in Z/2^m, i.e. |b| (and also b) divides a.
     502          Thus rr = 0.
     503      (b) g <> 1 and g divides a
     504          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
     505      (c) g <> 1 and g does not divide a
     506          Let's denote the division with remainder of a by g as follows:
     507          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
     508          fulfills (1) and (2), i.e. rr := t is the correct result. Hence
     509          in this third case, rr is the remainder of division of a by g in Z.
     510    This algorithm is the same as for the case Z/n, except that we may
     511    compute the gcd of |b| and 2^m "by hand": We just extract the highest
     512    power of 2 (<= 2^m) that is contained in b.
     513  */
     514  assume((unsigned long) b != 0);
     515  unsigned long g = 1;
     516  unsigned long b_div = (unsigned long) b;
     517
     518  /*
     519   * b_div is unsigned, so that (b_div < 0) evaluates false at compile-time
     520   *
     521  if (b_div < 0) b_div = -b_div; // b_div now represents |b|, BUT b_div is unsigned!
     522  */
     523
     524  unsigned long rr = 0;
     525  while ((g < r->mod2mMask ) && (b_div > 0) && (b_div % 2 == 0))
     526  {
     527    b_div = b_div >> 1;
     528    g = g << 1;
     529  } // g is now the gcd of 2^m and |b|
     530
     531  if (g != 1) rr = (unsigned long)a % g;
     532  return (number)rr;
     533}
     534
     535static number nr2mIntDiv(number a, number b, const coeffs r)
     536{
     537  if ((unsigned long)a == 0)
     538  {
     539    if ((unsigned long)b == 0)
     540      return (number)1;
     541    if ((unsigned long)b == 1)
     542      return (number)0;
     543    unsigned long c = r->mod2mMask + 1;
     544    if (c != 0) /* i.e., if no overflow */
     545      return (number)(c / (unsigned long)b);
     546    else
     547    {
     548      /* overflow: c = 2^32 resp. 2^64, depending on platform */
     549      mpz_ptr cc = (mpz_ptr)omAlloc(sizeof(mpz_t));
     550      mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
     551      mpz_div_ui(cc, cc, (unsigned long)(unsigned long)b);
     552      unsigned long s = mpz_get_ui(cc);
     553      mpz_clear(cc); omFree((ADDRESS)cc);
     554      return (number)(unsigned long)s;
     555    }
     556  }
     557  else
     558  {
     559    if ((unsigned long)b == 0)
     560      return (number)0;
     561    return (number)((unsigned long) a / (unsigned long) b);
     562  }
     563}
     564
     565static number nr2mAnn(number b, const coeffs r)
     566{
     567  if ((unsigned long)b == 0)
     568    return NULL;
     569  if ((unsigned long)b == 1)
     570    return NULL;
     571  unsigned long c = r->mod2mMask + 1;
     572  if (c != 0) /* i.e., if no overflow */
     573    return (number)(c / (unsigned long)b);
     574  else
     575  {
     576    /* overflow: c = 2^32 resp. 2^64, depending on platform */
     577    mpz_ptr cc = (mpz_ptr)omAlloc(sizeof(mpz_t));
     578    mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
     579    mpz_div_ui(cc, cc, (unsigned long)(unsigned long)b);
     580    unsigned long s = mpz_get_ui(cc);
     581    mpz_clear(cc); omFree((ADDRESS)cc);
     582    return (number)(unsigned long)s;
     583  }
     584}
     585
     586static number nr2mNeg(number c, const coeffs r)
     587{
     588  if ((unsigned long)c == 0) return c;
     589  return nr2mNegM(c, r);
     590}
     591
     592static number nr2mMapMachineInt(number from, const coeffs /*src*/, const coeffs dst)
     593{
     594  unsigned long i = ((unsigned long)from) % dst->mod2mMask ;
     595  return (number)i;
     596}
     597
     598static number nr2mMapProject(number from, const coeffs /*src*/, const coeffs dst)
     599{
     600  unsigned long i = ((unsigned long)from) % (dst->mod2mMask + 1);
     601  return (number)i;
     602}
     603
     604number nr2mMapZp(number from, const coeffs /*src*/, const coeffs dst)
     605{
     606  unsigned long j = (unsigned long)1;
     607  long ii = (long)from;
     608  if (ii < 0) { j = dst->mod2mMask; ii = -ii; }
     609  unsigned long i = (unsigned long)ii;
     610  i = i & dst->mod2mMask;
     611  /* now we have: from = j * i mod 2^m */
     612  return (number)nr2mMult((number)i, (number)j, dst);
     613}
     614
     615static number nr2mMapGMP(number from, const coeffs /*src*/, const coeffs dst)
     616{
     617  mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
     618  mpz_init(erg);
     619  mpz_ptr k = (mpz_ptr)omAlloc(sizeof(mpz_t));
     620  mpz_init_set_ui(k, dst->mod2mMask);
     621
     622  mpz_and(erg, (mpz_ptr)from, k);
     623  number res = (number) mpz_get_ui(erg);
     624
     625  mpz_clear(erg); omFree((ADDRESS)erg);
     626  mpz_clear(k);   omFree((ADDRESS)k);
     627
     628  return (number)res;
     629}
     630
     631static number nr2mMapQ(number from, const coeffs src, const coeffs dst)
     632{
     633  mpz_ptr gmp = (mpz_ptr)omAllocBin(gmp_nrz_bin);
     634  mpz_init(gmp);
     635  nlGMP(from, (number)gmp, src); // FIXME? TODO? // extern void   nlGMP(number &i, number n, const coeffs r); // to be replaced with n_MPZ(erg, from, src); // ?
     636  number res=nr2mMapGMP((number)gmp,src,dst);
     637  mpz_clear(gmp); omFree((ADDRESS)gmp);
     638  return res;
     639}
     640
     641static number nr2mMapZ(number from, const coeffs src, const coeffs dst)
     642{
     643  if (SR_HDL(from) & SR_INT)
     644  {
     645    long f_i=SR_TO_INT(from);
     646    return nr2mInit(f_i,dst);
     647  }
     648  return nr2mMapGMP(from,src,dst);
     649}
     650
     651static nMapFunc nr2mSetMap(const coeffs src, const coeffs dst)
     652{
     653  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
     654     && (src->mod2mMask == dst->mod2mMask))
     655  {
     656    return ndCopyMap;
     657  }
     658  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
     659     && (src->mod2mMask < dst->mod2mMask))
     660  { /* i.e. map an integer mod 2^s into Z mod 2^t, where t < s */
     661    return nr2mMapMachineInt;
     662  }
     663  if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
     664     && (src->mod2mMask > dst->mod2mMask))
     665  { /* i.e. map an integer mod 2^s into Z mod 2^t, where t > s */
     666    // to be done
     667    return nr2mMapProject;
     668  }
     669  if ((src->rep==n_rep_gmp) && nCoeff_is_Ring_Z(src))
     670  {
     671    return nr2mMapGMP;
     672  }
     673  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
     674  {
     675    return nr2mMapZ;
     676  }
     677  if ((src->rep==n_rep_gap_rat) && (nCoeff_is_Q(src)||nCoeff_is_Ring_Z(src)))
     678  {
     679    return nr2mMapQ;
     680  }
     681  if ((src->rep==n_rep_int) && nCoeff_is_Zp(src) && (src->ch == 2))
     682  {
     683    return nr2mMapZp;
     684  }
     685  if ((src->rep==n_rep_gmp) &&
     686  (nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src)))
     687  {
     688    if (mpz_divisible_2exp_p(src->modNumber,dst->modExponent))
     689      return nr2mMapGMP;
     690  }
     691  return NULL;      // default
     692}
     693
     694/*
     695 * set the exponent
     696 */
     697
     698static void nr2mSetExp(int m, coeffs r)
     699{
     700  if (m > 1)
     701  {
     702    /* we want mod2mMask to be the bit pattern
     703       '111..1' consisting of m one's: */
     704    r->modExponent= m;
     705    r->mod2mMask = 1;
     706    for (int i = 1; i < m; i++) r->mod2mMask = (r->mod2mMask << 1) + 1;
     707  }
     708  else
     709  {
     710    r->modExponent= 2;
     711    /* code unexpectedly called with m = 1; we continue with m = 2: */
     712    r->mod2mMask = 3; /* i.e., '11' in binary representation */
     713  }
     714}
     715
     716static void nr2mInitExp(int m, coeffs r)
     717{
     718  nr2mSetExp(m, r);
     719  if (m < 2)
     720    WarnS("nr2mInitExp unexpectedly called with m = 1 (we continue with Z/2^2");
     721}
     722
     723#ifdef LDEBUG
     724static BOOLEAN nr2mDBTest (number a, const char *, const int, const coeffs r)
     725{
     726  //if ((unsigned long)a < 0) return FALSE; // is unsigned!
     727  if (((unsigned long)a & r->mod2mMask) != (unsigned long)a) return FALSE;
     728  return TRUE;
     729}
     730#endif
     731
     732static void nr2mWrite (number a, const coeffs r)
     733{
     734  long i = nr2mInt(a, r);
     735  StringAppend("%ld", i);
     736}
     737
     738static const char* nr2mEati(const char *s, int *i, const coeffs r)
     739{
     740
     741  if (((*s) >= '0') && ((*s) <= '9'))
     742  {
     743    (*i) = 0;
     744    do
     745    {
     746      (*i) *= 10;
     747      (*i) += *s++ - '0';
     748      if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & r->mod2mMask;
     749    }
     750    while (((*s) >= '0') && ((*s) <= '9'));
     751    (*i) = (*i) & r->mod2mMask;
     752  }
     753  else (*i) = 1;
     754  return s;
     755}
     756
     757static const char * nr2mRead (const char *s, number *a, const coeffs r)
     758{
     759  int z;
     760  int n=1;
     761
     762  s = nr2mEati(s, &z,r);
     763  if ((*s) == '/')
     764  {
     765    s++;
     766    s = nr2mEati(s, &n,r);
     767  }
     768  if (n == 1)
     769    *a = (number)(long)z;
     770  else
     771      *a = nr2mDiv((number)(long)z,(number)(long)n,r);
     772  return s;
     773}
     774
    161775/* for initializing function pointers */
    162776BOOLEAN nr2mInitChar (coeffs r, void* p)
    … …  
    213827  r->cfExtGcd      = nr2mExtGcd;
    214828  r->cfCoeffWrite  = nr2mCoeffWrite;
     829  r->cfCoeffName   = nr2mCoeffName;
    215830  r->cfQuot1       = nr2mQuot1;
    216831#ifdef LDEBUG
    … …  
    221836}
    222837
    223 /*
    224  * Multiply two numbers
    225  */
    226 number nr2mMult(number a, number b, const coeffs r)
    227 {
    228   if (((unsigned long)a == 0) || ((unsigned long)b == 0))
    229     return (number)0;
    230   else
    231     return nr2mMultM(a, b, r);
    232 }
    233 
    234 /*
    235  * Give the smallest k, such that a * x = k = b * y has a solution
    236  */
    237 number nr2mLcm(number a, number b, const coeffs)
    238 {
    239   unsigned long res = 0;
    240   if ((unsigned long)a == 0) a = (number) 1;
    241   if ((unsigned long)b == 0) b = (number) 1;
    242   while ((unsigned long)a % 2 == 0)
    243   {
    244     a = (number)((unsigned long)a / 2);
    245     if ((unsigned long)b % 2 == 0) b = (number)((unsigned long)b / 2);
    246     res++;
    247   }
    248   while ((unsigned long)b % 2 == 0)
    249   {
    250     b = (number)((unsigned long)b / 2);
    251     res++;
    252   }
    253   return (number)(1L << res);  // (2**res)
    254 }
    255 
    256 /*
    257  * Give the largest k, such that a = x * k, b = y * k has
    258  * a solution.
    259  */
    260 number nr2mGcd(number a, number b, const coeffs)
    261 {
    262   unsigned long res = 0;
    263   if ((unsigned long)a == 0 && (unsigned long)b == 0) return (number)1;
    264   while ((unsigned long)a % 2 == 0 && (unsigned long)b % 2 == 0)
    265   {
    266     a = (number)((unsigned long)a / 2);
    267     b = (number)((unsigned long)b / 2);
    268     res++;
    269   }
    270 //  if ((unsigned long)b % 2 == 0)
    271 //  {
    272 //    return (number)((1L << res)); // * (unsigned long) a);  // (2**res)*a    a is a unit
    273 //  }
    274 //  else
    275 //  {
    276     return (number)((1L << res)); // * (unsigned long) b);  // (2**res)*b    b is a unit
    277 //  }
    278 }
    279 
    280 /*
    281  * Give the largest k, such that a = x * k, b = y * k has
    282  * a solution.
    283  */
    284 number nr2mExtGcd(number a, number b, number *s, number *t, const coeffs r)
    285 {
    286   unsigned long res = 0;
    287   if ((unsigned long)a == 0 && (unsigned long)b == 0) return (number)1;
    288   while ((unsigned long)a % 2 == 0 && (unsigned long)b % 2 == 0)
    289   {
    290     a = (number)((unsigned long)a / 2);
    291     b = (number)((unsigned long)b / 2);
    292     res++;
    293   }
    294   if ((unsigned long)b % 2 == 0)
    295   {
    296     *t = NULL;
    297     *s = nr2mInvers(a,r);
    298     return (number)((1L << res)); // * (unsigned long) a);  // (2**res)*a    a is a unit
    299   }
    300   else
    301   {
    302     *s = NULL;
    303     *t = nr2mInvers(b,r);
    304     return (number)((1L << res)); // * (unsigned long) b);  // (2**res)*b    b is a unit
    305   }
    306 }
    307 
    308 void nr2mPower(number a, int i, number * result, const coeffs r)
    309 {
    310   if (i == 0)
    311   {
    312     *(unsigned long *)result = 1;
    313   }
    314   else if (i == 1)
    315   {
    316     *result = a;
    317   }
    318   else
    319   {
    320     nr2mPower(a, i-1, result, r);
    321     *result = nr2mMultM(a, *result, r);
    322   }
    323 }
    324 
    325 /*
    326  * create a number from int
    327  */
    328 number nr2mInit(long i, const coeffs r)
    329 {
    330   if (i == 0) return (number)(unsigned long)i;
    331 
    332   long ii = i;
    333   unsigned long j = (unsigned long)1;
    334   if (ii < 0) { j = r->mod2mMask; ii = -ii; }
    335   unsigned long k = (unsigned long)ii;
    336   k = k & r->mod2mMask;
    337   /* now we have: i = j * k mod 2^m */
    338   return (number)nr2mMult((number)j, (number)k, r);
    339 }
    340 
    341 /*
    342  * convert a number to an int in ]-k/2 .. k/2],
    343  * where k = 2^m; i.e., an int in ]-2^(m-1) .. 2^(m-1)];
    344  */
    345 long nr2mInt(number &n, const coeffs r)
    346 {
    347   unsigned long nn = (unsigned long)(unsigned long)n & r->mod2mMask;
    348   unsigned long l = r->mod2mMask >> 1; l++; /* now: l = 2^(m-1) */
    349   if ((unsigned long)nn > l)
    350     return (long)((unsigned long)nn - r->mod2mMask - 1);
    351   else
    352     return (long)((unsigned long)nn);
    353 }
    354 
    355 number nr2mAdd(number a, number b, const coeffs r)
    356 {
    357   return nr2mAddM(a, b, r);
    358 }
    359 
    360 number nr2mSub(number a, number b, const coeffs r)
    361 {
    362   return nr2mSubM(a, b, r);
    363 }
    364 
    365 BOOLEAN nr2mIsUnit(number a, const coeffs)
    366 {
    367   return ((unsigned long)a % 2 == 1);
    368 }
    369 
    370 number nr2mGetUnit(number k, const coeffs)
    371 {
    372   if (k == NULL) return (number)1;
    373   unsigned long erg = (unsigned long)k;
    374   while (erg % 2 == 0) erg = erg / 2;
    375   return (number)erg;
    376 }
    377 
    378 BOOLEAN nr2mIsZero(number a, const coeffs)
    379 {
    380   return 0 == (unsigned long)a;
    381 }
    382 
    383 BOOLEAN nr2mIsOne(number a, const coeffs)
    384 {
    385   return 1 == (unsigned long)a;
    386 }
    387 
    388 BOOLEAN nr2mIsMOne(number a, const coeffs r)
    389 {
    390   return ((r->mod2mMask  == (unsigned long)a) &&(1L!=(long)a))/*for char 2^1*/;
    391 }
    392 
    393 BOOLEAN nr2mEqual(number a, number b, const coeffs)
    394 {
    395   return (a == b);
    396 }
    397 
    398 BOOLEAN nr2mGreater(number a, number b, const coeffs r)
    399 {
    400   return nr2mDivBy(a, b,r);
    401 }
    402 
    403 /* Is 'a' divisible by 'b'? There are two cases:
    404    1) a = 0 mod 2^m; then TRUE iff b = 0 or b is a power of 2
    405    2) a, b <> 0; then TRUE iff b/gcd(a, b) is a unit mod 2^m */
    406 BOOLEAN nr2mDivBy (number a, number b, const coeffs r)
    407 {
    408   if (a == NULL)
    409   {
    410     unsigned long c = r->mod2mMask + 1;
    411     if (c != 0) /* i.e., if no overflow */
    412       return (c % (unsigned long)b) == 0;
    413     else
    414     {
    415       /* overflow: we need to check whether b
    416          is zero or a power of 2: */
    417       c = (unsigned long)b;
    418       while (c != 0)
    419       {
    420         if ((c % 2) != 0) return FALSE;
    421         c = c >> 1;
    422       }
    423       return TRUE;
    424     }
    425   }
    426   else
    427   {
    428     number n = nr2mGcd(a, b, r);
    429     n = nr2mDiv(b, n, r);
    430     return nr2mIsUnit(n, r);
    431   }
    432 }
    433 
    434 int nr2mDivComp(number as, number bs, const coeffs)
    435 {
    436   unsigned long a = (unsigned long)as;
    437   unsigned long b = (unsigned long)bs;
    438   assume(a != 0 && b != 0);
    439   while (a % 2 == 0 && b % 2 == 0)
    440   {
    441     a = a / 2;
    442     b = b / 2;
    443   }
    444   if (a % 2 == 0)
    445   {
    446     return -1;
    447   }
    448   else
    449   {
    450     if (b % 2 == 1)
    451     {
    452       return 2;
    453     }
    454     else
    455     {
    456       return 1;
    457     }
    458   }
    459 }
    460 
    461 /* TRUE iff 0 < k <= 2^m / 2 */
    462 BOOLEAN nr2mGreaterZero(number k, const coeffs r)
    463 {
    464   if ((unsigned long)k == 0) return FALSE;
    465   if ((unsigned long)k > ((r->mod2mMask >> 1) + 1)) return FALSE;
    466   return TRUE;
    467 }
    468 
    469 /* assumes that 'a' is odd, i.e., a unit in Z/2^m, and computes
    470    the extended gcd of 'a' and 2^m, in order to find some 's'
    471    and 't' such that a * s + 2^m * t = gcd(a, 2^m) = 1;
    472    this code will always find a positive 's' */
    473 void specialXGCD(unsigned long& s, unsigned long a, const coeffs r)
    474 {
    475   mpz_ptr u = (mpz_ptr)omAlloc(sizeof(mpz_t));
    476   mpz_init_set_ui(u, a);
    477   mpz_ptr u0 = (mpz_ptr)omAlloc(sizeof(mpz_t));
    478   mpz_init(u0);
    479   mpz_ptr u1 = (mpz_ptr)omAlloc(sizeof(mpz_t));
    480   mpz_init_set_ui(u1, 1);
    481   mpz_ptr u2 = (mpz_ptr)omAlloc(sizeof(mpz_t));
    482   mpz_init(u2);
    483   mpz_ptr v = (mpz_ptr)omAlloc(sizeof(mpz_t));
    484   mpz_init_set_ui(v, r->mod2mMask);
    485   mpz_add_ui(v, v, 1); /* now: v = 2^m */
    486   mpz_ptr v0 = (mpz_ptr)omAlloc(sizeof(mpz_t));
    487   mpz_init(v0);
    488   mpz_ptr v1 = (mpz_ptr)omAlloc(sizeof(mpz_t));
    489   mpz_init(v1);
    490   mpz_ptr v2 = (mpz_ptr)omAlloc(sizeof(mpz_t));
    491   mpz_init_set_ui(v2, 1);
    492   mpz_ptr q = (mpz_ptr)omAlloc(sizeof(mpz_t));
    493   mpz_init(q);
    494   mpz_ptr rr = (mpz_ptr)omAlloc(sizeof(mpz_t));
    495   mpz_init(rr);
    496 
    497   while (mpz_cmp_ui(v, 0) != 0) /* i.e., while v != 0 */
    498   {
    499     mpz_div(q, u, v);
    500     mpz_mod(rr, u, v);
    501     mpz_set(u, v);
    502     mpz_set(v, rr);
    503     mpz_set(u0, u2);
    504     mpz_set(v0, v2);
    505     mpz_mul(u2, u2, q); mpz_sub(u2, u1, u2); /* u2 = u1 - q * u2 */
    506     mpz_mul(v2, v2, q); mpz_sub(v2, v1, v2); /* v2 = v1 - q * v2 */
    507     mpz_set(u1, u0);
    508     mpz_set(v1, v0);
    509   }
    510 
    511   while (mpz_cmp_ui(u1, 0) < 0) /* i.e., while u1 < 0 */
    512   {
    513     /* we add 2^m = (2^m - 1) + 1 to u1: */
    514     mpz_add_ui(u1, u1, r->mod2mMask);
    515     mpz_add_ui(u1, u1, 1);
    516   }
    517   s = mpz_get_ui(u1); /* now: 0 <= s <= 2^m - 1 */
    518 
    519   mpz_clear(u);  omFree((ADDRESS)u);
    520   mpz_clear(u0); omFree((ADDRESS)u0);
    521   mpz_clear(u1); omFree((ADDRESS)u1);
    522   mpz_clear(u2); omFree((ADDRESS)u2);
    523   mpz_clear(v);  omFree((ADDRESS)v);
    524   mpz_clear(v0); omFree((ADDRESS)v0);
    525   mpz_clear(v1); omFree((ADDRESS)v1);
    526   mpz_clear(v2); omFree((ADDRESS)v2);
    527   mpz_clear(q); omFree((ADDRESS)q);
    528   mpz_clear(rr); omFree((ADDRESS)rr);
    529 }
    530 
    531 unsigned long InvMod(unsigned long a, const coeffs r)
    532 {
    533   assume((unsigned long)a % 2 != 0);
    534   unsigned long s;
    535   specialXGCD(s, a, r);
    536   return s;
    537 }
    538 //#endif
    539 
    540 inline number nr2mInversM(number c, const coeffs r)
    541 {
    542   assume((unsigned long)c % 2 != 0);
    543   // Table !!!
    544   unsigned long inv;
    545   inv = InvMod((unsigned long)c,r);
    546   return (number)inv;
    547 }
    548 
    549 number nr2mDiv(number a, number b, const coeffs r)
    550 {
    551   if ((unsigned long)a == 0) return (number)0;
    552   else if ((unsigned long)b % 2 == 0)
    553   {
    554     if ((unsigned long)b != 0)
    555     {
    556       while (((unsigned long)b % 2 == 0) && ((unsigned long)a % 2 == 0))
    557       {
    558         a = (number)((unsigned long)a / 2);
    559         b = (number)((unsigned long)b / 2);
    560       }
    561     }
    562     if ((unsigned long)b % 2 == 0)
    563     {
    564       WerrorS("Division not possible, even by cancelling zero divisors.");
    565       WerrorS("Result is integer division without remainder.");
    566       return (number) ((unsigned long) a / (unsigned long) b);
    567     }
    568   }
    569   return (number)nr2mMult(a, nr2mInversM(b,r),r);
    570 }
    571 
    572 number nr2mMod(number a, number b, const coeffs r)
    573 {
    574   /*
    575     We need to return the number rr which is uniquely determined by the
    576     following two properties:
    577       (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
    578       (2) There exists some k in the integers Z such that a = k * b + rr.
    579     Consider g := gcd(2^m, |b|). Note that then |b|/g is a unit in Z/2^m.
    580     Now, there are three cases:
    581       (a) g = 1
    582           Then |b| is a unit in Z/2^m, i.e. |b| (and also b) divides a.
    583           Thus rr = 0.
    584       (b) g <> 1 and g divides a
    585           Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
    586       (c) g <> 1 and g does not divide a
    587           Let's denote the division with remainder of a by g as follows:
    588           a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
    589           fulfills (1) and (2), i.e. rr := t is the correct result. Hence
    590           in this third case, rr is the remainder of division of a by g in Z.
    591     This algorithm is the same as for the case Z/n, except that we may
    592     compute the gcd of |b| and 2^m "by hand": We just extract the highest
    593     power of 2 (<= 2^m) that is contained in b.
    594   */
    595   assume((unsigned long) b != 0);
    596   unsigned long g = 1;
    597   unsigned long b_div = (unsigned long) b;
    598 
    599   /*
    600    * b_div is unsigned, so that (b_div < 0) evaluates false at compile-time
    601    *
    602   if (b_div < 0) b_div = -b_div; // b_div now represents |b|, BUT b_div is unsigned!
    603   */
    604 
    605   unsigned long rr = 0;
    606   while ((g < r->mod2mMask ) && (b_div > 0) && (b_div % 2 == 0))
    607   {
    608     b_div = b_div >> 1;
    609     g = g << 1;
    610   } // g is now the gcd of 2^m and |b|
    611 
    612   if (g != 1) rr = (unsigned long)a % g;
    613   return (number)rr;
    614 }
    615 
    616 number nr2mIntDiv(number a, number b, const coeffs r)
    617 {
    618   if ((unsigned long)a == 0)
    619   {
    620     if ((unsigned long)b == 0)
    621       return (number)1;
    622     if ((unsigned long)b == 1)
    623       return (number)0;
    624     unsigned long c = r->mod2mMask + 1;
    625     if (c != 0) /* i.e., if no overflow */
    626       return (number)(c / (unsigned long)b);
    627     else
    628     {
    629       /* overflow: c = 2^32 resp. 2^64, depending on platform */
    630       mpz_ptr cc = (mpz_ptr)omAlloc(sizeof(mpz_t));
    631       mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
    632       mpz_div_ui(cc, cc, (unsigned long)(unsigned long)b);
    633       unsigned long s = mpz_get_ui(cc);
    634       mpz_clear(cc); omFree((ADDRESS)cc);
    635       return (number)(unsigned long)s;
    636     }
    637   }
    638   else
    639   {
    640     if ((unsigned long)b == 0)
    641       return (number)0;
    642     return (number)((unsigned long) a / (unsigned long) b);
    643   }
    644 }
    645 
    646 static number nr2mAnn(number b, const coeffs r)
    647 {
    648   if ((unsigned long)b == 0)
    649     return NULL;
    650   if ((unsigned long)b == 1)
    651     return NULL;
    652   unsigned long c = r->mod2mMask + 1;
    653   if (c != 0) /* i.e., if no overflow */
    654     return (number)(c / (unsigned long)b);
    655   else
    656   {
    657     /* overflow: c = 2^32 resp. 2^64, depending on platform */
    658     mpz_ptr cc = (mpz_ptr)omAlloc(sizeof(mpz_t));
    659     mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
    660     mpz_div_ui(cc, cc, (unsigned long)(unsigned long)b);
    661     unsigned long s = mpz_get_ui(cc);
    662     mpz_clear(cc); omFree((ADDRESS)cc);
    663     return (number)(unsigned long)s;
    664   }
    665 }
    666 
    667 number nr2mInvers(number c, const coeffs r)
    668 {
    669   if ((unsigned long)c % 2 == 0)
    670   {
    671     WerrorS("division by zero divisor");
    672     return (number)0;
    673   }
    674   return nr2mInversM(c, r);
    675 }
    676 
    677 number nr2mNeg(number c, const coeffs r)
    678 {
    679   if ((unsigned long)c == 0) return c;
    680   return nr2mNegM(c, r);
    681 }
    682 
    683 number nr2mMapMachineInt(number from, const coeffs /*src*/, const coeffs dst)
    684 {
    685   unsigned long i = ((unsigned long)from) % dst->mod2mMask ;
    686   return (number)i;
    687 }
    688 
    689 number nr2mMapProject(number from, const coeffs /*src*/, const coeffs dst)
    690 {
    691   unsigned long i = ((unsigned long)from) % (dst->mod2mMask + 1);
    692   return (number)i;
    693 }
    694 
    695 number nr2mMapZp(number from, const coeffs /*src*/, const coeffs dst)
    696 {
    697   unsigned long j = (unsigned long)1;
    698   long ii = (long)from;
    699   if (ii < 0) { j = dst->mod2mMask; ii = -ii; }
    700   unsigned long i = (unsigned long)ii;
    701   i = i & dst->mod2mMask;
    702   /* now we have: from = j * i mod 2^m */
    703   return (number)nr2mMult((number)i, (number)j, dst);
    704 }
    705 
    706 number nr2mMapGMP(number from, const coeffs /*src*/, const coeffs dst)
    707 {
    708   mpz_ptr erg = (mpz_ptr)omAllocBin(gmp_nrz_bin);
    709   mpz_init(erg);
    710   mpz_ptr k = (mpz_ptr)omAlloc(sizeof(mpz_t));
    711   mpz_init_set_ui(k, dst->mod2mMask);
    712 
    713   mpz_and(erg, (mpz_ptr)from, k);
    714   number res = (number) mpz_get_ui(erg);
    715 
    716   mpz_clear(erg); omFree((ADDRESS)erg);
    717   mpz_clear(k);   omFree((ADDRESS)k);
    718 
    719   return (number)res;
    720 }
    721 
    722 number nr2mMapQ(number from, const coeffs src, const coeffs dst)
    723 {
    724   mpz_ptr gmp = (mpz_ptr)omAllocBin(gmp_nrz_bin);
    725   mpz_init(gmp);
    726   nlGMP(from, (number)gmp, src); // FIXME? TODO? // extern void   nlGMP(number &i, number n, const coeffs r); // to be replaced with n_MPZ(erg, from, src); // ?
    727   number res=nr2mMapGMP((number)gmp,src,dst);
    728   mpz_clear(gmp); omFree((ADDRESS)gmp);
    729   return res;
    730 }
    731 
    732 number nr2mMapZ(number from, const coeffs src, const coeffs dst)
    733 {
    734   if (SR_HDL(from) & SR_INT)
    735   {
    736     long f_i=SR_TO_INT(from);
    737     return nr2mInit(f_i,dst);
    738   }
    739   return nr2mMapGMP(from,src,dst);
    740 }
    741 
    742 nMapFunc nr2mSetMap(const coeffs src, const coeffs dst)
    743 {
    744   if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
    745      && (src->mod2mMask == dst->mod2mMask))
    746   {
    747     return ndCopyMap;
    748   }
    749   if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
    750      && (src->mod2mMask < dst->mod2mMask))
    751   { /* i.e. map an integer mod 2^s into Z mod 2^t, where t < s */
    752     return nr2mMapMachineInt;
    753   }
    754   if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)
    755      && (src->mod2mMask > dst->mod2mMask))
    756   { /* i.e. map an integer mod 2^s into Z mod 2^t, where t > s */
    757     // to be done
    758     return nr2mMapProject;
    759   }
    760   if ((src->rep==n_rep_gmp) && nCoeff_is_Ring_Z(src))
    761   {
    762     return nr2mMapGMP;
    763   }
    764   if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
    765   {
    766     return nr2mMapZ;
    767   }
    768   if ((src->rep==n_rep_gap_rat) && (nCoeff_is_Q(src)||nCoeff_is_Ring_Z(src)))
    769   {
    770     return nr2mMapQ;
    771   }
    772   if ((src->rep==n_rep_int) && nCoeff_is_Zp(src) && (src->ch == 2))
    773   {
    774     return nr2mMapZp;
    775   }
    776   if ((src->rep==n_rep_gmp) &&
    777   (nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src)))
    778   {
    779     if (mpz_divisible_2exp_p(src->modNumber,dst->modExponent))
    780       return nr2mMapGMP;
    781   }
    782   return NULL;      // default
    783 }
    784 
    785 /*
    786  * set the exponent
    787  */
    788 
    789 void nr2mSetExp(int m, coeffs r)
    790 {
    791   if (m > 1)
    792   {
    793     /* we want mod2mMask to be the bit pattern
    794        '111..1' consisting of m one's: */
    795     r->modExponent= m;
    796     r->mod2mMask = 1;
    797     for (int i = 1; i < m; i++) r->mod2mMask = (r->mod2mMask << 1) + 1;
    798   }
    799   else
    800   {
    801     r->modExponent= 2;
    802     /* code unexpectedly called with m = 1; we continue with m = 2: */
    803     r->mod2mMask = 3; /* i.e., '11' in binary representation */
    804   }
    805 }
    806 
    807 void nr2mInitExp(int m, coeffs r)
    808 {
    809   nr2mSetExp(m, r);
    810   if (m < 2)
    811     WarnS("nr2mInitExp unexpectedly called with m = 1 (we continue with Z/2^2");
    812 }
    813 
    814 #ifdef LDEBUG
    815 BOOLEAN nr2mDBTest (number a, const char *, const int, const coeffs r)
    816 {
    817   //if ((unsigned long)a < 0) return FALSE; // is unsigned!
    818   if (((unsigned long)a & r->mod2mMask) != (unsigned long)a) return FALSE;
    819   return TRUE;
    820 }
    821 #endif
    822 
    823 void nr2mWrite (number a, const coeffs r)
    824 {
    825   long i = nr2mInt(a, r);
    826   StringAppend("%ld", i);
    827 }
    828 
    829 static const char* nr2mEati(const char *s, int *i, const coeffs r)
    830 {
    831 
    832   if (((*s) >= '0') && ((*s) <= '9'))
    833   {
    834     (*i) = 0;
    835     do
    836     {
    837       (*i) *= 10;
    838       (*i) += *s++ - '0';
    839       if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & r->mod2mMask;
    840     }
    841     while (((*s) >= '0') && ((*s) <= '9'));
    842     (*i) = (*i) & r->mod2mMask;
    843   }
    844   else (*i) = 1;
    845   return s;
    846 }
    847 
    848 const char * nr2mRead (const char *s, number *a, const coeffs r)
    849 {
    850   int z;
    851   int n=1;
    852 
    853   s = nr2mEati(s, &z,r);
    854   if ((*s) == '/')
    855   {
    856     s++;
    857     s = nr2mEati(s, &n,r);
    858   }
    859   if (n == 1)
    860     *a = (number)(long)z;
    861   else
    862       *a = nr2mDiv((number)(long)z,(number)(long)n,r);
    863   return s;
    864 }
    865838#endif
    866839/* #ifdef HAVE_RINGS */
Note: See TracChangeset for help on using the changeset viewer.