[35b1d7] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[85e68dd] | 4 | /* $Id: rmodulo2m.cc,v 1.17 2008-03-19 17:44:11 Singular Exp $ */ |
---|
[35b1d7] | 5 | /* |
---|
| 6 | * ABSTRACT: numbers modulo 2^m |
---|
| 7 | */ |
---|
| 8 | |
---|
| 9 | #include <string.h> |
---|
| 10 | #include "mod2.h" |
---|
[2a329d] | 11 | |
---|
| 12 | #ifdef HAVE_RING2TOM |
---|
[35b1d7] | 13 | #include <mylimits.h> |
---|
| 14 | #include "structs.h" |
---|
| 15 | #include "febase.h" |
---|
| 16 | #include "omalloc.h" |
---|
| 17 | #include "numbers.h" |
---|
| 18 | #include "longrat.h" |
---|
| 19 | #include "mpr_complex.h" |
---|
| 20 | #include "ring.h" |
---|
| 21 | #include "rmodulo2m.h" |
---|
| 22 | |
---|
| 23 | int nr2mExp; |
---|
[994445] | 24 | NATNUMBER nr2mModul; |
---|
[35b1d7] | 25 | |
---|
| 26 | /* |
---|
| 27 | * Multiply two numbers |
---|
| 28 | */ |
---|
| 29 | number nr2mMult (number a,number b) |
---|
| 30 | { |
---|
[994445] | 31 | if (((NATNUMBER)a == 0) || ((NATNUMBER)b == 0)) |
---|
[35b1d7] | 32 | return (number)0; |
---|
| 33 | else |
---|
| 34 | return nr2mMultM(a,b); |
---|
| 35 | } |
---|
| 36 | |
---|
[f92547] | 37 | /* |
---|
| 38 | * Give the smallest non unit k, such that a * x = k = b * y has a solution |
---|
| 39 | */ |
---|
| 40 | number nr2mLcm (number a,number b,ring r) |
---|
| 41 | { |
---|
[994445] | 42 | NATNUMBER res = 0; |
---|
| 43 | if ((NATNUMBER) a == 0) a = (number) 1; |
---|
| 44 | if ((NATNUMBER) b == 0) b = (number) 1; |
---|
| 45 | while ((NATNUMBER) a % 2 == 0) |
---|
[f92547] | 46 | { |
---|
[994445] | 47 | a = (number) ((NATNUMBER) a / 2); |
---|
| 48 | if ((NATNUMBER) b % 2 == 0) b = (number) ((NATNUMBER) b / 2); |
---|
[f92547] | 49 | res++; |
---|
| 50 | } |
---|
[994445] | 51 | while ((NATNUMBER) b % 2 == 0) |
---|
[f92547] | 52 | { |
---|
[994445] | 53 | b = (number) ((NATNUMBER) b / 2); |
---|
[f92547] | 54 | res++; |
---|
| 55 | } |
---|
| 56 | return (number) (1L << res); // (2**res) |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | /* |
---|
| 60 | * Give the largest non unit k, such that a = x * k, b = y * k has |
---|
[b429c16] | 61 | * a solution. |
---|
[f92547] | 62 | */ |
---|
| 63 | number nr2mGcd (number a,number b,ring r) |
---|
| 64 | { |
---|
[994445] | 65 | NATNUMBER res = 0; |
---|
| 66 | if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1; |
---|
| 67 | while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0) |
---|
[f92547] | 68 | { |
---|
[994445] | 69 | a = (number) ((NATNUMBER) a / 2); |
---|
| 70 | b = (number) ((NATNUMBER) b / 2); |
---|
[f92547] | 71 | res++; |
---|
| 72 | } |
---|
[206e158] | 73 | // if ((NATNUMBER) b % 2 == 0) |
---|
| 74 | // { |
---|
| 75 | // return (number) ((1L << res));// * (NATNUMBER) a); // (2**res)*a a ist Einheit |
---|
| 76 | // } |
---|
| 77 | // else |
---|
| 78 | // { |
---|
[994445] | 79 | return (number) ((1L << res));// * (NATNUMBER) b); // (2**res)*b b ist Einheit |
---|
[206e158] | 80 | // } |
---|
[f92547] | 81 | } |
---|
| 82 | |
---|
[1e579c6] | 83 | /* |
---|
| 84 | * Give the largest non unit k, such that a = x * k, b = y * k has |
---|
| 85 | * a solution. |
---|
| 86 | */ |
---|
| 87 | number nr2mExtGcd (number a, number b, number *s, number *t) |
---|
| 88 | { |
---|
| 89 | NATNUMBER res = 0; |
---|
| 90 | if ((NATNUMBER) a == 0 && (NATNUMBER) b == 0) return (number) 1; |
---|
| 91 | while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0) |
---|
| 92 | { |
---|
| 93 | a = (number) ((NATNUMBER) a / 2); |
---|
| 94 | b = (number) ((NATNUMBER) b / 2); |
---|
| 95 | res++; |
---|
| 96 | } |
---|
| 97 | if ((NATNUMBER) b % 2 == 0) |
---|
| 98 | { |
---|
| 99 | *t = NULL; |
---|
| 100 | *s = nr2mInvers(a); |
---|
| 101 | return (number) ((1L << res));// * (NATNUMBER) a); // (2**res)*a a ist Einheit |
---|
| 102 | } |
---|
| 103 | else |
---|
| 104 | { |
---|
| 105 | *s = NULL; |
---|
| 106 | *t = nr2mInvers(b); |
---|
| 107 | return (number) ((1L << res));// * (NATNUMBER) b); // (2**res)*b b ist Einheit |
---|
| 108 | } |
---|
| 109 | } |
---|
| 110 | |
---|
[35b1d7] | 111 | void nr2mPower (number a, int i, number * result) |
---|
| 112 | { |
---|
| 113 | if (i==0) |
---|
| 114 | { |
---|
| 115 | //npInit(1,result); |
---|
[994445] | 116 | *(NATNUMBER *)result = 1; |
---|
[35b1d7] | 117 | } |
---|
| 118 | else if (i==1) |
---|
| 119 | { |
---|
| 120 | *result = a; |
---|
| 121 | } |
---|
| 122 | else |
---|
| 123 | { |
---|
| 124 | nr2mPower(a,i-1,result); |
---|
| 125 | *result = nr2mMultM(a,*result); |
---|
| 126 | } |
---|
| 127 | } |
---|
| 128 | |
---|
| 129 | /* |
---|
| 130 | * create a number from int |
---|
| 131 | */ |
---|
| 132 | number nr2mInit (int i) |
---|
| 133 | { |
---|
[d3b2eb] | 134 | long ii = i; |
---|
[35b1d7] | 135 | while (ii < 0) ii += nr2mModul; |
---|
| 136 | while ((ii>1) && (ii >= nr2mModul)) ii -= nr2mModul; |
---|
| 137 | return (number) ii; |
---|
| 138 | } |
---|
| 139 | |
---|
| 140 | /* |
---|
| 141 | * convert a number to int (-p/2 .. p/2) |
---|
| 142 | */ |
---|
| 143 | int nr2mInt(number &n) |
---|
| 144 | { |
---|
[994445] | 145 | if ((NATNUMBER)n > (nr2mModul >>1)) return (int)((NATNUMBER)n - nr2mModul); |
---|
| 146 | else return (int)((NATNUMBER)n); |
---|
[35b1d7] | 147 | } |
---|
| 148 | |
---|
| 149 | number nr2mAdd (number a, number b) |
---|
| 150 | { |
---|
| 151 | return nr2mAddM(a,b); |
---|
| 152 | } |
---|
| 153 | |
---|
| 154 | number nr2mSub (number a, number b) |
---|
| 155 | { |
---|
| 156 | return nr2mSubM(a,b); |
---|
| 157 | } |
---|
| 158 | |
---|
[1e579c6] | 159 | BOOLEAN nr2mIsUnit (number a) |
---|
| 160 | { |
---|
| 161 | return ((NATNUMBER) a % 2 == 1); |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | number nr2mGetUnit (number k) |
---|
| 165 | { |
---|
| 166 | if (k == NULL) |
---|
| 167 | return (number) 1; |
---|
| 168 | NATNUMBER tmp = (NATNUMBER) k; |
---|
| 169 | while (tmp % 2 == 0) |
---|
| 170 | tmp = tmp / 2; |
---|
| 171 | return (number) tmp; |
---|
| 172 | } |
---|
| 173 | |
---|
[35b1d7] | 174 | BOOLEAN nr2mIsZero (number a) |
---|
| 175 | { |
---|
[994445] | 176 | return 0 == (NATNUMBER)a; |
---|
[35b1d7] | 177 | } |
---|
| 178 | |
---|
| 179 | BOOLEAN nr2mIsOne (number a) |
---|
| 180 | { |
---|
[994445] | 181 | return 1 == (NATNUMBER)a; |
---|
[35b1d7] | 182 | } |
---|
| 183 | |
---|
| 184 | BOOLEAN nr2mIsMOne (number a) |
---|
| 185 | { |
---|
[b18621] | 186 | return (nr2mModul == (NATNUMBER)a + 1) && (nr2mModul != 2); |
---|
[35b1d7] | 187 | } |
---|
| 188 | |
---|
| 189 | BOOLEAN nr2mEqual (number a,number b) |
---|
| 190 | { |
---|
| 191 | return nr2mEqualM(a,b); |
---|
| 192 | } |
---|
| 193 | |
---|
| 194 | BOOLEAN nr2mGreater (number a,number b) |
---|
[009d80] | 195 | { |
---|
| 196 | return nr2mDivBy(a, b); |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | BOOLEAN nr2mDivBy (number a,number b) |
---|
[35b1d7] | 200 | { |
---|
[093f30e] | 201 | if (a == NULL) |
---|
| 202 | return (nr2mModul % (NATNUMBER) b) == 0; |
---|
| 203 | else |
---|
| 204 | return ((NATNUMBER) a % (NATNUMBER) b) == 0; |
---|
[821a22] | 205 | /* |
---|
[994445] | 206 | if ((NATNUMBER) a == 0) return TRUE; |
---|
| 207 | if ((NATNUMBER) b == 0) return FALSE; |
---|
| 208 | while ((NATNUMBER) a % 2 == 0 && (NATNUMBER) b % 2 == 0) |
---|
[f92547] | 209 | { |
---|
[994445] | 210 | a = (number) ((NATNUMBER) a / 2); |
---|
| 211 | b = (number) ((NATNUMBER) b / 2); |
---|
[f92547] | 212 | } |
---|
[994445] | 213 | return ((NATNUMBER) b % 2 == 1); |
---|
[821a22] | 214 | */ |
---|
[35b1d7] | 215 | } |
---|
| 216 | |
---|
[206e158] | 217 | int nr2mComp(number as, number bs) |
---|
| 218 | { |
---|
| 219 | NATNUMBER a = (NATNUMBER) as; |
---|
| 220 | NATNUMBER b = (NATNUMBER) bs; |
---|
| 221 | assume(a != 0 && b != 0); |
---|
| 222 | while (a % 2 == 0 && b % 2 == 0) |
---|
| 223 | { |
---|
| 224 | a = a / 2; |
---|
| 225 | b = b / 2; |
---|
| 226 | } |
---|
| 227 | if (a % 2 == 0) |
---|
| 228 | { |
---|
| 229 | return -1; |
---|
| 230 | } |
---|
| 231 | else |
---|
| 232 | { |
---|
| 233 | if (b % 2 == 1) |
---|
| 234 | { |
---|
| 235 | return 0; |
---|
| 236 | } |
---|
| 237 | else |
---|
| 238 | { |
---|
| 239 | return 1; |
---|
| 240 | } |
---|
| 241 | } |
---|
| 242 | } |
---|
| 243 | |
---|
[35b1d7] | 244 | BOOLEAN nr2mGreaterZero (number k) |
---|
| 245 | { |
---|
[009d80] | 246 | return ((NATNUMBER) k !=0) && ((NATNUMBER) k <= (nr2mModul>>1)); |
---|
[35b1d7] | 247 | } |
---|
| 248 | |
---|
[cea6f3] | 249 | //#ifdef HAVE_DIV_MOD |
---|
[35b1d7] | 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 |
---|
| 301 | |
---|
[994445] | 302 | NATNUMBER InvMod(NATNUMBER a) |
---|
[35b1d7] | 303 | { |
---|
| 304 | long d, s, t; |
---|
| 305 | |
---|
| 306 | XGCD(d, s, t, a, nr2mModul); |
---|
| 307 | assume (d == 1); |
---|
| 308 | if (s < 0) |
---|
| 309 | return s + nr2mModul; |
---|
| 310 | else |
---|
| 311 | return s; |
---|
| 312 | } |
---|
[cea6f3] | 313 | //#endif |
---|
[35b1d7] | 314 | |
---|
| 315 | inline number nr2mInversM (number c) |
---|
| 316 | { |
---|
[cea6f3] | 317 | // Table !!! |
---|
[994445] | 318 | NATNUMBER inv; |
---|
| 319 | inv = InvMod((NATNUMBER)c); |
---|
[cea6f3] | 320 | return (number) inv; |
---|
[35b1d7] | 321 | } |
---|
| 322 | |
---|
| 323 | number nr2mDiv (number a,number b) |
---|
| 324 | { |
---|
[994445] | 325 | if ((NATNUMBER)a==0) |
---|
[35b1d7] | 326 | return (number)0; |
---|
[994445] | 327 | else if ((NATNUMBER)b%2==0) |
---|
[35b1d7] | 328 | { |
---|
[994445] | 329 | if ((NATNUMBER)b != 0) |
---|
[b429c16] | 330 | { |
---|
[994445] | 331 | while ((NATNUMBER) b%2 == 0 && (NATNUMBER) a%2 == 0) |
---|
[b429c16] | 332 | { |
---|
[994445] | 333 | a = (number) ((NATNUMBER) a / 2); |
---|
| 334 | b = (number) ((NATNUMBER) b / 2); |
---|
[b429c16] | 335 | } |
---|
| 336 | } |
---|
[994445] | 337 | if ((NATNUMBER) b%2 == 0) |
---|
[b429c16] | 338 | { |
---|
| 339 | WerrorS("div by zero divisor"); |
---|
| 340 | return (number)0; |
---|
| 341 | } |
---|
[35b1d7] | 342 | } |
---|
[b429c16] | 343 | return (number) nr2mMult(a, nr2mInversM(b)); |
---|
[35b1d7] | 344 | } |
---|
| 345 | |
---|
[f92547] | 346 | number nr2mIntDiv (number a,number b) |
---|
| 347 | { |
---|
[994445] | 348 | if ((NATNUMBER)a==0) |
---|
[f92547] | 349 | { |
---|
[206e158] | 350 | if ((NATNUMBER)b==0) |
---|
| 351 | return (number) 1; |
---|
[d681e8] | 352 | if ((NATNUMBER)b==1) |
---|
| 353 | return (number) 0; |
---|
[206e158] | 354 | return (number) (nr2mModul / (NATNUMBER) b); |
---|
[f92547] | 355 | } |
---|
| 356 | else |
---|
| 357 | { |
---|
[206e158] | 358 | if ((NATNUMBER)b==0) |
---|
| 359 | return (number) 0; |
---|
[994445] | 360 | return (number) ((NATNUMBER) a / (NATNUMBER) b); |
---|
[f92547] | 361 | } |
---|
| 362 | } |
---|
| 363 | |
---|
[35b1d7] | 364 | number nr2mInvers (number c) |
---|
| 365 | { |
---|
[994445] | 366 | if ((NATNUMBER)c%2==0) |
---|
[35b1d7] | 367 | { |
---|
| 368 | WerrorS("division by zero divisor"); |
---|
| 369 | return (number)0; |
---|
| 370 | } |
---|
| 371 | return nr2mInversM(c); |
---|
| 372 | } |
---|
| 373 | |
---|
| 374 | number nr2mNeg (number c) |
---|
| 375 | { |
---|
[994445] | 376 | if ((NATNUMBER)c==0) return c; |
---|
[35b1d7] | 377 | return nr2mNegM(c); |
---|
| 378 | } |
---|
| 379 | |
---|
| 380 | nMapFunc nr2mSetMap(ring src, ring dst) |
---|
| 381 | { |
---|
| 382 | return NULL; /* default */ |
---|
| 383 | } |
---|
| 384 | |
---|
| 385 | |
---|
| 386 | /* |
---|
| 387 | * set the exponent (allocate and init tables) (TODO) |
---|
| 388 | */ |
---|
| 389 | |
---|
| 390 | void nr2mSetExp(int m, ring r) |
---|
| 391 | { |
---|
| 392 | if (m>1) |
---|
| 393 | { |
---|
| 394 | nr2mExp = m; |
---|
| 395 | nr2mModul = 2; |
---|
| 396 | for (int i = 1; i < m; i++) { |
---|
| 397 | nr2mModul = nr2mModul * 2; |
---|
| 398 | } |
---|
| 399 | } |
---|
| 400 | else |
---|
| 401 | { |
---|
[093f30e] | 402 | nr2mExp=2; |
---|
| 403 | nr2mModul=4; |
---|
[35b1d7] | 404 | } |
---|
| 405 | } |
---|
| 406 | |
---|
| 407 | void nr2mInitExp(int m, ring r) |
---|
| 408 | { |
---|
[093f30e] | 409 | nr2mSetExp(m, r); |
---|
| 410 | if (m<2) WarnS("nInitExp failed: using Z/2^2"); |
---|
[35b1d7] | 411 | } |
---|
| 412 | |
---|
| 413 | #ifdef LDEBUG |
---|
[85e68dd] | 414 | BOOLEAN nr2mDBTest (number a, const char *f, const int l) |
---|
[35b1d7] | 415 | { |
---|
[994445] | 416 | if (((NATNUMBER)a<0) || ((NATNUMBER)a>nr2mModul)) |
---|
[35b1d7] | 417 | { |
---|
| 418 | return FALSE; |
---|
| 419 | } |
---|
| 420 | return TRUE; |
---|
| 421 | } |
---|
| 422 | #endif |
---|
| 423 | |
---|
| 424 | void nr2mWrite (number &a) |
---|
| 425 | { |
---|
[994445] | 426 | if ((NATNUMBER)a > (nr2mModul >>1)) StringAppend("-%d",(int)(nr2mModul-((NATNUMBER)a))); |
---|
| 427 | else StringAppend("%d",(int)((NATNUMBER)a)); |
---|
[35b1d7] | 428 | } |
---|
| 429 | |
---|
[85e68dd] | 430 | static const char* nr2mEati(const char *s, int *i) |
---|
[35b1d7] | 431 | { |
---|
| 432 | |
---|
| 433 | if (((*s) >= '0') && ((*s) <= '9')) |
---|
| 434 | { |
---|
| 435 | (*i) = 0; |
---|
| 436 | do |
---|
| 437 | { |
---|
| 438 | (*i) *= 10; |
---|
| 439 | (*i) += *s++ - '0'; |
---|
| 440 | if ((*i) >= (MAX_INT_VAL / 10)) (*i) = (*i) % nr2mModul; |
---|
| 441 | } |
---|
| 442 | while (((*s) >= '0') && ((*s) <= '9')); |
---|
| 443 | if ((*i) >= nr2mModul) (*i) = (*i) % nr2mModul; |
---|
| 444 | } |
---|
| 445 | else (*i) = 1; |
---|
| 446 | return s; |
---|
| 447 | } |
---|
| 448 | |
---|
[85e68dd] | 449 | const char * nr2mRead (const char *s, number *a) |
---|
[35b1d7] | 450 | { |
---|
| 451 | int z; |
---|
| 452 | int n=1; |
---|
| 453 | |
---|
| 454 | s = nr2mEati(s, &z); |
---|
| 455 | if ((*s) == '/') |
---|
| 456 | { |
---|
| 457 | s++; |
---|
| 458 | s = nr2mEati(s, &n); |
---|
| 459 | } |
---|
| 460 | if (n == 1) |
---|
| 461 | *a = (number)z; |
---|
[4f8867] | 462 | else |
---|
[35b1d7] | 463 | *a = nr2mDiv((number)z,(number)n); |
---|
| 464 | return s; |
---|
| 465 | } |
---|
[4f8867] | 466 | #endif |
---|