Changeset e90dfd6 in git


Ignore:
Timestamp:
Sep 6, 2010, 12:08:33 PM (14 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
e903733fa2ff9e4f7582f1faabf6cc2a6baea275
Parents:
5d594a96114bdb26100c2b8e2ae329a7529f0e45
git-author:
Frank Seelisch <seelisch@mathematik.uni-kl.de>2010-09-06 12:08:33+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 11:55:29+01:00
Message:
fixed svn merge bugs in Z/n and Z/2^m
Location:
coeffs
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • coeffs/coeffs.h

    • Property mode changed from 100644 to 100755
    r5d594a9 re90dfd6  
    190190
    191191#ifdef HAVE_RINGS
    192   int_number    ringflaga; /* Z/(ringflag^ringflagb)=Z/nrnModul*/
    193   unsigned long ringflagb;
    194   unsigned long nr2mModul;  /* Z/nr2mModul */
    195   int_number    nrnModul;
     192  /* The following members are for representing the ring Z/n,
     193     where n is not a prime. We distinguish three cases:
     194     1.) n has at least two distinct prime factors. Then
     195         modBase stores n, modExponent stores 1, modNumber
     196         stores n, and mod2mMask is not used;
     197     2.) n = p^k for some odd prime p and k > 1. Then
     198         modBase stores p, modExponent stores k, modNumber
     199         stores n, and mod2mMask is not used;
     200     3.) n = 2^k for some k > 1; moreover, 2^k - 1 fits in
     201         an unsigned long. Then modBase stores 2, modExponent
     202         stores k, modNumber is not used, and mod2mMask stores
     203         2^k - 1, i.e., the bit mask '111..1' of length k.
     204     4.) n = 2^k for some k > 1; but 2^k - 1 does not fit in
     205         an unsigned long. Then modBase stores 2, modExponent
     206         stores k, modNumber stores n, and mod2mMask is not
     207         used;
     208     Cases 1.), 2.), and 4.) are covered by the implementation
     209     in the files rmodulon.h and rmodulon.cc, whereas case 3.)
     210     is implemented in the files rmodulo2m.h and rmodulo2m.cc. */
     211  int_number    modBase;
     212  unsigned long modExponent;
     213  int_number    modNumber;
     214  unsigned long mod2mMask;
    196215#endif
    197216  int        ch;  /* characteristic, rInit */

  • coeffs/rmodulo2m.cc

    • Property mode changed from 100644 to 100755
    r5d594a9 re90dfd6  
    2424#include <string.h>
    2525
    26 int nr2mExp;
    27 
    2826extern omBin gmp_nrz_bin; /* init in rintegers*/
    2927
    … …  
    3432  nr2mInitExp((int)(long)(p), r);
    3533
    36   r->cfInit       = nr2mInit;
    37   r->cfCopy       = ndCopy;
    38   r->cfInt        = nr2mInt;
     34  r->cfInit        = nr2mInit;
     35  r->cfCopy        = ndCopy;
     36  r->cfInt         = nr2mInt;
    3937  r->cfAdd         = nr2mAdd;
    4038  r->cfSub         = nr2mSub;
    … …  
    5452  r->cfIsMOne      = nr2mIsMOne;
    5553  r->cfGreaterZero = nr2mGreaterZero;
    56   r->cfWrite      = nr2mWrite;
     54  r->cfWrite       = nr2mWrite;
    5755  r->cfRead        = nr2mRead;
    5856  r->cfPower       = nr2mPower;
    59   r->cfSetMap     = nr2mSetMap;
     57  r->cfSetMap      = nr2mSetMap;
    6058  r->cfNormalize   = ndNormalize;
    6159  r->cfLcm         = nr2mLcm;
    … …  
    7573 * Multiply two numbers
    7674 */
    77 number nr2mMult (number a, number b, const coeffs r)
     75number nr2mMult(number a, number b, const coeffs r)
    7876{
    7977  if (((NATNUMBER)a == 0) || ((NATNUMBER)b == 0))
    8078    return (number)0;
    8179  else
    82     return nr2mMultM(a,b,r);
    83 }
    84 
    85 /*
    86  * Give the smallest non unit k, such that a * x = k = b * y has a solution
     80    return nr2mMultM(a, b, r);
     81}
     82
     83/*
     84 * Give the smallest k, such that a * x = k = b * y has a solution
    8785 */
    88 number nr2mLcm (number a, number b, const coeffs r)
     86number nr2mLcm(number a, number b, const coeffs r)
    8987{
    9088  NATNUMBER res = 0;
    91   if ((NATNUMBER) a == 0) a = (number) 1;
    92   if ((NATNUMBER) b == 0) b = (number) 1;
    93   while ((NATNUMBER) a % 2 == 0)
    94   {
    95     a = (number) ((NATNUMBER) a / 2);
    96     if ((NATNUMBER) b % 2 == 0) b = (number) ((NATNUMBER) b / 2);
     89  if ((NATNUMBER)a == 0) a = (number) 1;
     90  if ((NATNUMBER)b == 0) b = (number) 1;
     91  while ((NATNUMBER)a % 2 == 0)
     92  {
     93    a = (number)((NATNUMBER)a / 2);
     94    if ((NATNUMBER)b % 2 == 0) b = (number)((NATNUMBER)b / 2);
    9795    res++;
    9896  }
    99   while ((NATNUMBER) b % 2 == 0)
    100   {
    101     b = (number) ((NATNUMBER) b / 2);
     97  while ((NATNUMBER)b % 2 == 0)
     98  {
     99    b = (number)((NATNUMBER)b / 2);
    102100    res++;
    103101  }
    104   return (number) (1L << res);  // (2**res)
    105 }
    106 
    107 /*
    108  * Give the largest non unit k, such that a = x * k, b = y * k has
     102  return (number)(1L << res);  // (2**res)
     103}
     104
     105/*
     106 * Give the largest k, such that a = x * k, b = y * k has
    109107 * a solution.
    110108 */
    111 number nr2mGcd (number a, number b, const coeffs r)
     109number nr2mGcd(number a, number b, const coeffs r)
    112110{
    113111  NATNUMBER res = 0;
    114   if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1;
    115   while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
    116   {
    117     a = (number) ((NATNUMBER) a / 2);
    118     b = (number) ((NATNUMBER) b / 2);
     112  if ((NATNUMBER)a == 0 && (NATNUMBER)b == 0) return (number)1;
     113  while ((NATNUMBER)a % 2 == 0 && (NATNUMBER)b % 2 == 0)
     114  {
     115    a = (number)((NATNUMBER)a / 2);
     116    b = (number)((NATNUMBER)b / 2);
    119117    res++;
    120118  }
    121 //  if ((NATNUMBER) b % 2 == 0)
     119//  if ((NATNUMBER)b % 2 == 0)
    122120//  {
    123 //    return (number) ((1L << res));// * (NATNUMBER) a);  // (2**res)*a    a ist Einheit
     121//    return (number)((1L << res)); // * (NATNUMBER) a);  // (2**res)*a    a is a unit
    124122//  }
    125123//  else
    126124//  {
    127     return (number) ((1L << res));// * (NATNUMBER) b);  // (2**res)*b    b ist Einheit
     125    return (number)((1L << res)); // * (NATNUMBER) b);  // (2**res)*b    b is a unit
    128126//  }
    129127}
    130128
    131129/*
    132  * Give the largest non unit k, such that a = x * k, b = y * k has
     130 * Give the largest k, such that a = x * k, b = y * k has
    133131 * a solution.
    134132 */
    135 number nr2mExtGcd (number a, number b, number *s, number *t, const coeffs r)
     133number nr2mExtGcd(number a, number b, number *s, number *t, const coeffs r)
    136134{
    137135  NATNUMBER res = 0;
    138   if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1;
    139   while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0)
    140   {
    141     a = (number) ((NATNUMBER) a / 2);
    142     b = (number) ((NATNUMBER) b / 2);
     136  if ((NATNUMBER)a == 0 && (NATNUMBER)b == 0) return (number)1;
     137  while ((NATNUMBER)a % 2 == 0 && (NATNUMBER)b % 2 == 0)
     138  {
     139    a = (number)((NATNUMBER)a / 2);
     140    b = (number)((NATNUMBER)b / 2);
    143141    res++;
    144142  }
    145   if ((NATNUMBER) b % 2 == 0)
     143  if ((NATNUMBER)b % 2 == 0)
    146144  {
    147145    *t = NULL;
    148146    *s = nr2mInvers(a,r);
    149     return (number) ((1L << res));// * (NATNUMBER) a);  // (2**res)*a    a ist Einheit
     147    return (number)((1L << res)); // * (NATNUMBER) a);  // (2**res)*a    a is a unit
    150148  }
    151149  else
    … …  
    153151    *s = NULL;
    154152    *t = nr2mInvers(b,r);
    155     return (number) ((1L << res));// * (NATNUMBER) b);  // (2**res)*b    b ist Einheit
    156   }
    157 }
    158 
    159 void nr2mPower (number a, int i, number * result, const coeffs r)
    160 {
    161   if (i==0)
     153    return (number)((1L << res)); // * (NATNUMBER) b);  // (2**res)*b    b is a unit
     154  }
     155}
     156
     157void nr2mPower(number a, int i, number * result, const coeffs r)
     158{
     159  if (i == 0)
    162160  {
    163161    *(NATNUMBER *)result = 1;
    164162  }
    165   else if (i==1)
     163  else if (i == 1)
    166164  {
    167165    *result = a;
    … …  
    169167  else
    170168  {
    171     nr2mPower(a,i-1,result,r);
    172     *result = nr2mMultM(a,*result,r);
     169    nr2mPower(a, i-1, result, r);
     170    *result = nr2mMultM(a, *result, r);
    173171  }
    174172}
    … …  
    177175 * create a number from int
    178176 */
    179 number nr2mInit (int i, const coeffs r)
     177number nr2mInit(int i, const coeffs r)
    180178{
    181179  if (i == 0) return (number)(NATNUMBER)i;
    … …  
    183181  long ii = i;
    184182  NATNUMBER j = (NATNUMBER)1;
    185   if (ii < 0) { j = r->nr2mModul; ii = -ii; }
     183  if (ii < 0) { j = r->mod2mMask; ii = -ii; }
    186184  NATNUMBER k = (NATNUMBER)ii;
    187   k = k & r->nr2mModul;
    188   /* now we have: from = j * k mod 2^m */
    189   return (number)nr2mMult((number)j, (number)k,r);
     185  k = k & r->mod2mMask;
     186  /* now we have: i = j * k mod 2^m */
     187  return (number)nr2mMult((number)j, (number)k, r);
    190188}
    191189
    … …  
    198196int nr2mInt(number &n, const coeffs r)
    199197{
    200   NATNUMBER nn = (unsigned long)(NATNUMBER)n & r->nr2mModul;
    201   unsigned long l = r->nr2mModul >> 1; l++; /* now: l = 2^(m-1) */
     198  NATNUMBER nn = (unsigned long)(NATNUMBER)n & r->mod2mMask;
     199  unsigned long l = r->mod2mMask >> 1; l++; /* now: l = 2^(m-1) */
    202200  if ((NATNUMBER)nn > l)
    203     return (int)((NATNUMBER)nn - r->nr2mModul - 1);
     201    return (int)((NATNUMBER)nn - r->mod2mMask - 1);
    204202  else
    205203    return (int)((NATNUMBER)nn);
    206204}
    207205
    208 number nr2mAdd (number a, number b, const coeffs r)
    209 {
    210   return nr2mAddM(a,b,r);
    211 }
    212 
    213 number nr2mSub (number a, number b, const coeffs r)
    214 {
    215   return nr2mSubM(a,b,r);
    216 }
    217 
    218 BOOLEAN nr2mIsUnit (number a, const coeffs r)
    219 {
    220   return ((NATNUMBER) a % 2 == 1);
    221 }
    222 
    223 number  nr2mGetUnit (number k, const coeffs r)
    224 {
    225   if (k == NULL)
    226     return (number) 1;
    227   NATNUMBER tmp = (NATNUMBER) k;
    228   while (tmp % 2 == 0)
    229     tmp = tmp / 2;
    230   return (number) tmp;
    231 }
    232 
    233 BOOLEAN nr2mIsZero (number a, const coeffs r)
     206number nr2mAdd(number a, number b, const coeffs r)
     207{
     208  return nr2mAddM(a, b, r);
     209}
     210
     211number nr2mSub(number a, number b, const coeffs r)
     212{
     213  return nr2mSubM(a, b, r);
     214}
     215
     216BOOLEAN nr2mIsUnit(number a, const coeffs r)
     217{
     218  return ((NATNUMBER)a % 2 == 1);
     219}
     220
     221number nr2mGetUnit(number k, const coeffs r)
     222{
     223  if (k == NULL) return (number)1;
     224  NATNUMBER erg = (NATNUMBER)k;
     225  while (erg % 2 == 0) erg = erg / 2;
     226  return (number)erg;
     227}
     228
     229BOOLEAN nr2mIsZero(number a, const coeffs r)
    234230{
    235231  return 0 == (NATNUMBER)a;
    236232}
    237233
    238 BOOLEAN nr2mIsOne (number a, const coeffs r)
     234BOOLEAN nr2mIsOne(number a, const coeffs r)
    239235{
    240236  return 1 == (NATNUMBER)a;
    241237}
    242238
    243 BOOLEAN nr2mIsMOne (number a, const coeffs r)
    244 {
    245   return (r->nr2mModul  == (NATNUMBER)a)
    246         && (r->nr2mModul != 2);
    247 }
    248 
    249 BOOLEAN nr2mEqual (number a, number b, const coeffs r)
    250 {
    251   return a==b;
    252 }
    253 
    254 BOOLEAN nr2mGreater (number a, number b, const coeffs r)
     239BOOLEAN nr2mIsMOne(number a, const coeffs r)
     240{
     241  return (r->mod2mMask  == (NATNUMBER)a);
     242}
     243
     244BOOLEAN nr2mEqual(number a, number b, const coeffs r)
     245{
     246  return a == b;
     247}
     248
     249BOOLEAN nr2mGreater(number a, number b, const coeffs r)
    255250{
    256251  return nr2mDivBy(a, b,r);
    … …  
    265260  if (a == NULL)
    266261  {
    267     NATNUMBER c = r->nr2mModul + 1;
     262    NATNUMBER c = r->mod2mMask + 1;
    268263    if (c != 0) /* i.e., if no overflow */
    269264      return (c % (NATNUMBER)b) == 0;
    … …  
    291286int nr2mDivComp(number as, number bs, const coeffs r)
    292287{
    293   NATNUMBER a = (NATNUMBER) as;
    294   NATNUMBER b = (NATNUMBER) bs;
     288  NATNUMBER a = (NATNUMBER)as;
     289  NATNUMBER b = (NATNUMBER)bs;
    295290  assume(a != 0 && b != 0);
    296291  while (a % 2 == 0 && b % 2 == 0)
    … …  
    317312
    318313/* TRUE iff 0 < k <= 2^m / 2 */
    319 BOOLEAN nr2mGreaterZero (number k, const coeffs r)
     314BOOLEAN nr2mGreaterZero(number k, const coeffs r)
    320315{
    321316  if ((NATNUMBER)k == 0) return FALSE;
    322   if ((NATNUMBER)k > ((r->nr2mModul >> 1) + 1)) return FALSE;
     317  if ((NATNUMBER)k > ((r->mod2mMask >> 1) + 1)) return FALSE;
    323318  return TRUE;
    324319}
    … …  
    328323   and 't' such that a * s + 2^m * t = gcd(a, 2^m) = 1;
    329324   this code will always find a positive 's' */
    330 void specialXGCD(unsigned long& s, unsigned long a, const coeffs R)
     325void specialXGCD(unsigned long& s, unsigned long a, const coeffs r)
    331326{
    332327  int_number u = (int_number)omAlloc(sizeof(mpz_t));
    … …  
    339334  mpz_init(u2);
    340335  int_number v = (int_number)omAlloc(sizeof(mpz_t));
    341   mpz_init_set_ui(v, R->nr2mModul);
     336  mpz_init_set_ui(v, r->mod2mMask);
    342337  mpz_add_ui(v, v, 1); /* now: v = 2^m */
    343338  int_number v0 = (int_number)omAlloc(sizeof(mpz_t));
    … …  
    349344  int_number q = (int_number)omAlloc(sizeof(mpz_t));
    350345  mpz_init(q);
    351   int_number r = (int_number)omAlloc(sizeof(mpz_t));
    352   mpz_init(r);
     346  int_number rr = (int_number)omAlloc(sizeof(mpz_t));
     347  mpz_init(rr);
    353348
    354349  while (mpz_cmp_ui(v, 0) != 0) /* i.e., while v != 0 */
    355350  {
    356351    mpz_div(q, u, v);
    357     mpz_mod(r, u, v);
     352    mpz_mod(rr, u, v);
    358353    mpz_set(u, v);
    359     mpz_set(v, r);
     354    mpz_set(v, rr);
    360355    mpz_set(u0, u2);
    361356    mpz_set(v0, v2);
    … …  
    369364  {
    370365    /* we add 2^m = (2^m - 1) + 1 to u1: */
    371     mpz_add_ui(u1, u1, R->nr2mModul);
     366    mpz_add_ui(u1, u1, r->mod2mMask);
    372367    mpz_add_ui(u1, u1, 1);
    373368  }
    … …  
    383378  mpz_clear(v2); omFree((ADDRESS)v2);
    384379  mpz_clear(q); omFree((ADDRESS)q);
    385   mpz_clear(r); omFree((ADDRESS)r);
     380  mpz_clear(rr); omFree((ADDRESS)rr);
    386381}
    387382
    … …  
    395390//#endif
    396391
    397 inline number nr2mInversM (number c, const coeffs r)
     392inline number nr2mInversM(number c, const coeffs r)
    398393{
    399394  assume((NATNUMBER)c % 2 != 0);
    … …  
    401396  NATNUMBER inv;
    402397  inv = InvMod((NATNUMBER)c,r);
    403   return (number) inv;
    404 }
    405 
    406 number nr2mDiv (number a, number b, const coeffs r)
    407 {
    408   if ((NATNUMBER)a==0)
    409     return (number)0;
    410   else if ((NATNUMBER)b%2==0)
     398  return (number)inv;
     399}
     400
     401number nr2mDiv(number a, number b, const coeffs r)
     402{
     403  if ((NATNUMBER)a == 0) return (number)0;
     404  else if ((NATNUMBER)b % 2 == 0)
    411405  {
    412406    if ((NATNUMBER)b != 0)
    413407    {
    414       while ((NATNUMBER) b%2 == 0 && (NATNUMBER) a%2 == 0)
     408      while (((NATNUMBER)b % 2 == 0) && ((NATNUMBER)a % 2 == 0))
    415409      {
    416         a = (number) ((NATNUMBER) a / 2);
    417         b = (number) ((NATNUMBER) b / 2);
     410        a = (number)((NATNUMBER)a / 2);
     411        b = (number)((NATNUMBER)b / 2);
    418412      }
    419413    }
    420     if ((NATNUMBER) b%2 == 0)
     414    if ((NATNUMBER)b % 2 == 0)
    421415    {
    422416      WerrorS("Division not possible, even by cancelling zero divisors.");
    … …  
    425419    }
    426420  }
    427   return (number) nr2mMult(a, nr2mInversM(b,r),r);
    428 }
    429 
    430 number nr2mMod (number a, number b, const coeffs R)
     421  return (number)nr2mMult(a, nr2mInversM(b,r),r);
     422}
     423
     424number nr2mMod(number a, number b, const coeffs r)
    431425{
    432426  /*
    433     We need to return the number r which is uniquely determined by the
     427    We need to return the number rr which is uniquely determined by the
    434428    following two properties:
    435       (1) 0 <= r < |b| (with respect to '<' and '<=' performed in Z x Z)
    436       (2) There exists some k in the integers Z such that a = k * b + r.
     429      (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
     430      (2) There exists some k in the integers Z such that a = k * b + rr.
    437431    Consider g := gcd(2^m, |b|). Note that then |b|/g is a unit in Z/2^m.
    438432    Now, there are three cases:
    439433      (a) g = 1
    440434          Then |b| is a unit in Z/2^m, i.e. |b| (and also b) divides a.
    441           Thus r = 0.
     435          Thus rr = 0.
    442436      (b) g <> 1 and g divides a
    443           Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again r = 0.
     437          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
    444438      (c) g <> 1 and g does not divide a
    445439          Let's denote the division with remainder of a by g as follows:
    446440          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
    447           fulfills (1) and (2), i.e. r := t is the correct result. Hence
    448           in this third case, r is the remainder of division of a by g in Z.
     441          fulfills (1) and (2), i.e. rr := t is the correct result. Hence
     442          in this third case, rr is the remainder of division of a by g in Z.
    449443    This algorithm is the same as for the case Z/n, except that we may
    450444    compute the gcd of |b| and 2^m "by hand": We just extract the highest
    … …  
    455449  NATNUMBER b_div = (NATNUMBER)b;
    456450  if (b_div < 0) b_div = -b_div; // b_div now represents |b|
    457   NATNUMBER r = 0;
    458   while ((g < R->nr2mModul ) && (b_div > 0) && (b_div % 2 == 0))
     451  NATNUMBER rr = 0;
     452  while ((g < r->mod2mMask ) && (b_div > 0) && (b_div % 2 == 0))
    459453  {
    460454    b_div = b_div >> 1;
    … …  
    462456  } // g is now the gcd of 2^m and |b|
    463457
    464   if (g != 1) r = (NATNUMBER)a % g;
    465   return (number)r;
    466 }
    467 
    468 number nr2mIntDiv (number a, number b, const coeffs r)
     458  if (g != 1) rr = (NATNUMBER)a % g;
     459  return (number)rr;
     460}
     461
     462number nr2mIntDiv(number a, number b, const coeffs r)
    469463{
    470464  if ((NATNUMBER)a == 0)
    … …  
    474468    if ((NATNUMBER)b == 1)
    475469      return (number)0;
    476     NATNUMBER c = r->nr2mModul + 1;
     470    NATNUMBER c = r->mod2mMask + 1;
    477471    if (c != 0) /* i.e., if no overflow */
    478472      return (number)(c / (NATNUMBER)b);
    … …  
    481475      /* overflow: c = 2^32 resp. 2^64, depending on platform */
    482476      int_number cc = (int_number)omAlloc(sizeof(mpz_t));
    483       mpz_init_set_ui(cc, r->nr2mModul); mpz_add_ui(cc, cc, 1);
     477      mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1);
    484478      mpz_div_ui(cc, cc, (unsigned long)(NATNUMBER)b);
    485479      unsigned long s = mpz_get_ui(cc);
    … …  
    496490}
    497491
    498 number  nr2mInvers (number c, const coeffs r)
    499 {
    500   if ((NATNUMBER)c%2==0)
     492number nr2mInvers(number c, const coeffs r)
     493{
     494  if ((NATNUMBER)c % 2 == 0)
    501495  {
    502496    WerrorS("division by zero divisor");
    503497    return (number)0;
    504498  }
    505   return nr2mInversM(c,r);
    506 }
    507 
    508 number nr2mNeg (number c, const coeffs r)
    509 {
    510   if ((NATNUMBER)c==0) return c;
    511   return nr2mNegM(c,r);
     499  return nr2mInversM(c, r);
     500}
     501
     502number nr2mNeg(number c, const coeffs r)
     503{
     504  if ((NATNUMBER)c == 0) return c;
     505  return nr2mNegM(c, r);
    512506}
    513507
    514508number nr2mMapMachineInt(number from, const coeffs src, const coeffs dst)
    515509{
    516   NATNUMBER i = ((NATNUMBER) from) % dst->nr2mModul ;
    517   return (number) i;
     510  NATNUMBER i = ((NATNUMBER)from) % dst->mod2mMask ;
     511  return (number)i;
    518512}
    519513
    … …  
    521515{
    522516  NATNUMBER j = (NATNUMBER)1;
    523   long ii = (long) from;
    524   if (ii < 0) { j = dst->nr2mModul; ii = -ii; }
     517  long ii = (long)from;
     518  if (ii < 0) { j = dst->mod2mMask; ii = -ii; }
    525519  NATNUMBER i = (NATNUMBER)ii;
    526   i = i & dst->nr2mModul;
     520  i = i & dst->mod2mMask;
    527521  /* now we have: from = j * i mod 2^m */
    528522  return (number)nr2mMult((number)i, (number)j, dst);
    … …  
    531525number nr2mMapQ(number from, const coeffs src, const coeffs dst)
    532526{
    533   int_number erg = (int_number)  omAllocBin(gmp_nrz_bin);
     527  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
    534528  mpz_init(erg);
    535   int_number k = (int_number) omAlloc(sizeof(mpz_t));
    536   mpz_init_set_ui(k, dst->nr2mModul);
     529  int_number k = (int_number)omAlloc(sizeof(mpz_t));
     530  mpz_init_set_ui(k, dst->mod2mMask);
    537531
    538532  nlGMP(from, (number)erg, dst);
    … …  
    543537  mpz_clear(k);   omFree((ADDRESS)k);
    544538
    545   return (number) res;
     539  return (number)res;
    546540}
    547541
    548542number nr2mMapGMP(number from, const coeffs src, const coeffs dst)
    549543{
    550   int_number erg = (int_number)  omAllocBin(gmp_nrz_bin);
     544  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
    551545  mpz_init(erg);
    552   int_number k = (int_number) omAlloc(sizeof(mpz_t));
    553   mpz_init_set_ui(k, dst->nr2mModul);
     546  int_number k = (int_number)omAlloc(sizeof(mpz_t));
     547  mpz_init_set_ui(k, dst->mod2mMask);
    554548
    555549  mpz_and(erg, (int_number)from, k);
    … …  
    559553  mpz_clear(k);   omFree((ADDRESS)k);
    560554
    561   return (number) res;
     555  return (number)res;
    562556}
    563557
    … …  
    565559{
    566560  if (nCoeff_is_Ring_2toM(src)
    567      && (src->ringflagb == dst->ringflagb))
     561     && (src->mod2mMask == dst->mod2mMask))
    568562  {
    569563    return ndCopyMap;
    … …  
    587581    return nr2mMapQ;
    588582  }
    589   if (nCoeff_is_Zp(src)
    590      && (src->ch == 2)
    591      && (dst->ringflagb == 1))
     583  if (nCoeff_is_Zp(src) && (src->ch == 2))
    592584  {
    593585    return nr2mMapZp;
    … …  
    596588  {
    597589    // Computing the n of Z/n
    598     int_number modul = (int_number)  omAllocBin(gmp_nrz_bin);
    599     mpz_init(modul);
    600     mpz_set(modul, src->ringflaga);
    601     mpz_pow_ui(modul, modul, src->ringflagb);
    602     if (mpz_divisible_2exp_p(modul, dst->ringflagb))
    603     {
    604       mpz_clear(modul);
    605       omFree((void *) modul);
     590    int_number modul = (int_number)omAllocBin(gmp_nrz_bin);
     591    mpz_init_set(modul, src->modNumber);
     592    int_number twoToTheK = (int_number)omAllocBin(gmp_nrz_bin);
     593    mpz_init_set_ui(twoToTheK, src->mod2mMask);
     594    mpz_add_ui(twoToTheK, twoToTheK, 1);
     595    if (mpz_divisible_p(modul, twoToTheK))
     596    {
     597      mpz_clear(modul);     omFree((void *)modul);
     598      mpz_clear(twoToTheK); omFree((void *)twoToTheK);
    606599      return nr2mMapGMP;
    607600    }
    608     mpz_clear(modul);
    609     omFree((void *) modul);
     601    mpz_clear(modul);     omFree((void *) modul);
     602    mpz_clear(twoToTheK); omFree((void *)twoToTheK);
    610603  }
    611604  return NULL;      // default
    … …  
    613606
    614607/*
    615  * set the exponent (allocate and init tables) (TODO)
     608 * set the exponent
    616609 */
    617610
    … …  
    620613  if (m > 1)
    621614  {
    622     nr2mExp = m;
    623     /* we want nr2mModul to be the bit pattern '11..1' consisting
    624        of m one's: */
    625     r->nr2mModul = 1;
    626     for (int i = 1; i < m; i++) r->nr2mModul = (r->nr2mModul * 2) + 1;
    627   }
    628   else
    629   {
    630     nr2mExp = 2;
    631     r->nr2mModul = 3; /* i.e., '11' in binary representation */
     615    /* we want mod2mMask to be the bit pattern
     616       '111..1' consisting of m one's: */
     617    r->mod2mMask = 1;
     618    for (int i = 1; i < m; i++) r->mod2mMask = (r->mod2mMask << 1) + 1;
     619  }
     620  else
     621  {
     622    /* code unexpectedly called with m = 1; we go on with m = 2: */
     623    r->mod2mMask = 3; /* i.e., '11' in binary representation */
    632624  }
    633625}
    … …  
    636628{
    637629  nr2mSetExp(m, r);
    638   if (m < 2) WarnS("nr2mInitExp failed: we go on with Z/2^2");
     630  if (m < 2)
     631    WarnS("nr2mInitExp unexpectedly called with m = 1 (we go on with Z/2^2");
    639632}
    640633
    … …  
    643636{
    644637  if ((NATNUMBER)a < 0) return FALSE;
    645   if (((NATNUMBER)a & r->nr2mModul) != (NATNUMBER)a) return FALSE;
     638  if (((NATNUMBER)a & r->mod2mMask) != (NATNUMBER)a) return FALSE;
    646639  return TRUE;
    647640}
    … …  
    664657      (*i) *= 10;
    665658      (*i) += *s++ - '0';
    666       if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & r->nr2mModul;
     659      if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & r->mod2mMask;
    667660    }
    668661    while (((*s) >= '0') && ((*s) <= '9'));
    669     (*i) = (*i) & r->nr2mModul;
     662    (*i) = (*i) & r->mod2mMask;
    670663  }
    671664  else (*i) = 1;

  • coeffs/rmodulo2m.h

    • Property mode changed from 100644 to 100755
    r5d594a9 re90dfd6  
    66/* $Id$ */
    77/*
    8 * ABSTRACT: numbers modulo 2^m
     8* ABSTRACT: numbers modulo 2^m such that 2^m - 1
     9*           fits in an unsigned long
    910*/
    1011#ifdef HAVE_RINGS
    … …  
    1617#endif
    1718
    18 
    1919BOOLEAN nr2mInitChar    (coeffs r, void*);
    20 
    21 extern int nr2mExp;
    22 extern NATNUMBER nr2mModul; /* for storing 2^m - 1, i.e., the
    23                                bit pattern '11..1' of length m */
    2420
    2521number  nr2mCopy        (number a, const coeffs r);
    … …  
    6258{
    6359  return (number)
    64     ((((NATNUMBER) a) * ((NATNUMBER) b)) & ((NATNUMBER) r->nr2mModul));
     60    ((((NATNUMBER) a) * ((NATNUMBER) b)) & ((NATNUMBER)r->mod2mMask));
    6561}
    6662
    … …  
    6864{
    6965  return (number)
    70     ((((NATNUMBER) a) + ((NATNUMBER) b)) & ((NATNUMBER) r->nr2mModul));
     66    ((((NATNUMBER) a) + ((NATNUMBER) b)) & ((NATNUMBER)r->mod2mMask));
    7167}
    7268
    7369static inline number nr2mSubM(number a, number b, const coeffs r)
    7470{
    75   return (number)((NATNUMBER)a<(NATNUMBER)b ?
    76                        r->nr2mModul-(NATNUMBER)b+(NATNUMBER)a + 1: (NATNUMBER)a-(NATNUMBER)b);
     71  return (number)((NATNUMBER)a < (NATNUMBER)b ?
     72                       r->mod2mMask - (NATNUMBER)b + (NATNUMBER)a + 1:
     73                       (NATNUMBER)a - (NATNUMBER)b);
    7774}
    7875
    79 #define nr2mNegM(A,r) (number)(r->nr2mModul-(NATNUMBER)(A) + 1)
     76#define nr2mNegM(A,r) (number)((r->mod2mMask - (NATNUMBER)(A) + 1) & r->mod2mMask)
    8077#define nr2mEqualM(A,B)  ((A)==(B))
    8178

  • coeffs/rmodulon.cc

    • Property mode changed from 100644 to 100755
    r5d594a9 re90dfd6  
    2626extern omBin gmp_nrz_bin;
    2727
    28 int_number nrnMinusOne = NULL;
    29 unsigned long nrnExponent = 0;
    30 
    3128/* for initializing function pointers */
    3229BOOLEAN nrnInitChar (coeffs r, void* p)
    … …  
    3431 
    3532  nrnInitExp((int)(long)(p), r);
    36  
    37      r->cfInit       = nrnInit;
    38      r->cfDelete     = nrnDelete;
    39      r->cfCopy       = nrnCopy;
    40      r->cfSize        = nrnSize;
    41      r->cfInt        = nrnInt;
    42      r->cfAdd         = nrnAdd;
    43      r->cfSub         = nrnSub;
    44      r->cfMult        = nrnMult;
    45      r->cfDiv         = nrnDiv;
    46      r->cfIntDiv      = nrnIntDiv;
    47      r->cfIntMod      = nrnMod;
    48      r->cfExactDiv    = nrnDiv;
    49      r->cfNeg         = nrnNeg;
    50      r->cfInvers      = nrnInvers;
    51      r->cfDivBy       = nrnDivBy;
    52      r->cfDivComp     = nrnDivComp;
    53      r->cfGreater     = nrnGreater;
    54      r->cfEqual       = nrnEqual;
    55      r->cfIsZero      = nrnIsZero;
    56      r->cfIsOne       = nrnIsOne;
    57      r->cfIsMOne      = nrnIsMOne;
    58      r->cfGreaterZero = nrnGreaterZero;
    59      r->cfWrite      = nrnWrite;
    60      r->cfRead        = nrnRead;
    61      r->cfPower       = nrnPower;
    62      r->cfSetMap     = nrnSetMap;
    63      r->cfNormalize   = ndNormalize;
    64      r->cfLcm         = nrnLcm;
    65      r->cfGcd         = nrnGcd;
    66      r->cfIsUnit      = nrnIsUnit;
    67      r->cfGetUnit     = nrnGetUnit;
    68      r->cfExtGcd      = nrnExtGcd;
    69      r->cfName        = ndName;
     33
     34  r->cfInit        = nrnInit;
     35  r->cfDelete      = nrnDelete;
     36  r->cfCopy        = nrnCopy;
     37  r->cfSize        = nrnSize;
     38  r->cfInt         = nrnInt;
     39  r->cfAdd         = nrnAdd;
     40  r->cfSub         = nrnSub;
     41  r->cfMult        = nrnMult;
     42  r->cfDiv         = nrnDiv;
     43  r->cfIntDiv      = nrnIntDiv;
     44  r->cfIntMod      = nrnMod;
     45  r->cfExactDiv    = nrnDiv;
     46  r->cfNeg         = nrnNeg;
     47  r->cfInvers      = nrnInvers;
     48  r->cfDivBy       = nrnDivBy;
     49  r->cfDivComp     = nrnDivComp;
     50  r->cfGreater     = nrnGreater;
     51  r->cfEqual       = nrnEqual;
     52  r->cfIsZero      = nrnIsZero;
     53  r->cfIsOne       = nrnIsOne;
     54  r->cfIsMOne      = nrnIsMOne;
     55  r->cfGreaterZero = nrnGreaterZero;
     56  r->cfWrite       = nrnWrite;
     57  r->cfRead        = nrnRead;
     58  r->cfPower       = nrnPower;
     59  r->cfSetMap      = nrnSetMap;
     60  r->cfNormalize   = ndNormalize;
     61  r->cfLcm         = nrnLcm;
     62  r->cfGcd         = nrnGcd;
     63  r->cfIsUnit      = nrnIsUnit;
     64  r->cfGetUnit     = nrnGetUnit;
     65  r->cfExtGcd      = nrnExtGcd;
     66  r->cfName        = ndName;
    7067#ifdef LDEBUG
    71      r->cfDBTest      = nrnDBTest;
     68  r->cfDBTest      = nrnDBTest;
    7269#endif
    7370  return FALSE;
    … …  
    7774 * create a number from int
    7875 */
    79 number nrnInit (int i, const coeffs r)
     76number nrnInit(int i, const coeffs r)
    8077{
    8178  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    8279  mpz_init_set_si(erg, i);
    83   mpz_mod(erg, erg, r->nrnModul);
     80  mpz_mod(erg, erg, r->modNumber);
    8481  return (number) erg;
    8582}
    … …  
    111108int nrnInt(number &n, const coeffs r)
    112109{
    113   return (int) mpz_get_si( (int_number) n);
     110  return (int)mpz_get_si((int_number) n);
    114111}
    115112
    … …  
    117114 * Multiply two numbers
    118115 */
    119 number nrnMult (number a, number b, const coeffs r)
    120 {
    121   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    122   mpz_init(erg);
    123   mpz_mul(erg, (int_number) a, (int_number) b);
    124   mpz_mod(erg, erg, r->nrnModul);
     116number nrnMult(number a, number b, const coeffs r)
     117{
     118  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     119  mpz_init(erg);
     120  mpz_mul(erg, (int_number)a, (int_number) b);
     121  mpz_mod(erg, erg, r->modNumber);
    125122  return (number) erg;
    126123}
    127124
    128 void nrnPower (number a, int i, number * result, const coeffs r)
    129 {
    130   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    131   mpz_init(erg);
    132   mpz_powm_ui(erg, (int_number) a, i, r->nrnModul);
     125void nrnPower(number a, int i, number * result, const coeffs r)
     126{
     127  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     128  mpz_init(erg);
     129  mpz_powm_ui(erg, (int_number)a, i, r->modNumber);
    133130  *result = (number) erg;
    134131}
    135132
    136 number nrnAdd (number a, number b, const coeffs r)
    137 {
    138   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    139   mpz_init(erg);
    140   mpz_add(erg, (int_number) a, (int_number) b);
    141   mpz_mod(erg, erg, r->nrnModul);
     133number nrnAdd(number a, number b, const coeffs r)
     134{
     135  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     136  mpz_init(erg);
     137  mpz_add(erg, (int_number)a, (int_number) b);
     138  mpz_mod(erg, erg, r->modNumber);
    142139  return (number) erg;
    143140}
    144141
    145 number nrnSub (number a, number b, const coeffs r)
    146 {
    147   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    148   mpz_init(erg);
    149   mpz_sub(erg, (int_number) a, (int_number) b);
    150   mpz_mod(erg, erg, r->nrnModul);
     142number nrnSub(number a, number b, const coeffs r)
     143{
     144  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     145  mpz_init(erg);
     146  mpz_sub(erg, (int_number)a, (int_number) b);
     147  mpz_mod(erg, erg, r->modNumber);
    151148  return (number) erg;
    152149}
    153150
    154 number nrnNeg (number c, const coeffs r)
    155 {
    156 // nNeg inplace !!!
    157   mpz_sub((int_number) c, r->nrnModul, (int_number) c);
     151number nrnNeg(number c, const coeffs r)
     152{
     153  // Attention: This method operates in-place.
     154  mpz_sub((int_number)c, r->modNumber, (int_number)c);
    158155  return c;
    159156}
    160157
    161 number  nrnInvers (number c, const coeffs r)
    162 {
    163   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    164   mpz_init(erg);
    165   mpz_invert(erg, (int_number) c, r->nrnModul);
     158number nrnInvers(number c, const coeffs r)
     159{
     160  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     161  mpz_init(erg);
     162  mpz_invert(erg, (int_number)c, r->modNumber);
    166163  return (number) erg;
    167164}
    168165
    169166/*
    170  * Give the smallest non unit k, such that a * x = k = b * y has a solution
    171  * TODO: lcm(gcd,gcd) besser als gcd(lcm) ?
    172  */
    173 number nrnLcm (number a, number b, const coeffs r)
     167 * Give the smallest k, such that a * x = k = b * y has a solution
     168 * TODO: lcm(gcd,gcd) better than gcd(lcm) ?
     169 */
     170number nrnLcm(number a, number b, const coeffs r)
    174171{
    175172  number erg = nrnGcd(NULL, a, r);
    176173  number tmp = nrnGcd(NULL, b, r);
    177   mpz_lcm((int_number) erg, (int_number) erg, (int_number) tmp);
     174  mpz_lcm((int_number)erg, (int_number)erg, (int_number)tmp);
    178175  nrnDelete(&tmp, NULL);
    179   return (number) erg;
    180 }
    181 
    182 /*
    183  * Give the largest non unit k, such that a = x * k, b = y * k has
     176  return (number)erg;
     177}
     178
     179/*
     180 * Give the largest k, such that a = x * k, b = y * k has
    184181 * a solution.
    185182 */
    186 number nrnGcd (number a, number b, const coeffs r)
     183number nrnGcd(number a, number b, const coeffs r)
    187184{
    188185  if ((a == NULL) && (b == NULL)) return nrnInit(0,r);
     186  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     187  mpz_init_set(erg, r->modNumber);
     188  if (a != NULL) mpz_gcd(erg, erg, (int_number)a);
     189  if (b != NULL) mpz_gcd(erg, erg, (int_number)b);
     190  return (number)erg;
     191}
     192
     193/* Not needed any more, but may have room for improvement
     194   number nrnGcd3(number a,number b, number c,ring r)
     195{
    189196  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    190   mpz_init_set(erg, r->nrnModul);
    191   if (a != NULL) mpz_gcd(erg, erg, (int_number) a);
    192   if (b != NULL) mpz_gcd(erg, erg, (int_number) b);
    193   return (number) erg;
    194 }
    195 
    196 /* Not needed any more, but may have room for improvement
    197 number nrnGcd3 (number a,number b, number c,ring r)
    198 {
    199   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    200   mpz_init(erg);
    201   if (a == NULL) a = (number) r->nrnModul;
    202   if (b == NULL) b = (number) r->nrnModul;
    203   if (c == NULL) c = (number) r->nrnModul;
    204   mpz_gcd(erg, (int_number) a, (int_number) b);
    205   mpz_gcd(erg, erg, (int_number) c);
    206   mpz_gcd(erg, erg, r->nrnModul);
    207   return (number) erg;
     197  mpz_init(erg);
     198  if (a == NULL) a = (number)r->modNumber;
     199  if (b == NULL) b = (number)r->modNumber;
     200  if (c == NULL) c = (number)r->modNumber;
     201  mpz_gcd(erg, (int_number)a, (int_number)b);
     202  mpz_gcd(erg, erg, (int_number)c);
     203  mpz_gcd(erg, erg, r->modNumber);
     204  return (number)erg;
    208205}
    209206*/
    210207
    211208/*
    212  * Give the largest non unit k, such that a = x * k, b = y * k has
     209 * Give the largest k, such that a = x * k, b = y * k has
    213210 * a solution and r, s, s.t. k = s*a + t*b
    214211 */
    215 number  nrnExtGcd (number a, number b, number *s, number *t, const coeffs r)
    216 {
    217   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    218   int_number bs = (int_number) omAllocBin(gmp_nrz_bin);
    219   int_number bt = (int_number) omAllocBin(gmp_nrz_bin);
     212number nrnExtGcd(number a, number b, number *s, number *t, const coeffs r)
     213{
     214  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     215  int_number bs  = (int_number)omAllocBin(gmp_nrz_bin);
     216  int_number bt  = (int_number)omAllocBin(gmp_nrz_bin);
    220217  mpz_init(erg);
    221218  mpz_init(bs);
    222219  mpz_init(bt);
    223   mpz_gcdext(erg, bs, bt, (int_number) a, (int_number) b);
    224   mpz_mod(bs, bs, r->nrnModul);
    225   mpz_mod(bt, bt, r->nrnModul);
    226   *s = (number) bs;
    227   *t = (number) bt;
    228   return (number) erg;
    229 }
    230 
    231 BOOLEAN nrnIsZero (number a, const coeffs r)
     220  mpz_gcdext(erg, bs, bt, (int_number)a, (int_number)b);
     221  mpz_mod(bs, bs, r->modNumber);
     222  mpz_mod(bt, bt, r->modNumber);
     223  *s = (number)bs;
     224  *t = (number)bt;
     225  return (number)erg;
     226}
     227
     228BOOLEAN nrnIsZero(number a, const coeffs r)
    232229{
    233230#ifdef LDEBUG
    234231  if (a == NULL) return FALSE;
    235232#endif
    236   return 0 == mpz_cmpabs_ui((int_number) a, 0);
    237 }
    238 
    239 BOOLEAN nrnIsOne (number a, const coeffs r)
     233  return 0 == mpz_cmpabs_ui((int_number)a, 0);
     234}
     235
     236BOOLEAN nrnIsOne(number a, const coeffs r)
    240237{
    241238#ifdef LDEBUG
    242239  if (a == NULL) return FALSE;
    243240#endif
    244   return 0 == mpz_cmp_si((int_number) a, 1);
    245 }
    246 
    247 BOOLEAN nrnIsMOne (number a, const coeffs r)
     241  return 0 == mpz_cmp_si((int_number)a, 1);
     242}
     243
     244BOOLEAN nrnIsMOne(number a, const coeffs r)
    248245{
    249246#ifdef LDEBUG
    250247  if (a == NULL) return FALSE;
    251248#endif
    252   return 0 == mpz_cmp((int_number) a, nrnMinusOne);
    253 }
    254 
    255 BOOLEAN nrnEqual (number a, number b, const coeffs r)
    256 {
    257   return 0 == mpz_cmp((int_number) a, (int_number) b);
    258 }
    259 
    260 BOOLEAN nrnGreater (number a, number b, const coeffs r)
    261 {
    262   return 0 < mpz_cmp((int_number) a, (int_number) b);
    263 }
    264 
    265 BOOLEAN nrnGreaterZero (number k, const coeffs r)
    266 {
    267   return 0 < mpz_cmp_si((int_number) k, 0);
    268 }
    269 
    270 BOOLEAN nrnIsUnit (number a, const coeffs r)
    271 {
    272   number tmp = nrnGcd(a, (number) r->nrnModul, r);
    273   bool res = nrnIsOne(tmp,r);
     249  mpz_t t; mpz_init_set(t, (int_number)a);
     250  mpz_add_ui(t, t, 1);
     251  bool erg = (0 == mpz_cmp(t, r->modNumber));
     252  mpz_clear(t);
     253  return erg;
     254}
     255
     256BOOLEAN nrnEqual(number a, number b, const coeffs r)
     257{
     258  return 0 == mpz_cmp((int_number)a, (int_number)b);
     259}
     260
     261BOOLEAN nrnGreater(number a, number b, const coeffs r)
     262{
     263  return 0 < mpz_cmp((int_number)a, (int_number)b);
     264}
     265
     266BOOLEAN nrnGreaterZero(number k, const coeffs r)
     267{
     268  return 0 < mpz_cmp_si((int_number)k, 0);
     269}
     270
     271BOOLEAN nrnIsUnit(number a, const coeffs r)
     272{
     273  number tmp = nrnGcd(a, (number)r->modNumber, r);
     274  bool res = nrnIsOne(tmp, r);
    274275  nrnDelete(&tmp, NULL);
    275276  return res;
    276277}
    277278
    278 number  nrnGetUnit (number k, const coeffs r)
    279 {
    280   if (mpz_divisible_p(r->nrnModul, (int_number) k)) return nrnInit(1,r);
    281 
    282   int_number unit = (int_number) nrnGcd(k, 0, r);
    283   mpz_tdiv_q(unit, (int_number) k, unit);
    284   int_number gcd = (int_number) nrnGcd((number) unit, 0, r);
    285   if (!nrnIsOne((number) gcd,r))
     279number nrnGetUnit(number k, const coeffs r)
     280{
     281  if (mpz_divisible_p(r->modNumber, (int_number)k)) return nrnInit(1,r);
     282
     283  int_number unit = (int_number)nrnGcd(k, 0, r);
     284  mpz_tdiv_q(unit, (int_number)k, unit);
     285  int_number gcd = (int_number)nrnGcd((number)unit, 0, r);
     286  if (!nrnIsOne((number)gcd,r))
    286287  {
    287288    int_number ctmp;
    … …  
    298299      // tmp := tmp * unit
    299300      mpz_mul(tmp, tmp, unit);
    300       mpz_mod(tmp, tmp, r->nrnModul);
     301      mpz_mod(tmp, tmp, r->modNumber);
    301302      // gcd_new := gcd(tmp, 0)
    302       mpz_gcd(gcd_new, tmp, r->nrnModul);
     303      mpz_gcd(gcd_new, tmp, r->modNumber);
    303304    }
    304     // unit := unit + nrnModul / gcd_new
    305     mpz_tdiv_q(tmp, r->nrnModul, gcd_new);
     305    // unit := unit + modNumber / gcd_new
     306    mpz_tdiv_q(tmp, r->modNumber, gcd_new);
    306307    mpz_add(unit, unit, tmp);
    307     mpz_mod(unit, unit, r->nrnModul);
     308    mpz_mod(unit, unit, r->modNumber);
    308309    nrnDelete((number*) &gcd_new, NULL);
    309310    nrnDelete((number*) &tmp, NULL);
    310311  }
    311312  nrnDelete((number*) &gcd, NULL);
    312   return (number) unit;
    313 }
    314 
    315 BOOLEAN nrnDivBy (number a, number b, const coeffs r)
     313  return (number)unit;
     314}
     315
     316BOOLEAN nrnDivBy(number a, number b, const coeffs r)
    316317{
    317318  if (a == NULL)
    318     return mpz_divisible_p(r->nrnModul, (int_number) b);
     319    return mpz_divisible_p(r->modNumber, (int_number)b);
    319320  else
    320321  { /* b divides a iff b/gcd(a, b) is a unit in the given ring: */
    … …  
    332333  if (mpz_divisible_p((int_number) a, (int_number) b)) return -1;
    333334  if (mpz_divisible_p((int_number) b, (int_number) a)) return 1;
    334   return 2;
    335 }
    336 
    337 number nrnDiv (number a, number b, const coeffs r)
    338 {
    339   if (a == NULL) a = (number) r->nrnModul;
    340   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    341   mpz_init(erg);
    342   if (mpz_divisible_p((int_number) a, (int_number) b))
    343   {
    344     mpz_divexact(erg, (int_number) a, (int_number) b);
    345     return (number) erg;
     335  return 0;
     336}
     337
     338number nrnDiv(number a, number b, const coeffs r)
     339{
     340  if (a == NULL) a = (number)r->modNumber;
     341  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     342  mpz_init(erg);
     343  if (mpz_divisible_p((int_number)a, (int_number)b))
     344  {
     345    mpz_divexact(erg, (int_number)a, (int_number)b);
     346    return (number)erg;
    346347  }
    347348  else
    348349  {
    349     int_number gcd = (int_number) nrnGcd(a, b, r);
    350     mpz_divexact(erg, (int_number) b, gcd);
    351     if (!nrnIsUnit((number) erg,r))
     350    int_number gcd = (int_number)nrnGcd(a, b, r);
     351    mpz_divexact(erg, (int_number)b, gcd);
     352    if (!nrnIsUnit((number)erg, r))
    352353    {
    353354      WerrorS("Division not possible, even by cancelling zero divisors.");
    … …  
    355356      mpz_tdiv_q(erg, (int_number) a, (int_number) b);
    356357      nrnDelete((number*) &gcd, NULL);
    357       return (number) erg;
     358      return (number)erg;
    358359    }
    359360    // a / gcd(a,b) * [b / gcd (a,b)]^(-1)
    360     int_number tmp = (int_number) nrnInvers((number) erg,r);
    361     mpz_divexact(erg, (int_number) a, gcd);
     361    int_number tmp = (int_number)nrnInvers((number) erg,r);
     362    mpz_divexact(erg, (int_number)a, gcd);
    362363    mpz_mul(erg, erg, tmp);
    363364    nrnDelete((number*) &gcd, NULL);
    364365    nrnDelete((number*) &tmp, NULL);
    365     mpz_mod(erg, erg, r->nrnModul);
    366     return (number) erg;
    367   }
    368 }
    369 
    370 number nrnMod (number a, number b, const coeffs R)
     366    mpz_mod(erg, erg, r->modNumber);
     367    return (number)erg;
     368  }
     369}
     370
     371number nrnMod(number a, number b, const coeffs r)
    371372{
    372373  /*
    373     We need to return the number r which is uniquely determined by the
     374    We need to return the number rr which is uniquely determined by the
    374375    following two properties:
    375       (1) 0 <= r < |b| (with respect to '<' and '<=' performed in Z x Z)
    376       (2) There exists some k in the integers Z such that a = k * b + r.
     376      (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z)
     377      (2) There exists some k in the integers Z such that a = k * b + rr.
    377378    Consider g := gcd(n, |b|). Note that then |b|/g is a unit in Z/n.
    378379    Now, there are three cases:
    379380      (a) g = 1
    380381          Then |b| is a unit in Z/n, i.e. |b| (and also b) divides a.
    381           Thus r = 0.
     382          Thus rr = 0.
    382383      (b) g <> 1 and g divides a
    383           Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again r = 0.
     384          Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0.
    384385      (c) g <> 1 and g does not divide a
    385386          Then denote the division with remainder of a by g as this:
    386387          a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b|
    387           fulfills (1) and (2), i.e. r := t is the correct result. Hence
    388           in this third case, r is the remainder of division of a by g in Z.
     388          fulfills (1) and (2), i.e. rr := t is the correct result. Hence
     389          in this third case, rr is the remainder of division of a by g in Z.
    389390     Remark: according to mpz_mod: a,b are always non-negative
    390391  */
    391   int_number g = (int_number) omAllocBin(gmp_nrz_bin);
    392   int_number r = (int_number) omAllocBin(gmp_nrz_bin);
     392  int_number g = (int_number)omAllocBin(gmp_nrz_bin);
     393  int_number rr = (int_number)omAllocBin(gmp_nrz_bin);
    393394  mpz_init(g);
    394   mpz_init_set_si(r,(long)0);
    395   mpz_gcd(g, (int_number) R->nrnModul, (int_number)b); // g is now as above
    396   if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(r, (int_number)a, g); // the case g <> 1
     395  mpz_init_set_si(rr, 0);
     396  mpz_gcd(g, (int_number)r->modNumber, (int_number)b); // g is now as above
     397  if (mpz_cmp_si(g, (long)1) != 0) mpz_mod(rr, (int_number)a, g); // the case g <> 1
    397398  mpz_clear(g);
    398399  omFreeBin(g, gmp_nrz_bin);
    399   return (number)r;
    400 }
    401 
    402 number nrnIntDiv (number a, number b, const coeffs r)
    403 {
    404   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    405   mpz_init(erg);
    406   if (a == NULL) a = (number) r->nrnModul;
    407   mpz_tdiv_q(erg, (int_number) a, (int_number) b);
    408   return (number) erg;
     400  return (number)rr;
     401}
     402
     403number nrnIntDiv(number a, number b, const coeffs r)
     404{
     405  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     406  mpz_init(erg);
     407  if (a == NULL) a = (number)r->modNumber;
     408  mpz_tdiv_q(erg, (int_number)a, (int_number)b);
     409  return (number)erg;
    409410}
    410411
    … …  
    422423number nrnMap2toM(number from, const coeffs src, const coeffs dst)
    423424{
    424   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    425   mpz_init(erg);
    426   mpz_mul_ui(erg, nrnMapCoef, (NATNUMBER) from);
    427   mpz_mod(erg, erg, dst->nrnModul);
    428   return (number) erg;
     425  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     426  mpz_init(erg);
     427  mpz_mul_ui(erg, nrnMapCoef, (NATNUMBER)from);
     428  mpz_mod(erg, erg, dst->modNumber);
     429  return (number)erg;
    429430}
    430431
    431432number nrnMapZp(number from, const coeffs src, const coeffs dst)
    432433{
    433   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
     434  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
    434435  mpz_init(erg);
    435436  // TODO: use npInt(...)
    436   mpz_mul_si(erg, nrnMapCoef, (NATNUMBER) from);
    437   mpz_mod(erg, erg, dst->nrnModul);
    438   return (number) erg;
     437  mpz_mul_si(erg, nrnMapCoef, (NATNUMBER)from);
     438  mpz_mod(erg, erg, dst->modNumber);
     439  return (number)erg;
    439440}
    440441
    441442number nrnMapGMP(number from, const coeffs src, const coeffs dst)
    442443{
    443   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    444   mpz_init(erg);
    445   mpz_mod(erg, (int_number) from, dst->nrnModul);
    446   return (number) erg;
     444  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     445  mpz_init(erg);
     446  mpz_mod(erg, (int_number)from, dst->modNumber);
     447  return (number)erg;
    447448}
    448449
    449450number nrnMapQ(number from, const coeffs src, const coeffs dst)
    450451{
    451   int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    452   mpz_init(erg);
    453   nlGMP(from, (number) erg, src);
    454   mpz_mod(erg, erg, src->nrnModul);
    455   return (number) erg;
     452  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     453  mpz_init(erg);
     454  nlGMP(from, (number)erg, src);
     455  mpz_mod(erg, erg, src->modNumber);
     456  return (number)erg;
    456457}
    457458
    … …  
    471472  {
    472473    if (   (src->ringtype > 0)
    473         && (mpz_cmp(src->ringflaga, dst->ringflaga) == 0)
    474         && (src->ringflagb == dst->ringflagb)) return nrnCopy;
     474        && (mpz_cmp(src->modBase, dst->modBase) == 0)
     475        && (src->modExponent == dst->modExponent)) return nrnMapGMP;
    475476    else
    476477    {
    … …  
    484485      {
    485486        mpz_init(nrnMapModul);
    486         mpz_set(nrnMapModul, src->ringflaga);
    487         mpz_pow_ui(nrnMapModul, nrnMapModul, src->ringflagb);
     487        mpz_set(nrnMapModul, src->modNumber);
    488488      }
    489489      // nrnMapCoef = 1 in dst       if dst is a subring of src
    … …  
    494494        mpz_init(nrnMapCoef);
    495495      }
    496       if (mpz_divisible_p(nrnMapModul, dst->nrnModul))
     496      if (mpz_divisible_p(nrnMapModul, dst->modNumber))
    497497      {
    498498        mpz_set_si(nrnMapCoef, 1);
    … …  
    501501      if (nrnDivBy(NULL, (number) nrnMapModul,dst))
    502502      {
    503         mpz_divexact(nrnMapCoef, dst->nrnModul, nrnMapModul);
    504         int_number tmp = dst->nrnModul;
    505         dst->nrnModul = nrnMapModul;
     503        mpz_divexact(nrnMapCoef, dst->modNumber, nrnMapModul);
     504        int_number tmp = dst->modNumber;
     505        dst->modNumber = nrnMapModul;
    506506        if (!nrnIsUnit((number) nrnMapCoef,dst))
    507507        {
    508           dst->nrnModul = tmp;
     508          dst->modNumber = tmp;
    509509          nrnDelete((number*) &nrnMapModul, dst);
    510510          return NULL;
    511511        }
    512512        int_number inv = (int_number) nrnInvers((number) nrnMapCoef,dst);
    513         dst->nrnModul = tmp;
     513        dst->modNumber = tmp;
    514514        mpz_mul(nrnMapCoef, nrnMapCoef, inv);
    515         mpz_mod(nrnMapCoef, nrnMapCoef, dst->nrnModul);
     515        mpz_mod(nrnMapCoef, nrnMapCoef, dst->modNumber);
    516516        nrnDelete((number*) &inv, dst);
    517517      }
    … …  
    539539void nrnSetExp(int m, coeffs r)
    540540{
    541   if ((r->nrnModul != NULL) && (mpz_cmp(r->nrnModul, r->ringflaga) == 0) && (nrnExponent == r->ringflagb)) return;
    542 
    543   nrnExponent = r->ringflagb;
    544   if (r->nrnModul == NULL)
    545   {
    546     r->nrnModul = (int_number) omAllocBin(gmp_nrz_bin);
    547     mpz_init(r->nrnModul);
    548     nrnMinusOne = (int_number) omAllocBin(gmp_nrz_bin);
    549     mpz_init(nrnMinusOne);
    550   }
    551   // BUG:  r->ringflaga is undefined!
    552   mpz_set(r->nrnModul, r->ringflaga);
    553   mpz_pow_ui(r->nrnModul, r->nrnModul, nrnExponent);
    554   mpz_sub_ui(nrnMinusOne, r->nrnModul, 1);
    555 }
    556 
     541  /* clean up former stuff */
     542  if (r->modBase   != NULL) mpz_clear(r->modBase);
     543  if (r->modNumber != NULL) mpz_clear(r->modNumber);
     544
     545  /* this is Z/m = Z/(m^1), hence set modBase = m, modExponent = 1: */
     546  r->modBase = (int_number)omAllocBin(gmp_nrz_bin);
     547  mpz_init(r->modBase);
     548  mpz_set_ui(r->modBase, (unsigned long)m);
     549  r->modExponent = 1;
     550  r->modNumber = (int_number)omAllocBin(gmp_nrz_bin);
     551  mpz_init(r->modNumber);
     552  mpz_set(r->modNumber, r->modBase);
     553  /* mpz_pow_ui(r->modNumber, r->modNumber, r->modExponent); */
     554}
     555
     556/* We expect this ring to be Z/m for some m > 2 which is not a prime. */
    557557void nrnInitExp(int m, coeffs r)
    558558{
     559  if (m <= 2) WarnS("nrnInitExp failed (m in Z/m too small)");
    559560  nrnSetExp(m, r);
    560 
    561   if (mpz_cmp_ui(r->nrnModul, 2) <= 0) // ???
    562   {
    563     WarnS("nrnInitExp failed");
    564   }
    565  
    566   r->ch = m; // ???
    567561}
    568562
    … …  
    571565{
    572566  if (a==NULL) return TRUE;
    573   if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, r->nrnModul) > 0) )
     567  if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, r->modNumber) > 0) )
    574568  {
    575569    return FALSE;
    … …  
    612606    s = nlCPEatLongC((char *)s, z);
    613607  }
    614   mpz_mod(z, z, r->nrnModul);
     608  mpz_mod(z, z, r->modNumber);
    615609  *a = (number) z;
    616610  return s;
Note: See TracChangeset for help on using the changeset viewer.