[35aab3] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[f7aa395] | 4 | /* $Id: modulop.cc 14402 2011-09-29 17:16:19Z hannes $ */ |
---|
[35aab3] | 5 | /* |
---|
| 6 | * ABSTRACT: numbers modulo p (<=32003) |
---|
| 7 | */ |
---|
| 8 | |
---|
| 9 | #include <string.h> |
---|
[599326] | 10 | #include <kernel/mod2.h> |
---|
[b1dfaf] | 11 | #include <omalloc/mylimits.h> |
---|
[599326] | 12 | #include <kernel/structs.h> |
---|
| 13 | #include <kernel/febase.h> |
---|
[b1dfaf] | 14 | #include <omalloc/omalloc.h> |
---|
[599326] | 15 | #include <kernel/numbers.h> |
---|
| 16 | #include <kernel/longrat.h> |
---|
| 17 | #include <kernel/mpr_complex.h> |
---|
| 18 | #include <kernel/ring.h> |
---|
[894f5b1] | 19 | #ifdef HAVE_RINGS |
---|
[c1526f] | 20 | #include <si_gmp.h> |
---|
[894f5b1] | 21 | #endif |
---|
[599326] | 22 | #include <kernel/modulop.h> |
---|
[35aab3] | 23 | |
---|
| 24 | long npPrimeM=0; |
---|
| 25 | int npGen=0; |
---|
| 26 | long npPminus1M=0; |
---|
| 27 | long npMapPrime; |
---|
| 28 | |
---|
| 29 | #ifdef HAVE_DIV_MOD |
---|
[900802] | 30 | unsigned short *npInvTable=NULL; |
---|
[35aab3] | 31 | #endif |
---|
| 32 | |
---|
| 33 | #if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD) |
---|
[900802] | 34 | unsigned short *npExpTable=NULL; |
---|
| 35 | unsigned short *npLogTable=NULL; |
---|
[35aab3] | 36 | #endif |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | BOOLEAN npGreaterZero (number k) |
---|
| 40 | { |
---|
[7447d8] | 41 | int h = (int)((long) k); |
---|
[35aab3] | 42 | return ((int)h !=0) && (h <= (npPrimeM>>1)); |
---|
| 43 | } |
---|
| 44 | |
---|
| 45 | //unsigned long npMultMod(unsigned long a, unsigned long b) |
---|
| 46 | //{ |
---|
| 47 | // unsigned long c = a*b; |
---|
| 48 | // c = c % npPrimeM; |
---|
| 49 | // assume(c == (unsigned long) npMultM((number) a, (number) b)); |
---|
| 50 | // return c; |
---|
| 51 | //} |
---|
| 52 | |
---|
| 53 | number npMult (number a,number b) |
---|
| 54 | { |
---|
| 55 | if (((long)a == 0) || ((long)b == 0)) |
---|
| 56 | return (number)0; |
---|
| 57 | else |
---|
| 58 | return npMultM(a,b); |
---|
| 59 | } |
---|
| 60 | |
---|
| 61 | /*2 |
---|
| 62 | * create a number from int |
---|
| 63 | */ |
---|
[8391d8] | 64 | number npInit (int i, const ring r) |
---|
[35aab3] | 65 | { |
---|
| 66 | long ii=i; |
---|
[f7aa395] | 67 | long p=(long)ABS(r->ch); |
---|
| 68 | while (ii < 0L) ii += p; |
---|
| 69 | while ((ii>1L) && (ii >= p)) ii -= p; |
---|
[35aab3] | 70 | return (number)ii; |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | /*2 |
---|
| 74 | * convert a number to int (-p/2 .. p/2) |
---|
| 75 | */ |
---|
[cf74cd6] | 76 | int npInt(number &n, const ring r) |
---|
[35aab3] | 77 | { |
---|
[03579d] | 78 | if ((long)n > (((long)r->ch) >>1)) return (int)((long)n -((long)r->ch)); |
---|
| 79 | else return (int)((long)n); |
---|
[35aab3] | 80 | } |
---|
| 81 | |
---|
| 82 | number npAdd (number a, number b) |
---|
| 83 | { |
---|
| 84 | return npAddM(a,b); |
---|
| 85 | } |
---|
| 86 | |
---|
| 87 | number npSub (number a, number b) |
---|
| 88 | { |
---|
| 89 | return npSubM(a,b); |
---|
| 90 | } |
---|
| 91 | |
---|
| 92 | BOOLEAN npIsZero (number a) |
---|
| 93 | { |
---|
| 94 | return 0 == (long)a; |
---|
| 95 | } |
---|
| 96 | |
---|
| 97 | BOOLEAN npIsOne (number a) |
---|
| 98 | { |
---|
| 99 | return 1 == (long)a; |
---|
| 100 | } |
---|
| 101 | |
---|
| 102 | BOOLEAN npIsMOne (number a) |
---|
| 103 | { |
---|
| 104 | return ((npPminus1M == (long)a)&&((long)1!=(long)a)); |
---|
| 105 | } |
---|
| 106 | |
---|
| 107 | #ifdef HAVE_DIV_MOD |
---|
[1fc18b] | 108 | #ifdef USE_NTL_XGCD |
---|
| 109 | //ifdef HAVE_NTL // in ntl.a |
---|
[35aab3] | 110 | //extern void XGCD(long& d, long& s, long& t, long a, long b); |
---|
| 111 | #include <NTL/ZZ.h> |
---|
| 112 | #ifdef NTL_CLIENT |
---|
| 113 | NTL_CLIENT |
---|
| 114 | #endif |
---|
[1fc18b] | 115 | #endif |
---|
[35aab3] | 116 | |
---|
[1fc18b] | 117 | long InvMod(long a) |
---|
| 118 | { |
---|
| 119 | long d, s, t; |
---|
[35aab3] | 120 | |
---|
[1fc18b] | 121 | #ifdef USE_NTL_XGCD |
---|
| 122 | XGCD(d, s, t, a, npPrimeM); |
---|
| 123 | assume (d == 1); |
---|
| 124 | #else |
---|
| 125 | long u, v, u0, v0, u1, v1, u2, v2, q, r; |
---|
[35aab3] | 126 | |
---|
[1fc18b] | 127 | assume(a>0); |
---|
| 128 | u1=1; u2=0; |
---|
| 129 | u = a; v = npPrimeM; |
---|
[35aab3] | 130 | |
---|
[1fc18b] | 131 | while (v != 0) |
---|
| 132 | { |
---|
[35aab3] | 133 | q = u / v; |
---|
| 134 | r = u % v; |
---|
| 135 | u = v; |
---|
| 136 | v = r; |
---|
| 137 | u0 = u2; |
---|
[1fc18b] | 138 | u2 = u1 - q*u2; |
---|
[35aab3] | 139 | u1 = u0; |
---|
| 140 | } |
---|
| 141 | |
---|
[1fc18b] | 142 | assume(u==1); |
---|
[35aab3] | 143 | s = u1; |
---|
| 144 | #endif |
---|
| 145 | if (s < 0) |
---|
| 146 | return s + npPrimeM; |
---|
| 147 | else |
---|
| 148 | return s; |
---|
| 149 | } |
---|
| 150 | #endif |
---|
| 151 | |
---|
[24cece] | 152 | static inline number npInversM (number c) |
---|
[35aab3] | 153 | { |
---|
| 154 | #ifndef HAVE_DIV_MOD |
---|
[0a78c3] | 155 | return (number)(long)npExpTable[npPminus1M - npLogTable[(long)c]]; |
---|
[35aab3] | 156 | #else |
---|
[0a78c3] | 157 | long inv=(long)npInvTable[(long)c]; |
---|
[35aab3] | 158 | if (inv==0) |
---|
| 159 | { |
---|
| 160 | inv=InvMod((long)c); |
---|
| 161 | npInvTable[(long)c]=inv; |
---|
| 162 | } |
---|
| 163 | return (number)inv; |
---|
| 164 | #endif |
---|
| 165 | } |
---|
| 166 | |
---|
| 167 | number npDiv (number a,number b) |
---|
| 168 | { |
---|
[6c15ec] | 169 | //#ifdef NV_OPS |
---|
| 170 | // if (npPrimeM>NV_MAX_PRIME) |
---|
| 171 | // return nvDiv(a,b); |
---|
| 172 | //#endif |
---|
[35aab3] | 173 | if ((long)a==0) |
---|
| 174 | return (number)0; |
---|
| 175 | #ifndef HAVE_DIV_MOD |
---|
[6c15ec] | 176 | if ((long)b==0) |
---|
[35aab3] | 177 | { |
---|
[b7e838] | 178 | WerrorS(nDivBy0); |
---|
[35aab3] | 179 | return (number)0; |
---|
| 180 | } |
---|
| 181 | else |
---|
| 182 | { |
---|
| 183 | int s = npLogTable[(long)a] - npLogTable[(long)b]; |
---|
| 184 | if (s < 0) |
---|
| 185 | s += npPminus1M; |
---|
[0a78c3] | 186 | return (number)(long)npExpTable[s]; |
---|
[35aab3] | 187 | } |
---|
| 188 | #else |
---|
| 189 | number inv=npInversM(b); |
---|
| 190 | return npMultM(a,inv); |
---|
| 191 | #endif |
---|
| 192 | } |
---|
| 193 | number npInvers (number c) |
---|
| 194 | { |
---|
| 195 | if ((long)c==0) |
---|
| 196 | { |
---|
| 197 | WerrorS("1/0"); |
---|
| 198 | return (number)0; |
---|
| 199 | } |
---|
| 200 | return npInversM(c); |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | number npNeg (number c) |
---|
| 204 | { |
---|
| 205 | if ((long)c==0) return c; |
---|
| 206 | return npNegM(c); |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | BOOLEAN npGreater (number a,number b) |
---|
| 210 | { |
---|
| 211 | //return (long)a != (long)b; |
---|
| 212 | return (long)a > (long)b; |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | BOOLEAN npEqual (number a,number b) |
---|
| 216 | { |
---|
| 217 | // return (long)a == (long)b; |
---|
| 218 | return npEqualM(a,b); |
---|
| 219 | } |
---|
| 220 | |
---|
[493225] | 221 | void npWrite (number &a, const ring r) |
---|
[35aab3] | 222 | { |
---|
[493225] | 223 | if ((long)a>(((long)r->ch) >>1)) StringAppend("-%d",(int)(((long)r->ch)-((long)a))); |
---|
| 224 | else StringAppend("%d",(int)((long)a)); |
---|
[35aab3] | 225 | } |
---|
| 226 | |
---|
| 227 | void npPower (number a, int i, number * result) |
---|
| 228 | { |
---|
| 229 | if (i==0) |
---|
| 230 | { |
---|
| 231 | //npInit(1,result); |
---|
| 232 | *(long *)result = 1; |
---|
| 233 | } |
---|
| 234 | else if (i==1) |
---|
| 235 | { |
---|
| 236 | *result = a; |
---|
| 237 | } |
---|
| 238 | else |
---|
| 239 | { |
---|
| 240 | npPower(a,i-1,result); |
---|
| 241 | *result = npMultM(a,*result); |
---|
| 242 | } |
---|
| 243 | } |
---|
| 244 | |
---|
[85e68dd] | 245 | static const char* npEati(const char *s, int *i) |
---|
[35aab3] | 246 | { |
---|
| 247 | |
---|
| 248 | if (((*s) >= '0') && ((*s) <= '9')) |
---|
| 249 | { |
---|
[923b64] | 250 | unsigned long ii=0L; |
---|
[35aab3] | 251 | do |
---|
| 252 | { |
---|
[923b64] | 253 | ii *= 10; |
---|
| 254 | ii += *s++ - '0'; |
---|
[ead697] | 255 | if (ii >= (MAX_INT_VAL / 10)) ii = ii % npPrimeM; |
---|
[35aab3] | 256 | } |
---|
| 257 | while (((*s) >= '0') && ((*s) <= '9')); |
---|
[923b64] | 258 | if (ii >= npPrimeM) ii = ii % npPrimeM; |
---|
| 259 | *i=(int)ii; |
---|
[35aab3] | 260 | } |
---|
| 261 | else (*i) = 1; |
---|
| 262 | return s; |
---|
| 263 | } |
---|
| 264 | |
---|
[85e68dd] | 265 | const char * npRead (const char *s, number *a) |
---|
[35aab3] | 266 | { |
---|
| 267 | int z; |
---|
| 268 | int n=1; |
---|
| 269 | |
---|
| 270 | s = npEati(s, &z); |
---|
| 271 | if ((*s) == '/') |
---|
| 272 | { |
---|
| 273 | s++; |
---|
| 274 | s = npEati(s, &n); |
---|
| 275 | } |
---|
| 276 | if (n == 1) |
---|
| 277 | *a = (number)z; |
---|
[6c15ec] | 278 | else |
---|
[100025] | 279 | { |
---|
| 280 | if ((z==0)&&(n==0)) WerrorS(nDivBy0); |
---|
[35aab3] | 281 | else |
---|
[100025] | 282 | { |
---|
| 283 | #ifdef NV_OPS |
---|
| 284 | if (npPrimeM>NV_MAX_PRIME) |
---|
| 285 | *a = nvDiv((number)z,(number)n); |
---|
| 286 | else |
---|
[6c15ec] | 287 | #endif |
---|
[100025] | 288 | *a = npDiv((number)z,(number)n); |
---|
| 289 | } |
---|
| 290 | } |
---|
[35aab3] | 291 | return s; |
---|
| 292 | } |
---|
| 293 | |
---|
| 294 | /*2 |
---|
| 295 | * the last used charcteristic |
---|
| 296 | */ |
---|
| 297 | //int npGetChar() |
---|
| 298 | //{ |
---|
| 299 | // return npPrimeM; |
---|
| 300 | //} |
---|
| 301 | |
---|
| 302 | /*2 |
---|
| 303 | * set the charcteristic (allocate and init tables) |
---|
| 304 | */ |
---|
| 305 | |
---|
| 306 | void npSetChar(int c, ring r) |
---|
| 307 | { |
---|
[e8a519] | 308 | |
---|
[35aab3] | 309 | // if (c==npPrimeM) return; |
---|
| 310 | if ((c>1) || (c<(-1))) |
---|
| 311 | { |
---|
| 312 | if (c>1) npPrimeM = c; |
---|
| 313 | else npPrimeM = -c; |
---|
| 314 | npPminus1M = npPrimeM - 1; |
---|
[5bca73] | 315 | #ifdef NV_OPS |
---|
| 316 | if (r->cf->npPrimeM >NV_MAX_PRIME) return; |
---|
| 317 | #endif |
---|
[35aab3] | 318 | #ifdef HAVE_DIV_MOD |
---|
| 319 | npInvTable=r->cf->npInvTable; |
---|
| 320 | #endif |
---|
| 321 | #if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD) |
---|
| 322 | npExpTable=r->cf->npExpTable; |
---|
| 323 | npLogTable=r->cf->npLogTable; |
---|
| 324 | npGen = npExpTable[1]; |
---|
| 325 | #endif |
---|
| 326 | } |
---|
| 327 | else |
---|
| 328 | { |
---|
| 329 | npPrimeM=0; |
---|
| 330 | #ifdef HAVE_DIV_MOD |
---|
| 331 | npInvTable=NULL; |
---|
| 332 | #endif |
---|
| 333 | #if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD) |
---|
| 334 | npExpTable=NULL; |
---|
| 335 | npLogTable=NULL; |
---|
| 336 | #endif |
---|
| 337 | } |
---|
| 338 | } |
---|
| 339 | |
---|
| 340 | void npInitChar(int c, ring r) |
---|
| 341 | { |
---|
| 342 | int i, w; |
---|
| 343 | |
---|
| 344 | if ((c>1) || (c<(-1))) |
---|
| 345 | { |
---|
| 346 | if (c>1) r->cf->npPrimeM = c; |
---|
| 347 | else r->cf->npPrimeM = -c; |
---|
| 348 | r->cf->npPminus1M = r->cf->npPrimeM - 1; |
---|
| 349 | #ifdef NV_OPS |
---|
| 350 | if (r->cf->npPrimeM <=NV_MAX_PRIME) |
---|
| 351 | #endif |
---|
| 352 | { |
---|
| 353 | #if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD) |
---|
[900802] | 354 | r->cf->npExpTable=(unsigned short *)omAlloc( r->cf->npPrimeM*sizeof(unsigned short) ); |
---|
| 355 | r->cf->npLogTable=(unsigned short *)omAlloc( r->cf->npPrimeM*sizeof(unsigned short) ); |
---|
[35aab3] | 356 | r->cf->npExpTable[0] = 1; |
---|
| 357 | r->cf->npLogTable[0] = 0; |
---|
| 358 | if (r->cf->npPrimeM > 2) |
---|
| 359 | { |
---|
| 360 | w = 1; |
---|
| 361 | loop |
---|
| 362 | { |
---|
| 363 | r->cf->npLogTable[1] = 0; |
---|
| 364 | w++; |
---|
| 365 | i = 0; |
---|
| 366 | loop |
---|
| 367 | { |
---|
| 368 | i++; |
---|
| 369 | r->cf->npExpTable[i] =(int)(((long)w * (long)r->cf->npExpTable[i-1]) |
---|
| 370 | % r->cf->npPrimeM); |
---|
| 371 | r->cf->npLogTable[r->cf->npExpTable[i]] = i; |
---|
| 372 | if (/*(i == npPrimeM - 1 ) ||*/ (r->cf->npExpTable[i] == 1)) |
---|
| 373 | break; |
---|
| 374 | } |
---|
| 375 | if (i == r->cf->npPrimeM - 1) |
---|
| 376 | break; |
---|
| 377 | } |
---|
| 378 | } |
---|
| 379 | else |
---|
| 380 | { |
---|
| 381 | r->cf->npExpTable[1] = 1; |
---|
| 382 | r->cf->npLogTable[1] = 0; |
---|
| 383 | } |
---|
| 384 | #endif |
---|
| 385 | #ifdef HAVE_DIV_MOD |
---|
[900802] | 386 | r->cf->npInvTable=(unsigned short*)omAlloc0( r->cf->npPrimeM*sizeof(unsigned short) ); |
---|
[35aab3] | 387 | #endif |
---|
| 388 | } |
---|
| 389 | } |
---|
| 390 | else |
---|
| 391 | { |
---|
| 392 | WarnS("nInitChar failed"); |
---|
| 393 | } |
---|
| 394 | } |
---|
| 395 | |
---|
| 396 | #ifdef LDEBUG |
---|
[85e68dd] | 397 | BOOLEAN npDBTest (number a, const char *f, const int l) |
---|
[35aab3] | 398 | { |
---|
| 399 | if (((long)a<0) || ((long)a>npPrimeM)) |
---|
| 400 | { |
---|
[f1e33bb] | 401 | Print("wrong mod p number %ld at %s,%d\n",(long)a,f,l); |
---|
[35aab3] | 402 | return FALSE; |
---|
| 403 | } |
---|
| 404 | return TRUE; |
---|
| 405 | } |
---|
| 406 | #endif |
---|
| 407 | |
---|
| 408 | number npMap0(number from) |
---|
| 409 | { |
---|
[8391d8] | 410 | return npInit(nlModP(from,npPrimeM),currRing); |
---|
[35aab3] | 411 | } |
---|
| 412 | |
---|
| 413 | number npMapP(number from) |
---|
| 414 | { |
---|
| 415 | long i = (long)from; |
---|
| 416 | if (i>npMapPrime/2) |
---|
| 417 | { |
---|
| 418 | i-=npMapPrime; |
---|
| 419 | while (i < 0) i+=npPrimeM; |
---|
| 420 | } |
---|
| 421 | i%=npPrimeM; |
---|
| 422 | return (number)i; |
---|
| 423 | } |
---|
| 424 | |
---|
| 425 | static number npMapLongR(number from) |
---|
| 426 | { |
---|
| 427 | gmp_float *ff=(gmp_float*)from; |
---|
| 428 | mpf_t *f=ff->_mpfp(); |
---|
| 429 | number res; |
---|
[a604c3] | 430 | mpz_ptr dest,ndest; |
---|
[35aab3] | 431 | int size,i; |
---|
[7447d8] | 432 | int e,al,bl,in; |
---|
| 433 | long iz; |
---|
[35aab3] | 434 | mp_ptr qp,dd,nn; |
---|
| 435 | |
---|
| 436 | size = (*f)[0]._mp_size; |
---|
| 437 | if (size == 0) |
---|
[8391d8] | 438 | return npInit(0,currRing); |
---|
[35aab3] | 439 | if(size<0) |
---|
| 440 | size = -size; |
---|
| 441 | |
---|
| 442 | qp = (*f)[0]._mp_d; |
---|
| 443 | while(qp[0]==0) |
---|
| 444 | { |
---|
| 445 | qp++; |
---|
| 446 | size--; |
---|
| 447 | } |
---|
| 448 | |
---|
| 449 | if(npPrimeM>2) |
---|
| 450 | e=(*f)[0]._mp_exp-size; |
---|
| 451 | else |
---|
| 452 | e=0; |
---|
[896561] | 453 | res = ALLOC_RNUMBER(); |
---|
[35aab3] | 454 | #if defined(LDEBUG) |
---|
| 455 | res->debug=123456; |
---|
| 456 | #endif |
---|
[a604c3] | 457 | dest = res->z; |
---|
[35aab3] | 458 | |
---|
| 459 | if (e<0) |
---|
| 460 | { |
---|
| 461 | al = dest->_mp_size = size; |
---|
| 462 | if (al<2) al = 2; |
---|
| 463 | dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al); |
---|
| 464 | for (i=0;i<size;i++) dd[i] = qp[i]; |
---|
| 465 | bl = 1-e; |
---|
| 466 | nn = (mp_ptr)omAlloc(sizeof(mp_limb_t)*bl); |
---|
| 467 | nn[bl-1] = 1; |
---|
| 468 | for (i=bl-2;i>=0;i--) nn[i] = 0; |
---|
[a604c3] | 469 | ndest = res->n; |
---|
[35aab3] | 470 | ndest->_mp_d = nn; |
---|
| 471 | ndest->_mp_alloc = ndest->_mp_size = bl; |
---|
| 472 | res->s = 0; |
---|
[603aebf] | 473 | in=mpz_fdiv_ui(ndest,npPrimeM); |
---|
[35aab3] | 474 | mpz_clear(ndest); |
---|
| 475 | } |
---|
| 476 | else |
---|
| 477 | { |
---|
| 478 | al = dest->_mp_size = size+e; |
---|
| 479 | if (al<2) al = 2; |
---|
| 480 | dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al); |
---|
| 481 | for (i=0;i<size;i++) dd[i+e] = qp[i]; |
---|
| 482 | for (i=0;i<e;i++) dd[i] = 0; |
---|
| 483 | res->s = 3; |
---|
| 484 | } |
---|
| 485 | |
---|
| 486 | dest->_mp_d = dd; |
---|
| 487 | dest->_mp_alloc = al; |
---|
[603aebf] | 488 | iz=mpz_fdiv_ui(dest,npPrimeM); |
---|
[35aab3] | 489 | mpz_clear(dest); |
---|
| 490 | if(res->s==0) |
---|
[7447d8] | 491 | iz=(long)npDiv((number)iz,(number)in); |
---|
[896561] | 492 | FREE_RNUMBER(res); |
---|
[35aab3] | 493 | return (number)iz; |
---|
| 494 | } |
---|
| 495 | |
---|
[894f5b1] | 496 | #ifdef HAVE_RINGS |
---|
| 497 | /*2 |
---|
| 498 | * convert from a GMP integer |
---|
| 499 | */ |
---|
| 500 | number npMapGMP(number from) |
---|
| 501 | { |
---|
[a604c3] | 502 | int_number erg = (int_number) omAlloc(sizeof(mpz_t)); // evtl. spaeter mit bin |
---|
[894f5b1] | 503 | mpz_init(erg); |
---|
| 504 | |
---|
| 505 | mpz_mod_ui(erg, (int_number) from, npPrimeM); |
---|
| 506 | number r = (number) mpz_get_si(erg); |
---|
| 507 | |
---|
| 508 | mpz_clear(erg); |
---|
| 509 | omFree((ADDRESS) erg); |
---|
| 510 | return (number) r; |
---|
| 511 | } |
---|
| 512 | |
---|
| 513 | /*2 |
---|
| 514 | * convert from an machine long |
---|
| 515 | */ |
---|
| 516 | number npMapMachineInt(number from) |
---|
| 517 | { |
---|
| 518 | long i = (long) (((unsigned long) from) % npPrimeM); |
---|
| 519 | return (number) i; |
---|
| 520 | } |
---|
| 521 | #endif |
---|
| 522 | |
---|
[208e0c] | 523 | nMapFunc npSetMap(const ring src, const ring dst) |
---|
[35aab3] | 524 | { |
---|
[894f5b1] | 525 | #ifdef HAVE_RINGS |
---|
| 526 | if (rField_is_Ring_2toM(src)) |
---|
| 527 | { |
---|
| 528 | return npMapMachineInt; |
---|
| 529 | } |
---|
| 530 | if (rField_is_Ring_Z(src) || rField_is_Ring_PtoM(src) || rField_is_Ring_ModN(src)) |
---|
| 531 | { |
---|
| 532 | return npMapGMP; |
---|
| 533 | } |
---|
| 534 | #endif |
---|
[35aab3] | 535 | if (rField_is_Q(src)) |
---|
| 536 | { |
---|
| 537 | return npMap0; |
---|
| 538 | } |
---|
| 539 | if ( rField_is_Zp(src) ) |
---|
| 540 | { |
---|
| 541 | if (rChar(src) == rChar(dst)) |
---|
| 542 | { |
---|
| 543 | return ndCopy; |
---|
| 544 | } |
---|
| 545 | else |
---|
| 546 | { |
---|
| 547 | npMapPrime=rChar(src); |
---|
| 548 | return npMapP; |
---|
| 549 | } |
---|
| 550 | } |
---|
| 551 | if (rField_is_long_R(src)) |
---|
| 552 | { |
---|
| 553 | return npMapLongR; |
---|
| 554 | } |
---|
| 555 | return NULL; /* default */ |
---|
| 556 | } |
---|
| 557 | |
---|
| 558 | // ----------------------------------------------------------- |
---|
| 559 | // operation for very large primes (32003< p < 2^31-1) |
---|
| 560 | // ---------------------------------------------------------- |
---|
| 561 | #ifdef NV_OPS |
---|
| 562 | |
---|
| 563 | number nvMult (number a,number b) |
---|
| 564 | { |
---|
[2e1b74] | 565 | //if (((long)a == 0) || ((long)b == 0)) |
---|
| 566 | // return (number)0; |
---|
| 567 | //else |
---|
[35aab3] | 568 | return nvMultM(a,b); |
---|
| 569 | } |
---|
| 570 | |
---|
[2e1b74] | 571 | void nvInpMult(number &a, number b, const ring r) |
---|
| 572 | { |
---|
| 573 | number n=nvMultM(a,b); |
---|
| 574 | a=n; |
---|
| 575 | } |
---|
| 576 | |
---|
| 577 | |
---|
[35aab3] | 578 | long nvInvMod(long a) |
---|
| 579 | { |
---|
| 580 | long s, t; |
---|
| 581 | |
---|
| 582 | long u, v, u0, v0, u1, v1, u2, v2, q, r; |
---|
| 583 | |
---|
| 584 | u1=1; v1=0; |
---|
| 585 | u2=0; v2=1; |
---|
| 586 | u = a; v = npPrimeM; |
---|
| 587 | |
---|
[9c8927] | 588 | while (v != 0) |
---|
| 589 | { |
---|
[35aab3] | 590 | q = u / v; |
---|
| 591 | r = u % v; |
---|
| 592 | u = v; |
---|
| 593 | v = r; |
---|
| 594 | u0 = u2; |
---|
| 595 | v0 = v2; |
---|
| 596 | u2 = u1 - q*u2; |
---|
| 597 | v2 = v1- q*v2; |
---|
| 598 | u1 = u0; |
---|
| 599 | v1 = v0; |
---|
| 600 | } |
---|
| 601 | |
---|
| 602 | s = u1; |
---|
| 603 | //t = v1; |
---|
| 604 | if (s < 0) |
---|
| 605 | return s + npPrimeM; |
---|
| 606 | else |
---|
| 607 | return s; |
---|
| 608 | } |
---|
| 609 | |
---|
[24cece] | 610 | static inline number nvInversM (number c) |
---|
[35aab3] | 611 | { |
---|
| 612 | long inv=nvInvMod((long)c); |
---|
| 613 | return (number)inv; |
---|
| 614 | } |
---|
| 615 | |
---|
| 616 | number nvDiv (number a,number b) |
---|
| 617 | { |
---|
| 618 | if ((long)a==0) |
---|
| 619 | return (number)0; |
---|
| 620 | else if ((long)b==0) |
---|
| 621 | { |
---|
[b7e838] | 622 | WerrorS(nDivBy0); |
---|
[35aab3] | 623 | return (number)0; |
---|
| 624 | } |
---|
| 625 | else |
---|
| 626 | { |
---|
| 627 | number inv=nvInversM(b); |
---|
| 628 | return nvMultM(a,inv); |
---|
| 629 | } |
---|
| 630 | } |
---|
| 631 | number nvInvers (number c) |
---|
| 632 | { |
---|
| 633 | if ((long)c==0) |
---|
| 634 | { |
---|
[b7e838] | 635 | WerrorS(nDivBy0); |
---|
[35aab3] | 636 | return (number)0; |
---|
| 637 | } |
---|
| 638 | return nvInversM(c); |
---|
| 639 | } |
---|
[9c8927] | 640 | void nvPower (number a, int i, number * result) |
---|
| 641 | { |
---|
| 642 | if (i==0) |
---|
| 643 | { |
---|
| 644 | //npInit(1,result); |
---|
| 645 | *(long *)result = 1; |
---|
| 646 | } |
---|
| 647 | else if (i==1) |
---|
| 648 | { |
---|
| 649 | *result = a; |
---|
| 650 | } |
---|
| 651 | else |
---|
| 652 | { |
---|
| 653 | nvPower(a,i-1,result); |
---|
| 654 | *result = nvMultM(a,*result); |
---|
| 655 | } |
---|
| 656 | } |
---|
[35aab3] | 657 | #endif |
---|