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