[35b1d7] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
| 4 | /* |
---|
| 5 | * ABSTRACT: numbers modulo 2^m |
---|
| 6 | */ |
---|
| 7 | |
---|
[691dba] | 8 | #include "config.h" |
---|
[18cb65] | 9 | #include <misc/auxiliary.h> |
---|
[2a329d] | 10 | |
---|
[c90b43] | 11 | #ifdef HAVE_RINGS |
---|
[f1c465f] | 12 | |
---|
[18cb65] | 13 | #include <misc/mylimits.h> |
---|
[2d805a] | 14 | #include <coeffs/coeffs.h> |
---|
[31213a4] | 15 | #include <reporter/reporter.h> |
---|
| 16 | #include <omalloc/omalloc.h> |
---|
[2d805a] | 17 | #include <coeffs/numbers.h> |
---|
| 18 | #include <coeffs/longrat.h> |
---|
| 19 | #include <coeffs/mpr_complex.h> |
---|
| 20 | #include <coeffs/rmodulo2m.h> |
---|
[e3b233] | 21 | #include "si_gmp.h" |
---|
[894f5b1] | 22 | |
---|
[f1c465f] | 23 | #include <string.h> |
---|
| 24 | |
---|
[73a9ffb] | 25 | /// Our Type! |
---|
| 26 | static const n_coeffType ID = n_Z2m; |
---|
| 27 | |
---|
[925a43c] | 28 | extern omBin gmp_nrz_bin; /* init in rintegers*/ |
---|
| 29 | |
---|
[03f7b5] | 30 | void nr2mCoeffWrite (const coeffs r, BOOLEAN /*details*/) |
---|
[7a8011] | 31 | { |
---|
[613794] | 32 | Print("// Z/2^%lu\n", r->modExponent); |
---|
[7a8011] | 33 | } |
---|
| 34 | |
---|
[613794] | 35 | BOOLEAN nr2mCoeffIsEqual(const coeffs r, n_coeffType n, void * p) |
---|
| 36 | { |
---|
| 37 | if (n==n_Z2m) |
---|
| 38 | { |
---|
| 39 | int m=(int)(long)(p); |
---|
| 40 | unsigned long mm=r->mod2mMask; |
---|
| 41 | if ((mm>>m)==1L) return TRUE; |
---|
| 42 | } |
---|
| 43 | return FALSE; |
---|
| 44 | } |
---|
[14b11bb] | 45 | /* for initializing function pointers */ |
---|
[1cce47] | 46 | BOOLEAN nr2mInitChar (coeffs r, void* p) |
---|
[14b11bb] | 47 | { |
---|
[73a9ffb] | 48 | assume( getCoeffType(r) == ID ); |
---|
[1112b76] | 49 | nr2mInitExp((int)(long)(p), r); |
---|
[613794] | 50 | r->cfKillChar = ndKillChar; /* dummy*/ |
---|
| 51 | r->nCoeffIsEqual = nr2mCoeffIsEqual; |
---|
[1112b76] | 52 | |
---|
[f0797c] | 53 | r->ringtype = 1; |
---|
[9bb5457] | 54 | |
---|
[73a9ffb] | 55 | /* next cast may yield an overflow as mod2mMask is an unsigned long */ |
---|
| 56 | r->ch = (int)r->mod2mMask + 1; |
---|
[f0797c] | 57 | |
---|
[e90dfd6] | 58 | r->cfInit = nr2mInit; |
---|
| 59 | r->cfCopy = ndCopy; |
---|
| 60 | r->cfInt = nr2mInt; |
---|
[5d594a9] | 61 | r->cfAdd = nr2mAdd; |
---|
| 62 | r->cfSub = nr2mSub; |
---|
| 63 | r->cfMult = nr2mMult; |
---|
| 64 | r->cfDiv = nr2mDiv; |
---|
| 65 | r->cfIntDiv = nr2mIntDiv; |
---|
| 66 | r->cfIntMod = nr2mMod; |
---|
| 67 | r->cfExactDiv = nr2mDiv; |
---|
| 68 | r->cfNeg = nr2mNeg; |
---|
| 69 | r->cfInvers = nr2mInvers; |
---|
| 70 | r->cfDivBy = nr2mDivBy; |
---|
| 71 | r->cfDivComp = nr2mDivComp; |
---|
| 72 | r->cfGreater = nr2mGreater; |
---|
| 73 | r->cfEqual = nr2mEqual; |
---|
| 74 | r->cfIsZero = nr2mIsZero; |
---|
| 75 | r->cfIsOne = nr2mIsOne; |
---|
| 76 | r->cfIsMOne = nr2mIsMOne; |
---|
| 77 | r->cfGreaterZero = nr2mGreaterZero; |
---|
[ce1f78] | 78 | r->cfWriteLong = nr2mWrite; |
---|
[5d594a9] | 79 | r->cfRead = nr2mRead; |
---|
| 80 | r->cfPower = nr2mPower; |
---|
[e90dfd6] | 81 | r->cfSetMap = nr2mSetMap; |
---|
[5d594a9] | 82 | r->cfNormalize = ndNormalize; |
---|
| 83 | r->cfLcm = nr2mLcm; |
---|
| 84 | r->cfGcd = nr2mGcd; |
---|
| 85 | r->cfIsUnit = nr2mIsUnit; |
---|
| 86 | r->cfGetUnit = nr2mGetUnit; |
---|
| 87 | r->cfExtGcd = nr2mExtGcd; |
---|
| 88 | r->cfName = ndName; |
---|
[7a8011] | 89 | r->cfCoeffWrite = nr2mCoeffWrite; |
---|
[14b11bb] | 90 | #ifdef LDEBUG |
---|
[5d594a9] | 91 | r->cfDBTest = nr2mDBTest; |
---|
[14b11bb] | 92 | #endif |
---|
[5d594a9] | 93 | r->has_simple_Alloc=TRUE; |
---|
| 94 | return FALSE; |
---|
[14b11bb] | 95 | } |
---|
| 96 | |
---|
[35b1d7] | 97 | /* |
---|
| 98 | * Multiply two numbers |
---|
| 99 | */ |
---|
[e90dfd6] | 100 | number nr2mMult(number a, number b, const coeffs r) |
---|
[35b1d7] | 101 | { |
---|
[994445] | 102 | if (((NATNUMBER)a == 0) || ((NATNUMBER)b == 0)) |
---|
[35b1d7] | 103 | return (number)0; |
---|
| 104 | else |
---|
[e90dfd6] | 105 | return nr2mMultM(a, b, r); |
---|
[35b1d7] | 106 | } |
---|
| 107 | |
---|
[f92547] | 108 | /* |
---|
[e90dfd6] | 109 | * Give the smallest k, such that a * x = k = b * y has a solution |
---|
[f92547] | 110 | */ |
---|
[9bb5457] | 111 | number nr2mLcm(number a, number b, const coeffs) |
---|
[f92547] | 112 | { |
---|
[994445] | 113 | NATNUMBER res = 0; |
---|
[e90dfd6] | 114 | if ((NATNUMBER)a == 0) a = (number) 1; |
---|
| 115 | if ((NATNUMBER)b == 0) b = (number) 1; |
---|
| 116 | while ((NATNUMBER)a % 2 == 0) |
---|
[f92547] | 117 | { |
---|
[e90dfd6] | 118 | a = (number)((NATNUMBER)a / 2); |
---|
| 119 | if ((NATNUMBER)b % 2 == 0) b = (number)((NATNUMBER)b / 2); |
---|
[f92547] | 120 | res++; |
---|
| 121 | } |
---|
[e90dfd6] | 122 | while ((NATNUMBER)b % 2 == 0) |
---|
[f92547] | 123 | { |
---|
[e90dfd6] | 124 | b = (number)((NATNUMBER)b / 2); |
---|
[f92547] | 125 | res++; |
---|
| 126 | } |
---|
[e90dfd6] | 127 | return (number)(1L << res); // (2**res) |
---|
[f92547] | 128 | } |
---|
| 129 | |
---|
| 130 | /* |
---|
[e90dfd6] | 131 | * Give the largest k, such that a = x * k, b = y * k has |
---|
[b429c16] | 132 | * a solution. |
---|
[f92547] | 133 | */ |
---|
[9bb5457] | 134 | number nr2mGcd(number a, number b, const coeffs) |
---|
[f92547] | 135 | { |
---|
[994445] | 136 | NATNUMBER res = 0; |
---|
[e90dfd6] | 137 | if ((NATNUMBER)a == 0 && (NATNUMBER)b == 0) return (number)1; |
---|
| 138 | while ((NATNUMBER)a % 2 == 0 && (NATNUMBER)b % 2 == 0) |
---|
[f92547] | 139 | { |
---|
[e90dfd6] | 140 | a = (number)((NATNUMBER)a / 2); |
---|
| 141 | b = (number)((NATNUMBER)b / 2); |
---|
[f92547] | 142 | res++; |
---|
| 143 | } |
---|
[e90dfd6] | 144 | // if ((NATNUMBER)b % 2 == 0) |
---|
[206e158] | 145 | // { |
---|
[e90dfd6] | 146 | // return (number)((1L << res)); // * (NATNUMBER) a); // (2**res)*a a is a unit |
---|
[206e158] | 147 | // } |
---|
| 148 | // else |
---|
| 149 | // { |
---|
[e90dfd6] | 150 | return (number)((1L << res)); // * (NATNUMBER) b); // (2**res)*b b is a unit |
---|
[206e158] | 151 | // } |
---|
[f92547] | 152 | } |
---|
| 153 | |
---|
[1e579c6] | 154 | /* |
---|
[e90dfd6] | 155 | * Give the largest k, such that a = x * k, b = y * k has |
---|
[1e579c6] | 156 | * a solution. |
---|
| 157 | */ |
---|
[e90dfd6] | 158 | number nr2mExtGcd(number a, number b, number *s, number *t, const coeffs r) |
---|
[1e579c6] | 159 | { |
---|
| 160 | NATNUMBER res = 0; |
---|
[e90dfd6] | 161 | if ((NATNUMBER)a == 0 && (NATNUMBER)b == 0) return (number)1; |
---|
| 162 | while ((NATNUMBER)a % 2 == 0 && (NATNUMBER)b % 2 == 0) |
---|
[1e579c6] | 163 | { |
---|
[e90dfd6] | 164 | a = (number)((NATNUMBER)a / 2); |
---|
| 165 | b = (number)((NATNUMBER)b / 2); |
---|
[1e579c6] | 166 | res++; |
---|
| 167 | } |
---|
[e90dfd6] | 168 | if ((NATNUMBER)b % 2 == 0) |
---|
[1e579c6] | 169 | { |
---|
| 170 | *t = NULL; |
---|
[925a43c] | 171 | *s = nr2mInvers(a,r); |
---|
[e90dfd6] | 172 | return (number)((1L << res)); // * (NATNUMBER) a); // (2**res)*a a is a unit |
---|
[1e579c6] | 173 | } |
---|
| 174 | else |
---|
| 175 | { |
---|
| 176 | *s = NULL; |
---|
[925a43c] | 177 | *t = nr2mInvers(b,r); |
---|
[e90dfd6] | 178 | return (number)((1L << res)); // * (NATNUMBER) b); // (2**res)*b b is a unit |
---|
[1e579c6] | 179 | } |
---|
| 180 | } |
---|
| 181 | |
---|
[e90dfd6] | 182 | void nr2mPower(number a, int i, number * result, const coeffs r) |
---|
[35b1d7] | 183 | { |
---|
[e90dfd6] | 184 | if (i == 0) |
---|
[35b1d7] | 185 | { |
---|
[994445] | 186 | *(NATNUMBER *)result = 1; |
---|
[35b1d7] | 187 | } |
---|
[e90dfd6] | 188 | else if (i == 1) |
---|
[35b1d7] | 189 | { |
---|
| 190 | *result = a; |
---|
| 191 | } |
---|
| 192 | else |
---|
| 193 | { |
---|
[e90dfd6] | 194 | nr2mPower(a, i-1, result, r); |
---|
| 195 | *result = nr2mMultM(a, *result, r); |
---|
[35b1d7] | 196 | } |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | /* |
---|
| 200 | * create a number from int |
---|
| 201 | */ |
---|
[2f3764] | 202 | number nr2mInit(long i, const coeffs r) |
---|
[35b1d7] | 203 | { |
---|
[76e501] | 204 | if (i == 0) return (number)(NATNUMBER)i; |
---|
| 205 | |
---|
[d3b2eb] | 206 | long ii = i; |
---|
[76e501] | 207 | NATNUMBER j = (NATNUMBER)1; |
---|
[e90dfd6] | 208 | if (ii < 0) { j = r->mod2mMask; ii = -ii; } |
---|
[76e501] | 209 | NATNUMBER k = (NATNUMBER)ii; |
---|
[e90dfd6] | 210 | k = k & r->mod2mMask; |
---|
| 211 | /* now we have: i = j * k mod 2^m */ |
---|
| 212 | return (number)nr2mMult((number)j, (number)k, r); |
---|
[35b1d7] | 213 | } |
---|
| 214 | |
---|
| 215 | /* |
---|
[76e501] | 216 | * convert a number to an int in ]-k/2 .. k/2], |
---|
| 217 | * where k = 2^m; i.e., an int in ]-2^(m-1) .. 2^(m-1)]; |
---|
| 218 | * note that the code computes a long which will then |
---|
| 219 | * automatically casted to int |
---|
[35b1d7] | 220 | */ |
---|
[aa2bcca] | 221 | static long nr2mLong(number &n, const coeffs r) |
---|
[35b1d7] | 222 | { |
---|
[e90dfd6] | 223 | NATNUMBER nn = (unsigned long)(NATNUMBER)n & r->mod2mMask; |
---|
| 224 | unsigned long l = r->mod2mMask >> 1; l++; /* now: l = 2^(m-1) */ |
---|
[4fb4f3] | 225 | if ((NATNUMBER)nn > l) |
---|
[aa2bcca] | 226 | return (long)((NATNUMBER)nn - r->mod2mMask - 1); |
---|
[76e501] | 227 | else |
---|
[aa2bcca] | 228 | return (long)((NATNUMBER)nn); |
---|
| 229 | } |
---|
| 230 | int nr2mInt(number &n, const coeffs r) |
---|
| 231 | { |
---|
| 232 | return (int)nr2mLong(n,r); |
---|
[35b1d7] | 233 | } |
---|
| 234 | |
---|
[e90dfd6] | 235 | number nr2mAdd(number a, number b, const coeffs r) |
---|
[35b1d7] | 236 | { |
---|
[e90dfd6] | 237 | return nr2mAddM(a, b, r); |
---|
[35b1d7] | 238 | } |
---|
| 239 | |
---|
[e90dfd6] | 240 | number nr2mSub(number a, number b, const coeffs r) |
---|
[35b1d7] | 241 | { |
---|
[e90dfd6] | 242 | return nr2mSubM(a, b, r); |
---|
[35b1d7] | 243 | } |
---|
| 244 | |
---|
[9bb5457] | 245 | BOOLEAN nr2mIsUnit(number a, const coeffs) |
---|
[1e579c6] | 246 | { |
---|
[e90dfd6] | 247 | return ((NATNUMBER)a % 2 == 1); |
---|
[1e579c6] | 248 | } |
---|
| 249 | |
---|
[9bb5457] | 250 | number nr2mGetUnit(number k, const coeffs) |
---|
[1e579c6] | 251 | { |
---|
[e90dfd6] | 252 | if (k == NULL) return (number)1; |
---|
| 253 | NATNUMBER erg = (NATNUMBER)k; |
---|
| 254 | while (erg % 2 == 0) erg = erg / 2; |
---|
| 255 | return (number)erg; |
---|
[1e579c6] | 256 | } |
---|
| 257 | |
---|
[9bb5457] | 258 | BOOLEAN nr2mIsZero(number a, const coeffs) |
---|
[35b1d7] | 259 | { |
---|
[994445] | 260 | return 0 == (NATNUMBER)a; |
---|
[35b1d7] | 261 | } |
---|
| 262 | |
---|
[9bb5457] | 263 | BOOLEAN nr2mIsOne(number a, const coeffs) |
---|
[35b1d7] | 264 | { |
---|
[994445] | 265 | return 1 == (NATNUMBER)a; |
---|
[35b1d7] | 266 | } |
---|
| 267 | |
---|
[e90dfd6] | 268 | BOOLEAN nr2mIsMOne(number a, const coeffs r) |
---|
[35b1d7] | 269 | { |
---|
[e90dfd6] | 270 | return (r->mod2mMask == (NATNUMBER)a); |
---|
[35b1d7] | 271 | } |
---|
| 272 | |
---|
[9bb5457] | 273 | BOOLEAN nr2mEqual(number a, number b, const coeffs) |
---|
[35b1d7] | 274 | { |
---|
[e5422d] | 275 | return (a == b); |
---|
[35b1d7] | 276 | } |
---|
| 277 | |
---|
[e90dfd6] | 278 | BOOLEAN nr2mGreater(number a, number b, const coeffs r) |
---|
[009d80] | 279 | { |
---|
[925a43c] | 280 | return nr2mDivBy(a, b,r); |
---|
[009d80] | 281 | } |
---|
| 282 | |
---|
[44d5ad] | 283 | /* Is 'a' divisible by 'b'? There are two cases: |
---|
[9cd697] | 284 | 1) a = 0 mod 2^m; then TRUE iff b = 0 or b is a power of 2 |
---|
[44d5ad] | 285 | 2) a, b <> 0; then TRUE iff b/gcd(a, b) is a unit mod 2^m */ |
---|
[14b11bb] | 286 | BOOLEAN nr2mDivBy (number a, number b, const coeffs r) |
---|
[35b1d7] | 287 | { |
---|
[4fb4f3] | 288 | if (a == NULL) |
---|
| 289 | { |
---|
[e90dfd6] | 290 | NATNUMBER c = r->mod2mMask + 1; |
---|
[4fb4f3] | 291 | if (c != 0) /* i.e., if no overflow */ |
---|
| 292 | return (c % (NATNUMBER)b) == 0; |
---|
| 293 | else |
---|
| 294 | { |
---|
| 295 | /* overflow: we need to check whether b |
---|
[9cd697] | 296 | is zero or a power of 2: */ |
---|
[4fb4f3] | 297 | c = (NATNUMBER)b; |
---|
| 298 | while (c != 0) |
---|
| 299 | { |
---|
| 300 | if ((c % 2) != 0) return FALSE; |
---|
| 301 | c = c >> 1; |
---|
| 302 | } |
---|
| 303 | return TRUE; |
---|
| 304 | } |
---|
| 305 | } |
---|
| 306 | else |
---|
[9cd697] | 307 | { |
---|
[14b11bb] | 308 | number n = nr2mGcd(a, b, r); |
---|
| 309 | n = nr2mDiv(b, n, r); |
---|
| 310 | return nr2mIsUnit(n, r); |
---|
[9cd697] | 311 | } |
---|
[35b1d7] | 312 | } |
---|
| 313 | |
---|
[9bb5457] | 314 | int nr2mDivComp(number as, number bs, const coeffs) |
---|
[206e158] | 315 | { |
---|
[e90dfd6] | 316 | NATNUMBER a = (NATNUMBER)as; |
---|
| 317 | NATNUMBER b = (NATNUMBER)bs; |
---|
[206e158] | 318 | assume(a != 0 && b != 0); |
---|
| 319 | while (a % 2 == 0 && b % 2 == 0) |
---|
| 320 | { |
---|
| 321 | a = a / 2; |
---|
| 322 | b = b / 2; |
---|
| 323 | } |
---|
| 324 | if (a % 2 == 0) |
---|
| 325 | { |
---|
| 326 | return -1; |
---|
| 327 | } |
---|
| 328 | else |
---|
| 329 | { |
---|
| 330 | if (b % 2 == 1) |
---|
| 331 | { |
---|
[91d286] | 332 | return 2; |
---|
[206e158] | 333 | } |
---|
| 334 | else |
---|
| 335 | { |
---|
| 336 | return 1; |
---|
| 337 | } |
---|
| 338 | } |
---|
| 339 | } |
---|
| 340 | |
---|
[76e501] | 341 | /* TRUE iff 0 < k <= 2^m / 2 */ |
---|
[e90dfd6] | 342 | BOOLEAN nr2mGreaterZero(number k, const coeffs r) |
---|
[35b1d7] | 343 | { |
---|
[76e501] | 344 | if ((NATNUMBER)k == 0) return FALSE; |
---|
[e90dfd6] | 345 | if ((NATNUMBER)k > ((r->mod2mMask >> 1) + 1)) return FALSE; |
---|
[76e501] | 346 | return TRUE; |
---|
[35b1d7] | 347 | } |
---|
| 348 | |
---|
[76e501] | 349 | /* assumes that 'a' is odd, i.e., a unit in Z/2^m, and computes |
---|
| 350 | the extended gcd of 'a' and 2^m, in order to find some 's' |
---|
| 351 | and 't' such that a * s + 2^m * t = gcd(a, 2^m) = 1; |
---|
| 352 | this code will always find a positive 's' */ |
---|
[e90dfd6] | 353 | void specialXGCD(unsigned long& s, unsigned long a, const coeffs r) |
---|
[76e501] | 354 | { |
---|
| 355 | int_number u = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 356 | mpz_init_set_ui(u, a); |
---|
| 357 | int_number u0 = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 358 | mpz_init(u0); |
---|
| 359 | int_number u1 = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 360 | mpz_init_set_ui(u1, 1); |
---|
| 361 | int_number u2 = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 362 | mpz_init(u2); |
---|
| 363 | int_number v = (int_number)omAlloc(sizeof(mpz_t)); |
---|
[e90dfd6] | 364 | mpz_init_set_ui(v, r->mod2mMask); |
---|
[76e501] | 365 | mpz_add_ui(v, v, 1); /* now: v = 2^m */ |
---|
| 366 | int_number v0 = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 367 | mpz_init(v0); |
---|
| 368 | int_number v1 = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 369 | mpz_init(v1); |
---|
| 370 | int_number v2 = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 371 | mpz_init_set_ui(v2, 1); |
---|
| 372 | int_number q = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 373 | mpz_init(q); |
---|
[e90dfd6] | 374 | int_number rr = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 375 | mpz_init(rr); |
---|
[76e501] | 376 | |
---|
| 377 | while (mpz_cmp_ui(v, 0) != 0) /* i.e., while v != 0 */ |
---|
| 378 | { |
---|
| 379 | mpz_div(q, u, v); |
---|
[e90dfd6] | 380 | mpz_mod(rr, u, v); |
---|
[76e501] | 381 | mpz_set(u, v); |
---|
[e90dfd6] | 382 | mpz_set(v, rr); |
---|
[76e501] | 383 | mpz_set(u0, u2); |
---|
| 384 | mpz_set(v0, v2); |
---|
| 385 | mpz_mul(u2, u2, q); mpz_sub(u2, u1, u2); /* u2 = u1 - q * u2 */ |
---|
| 386 | mpz_mul(v2, v2, q); mpz_sub(v2, v1, v2); /* v2 = v1 - q * v2 */ |
---|
| 387 | mpz_set(u1, u0); |
---|
| 388 | mpz_set(v1, v0); |
---|
| 389 | } |
---|
| 390 | |
---|
| 391 | while (mpz_cmp_ui(u1, 0) < 0) /* i.e., while u1 < 0 */ |
---|
| 392 | { |
---|
| 393 | /* we add 2^m = (2^m - 1) + 1 to u1: */ |
---|
[e90dfd6] | 394 | mpz_add_ui(u1, u1, r->mod2mMask); |
---|
[76e501] | 395 | mpz_add_ui(u1, u1, 1); |
---|
| 396 | } |
---|
| 397 | s = mpz_get_ui(u1); /* now: 0 <= s <= 2^m - 1 */ |
---|
| 398 | |
---|
| 399 | mpz_clear(u); omFree((ADDRESS)u); |
---|
| 400 | mpz_clear(u0); omFree((ADDRESS)u0); |
---|
| 401 | mpz_clear(u1); omFree((ADDRESS)u1); |
---|
| 402 | mpz_clear(u2); omFree((ADDRESS)u2); |
---|
| 403 | mpz_clear(v); omFree((ADDRESS)v); |
---|
| 404 | mpz_clear(v0); omFree((ADDRESS)v0); |
---|
| 405 | mpz_clear(v1); omFree((ADDRESS)v1); |
---|
| 406 | mpz_clear(v2); omFree((ADDRESS)v2); |
---|
| 407 | mpz_clear(q); omFree((ADDRESS)q); |
---|
[e90dfd6] | 408 | mpz_clear(rr); omFree((ADDRESS)rr); |
---|
[35b1d7] | 409 | } |
---|
| 410 | |
---|
[14b11bb] | 411 | NATNUMBER InvMod(NATNUMBER a, const coeffs r) |
---|
[35b1d7] | 412 | { |
---|
[76e501] | 413 | assume((NATNUMBER)a % 2 != 0); |
---|
| 414 | unsigned long s; |
---|
[14b11bb] | 415 | specialXGCD(s, a, r); |
---|
[76e501] | 416 | return s; |
---|
[35b1d7] | 417 | } |
---|
[cea6f3] | 418 | //#endif |
---|
[35b1d7] | 419 | |
---|
[e90dfd6] | 420 | inline number nr2mInversM(number c, const coeffs r) |
---|
[35b1d7] | 421 | { |
---|
[76e501] | 422 | assume((NATNUMBER)c % 2 != 0); |
---|
[925a43c] | 423 | // Table !!! |
---|
| 424 | NATNUMBER inv; |
---|
| 425 | inv = InvMod((NATNUMBER)c,r); |
---|
[e90dfd6] | 426 | return (number)inv; |
---|
[35b1d7] | 427 | } |
---|
| 428 | |
---|
[e90dfd6] | 429 | number nr2mDiv(number a, number b, const coeffs r) |
---|
[35b1d7] | 430 | { |
---|
[e90dfd6] | 431 | if ((NATNUMBER)a == 0) return (number)0; |
---|
| 432 | else if ((NATNUMBER)b % 2 == 0) |
---|
[35b1d7] | 433 | { |
---|
[994445] | 434 | if ((NATNUMBER)b != 0) |
---|
[b429c16] | 435 | { |
---|
[e90dfd6] | 436 | while (((NATNUMBER)b % 2 == 0) && ((NATNUMBER)a % 2 == 0)) |
---|
[b429c16] | 437 | { |
---|
[e90dfd6] | 438 | a = (number)((NATNUMBER)a / 2); |
---|
| 439 | b = (number)((NATNUMBER)b / 2); |
---|
[b429c16] | 440 | } |
---|
| 441 | } |
---|
[e90dfd6] | 442 | if ((NATNUMBER)b % 2 == 0) |
---|
[b429c16] | 443 | { |
---|
[776bf3e] | 444 | WerrorS("Division not possible, even by cancelling zero divisors."); |
---|
| 445 | WerrorS("Result is integer division without remainder."); |
---|
[67dbdb] | 446 | return (number) ((NATNUMBER) a / (NATNUMBER) b); |
---|
[b429c16] | 447 | } |
---|
[35b1d7] | 448 | } |
---|
[e90dfd6] | 449 | return (number)nr2mMult(a, nr2mInversM(b,r),r); |
---|
[35b1d7] | 450 | } |
---|
| 451 | |
---|
[e90dfd6] | 452 | number nr2mMod(number a, number b, const coeffs r) |
---|
[6ea941] | 453 | { |
---|
| 454 | /* |
---|
[e90dfd6] | 455 | We need to return the number rr which is uniquely determined by the |
---|
[6ea941] | 456 | following two properties: |
---|
[e90dfd6] | 457 | (1) 0 <= rr < |b| (with respect to '<' and '<=' performed in Z x Z) |
---|
| 458 | (2) There exists some k in the integers Z such that a = k * b + rr. |
---|
[6ea941] | 459 | Consider g := gcd(2^m, |b|). Note that then |b|/g is a unit in Z/2^m. |
---|
| 460 | Now, there are three cases: |
---|
| 461 | (a) g = 1 |
---|
| 462 | Then |b| is a unit in Z/2^m, i.e. |b| (and also b) divides a. |
---|
[e90dfd6] | 463 | Thus rr = 0. |
---|
[6ea941] | 464 | (b) g <> 1 and g divides a |
---|
[e90dfd6] | 465 | Then a = (a/g) * (|b|/g)^(-1) * b (up to sign), i.e. again rr = 0. |
---|
[6ea941] | 466 | (c) g <> 1 and g does not divide a |
---|
| 467 | Let's denote the division with remainder of a by g as follows: |
---|
| 468 | a = s * g + t. Then t = a - s * g = a - s * (|b|/g)^(-1) * |b| |
---|
[e90dfd6] | 469 | fulfills (1) and (2), i.e. rr := t is the correct result. Hence |
---|
| 470 | in this third case, rr is the remainder of division of a by g in Z. |
---|
[6ea941] | 471 | This algorithm is the same as for the case Z/n, except that we may |
---|
| 472 | compute the gcd of |b| and 2^m "by hand": We just extract the highest |
---|
| 473 | power of 2 (<= 2^m) that is contained in b. |
---|
| 474 | */ |
---|
[9fc2688] | 475 | assume((NATNUMBER) b != 0); |
---|
[6ea941] | 476 | NATNUMBER g = 1; |
---|
[9fc2688] | 477 | NATNUMBER b_div = (NATNUMBER) b; |
---|
[aa2bcca] | 478 | |
---|
[9fc2688] | 479 | /* |
---|
| 480 | * b_div is unsigned, so that (b_div < 0) evaluates false at compile-time |
---|
| 481 | * |
---|
[9bb5457] | 482 | if (b_div < 0) b_div = -b_div; // b_div now represents |b|, BUT b_div is unsigned! |
---|
[9fc2688] | 483 | */ |
---|
| 484 | |
---|
[e90dfd6] | 485 | NATNUMBER rr = 0; |
---|
| 486 | while ((g < r->mod2mMask ) && (b_div > 0) && (b_div % 2 == 0)) |
---|
[6ea941] | 487 | { |
---|
| 488 | b_div = b_div >> 1; |
---|
| 489 | g = g << 1; |
---|
| 490 | } // g is now the gcd of 2^m and |b| |
---|
| 491 | |
---|
[e90dfd6] | 492 | if (g != 1) rr = (NATNUMBER)a % g; |
---|
| 493 | return (number)rr; |
---|
[6ea941] | 494 | } |
---|
| 495 | |
---|
[e90dfd6] | 496 | number nr2mIntDiv(number a, number b, const coeffs r) |
---|
[f92547] | 497 | { |
---|
[4fb4f3] | 498 | if ((NATNUMBER)a == 0) |
---|
| 499 | { |
---|
| 500 | if ((NATNUMBER)b == 0) |
---|
| 501 | return (number)1; |
---|
| 502 | if ((NATNUMBER)b == 1) |
---|
| 503 | return (number)0; |
---|
[e90dfd6] | 504 | NATNUMBER c = r->mod2mMask + 1; |
---|
[4fb4f3] | 505 | if (c != 0) /* i.e., if no overflow */ |
---|
| 506 | return (number)(c / (NATNUMBER)b); |
---|
| 507 | else |
---|
| 508 | { |
---|
| 509 | /* overflow: c = 2^32 resp. 2^64, depending on platform */ |
---|
| 510 | int_number cc = (int_number)omAlloc(sizeof(mpz_t)); |
---|
[e90dfd6] | 511 | mpz_init_set_ui(cc, r->mod2mMask); mpz_add_ui(cc, cc, 1); |
---|
[4fb4f3] | 512 | mpz_div_ui(cc, cc, (unsigned long)(NATNUMBER)b); |
---|
| 513 | unsigned long s = mpz_get_ui(cc); |
---|
| 514 | mpz_clear(cc); omFree((ADDRESS)cc); |
---|
| 515 | return (number)(NATNUMBER)s; |
---|
| 516 | } |
---|
| 517 | } |
---|
| 518 | else |
---|
| 519 | { |
---|
| 520 | if ((NATNUMBER)b == 0) |
---|
| 521 | return (number)0; |
---|
| 522 | return (number)((NATNUMBER) a / (NATNUMBER) b); |
---|
| 523 | } |
---|
[f92547] | 524 | } |
---|
| 525 | |
---|
[e90dfd6] | 526 | number nr2mInvers(number c, const coeffs r) |
---|
[35b1d7] | 527 | { |
---|
[e90dfd6] | 528 | if ((NATNUMBER)c % 2 == 0) |
---|
[35b1d7] | 529 | { |
---|
| 530 | WerrorS("division by zero divisor"); |
---|
| 531 | return (number)0; |
---|
| 532 | } |
---|
[e90dfd6] | 533 | return nr2mInversM(c, r); |
---|
[35b1d7] | 534 | } |
---|
| 535 | |
---|
[e90dfd6] | 536 | number nr2mNeg(number c, const coeffs r) |
---|
[35b1d7] | 537 | { |
---|
[e90dfd6] | 538 | if ((NATNUMBER)c == 0) return c; |
---|
| 539 | return nr2mNegM(c, r); |
---|
[35b1d7] | 540 | } |
---|
| 541 | |
---|
[9bb5457] | 542 | number nr2mMapMachineInt(number from, const coeffs /*src*/, const coeffs dst) |
---|
[35b1d7] | 543 | { |
---|
[e90dfd6] | 544 | NATNUMBER i = ((NATNUMBER)from) % dst->mod2mMask ; |
---|
| 545 | return (number)i; |
---|
[35b1d7] | 546 | } |
---|
| 547 | |
---|
[9bb5457] | 548 | number nr2mMapZp(number from, const coeffs /*src*/, const coeffs dst) |
---|
[0f93f5] | 549 | { |
---|
[76e501] | 550 | NATNUMBER j = (NATNUMBER)1; |
---|
[e90dfd6] | 551 | long ii = (long)from; |
---|
| 552 | if (ii < 0) { j = dst->mod2mMask; ii = -ii; } |
---|
[76e501] | 553 | NATNUMBER i = (NATNUMBER)ii; |
---|
[e90dfd6] | 554 | i = i & dst->mod2mMask; |
---|
[76e501] | 555 | /* now we have: from = j * i mod 2^m */ |
---|
[925a43c] | 556 | return (number)nr2mMult((number)i, (number)j, dst); |
---|
[894f5b1] | 557 | } |
---|
| 558 | |
---|
[9bb5457] | 559 | number nr2mMapQ(number from, const coeffs /*src*/, const coeffs dst) |
---|
[894f5b1] | 560 | { |
---|
[e90dfd6] | 561 | int_number erg = (int_number)omAllocBin(gmp_nrz_bin); |
---|
[894f5b1] | 562 | mpz_init(erg); |
---|
[e90dfd6] | 563 | int_number k = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 564 | mpz_init_set_ui(k, dst->mod2mMask); |
---|
[894f5b1] | 565 | |
---|
[4d1ae5] | 566 | nlGMP(from, (number)erg, dst); |
---|
[76e501] | 567 | mpz_and(erg, erg, k); |
---|
[14b11bb] | 568 | number res = (number)mpz_get_ui(erg); |
---|
[76e501] | 569 | |
---|
| 570 | mpz_clear(erg); omFree((ADDRESS)erg); |
---|
| 571 | mpz_clear(k); omFree((ADDRESS)k); |
---|
[894f5b1] | 572 | |
---|
[e90dfd6] | 573 | return (number)res; |
---|
[894f5b1] | 574 | } |
---|
| 575 | |
---|
[9bb5457] | 576 | number nr2mMapGMP(number from, const coeffs /*src*/, const coeffs dst) |
---|
[894f5b1] | 577 | { |
---|
[e90dfd6] | 578 | int_number erg = (int_number)omAllocBin(gmp_nrz_bin); |
---|
[894f5b1] | 579 | mpz_init(erg); |
---|
[e90dfd6] | 580 | int_number k = (int_number)omAlloc(sizeof(mpz_t)); |
---|
| 581 | mpz_init_set_ui(k, dst->mod2mMask); |
---|
[894f5b1] | 582 | |
---|
[76e501] | 583 | mpz_and(erg, (int_number)from, k); |
---|
[14b11bb] | 584 | number res = (number) mpz_get_ui(erg); |
---|
[894f5b1] | 585 | |
---|
[76e501] | 586 | mpz_clear(erg); omFree((ADDRESS)erg); |
---|
| 587 | mpz_clear(k); omFree((ADDRESS)k); |
---|
| 588 | |
---|
[e90dfd6] | 589 | return (number)res; |
---|
[894f5b1] | 590 | } |
---|
| 591 | |
---|
[925a43c] | 592 | nMapFunc nr2mSetMap(const coeffs src, const coeffs dst) |
---|
[894f5b1] | 593 | { |
---|
[1cce47] | 594 | if (nCoeff_is_Ring_2toM(src) |
---|
[e90dfd6] | 595 | && (src->mod2mMask == dst->mod2mMask)) |
---|
[894f5b1] | 596 | { |
---|
[925a43c] | 597 | return ndCopyMap; |
---|
[894f5b1] | 598 | } |
---|
[1cce47] | 599 | if (nCoeff_is_Ring_2toM(src) |
---|
[e25a99] | 600 | && (src->mod2mMask < dst->mod2mMask)) |
---|
[5beac5] | 601 | { /* i.e. map an integer mod 2^s into Z mod 2^t, where t < s */ |
---|
[e1375d] | 602 | return nr2mMapMachineInt; |
---|
| 603 | } |
---|
[1cce47] | 604 | if (nCoeff_is_Ring_2toM(src) |
---|
[e25a99] | 605 | && (src->mod2mMask > dst->mod2mMask)) |
---|
[5beac5] | 606 | { /* i.e. map an integer mod 2^s into Z mod 2^t, where t > s */ |
---|
| 607 | // to be done |
---|
| 608 | } |
---|
[1cce47] | 609 | if (nCoeff_is_Ring_Z(src)) |
---|
[894f5b1] | 610 | { |
---|
| 611 | return nr2mMapGMP; |
---|
| 612 | } |
---|
[1cce47] | 613 | if (nCoeff_is_Q(src)) |
---|
[894f5b1] | 614 | { |
---|
| 615 | return nr2mMapQ; |
---|
| 616 | } |
---|
[e90dfd6] | 617 | if (nCoeff_is_Zp(src) && (src->ch == 2)) |
---|
[894f5b1] | 618 | { |
---|
| 619 | return nr2mMapZp; |
---|
| 620 | } |
---|
[1cce47] | 621 | if (nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src)) |
---|
[894f5b1] | 622 | { |
---|
| 623 | // Computing the n of Z/n |
---|
[e90dfd6] | 624 | int_number modul = (int_number)omAllocBin(gmp_nrz_bin); |
---|
| 625 | mpz_init_set(modul, src->modNumber); |
---|
| 626 | int_number twoToTheK = (int_number)omAllocBin(gmp_nrz_bin); |
---|
| 627 | mpz_init_set_ui(twoToTheK, src->mod2mMask); |
---|
| 628 | mpz_add_ui(twoToTheK, twoToTheK, 1); |
---|
| 629 | if (mpz_divisible_p(modul, twoToTheK)) |
---|
[894f5b1] | 630 | { |
---|
[e90dfd6] | 631 | mpz_clear(modul); omFree((void *)modul); |
---|
| 632 | mpz_clear(twoToTheK); omFree((void *)twoToTheK); |
---|
[894f5b1] | 633 | return nr2mMapGMP; |
---|
| 634 | } |
---|
[e90dfd6] | 635 | mpz_clear(modul); omFree((void *) modul); |
---|
| 636 | mpz_clear(twoToTheK); omFree((void *)twoToTheK); |
---|
[894f5b1] | 637 | } |
---|
| 638 | return NULL; // default |
---|
| 639 | } |
---|
[35b1d7] | 640 | |
---|
| 641 | /* |
---|
[e90dfd6] | 642 | * set the exponent |
---|
[35b1d7] | 643 | */ |
---|
| 644 | |
---|
[1112b76] | 645 | void nr2mSetExp(int m, coeffs r) |
---|
[35b1d7] | 646 | { |
---|
[76e501] | 647 | if (m > 1) |
---|
[35b1d7] | 648 | { |
---|
[e90dfd6] | 649 | /* we want mod2mMask to be the bit pattern |
---|
| 650 | '111..1' consisting of m one's: */ |
---|
| 651 | r->mod2mMask = 1; |
---|
| 652 | for (int i = 1; i < m; i++) r->mod2mMask = (r->mod2mMask << 1) + 1; |
---|
[35b1d7] | 653 | } |
---|
| 654 | else |
---|
| 655 | { |
---|
[73a9ffb] | 656 | /* code unexpectedly called with m = 1; we continue with m = 2: */ |
---|
[e90dfd6] | 657 | r->mod2mMask = 3; /* i.e., '11' in binary representation */ |
---|
[35b1d7] | 658 | } |
---|
| 659 | } |
---|
| 660 | |
---|
[1112b76] | 661 | void nr2mInitExp(int m, coeffs r) |
---|
[35b1d7] | 662 | { |
---|
[093f30e] | 663 | nr2mSetExp(m, r); |
---|
[e90dfd6] | 664 | if (m < 2) |
---|
[73a9ffb] | 665 | WarnS("nr2mInitExp unexpectedly called with m = 1 (we continue with Z/2^2"); |
---|
[35b1d7] | 666 | } |
---|
| 667 | |
---|
| 668 | #ifdef LDEBUG |
---|
[9bb5457] | 669 | BOOLEAN nr2mDBTest (number a, const char *, const int, const coeffs r) |
---|
[35b1d7] | 670 | { |
---|
[9bb5457] | 671 | //if ((NATNUMBER)a < 0) return FALSE; // is unsigned! |
---|
[e90dfd6] | 672 | if (((NATNUMBER)a & r->mod2mMask) != (NATNUMBER)a) return FALSE; |
---|
[35b1d7] | 673 | return TRUE; |
---|
| 674 | } |
---|
| 675 | #endif |
---|
| 676 | |
---|
[14b11bb] | 677 | void nr2mWrite (number &a, const coeffs r) |
---|
[35b1d7] | 678 | { |
---|
[aa2bcca] | 679 | long i = nr2mLong(a, r); |
---|
| 680 | StringAppend("%ld", i); |
---|
[35b1d7] | 681 | } |
---|
| 682 | |
---|
[14b11bb] | 683 | static const char* nr2mEati(const char *s, int *i, const coeffs r) |
---|
[35b1d7] | 684 | { |
---|
| 685 | |
---|
| 686 | if (((*s) >= '0') && ((*s) <= '9')) |
---|
| 687 | { |
---|
| 688 | (*i) = 0; |
---|
| 689 | do |
---|
| 690 | { |
---|
| 691 | (*i) *= 10; |
---|
| 692 | (*i) += *s++ - '0'; |
---|
[e90dfd6] | 693 | if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & r->mod2mMask; |
---|
[35b1d7] | 694 | } |
---|
| 695 | while (((*s) >= '0') && ((*s) <= '9')); |
---|
[e90dfd6] | 696 | (*i) = (*i) & r->mod2mMask; |
---|
[35b1d7] | 697 | } |
---|
| 698 | else (*i) = 1; |
---|
| 699 | return s; |
---|
| 700 | } |
---|
| 701 | |
---|
[14b11bb] | 702 | const char * nr2mRead (const char *s, number *a, const coeffs r) |
---|
[35b1d7] | 703 | { |
---|
| 704 | int z; |
---|
| 705 | int n=1; |
---|
| 706 | |
---|
[925a43c] | 707 | s = nr2mEati(s, &z,r); |
---|
[35b1d7] | 708 | if ((*s) == '/') |
---|
| 709 | { |
---|
| 710 | s++; |
---|
[925a43c] | 711 | s = nr2mEati(s, &n,r); |
---|
[35b1d7] | 712 | } |
---|
| 713 | if (n == 1) |
---|
| 714 | *a = (number)z; |
---|
[4f8867] | 715 | else |
---|
[925a43c] | 716 | *a = nr2mDiv((number)z,(number)n,r); |
---|
[35b1d7] | 717 | return s; |
---|
| 718 | } |
---|
[4f8867] | 719 | #endif |
---|
[8d0331d] | 720 | /* #ifdef HAVE_RINGS */ |
---|