Changeset 76e501 in git for kernel/rmodulo2m.cc
 Timestamp:
 Jul 21, 2010, 10:37:19 AM (13 years ago)
 Branches:
 (u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
 Children:
 f0a801bdd99a5337cb60770f49d314447999474b
 Parents:
 cf21dd42f35f06862c90e2223eea0f33e4ff0fbc
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

kernel/rmodulo2m.cc
rcf21dd4 r76e501 27 27 * Multiply two numbers 28 28 */ 29 number nr2mMult (number a, number b)29 number nr2mMult (number a, number b) 30 30 { 31 31 if (((NATNUMBER)a == 0)  ((NATNUMBER)b == 0)) … … 131 131 number nr2mInit (int i, const ring r) 132 132 { 133 if (i == 0) return (number)(NATNUMBER)i; 134 133 135 long ii = i; 134 while (ii < 0) ii += r>nr2mModul ; 135 while ((ii>1) && (ii >= r>nr2mModul )) ii = r>nr2mModul ; 136 return (number) ii; 137 } 138 139 /* 140 * convert a number to int (p/2 .. p/2) 136 NATNUMBER j = (NATNUMBER)1; 137 if (ii < 0) { j = currRing>nr2mModul; ii = ii; } 138 NATNUMBER k = (NATNUMBER)ii; 139 k = k & currRing>nr2mModul; 140 /* now we have: from = j * k mod 2^m */ 141 return (number)nr2mMult((number)j, (number)k); 142 } 143 144 /* 145 * convert a number to an int in ]k/2 .. k/2], 146 * where k = 2^m; i.e., an int in ]2^(m1) .. 2^(m1)]; 147 * note that the code computes a long which will then 148 * automatically casted to int 141 149 */ 142 150 int nr2mInt(number &n, const ring r) 143 151 { 144 if ((NATNUMBER)n > (r>nr2mModul >>1)) return (int)((NATNUMBER)n  r>nr2mModul ); 145 else return (int)((NATNUMBER)n); 152 NATNUMBER nn = (unsigned long)(NATNUMBER)n & r>nr2mModul; 153 unsigned long l = r>nr2mModul >> 1; l++; 154 if (l == 0) 155 return (int)(signed long)(NATNUMBER)nn; 156 else if ((NATNUMBER)nn > l) 157 return (int)((NATNUMBER)nn  r>nr2mModul  1); 158 else 159 return (int)((NATNUMBER)nn); 146 160 } 147 161 … … 183 197 BOOLEAN nr2mIsMOne (number a) 184 198 { 185 return (currRing>nr2mModul == (NATNUMBER)a + 1) 186 && (currRing>nr2mModul != 2); 187 } 188 189 BOOLEAN nr2mEqual (number a,number b) 199 return (currRing>nr2mModul == (NATNUMBER)a); 200 } 201 202 BOOLEAN nr2mEqual (number a, number b) 190 203 { 191 204 return nr2mEqualM(a,b); 192 205 } 193 206 194 BOOLEAN nr2mGreater (number a, number b)207 BOOLEAN nr2mGreater (number a, number b) 195 208 { 196 209 return nr2mDivBy(a, b); 197 210 } 198 211 199 BOOLEAN nr2mDivBy (number a,number b) 200 { 201 if (a == NULL) 202 return (currRing>nr2mModul % (NATNUMBER) b) == 0; 203 else 204 return ((NATNUMBER) a % (NATNUMBER) b) == 0; 205 /* 206 if ((NATNUMBER) a == 0) return TRUE; 207 if ((NATNUMBER) b == 0) return FALSE; 208 while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0) 209 { 210 a = (number) ((NATNUMBER) a / 2); 211 b = (number) ((NATNUMBER) b / 2); 212 } 213 return ((NATNUMBER) b % 2 == 1); 214 */ 212 BOOLEAN nr2mDivBy (number a, number b) 213 { 214 if ((NATNUMBER)a == 0) return TRUE; 215 if ((NATNUMBER)b == 0) return FALSE; 216 return ((NATNUMBER)a % (NATNUMBER)b) == 0; 215 217 } 216 218 … … 242 244 } 243 245 246 /* TRUE iff 0 < k <= 2^m / 2 */ 244 247 BOOLEAN nr2mGreaterZero (number k) 245 248 { 246 return ((NATNUMBER) k !=0) && ((NATNUMBER) k <= (currRing>nr2mModul >>1)); 247 } 248 249 //#ifdef HAVE_DIV_MOD 250 #if 1 //ifdef HAVE_NTL // in ntl.a 251 //extern void XGCD(long& d, long& s, long& t, long a, long b); 252 #include <NTL/ZZ.h> 253 #ifdef NTL_CLIENT 254 NTL_CLIENT 255 #endif 256 #else 257 void XGCD(long& d, long& s, long& t, long a, long b) 258 { 259 long u, v, u0, v0, u1, v1, u2, v2, q, r; 260 261 long aneg = 0, bneg = 0; 262 263 if (a < 0) { 264 a = a; 265 aneg = 1; 266 } 267 268 if (b < 0) { 269 b = b; 270 bneg = 1; 271 } 272 273 u1=1; v1=0; 274 u2=0; v2=1; 275 u = a; v = b; 276 277 while (v != 0) { 278 q = u / v; 279 r = u % v; 280 u = v; 281 v = r; 282 u0 = u2; 283 v0 = v2; 284 u2 = u1  q*u2; 285 v2 = v1 q*v2; 286 u1 = u0; 287 v1 = v0; 288 } 289 290 if (aneg) 291 u1 = u1; 292 293 if (bneg) 294 v1 = v1; 295 296 d = u; 297 s = u1; 298 t = v1; 299 } 300 #endif 249 if ((NATNUMBER)k == 0) return FALSE; 250 if ((NATNUMBER)k > ((currRing>nr2mModul >> 1) + 1)) return FALSE; 251 return TRUE; 252 } 253 254 /* assumes that 'a' is odd, i.e., a unit in Z/2^m, and computes 255 the extended gcd of 'a' and 2^m, in order to find some 's' 256 and 't' such that a * s + 2^m * t = gcd(a, 2^m) = 1; 257 this code will always find a positive 's' */ 258 void specialXGCD(unsigned long& s, unsigned long a) 259 { 260 int_number u = (int_number)omAlloc(sizeof(mpz_t)); 261 mpz_init_set_ui(u, a); 262 int_number u0 = (int_number)omAlloc(sizeof(mpz_t)); 263 mpz_init(u0); 264 int_number u1 = (int_number)omAlloc(sizeof(mpz_t)); 265 mpz_init_set_ui(u1, 1); 266 int_number u2 = (int_number)omAlloc(sizeof(mpz_t)); 267 mpz_init(u2); 268 int_number v = (int_number)omAlloc(sizeof(mpz_t)); 269 mpz_init_set_ui(v, currRing>nr2mModul); 270 mpz_add_ui(v, v, 1); /* now: v = 2^m */ 271 int_number v0 = (int_number)omAlloc(sizeof(mpz_t)); 272 mpz_init(v0); 273 int_number v1 = (int_number)omAlloc(sizeof(mpz_t)); 274 mpz_init(v1); 275 int_number v2 = (int_number)omAlloc(sizeof(mpz_t)); 276 mpz_init_set_ui(v2, 1); 277 int_number q = (int_number)omAlloc(sizeof(mpz_t)); 278 mpz_init(q); 279 int_number r = (int_number)omAlloc(sizeof(mpz_t)); 280 mpz_init(r); 281 282 while (mpz_cmp_ui(v, 0) != 0) /* i.e., while v != 0 */ 283 { 284 mpz_div(q, u, v); 285 mpz_mod(r, u, v); 286 mpz_set(u, v); 287 mpz_set(v, r); 288 mpz_set(u0, u2); 289 mpz_set(v0, v2); 290 mpz_mul(u2, u2, q); mpz_sub(u2, u1, u2); /* u2 = u1  q * u2 */ 291 mpz_mul(v2, v2, q); mpz_sub(v2, v1, v2); /* v2 = v1  q * v2 */ 292 mpz_set(u1, u0); 293 mpz_set(v1, v0); 294 } 295 296 while (mpz_cmp_ui(u1, 0) < 0) /* i.e., while u1 < 0 */ 297 { 298 /* we add 2^m = (2^m  1) + 1 to u1: */ 299 mpz_add_ui(u1, u1, currRing>nr2mModul); 300 mpz_add_ui(u1, u1, 1); 301 } 302 s = mpz_get_ui(u1); /* now: 0 <= s <= 2^m  1 */ 303 304 mpz_clear(u); omFree((ADDRESS)u); 305 mpz_clear(u0); omFree((ADDRESS)u0); 306 mpz_clear(u1); omFree((ADDRESS)u1); 307 mpz_clear(u2); omFree((ADDRESS)u2); 308 mpz_clear(v); omFree((ADDRESS)v); 309 mpz_clear(v0); omFree((ADDRESS)v0); 310 mpz_clear(v1); omFree((ADDRESS)v1); 311 mpz_clear(v2); omFree((ADDRESS)v2); 312 mpz_clear(q); omFree((ADDRESS)q); 313 mpz_clear(r); omFree((ADDRESS)r); 314 } 301 315 302 316 NATNUMBER InvMod(NATNUMBER a) 303 317 { 304 long d, s, t; 305 306 XGCD(d, s, t, a, currRing>nr2mModul ); 307 assume (d == 1); 308 if (s < 0) 309 return s + currRing>nr2mModul ; 310 else 311 return s; 318 assume((NATNUMBER)a % 2 != 0); 319 unsigned long s; 320 specialXGCD(s, a); 321 return s; 312 322 } 313 323 //#endif … … 315 325 inline number nr2mInversM (number c) 316 326 { 317 // Table !!! 318 NATNUMBER inv; 319 inv = InvMod((NATNUMBER)c); 320 return (number) inv; 327 assume((NATNUMBER)c % 2 != 0); 328 return (number)InvMod((NATNUMBER)c); 321 329 } 322 330 … … 368 376 power of 2 (<= 2^m) that is contained in b. 369 377 */ 378 assume((NATNUMBER)b != 0); 370 379 NATNUMBER g = 1; 371 380 NATNUMBER b_div = (NATNUMBER)b; 372 if (b_div < 0) b_div =  381 if (b_div < 0) b_div = b_div; // b_div now represents b 373 382 NATNUMBER r = 0; 374 while ((g < currRing>nr2mModul 383 while ((g < currRing>nr2mModul) && (b_div > 0) && (b_div % 2 == 0)) 375 384 { 376 385 b_div = b_div >> 1; … … 384 393 number nr2mIntDiv (number a,number b) 385 394 { 386 if ((NATNUMBER)a==0) 387 { 388 if ((NATNUMBER)b==0) 389 return (number) 1; 390 if ((NATNUMBER)b==1) 391 return (number) 0; 392 return (number) (currRing>nr2mModul / (NATNUMBER) b); 393 } 394 else 395 { 396 if ((NATNUMBER)b==0) 397 return (number) 0; 398 return (number) ((NATNUMBER) a / (NATNUMBER) b); 399 } 395 assume((NATNUMBER)b != 0); 396 return (number) ((NATNUMBER) a / (NATNUMBER) b); 400 397 } 401 398 … … 418 415 number nr2mMapMachineInt(number from) 419 416 { 420 NATNUMBER i = ((NATNUMBER) from) % currRing>nr2mModul;417 NATNUMBER i = ((NATNUMBER) from) & currRing>nr2mModul; 421 418 return (number) i; 422 419 } … … 424 421 number nr2mMapZp(number from) 425 422 { 426 long ii = (long) from; 427 while (ii < 0) ii += currRing>nr2mModul ; 428 while ((ii>1) && (ii >= currRing>nr2mModul )) ii = currRing>nr2mModul ; 429 return (number) ii; 423 long ii = (long)from; 424 NATNUMBER j = (NATNUMBER)1; 425 if (ii < 0) { j = currRing>nr2mModul; ii = ii; } 426 NATNUMBER i = (NATNUMBER)ii; 427 i = i & currRing>nr2mModul; 428 /* now we have: from = j * i mod 2^m */ 429 return (number)nr2mMult((number)i, (number)j); 430 430 } 431 431 432 432 number nr2mMapQ(number from) 433 433 { 434 int_number erg = (int_number) omAlloc(sizeof(mpz_t)); // evtl. spaeter mit bin434 int_number erg = (int_number) omAlloc(sizeof(mpz_t)); 435 435 mpz_init(erg); 436 437 nlGMP(from, (number) erg); 438 mpz_mod_ui(erg, erg, currRing>nr2mModul ); 436 int_number k = (int_number) omAlloc(sizeof(mpz_t)); 437 mpz_init_set_ui(k, currRing>nr2mModul); 438 439 nlGMP(from, (number)erg); 440 mpz_and(erg, erg, k); 441 number r = (number)mpz_get_ui(erg); 442 443 mpz_clear(erg); omFree((ADDRESS)erg); 444 mpz_clear(k); omFree((ADDRESS)k); 445 446 return (number) r; 447 } 448 449 number nr2mMapGMP(number from) 450 { 451 int_number erg = (int_number) omAlloc(sizeof(mpz_t)); 452 mpz_init(erg); 453 int_number k = (int_number) omAlloc(sizeof(mpz_t)); 454 mpz_init_set_ui(k, currRing>nr2mModul); 455 456 mpz_and(erg, (int_number)from, k); 439 457 number r = (number) mpz_get_ui(erg); 440 458 441 mpz_clear(erg); 442 omFree((ADDRESS) erg); 443 return (number) r; 444 } 445 446 number nr2mMapGMP(number from) 447 { 448 int_number erg = (int_number) omAlloc(sizeof(mpz_t)); // evtl. spaeter mit bin 449 mpz_init(erg); 450 451 mpz_mod_ui(erg, (int_number) from, currRing>nr2mModul ); 452 number r = (number) mpz_get_ui(erg); 453 454 mpz_clear(erg); 455 omFree((ADDRESS) erg); 459 mpz_clear(erg); omFree((ADDRESS)erg); 460 mpz_clear(k); omFree((ADDRESS)k); 461 456 462 return (number) r; 457 463 } … … 503 509 void nr2mSetExp(int m, const ring r) 504 510 { 505 if (m >1)511 if (m > 1) 506 512 { 507 513 nr2mExp = m; 508 r>nr2mModul = 2; 509 for (int i = 1; i < m; i++) 510 { 511 r>nr2mModul = r>nr2mModul * 2; 512 } 513 } 514 else 515 { 516 nr2mExp=2; 517 r>nr2mModul =4; 514 /* we want nr2mModul to be the bit pattern '11..1' consisting 515 of m one's: */ 516 r>nr2mModul = 1; 517 for (int i = 1; i < m; i++) r>nr2mModul = (r>nr2mModul * 2) + 1; 518 } 519 else 520 { 521 nr2mExp = 2; 522 r>nr2mModul = 3; /* i.e., '11' in binary representation */ 518 523 } 519 524 } … … 522 527 { 523 528 nr2mSetExp(m, r); 524 if (m <2) WarnS("nInitExp failed: usingZ/2^2");529 if (m < 2) WarnS("nr2mInitExp failed: we go on with Z/2^2"); 525 530 } 526 531 … … 528 533 BOOLEAN nr2mDBTest (number a, const char *f, const int l) 529 534 { 530 if (((NATNUMBER)a<0)  ((NATNUMBER)a>currRing>nr2mModul )) 531 { 532 return FALSE; 533 } 535 if ((NATNUMBER)a < 0) return FALSE; 536 if (((NATNUMBER)a & currRing>nr2mModul) != (NATNUMBER)a) return FALSE; 534 537 return TRUE; 535 538 } … … 538 541 void nr2mWrite (number &a, const ring r) 539 542 { 540 if ((NATNUMBER)a > (r>nr2mModul >>1)) 541 StringAppend("%d",(int)(r>nr2mModul ((NATNUMBER)a))); 542 else StringAppend("%d",(int)((NATNUMBER)a)); 543 int i = nr2mInt(a, r); 544 StringAppend("%d", i); 543 545 } 544 546 … … 553 555 (*i) *= 10; 554 556 (*i) += *s++  '0'; 555 if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) % currRing>nr2mModul;557 if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) & currRing>nr2mModul; 556 558 } 557 559 while (((*s) >= '0') && ((*s) <= '9')); 558 (*i) = (*i) % currRing>nr2mModul;560 (*i) = (*i) & currRing>nr2mModul; 559 561 } 560 562 else (*i) = 1;
Note: See TracChangeset
for help on using the changeset viewer.