[275ecc] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[31e857] | 4 | /* $Id: rmodulon.cc,v 1.14 2008-02-01 10:31:44 wienand Exp $ */ |
---|
[275ecc] | 5 | /* |
---|
| 6 | * ABSTRACT: numbers modulo n |
---|
| 7 | */ |
---|
| 8 | |
---|
| 9 | #include <string.h> |
---|
| 10 | #include "mod2.h" |
---|
| 11 | #include <mylimits.h> |
---|
| 12 | #include "structs.h" |
---|
| 13 | #include "febase.h" |
---|
| 14 | #include "omalloc.h" |
---|
| 15 | #include "numbers.h" |
---|
| 16 | #include "longrat.h" |
---|
| 17 | #include "mpr_complex.h" |
---|
| 18 | #include "ring.h" |
---|
| 19 | #include "rmodulon.h" |
---|
[8e56ad] | 20 | #include "gmp.h" |
---|
[275ecc] | 21 | |
---|
| 22 | #ifdef HAVE_RINGMODN |
---|
| 23 | |
---|
[8e56ad] | 24 | typedef MP_INT *int_number; |
---|
[31e857] | 25 | omBin gmp_bin = omGetSpecBin(sizeof(MP_INT)); |
---|
[275ecc] | 26 | |
---|
[12ea9d] | 27 | int_number nrnModul = NULL; |
---|
| 28 | int_number nrnMinusOne = NULL; |
---|
| 29 | unsigned long nrnExponent = 0; |
---|
| 30 | unsigned long long nrnBase = 0; |
---|
[275ecc] | 31 | |
---|
[8e1c4e] | 32 | /* |
---|
| 33 | * create a number from int |
---|
| 34 | */ |
---|
| 35 | number nrnInit (int i) |
---|
| 36 | { |
---|
[31e857] | 37 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e1c4e] | 38 | mpz_init_set_si(erg, i); |
---|
| 39 | mpz_mod(erg, erg, nrnModul); |
---|
| 40 | return (number) erg; |
---|
| 41 | } |
---|
| 42 | |
---|
| 43 | void nrnDelete(number *a, const ring r) |
---|
| 44 | { |
---|
| 45 | mpz_clear((int_number) *a); |
---|
[31e857] | 46 | omFreeBin((ADDRESS) *a, gmp_bin); |
---|
[8e1c4e] | 47 | } |
---|
| 48 | |
---|
| 49 | /* |
---|
| 50 | * convert a number to int (-p/2 .. p/2) |
---|
| 51 | */ |
---|
| 52 | int nrnInt(number &n) |
---|
| 53 | { |
---|
| 54 | return (int) mpz_get_si( (int_number) &n); |
---|
| 55 | } |
---|
| 56 | |
---|
[275ecc] | 57 | /* |
---|
| 58 | * Multiply two numbers |
---|
| 59 | */ |
---|
[8e56ad] | 60 | number nrnMult (number a, number b) |
---|
[275ecc] | 61 | { |
---|
[31e857] | 62 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e56ad] | 63 | mpz_init(erg); |
---|
| 64 | mpz_mul(erg, (int_number) a, (int_number) b); |
---|
| 65 | mpz_mod(erg, erg, nrnModul); |
---|
| 66 | return (number) erg; |
---|
[275ecc] | 67 | } |
---|
| 68 | |
---|
[8e1c4e] | 69 | void nrnPower (number a, int i, number * result) |
---|
| 70 | { |
---|
[31e857] | 71 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e1c4e] | 72 | mpz_init(erg); |
---|
| 73 | mpz_powm_ui(erg, (int_number) a, i, nrnModul); |
---|
| 74 | *result = (number) erg; |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | number nrnAdd (number a, number b) |
---|
| 78 | { |
---|
[31e857] | 79 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e1c4e] | 80 | mpz_init(erg); |
---|
| 81 | mpz_add(erg, (int_number) a, (int_number) b); |
---|
| 82 | mpz_mod(erg, erg, nrnModul); |
---|
| 83 | return (number) erg; |
---|
| 84 | } |
---|
| 85 | |
---|
| 86 | number nrnSub (number a, number b) |
---|
| 87 | { |
---|
[31e857] | 88 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e1c4e] | 89 | mpz_init(erg); |
---|
| 90 | mpz_sub(erg, (int_number) a, (int_number) b); |
---|
| 91 | mpz_mod(erg, erg, nrnModul); |
---|
| 92 | return (number) erg; |
---|
| 93 | } |
---|
| 94 | |
---|
| 95 | number nrnNeg (number c) |
---|
| 96 | { |
---|
[31e857] | 97 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e1c4e] | 98 | mpz_init(erg); |
---|
| 99 | mpz_sub(erg, nrnModul, (int_number) c); |
---|
| 100 | return (number) erg; |
---|
| 101 | } |
---|
| 102 | |
---|
| 103 | number nrnInvers (number c) |
---|
| 104 | { |
---|
[31e857] | 105 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e1c4e] | 106 | mpz_init(erg); |
---|
| 107 | mpz_invert(erg, (int_number) c, nrnModul); |
---|
| 108 | return (number) erg; |
---|
| 109 | } |
---|
| 110 | |
---|
[275ecc] | 111 | /* |
---|
| 112 | * Give the smallest non unit k, such that a * x = k = b * y has a solution |
---|
[8e1c4e] | 113 | * TODO: lcm(gcd,gcd) besser als gcd(lcm) ? |
---|
[275ecc] | 114 | */ |
---|
| 115 | number nrnLcm (number a,number b,ring r) |
---|
| 116 | { |
---|
[8e1c4e] | 117 | number erg = nrnGcd(NULL, a, NULL); |
---|
| 118 | number tmp = nrnGcd(NULL, b, NULL); |
---|
| 119 | mpz_lcm((int_number) erg, (int_number) erg, (int_number) tmp); |
---|
| 120 | nrnDelete(&tmp, NULL); |
---|
[d681e8] | 121 | return (number) erg; |
---|
[275ecc] | 122 | } |
---|
| 123 | |
---|
| 124 | /* |
---|
| 125 | * Give the largest non unit k, such that a = x * k, b = y * k has |
---|
| 126 | * a solution. |
---|
| 127 | */ |
---|
| 128 | number nrnGcd (number a,number b,ring r) |
---|
| 129 | { |
---|
[31e857] | 130 | if ((a == NULL) && (b == NULL)) return nrnInit(0); |
---|
| 131 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
| 132 | mpz_init_set(erg, nrnModul); |
---|
| 133 | if (a != NULL) mpz_gcd(erg, erg, (int_number) a); |
---|
| 134 | if (b != NULL) mpz_gcd(erg, erg, (int_number) b); |
---|
[8e56ad] | 135 | return (number) erg; |
---|
| 136 | } |
---|
| 137 | |
---|
[8e1c4e] | 138 | /* Not needed any more, but may have room for improvement |
---|
[af378f7] | 139 | number nrnGcd3 (number a,number b, number c,ring r) |
---|
| 140 | { |
---|
[31e857] | 141 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[af378f7] | 142 | mpz_init(erg); |
---|
| 143 | if (a == NULL) a = (number) nrnModul; |
---|
| 144 | if (b == NULL) b = (number) nrnModul; |
---|
| 145 | if (c == NULL) c = (number) nrnModul; |
---|
| 146 | mpz_gcd(erg, (int_number) a, (int_number) b); |
---|
| 147 | mpz_gcd(erg, erg, (int_number) c); |
---|
[8e1c4e] | 148 | mpz_gcd(erg, erg, nrnModul); |
---|
[af378f7] | 149 | return (number) erg; |
---|
| 150 | } |
---|
[8e1c4e] | 151 | */ |
---|
[af378f7] | 152 | |
---|
[8e56ad] | 153 | /* |
---|
| 154 | * Give the largest non unit k, such that a = x * k, b = y * k has |
---|
| 155 | * a solution and r, s, s.t. k = s*a + t*b |
---|
| 156 | */ |
---|
| 157 | number nrnExtGcd (number a, number b, number *s, number *t) |
---|
| 158 | { |
---|
[31e857] | 159 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
| 160 | int_number bs = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
| 161 | int_number bt = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e56ad] | 162 | mpz_init(erg); |
---|
| 163 | mpz_init(bs); |
---|
| 164 | mpz_init(bt); |
---|
| 165 | mpz_gcdext(erg, bs, bt, (int_number) a, (int_number) b); |
---|
[8e1c4e] | 166 | mpz_mod(bs, bs, nrnModul); |
---|
| 167 | mpz_mod(bt, bt, nrnModul); |
---|
[8e56ad] | 168 | *s = (number) bs; |
---|
| 169 | *t = (number) bt; |
---|
| 170 | return (number) erg; |
---|
[275ecc] | 171 | } |
---|
| 172 | |
---|
[8e1c4e] | 173 | BOOLEAN nrnIsZero (number a) |
---|
[275ecc] | 174 | { |
---|
[8e1c4e] | 175 | return 0 == mpz_cmpabs_ui((int_number) a, 0); |
---|
[275ecc] | 176 | } |
---|
| 177 | |
---|
[8e1c4e] | 178 | BOOLEAN nrnIsOne (number a) |
---|
[275ecc] | 179 | { |
---|
[8e1c4e] | 180 | return 0 == mpz_cmp_si((int_number) a, 1); |
---|
[8e56ad] | 181 | } |
---|
| 182 | |
---|
[8e1c4e] | 183 | BOOLEAN nrnIsMOne (number a) |
---|
[8e56ad] | 184 | { |
---|
[8e1c4e] | 185 | return 0 == mpz_cmp((int_number) a, nrnMinusOne); |
---|
[275ecc] | 186 | } |
---|
| 187 | |
---|
[8e1c4e] | 188 | BOOLEAN nrnEqual (number a,number b) |
---|
[275ecc] | 189 | { |
---|
[8e1c4e] | 190 | return 0 == mpz_cmp((int_number) a, (int_number) b); |
---|
[275ecc] | 191 | } |
---|
| 192 | |
---|
[8e1c4e] | 193 | BOOLEAN nrnGreater (number a,number b) |
---|
[275ecc] | 194 | { |
---|
[8e1c4e] | 195 | return 0 < mpz_cmp((int_number) a, (int_number) b); |
---|
[275ecc] | 196 | } |
---|
| 197 | |
---|
[8e1c4e] | 198 | BOOLEAN nrnGreaterZero (number k) |
---|
[275ecc] | 199 | { |
---|
[8e1c4e] | 200 | return 0 <= mpz_cmp_si((int_number) k, 0); |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | BOOLEAN nrnIsUnit (number a) |
---|
| 204 | { |
---|
| 205 | number tmp = nrnGcd(a, (number) nrnModul, NULL); |
---|
| 206 | bool res = nrnIsOne(tmp); |
---|
| 207 | nrnDelete(&tmp, NULL); |
---|
| 208 | return res; |
---|
[275ecc] | 209 | } |
---|
| 210 | |
---|
[af378f7] | 211 | number nrnGetUnit (number k) |
---|
[1e579c6] | 212 | { |
---|
[31e857] | 213 | int_number unit = (int_number) nrnGcd(k, 0, currRing); |
---|
| 214 | mpz_tdiv_q(unit, (int_number) k, unit); |
---|
| 215 | int_number gcd = (int_number) nrnGcd((number) unit, 0, currRing); |
---|
| 216 | if (!nrnIsOne((number) gcd)) |
---|
[af378f7] | 217 | { |
---|
[31e857] | 218 | int_number ctmp; |
---|
| 219 | // tmp := unit^2 |
---|
| 220 | int_number tmp = (int_number) nrnMult((number) unit,(number) unit); |
---|
| 221 | // gcd_new := gcd(tmp, 0) |
---|
| 222 | int_number gcd_new = (int_number) nrnGcd((number) tmp, 0, currRing); |
---|
| 223 | while (!nrnEqual((number) gcd_new,(number) gcd)) |
---|
[af378f7] | 224 | { |
---|
[31e857] | 225 | // gcd := gcd_new |
---|
| 226 | ctmp = gcd; |
---|
[af378f7] | 227 | gcd = gcd_new; |
---|
[31e857] | 228 | gcd_new = ctmp; |
---|
| 229 | // tmp := tmp * unit |
---|
| 230 | mpz_mul(tmp, tmp, unit); |
---|
| 231 | mpz_mod(tmp, tmp, nrnModul); |
---|
| 232 | // gcd_new := gcd(tmp, 0) |
---|
| 233 | mpz_gcd(gcd_new, tmp, nrnModul); |
---|
[af378f7] | 234 | } |
---|
[31e857] | 235 | // unit := unit + nrnModul / gcd_new |
---|
| 236 | mpz_tdiv_q(tmp, nrnModul, gcd_new); |
---|
| 237 | mpz_add(unit, unit, tmp); |
---|
| 238 | mpz_mod(unit, unit, nrnModul); |
---|
| 239 | nrnDelete((number*) &gcd_new, NULL); |
---|
| 240 | nrnDelete((number*) &tmp, NULL); |
---|
[af378f7] | 241 | } |
---|
[31e857] | 242 | nrnDelete((number*) &gcd, NULL); |
---|
| 243 | return (number) unit; |
---|
[1e579c6] | 244 | } |
---|
| 245 | |
---|
[8e1c4e] | 246 | BOOLEAN nrnDivBy (number a,number b) |
---|
[275ecc] | 247 | { |
---|
[31e857] | 248 | number bs = nrnGcd(a, b, NULL); |
---|
| 249 | mpz_tdiv_q((int_number) bs, (int_number) b, (int_number) bs); |
---|
[8e1c4e] | 250 | bool res = nrnIsUnit(bs); |
---|
| 251 | nrnDelete(&bs, NULL); |
---|
[af378f7] | 252 | return res; |
---|
[275ecc] | 253 | } |
---|
| 254 | |
---|
[8e56ad] | 255 | int nrnComp(number a, number b) |
---|
[275ecc] | 256 | { |
---|
[8e56ad] | 257 | if (nrnEqual(a, b)) return 0; |
---|
[31e857] | 258 | if (mpz_divisible_p((int_number) a, (int_number) b)) return -1; |
---|
| 259 | if (mpz_divisible_p((int_number) b, (int_number) a)) return 1; |
---|
[8e56ad] | 260 | return 2; |
---|
[275ecc] | 261 | } |
---|
| 262 | |
---|
| 263 | number nrnDiv (number a,number b) |
---|
| 264 | { |
---|
[af378f7] | 265 | if (a == NULL) a = (number) nrnModul; |
---|
[31e857] | 266 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e56ad] | 267 | mpz_init(erg); |
---|
[af378f7] | 268 | if (mpz_divisible_p((int_number) a, (int_number) b)) |
---|
[275ecc] | 269 | { |
---|
[8e56ad] | 270 | mpz_divexact(erg, (int_number) a, (int_number) b); |
---|
| 271 | return (number) erg; |
---|
[275ecc] | 272 | } |
---|
| 273 | else |
---|
| 274 | { |
---|
[8e56ad] | 275 | int_number gcd = (int_number) nrnGcd(a, b, NULL); |
---|
| 276 | mpz_divexact(erg, (int_number) b, gcd); |
---|
| 277 | if (!nrnIsUnit((number) erg)) |
---|
| 278 | { |
---|
[12ea9d] | 279 | WarnS("Division not possible, even by cancelling zero divisors."); |
---|
[8e56ad] | 280 | WarnS("Result is zero."); |
---|
[12ea9d] | 281 | mpz_set_ui(erg, 0); |
---|
[31e857] | 282 | nrnDelete((number*) &gcd, NULL); |
---|
[12ea9d] | 283 | return (number) erg; |
---|
[8e56ad] | 284 | } |
---|
[31e857] | 285 | // a / gcd(a,b) * [b / gcd (a,b)]^(-1) |
---|
| 286 | int_number tmp = (int_number) nrnInvers((number) erg); |
---|
[8e56ad] | 287 | mpz_divexact(erg, (int_number) a, gcd); |
---|
[12ea9d] | 288 | mpz_mul(erg, erg, tmp); |
---|
[31e857] | 289 | nrnDelete((number*) &gcd, NULL); |
---|
| 290 | nrnDelete((number*) &tmp, NULL); |
---|
[8e56ad] | 291 | mpz_mod(erg, erg, nrnModul); |
---|
| 292 | return (number) erg; |
---|
[275ecc] | 293 | } |
---|
| 294 | } |
---|
| 295 | |
---|
| 296 | number nrnIntDiv (number a,number b) |
---|
| 297 | { |
---|
[31e857] | 298 | int_number erg = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e56ad] | 299 | mpz_init(erg); |
---|
[af378f7] | 300 | if (a == NULL) a = (number) nrnModul; |
---|
[8e56ad] | 301 | mpz_tdiv_q(erg, (int_number) a, (int_number) b); |
---|
| 302 | return (number) erg; |
---|
[275ecc] | 303 | } |
---|
| 304 | |
---|
| 305 | nMapFunc nrnSetMap(ring src, ring dst) |
---|
| 306 | { |
---|
| 307 | return NULL; /* default */ |
---|
| 308 | } |
---|
| 309 | |
---|
| 310 | /* |
---|
| 311 | * set the exponent (allocate and init tables) (TODO) |
---|
| 312 | */ |
---|
| 313 | |
---|
[12ea9d] | 314 | void mpz_set_ull(int_number res, unsigned long long xx) |
---|
| 315 | { |
---|
| 316 | unsigned long h = xx >> 32; |
---|
| 317 | mpz_set_ui (res, h); |
---|
| 318 | mpz_mul_2exp (res, res, 32); |
---|
| 319 | mpz_add_ui (res, res, (unsigned long) xx); |
---|
| 320 | } |
---|
| 321 | |
---|
[275ecc] | 322 | void nrnSetExp(int m, ring r) |
---|
| 323 | { |
---|
[12ea9d] | 324 | if ((nrnBase == r->ringflaga) && (nrnExponent == r->ringflagb)) return; |
---|
| 325 | nrnBase = r->ringflaga; |
---|
| 326 | nrnExponent = r->ringflagb; |
---|
| 327 | if (nrnModul != NULL) |
---|
| 328 | { |
---|
[31e857] | 329 | nrnDelete((number*) &nrnModul, NULL); |
---|
| 330 | nrnDelete((number*) &nrnMinusOne, NULL); |
---|
[12ea9d] | 331 | } |
---|
[31e857] | 332 | nrnModul = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[12ea9d] | 333 | mpz_init(nrnModul); |
---|
| 334 | mpz_set_ull(nrnModul, nrnBase); |
---|
| 335 | mpz_pow_ui(nrnModul, nrnModul, nrnExponent); |
---|
| 336 | |
---|
[31e857] | 337 | nrnMinusOne = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[8e56ad] | 338 | mpz_init(nrnMinusOne); |
---|
| 339 | mpz_sub_ui(nrnMinusOne, nrnModul, 1); |
---|
[275ecc] | 340 | } |
---|
| 341 | |
---|
| 342 | void nrnInitExp(int m, ring r) |
---|
| 343 | { |
---|
[12ea9d] | 344 | nrnSetExp(m, r); |
---|
| 345 | if (mpz_cmp_ui(nrnModul,2) <= 0) |
---|
[275ecc] | 346 | { |
---|
| 347 | WarnS("nInitChar failed"); |
---|
| 348 | } |
---|
| 349 | } |
---|
| 350 | |
---|
| 351 | #ifdef LDEBUG |
---|
| 352 | BOOLEAN nrnDBTest (number a, char *f, int l) |
---|
| 353 | { |
---|
[8e56ad] | 354 | if ( (mpz_cmp_si((int_number) a, 0) < 0) || (mpz_cmp((int_number) a, nrnModul) > 0) ) |
---|
[275ecc] | 355 | { |
---|
| 356 | return FALSE; |
---|
| 357 | } |
---|
| 358 | return TRUE; |
---|
| 359 | } |
---|
| 360 | #endif |
---|
| 361 | |
---|
| 362 | void nrnWrite (number &a) |
---|
| 363 | { |
---|
[8e56ad] | 364 | char *s,*z; |
---|
| 365 | if (a==NULL) |
---|
| 366 | { |
---|
| 367 | StringAppendS("o"); |
---|
| 368 | } |
---|
| 369 | else |
---|
| 370 | { |
---|
| 371 | int l=mpz_sizeinbase((int_number) a, 10); |
---|
| 372 | if (a->s<2) l=si_max(l,mpz_sizeinbase((int_number) a,10)); |
---|
| 373 | l+=2; |
---|
| 374 | s=(char*)omAlloc(l); |
---|
| 375 | z=mpz_get_str(s,10,(int_number) a); |
---|
| 376 | StringAppendS(z); |
---|
| 377 | omFreeSize((ADDRESS)s,l); |
---|
| 378 | } |
---|
[275ecc] | 379 | } |
---|
| 380 | |
---|
[8e56ad] | 381 | /*2 |
---|
| 382 | * extracts a long integer from s, returns the rest (COPY FROM longrat0.cc) |
---|
| 383 | */ |
---|
| 384 | char * nlCPEatLongC(char *s, MP_INT *i) |
---|
[275ecc] | 385 | { |
---|
[8e56ad] | 386 | char * start=s; |
---|
[af378f7] | 387 | if (!(*s >= '0' && *s <= '9')) |
---|
| 388 | { |
---|
| 389 | mpz_init_set_si(i, 1); |
---|
| 390 | return s; |
---|
| 391 | } |
---|
| 392 | mpz_init(i); |
---|
[8e56ad] | 393 | while (*s >= '0' && *s <= '9') s++; |
---|
| 394 | if (*s=='\0') |
---|
[275ecc] | 395 | { |
---|
[8e56ad] | 396 | mpz_set_str(i,start,10); |
---|
| 397 | } |
---|
| 398 | else |
---|
| 399 | { |
---|
| 400 | char c=*s; |
---|
| 401 | *s='\0'; |
---|
| 402 | mpz_set_str(i,start,10); |
---|
| 403 | *s=c; |
---|
[275ecc] | 404 | } |
---|
| 405 | return s; |
---|
| 406 | } |
---|
| 407 | |
---|
| 408 | char * nrnRead (char *s, number *a) |
---|
| 409 | { |
---|
[31e857] | 410 | int_number z = (int_number) omAllocBin(gmp_bin); // evtl. spaeter mit bin |
---|
[275ecc] | 411 | { |
---|
[8e56ad] | 412 | s = nlCPEatLongC(s, z); |
---|
[275ecc] | 413 | } |
---|
[8e56ad] | 414 | *a = (number) z; |
---|
[275ecc] | 415 | return s; |
---|
| 416 | } |
---|
| 417 | #endif |
---|