| 1 | #include <math.h> |
| 2 | #include <stdio.h> |
| 3 | #include <stdlib.h> |
| 4 | #include <string.h> |
| 5 | |
| 6 | #include <mLib/alloc.h> |
| 7 | #include <mLib/darray.h> |
| 8 | #include <mLib/dstr.h> |
| 9 | #include <mLib/quis.h> |
| 10 | #include <mLib/report.h> |
| 11 | |
| 12 | #include "fibrand.h" |
| 13 | #include "field.h" |
| 14 | #include "ec.h" |
| 15 | #include "gf.h" |
| 16 | #include "gfreduce.h" |
| 17 | #include "mp.h" |
| 18 | #include "mpbarrett.h" |
| 19 | #include "mpint.h" |
| 20 | #include "mprand.h" |
| 21 | #include "pgen.h" |
| 22 | #include "primetab.h" |
| 23 | |
| 24 | DA_DECL(mp_v, mp *); |
| 25 | |
| 26 | typedef struct fbe { |
| 27 | unsigned long p; |
| 28 | unsigned e; |
| 29 | mp *n; |
| 30 | } fbe; |
| 31 | DA_DECL(fbe_v, fbe); |
| 32 | |
| 33 | typedef struct fact { |
| 34 | mp *p; |
| 35 | unsigned e; |
| 36 | mp *n; |
| 37 | mp *r; |
| 38 | } fact; |
| 39 | DA_DECL(fact_v, fact); |
| 40 | |
| 41 | static grand *rng; |
| 42 | |
| 43 | /*----- Factoring ---------------------------------------------------------*/ |
| 44 | |
| 45 | static void mpv_free(mp_v *v) |
| 46 | { |
| 47 | size_t i; |
| 48 | |
| 49 | for (i = 0; i < DA_LEN(v); i++) |
| 50 | mp_drop(DA(v)[i]); |
| 51 | DA_DESTROY(v); |
| 52 | } |
| 53 | |
| 54 | static mp *smallfactors(mp *x, mp_v *v) |
| 55 | { |
| 56 | mp f; |
| 57 | mpw fw; |
| 58 | mp *y = MP_NEW, *z = MP_NEW; |
| 59 | unsigned i; |
| 60 | |
| 61 | MP_COPY(x); |
| 62 | mp_build(&f, &fw, &fw + 1); |
| 63 | #ifdef DEBUG |
| 64 | MP_PRINT("target", x); |
| 65 | #endif |
| 66 | for (i = 0; i < NPRIME; i++) { |
| 67 | again: |
| 68 | fw = primetab[i]; |
| 69 | mp_div(&y, &z, x, &f); |
| 70 | if (MP_ZEROP(z)) { |
| 71 | #ifdef DEBUG |
| 72 | MP_PRINT("factor", &f); |
| 73 | #endif |
| 74 | MP_DROP(x); |
| 75 | x = y; |
| 76 | y = MP_NEW; |
| 77 | #ifdef DEBUG |
| 78 | MP_PRINT("remainder", x); |
| 79 | #endif |
| 80 | DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i])); |
| 81 | goto again; |
| 82 | } |
| 83 | } |
| 84 | mp_drop(y); |
| 85 | mp_drop(z); |
| 86 | return (x); |
| 87 | } |
| 88 | |
| 89 | #ifdef POLLARDRHO |
| 90 | static mp *pollardrho(mp *x) |
| 91 | { |
| 92 | unsigned i; |
| 93 | mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW; |
| 94 | mpbarrett mb; |
| 95 | |
| 96 | mpbarrett_create(&mb, x); |
| 97 | y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y); |
| 98 | for (i = 0; i < 1024; i++) { |
| 99 | y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y); |
| 100 | z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z); |
| 101 | z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z); |
| 102 | t = mp_sub(t, y, z); |
| 103 | mp_gcd(&g, 0, 0, x, t); |
| 104 | if (!MP_EQ(g, MP_ONE)) |
| 105 | goto done; |
| 106 | } |
| 107 | MP_DROP(g); |
| 108 | g = 0; |
| 109 | |
| 110 | done: |
| 111 | mpbarrett_destroy(&mb); |
| 112 | MP_DROP(y); |
| 113 | MP_DROP(z); |
| 114 | MP_DROP(t); |
| 115 | #ifdef DEBUG |
| 116 | MP_PRINT("rho factor", g); |
| 117 | #endif |
| 118 | return (g); |
| 119 | } |
| 120 | #endif |
| 121 | |
| 122 | static void addfbe(fbe_v *v, mp *B, unsigned long p) |
| 123 | { |
| 124 | fbe e; |
| 125 | unsigned w; |
| 126 | mp *pp = mp_fromulong(MP_NEW, p); |
| 127 | mp *qq = MP_NEW, *rr = MP_NEW; |
| 128 | mp *t; |
| 129 | |
| 130 | w = 0; |
| 131 | qq = MP_COPY(pp); |
| 132 | for (;;) { |
| 133 | w++; |
| 134 | rr = mp_mul(rr, qq, pp); |
| 135 | if (MP_CMP(rr, >=, B)) break; |
| 136 | t = rr; rr = qq; qq = t; |
| 137 | } |
| 138 | e.p = p; |
| 139 | e.e = w; |
| 140 | e.n = qq; |
| 141 | MP_DROP(rr); |
| 142 | MP_DROP(pp); |
| 143 | DA_PUSH(v, e); |
| 144 | } |
| 145 | |
| 146 | static void makefactorbase(mp *x, fbe_v *v) |
| 147 | { |
| 148 | mp *B; |
| 149 | unsigned nb; |
| 150 | unsigned long max; |
| 151 | unsigned long i; |
| 152 | size_t j; |
| 153 | |
| 154 | nb = sqrt(mp_bits(x)) * 1.5; |
| 155 | B = mp_lsl(MP_NEW, MP_ONE, nb); |
| 156 | if (nb < 32) |
| 157 | max = 1 << nb; |
| 158 | else |
| 159 | max = 0xfffffff0; |
| 160 | |
| 161 | addfbe(v, B, 2); |
| 162 | for (i = 3; i < max; i++) { |
| 163 | for (j = 0; j < DA_LEN(v); j++) { |
| 164 | if (i % DA(v)[j].p == 0) |
| 165 | goto composite; |
| 166 | } |
| 167 | addfbe(v, B, i); |
| 168 | composite:; |
| 169 | } |
| 170 | |
| 171 | #ifdef DEBUG |
| 172 | for (j = 0; j < DA_LEN(v); j++) { |
| 173 | printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e); |
| 174 | mp_writefile(DA(v)[j].n, stdout, 10); |
| 175 | putchar('\n'); |
| 176 | } |
| 177 | #endif |
| 178 | } |
| 179 | |
| 180 | static void freefactorbase(fbe_v *v) |
| 181 | { |
| 182 | size_t i; |
| 183 | |
| 184 | for (i = 0; i < DA_LEN(v); i++) |
| 185 | mp_drop(DA(v)[i].n); |
| 186 | DA_DESTROY(v); |
| 187 | } |
| 188 | |
| 189 | static mp *ecm(mp *x) |
| 190 | { |
| 191 | fbe_v fb = DA_INIT; |
| 192 | field *f = field_prime(x); |
| 193 | ec_curve *c; |
| 194 | ec p; |
| 195 | mp *a = MP_NEW, *b = MP_NEW; |
| 196 | mp *g = MP_NEW; |
| 197 | unsigned i; |
| 198 | size_t j; |
| 199 | |
| 200 | makefactorbase(x, &fb); |
| 201 | for (;;) { |
| 202 | a = mprand_range(a, x, rng, 0); |
| 203 | b = mprand_range(b, x, rng, 0); |
| 204 | c = ec_primeproj(f, a, b); |
| 205 | for (i = 0; i < 1024; i++) { |
| 206 | EC_CREATE(&p); |
| 207 | p.x = mprand_range(p.x, x, rng, 0); |
| 208 | p.y = mprand_range(p.y, x, rng, 0); |
| 209 | p.z = MP_ONE; |
| 210 | for (j = 0; j < DA_LEN(&fb); j++) |
| 211 | ec_imul(c, &p, &p, DA(&fb)[j].n); |
| 212 | if (EC_ATINF(&p)) |
| 213 | continue; |
| 214 | mp_gcd(&g, 0, 0, x, p.z); |
| 215 | if (!MP_EQ(g, MP_ONE)) |
| 216 | goto done; |
| 217 | EC_DESTROY(&p); |
| 218 | } |
| 219 | ec_destroycurve(c); |
| 220 | } |
| 221 | |
| 222 | done: |
| 223 | MP_DROP(a); MP_DROP(b); |
| 224 | EC_DESTROY(&p); |
| 225 | ec_destroycurve(c); |
| 226 | F_DESTROY(f); |
| 227 | freefactorbase(&fb); |
| 228 | #ifdef DEBUG |
| 229 | MP_PRINT("ecm factor", g); |
| 230 | #endif |
| 231 | return (g); |
| 232 | } |
| 233 | |
| 234 | static void dofactor(mp *x, mp_v *v) |
| 235 | { |
| 236 | mp *f; |
| 237 | |
| 238 | x = smallfactors(x, v); |
| 239 | if (MP_EQ(x, MP_ONE)) return; |
| 240 | |
| 241 | #ifdef POLLARDRHO |
| 242 | for (;;) { |
| 243 | if (pgen_primep(x, rng)) |
| 244 | goto done; |
| 245 | f = pollardrho(x); |
| 246 | if (!f) break; |
| 247 | dofactor(f, v); |
| 248 | mp_div(&x, 0, x, f); |
| 249 | MP_DROP(f); |
| 250 | } |
| 251 | #endif |
| 252 | |
| 253 | for (;;) { |
| 254 | if (pgen_primep(x, rng)) |
| 255 | goto done; |
| 256 | f = ecm(x); |
| 257 | if (!f) break; |
| 258 | dofactor(f, v); |
| 259 | mp_div(&x, 0, x, f); |
| 260 | MP_DROP(f); |
| 261 | } |
| 262 | |
| 263 | done: |
| 264 | DA_PUSH(v, x); |
| 265 | } |
| 266 | |
| 267 | static int cmpmp(const void *a, const void *b) |
| 268 | { |
| 269 | mp *const *aa = a, *const *bb = b; |
| 270 | return (mp_cmp(*aa, *bb)); |
| 271 | } |
| 272 | |
| 273 | static void factor(mp *x, fact_v *v) |
| 274 | { |
| 275 | fact f; |
| 276 | mp_v fv = DA_INIT; |
| 277 | mp *e = MP_NEW; |
| 278 | size_t i; |
| 279 | |
| 280 | dofactor(x, &fv); |
| 281 | qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp); |
| 282 | for (i = 0; i < DA_LEN(&fv); ) { |
| 283 | f.p = MP_COPY(DA(&fv)[i]); |
| 284 | for (f.e = 1, i++; |
| 285 | i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]); |
| 286 | f.e++, i++) ; |
| 287 | e = mp_fromuint(e, f.e); |
| 288 | f.n = mp_exp(MP_NEW, f.p, e); |
| 289 | f.r = MP_NEW; |
| 290 | mp_div(&f.r, 0, x, f.p); |
| 291 | DA_PUSH(v, f); |
| 292 | } |
| 293 | mp_drop(e); |
| 294 | mpv_free(&fv); |
| 295 | } |
| 296 | |
| 297 | static void freefactors(fact_v *v) |
| 298 | { |
| 299 | size_t i; |
| 300 | |
| 301 | for (i = 0; i < DA_LEN(v); i++) { |
| 302 | mp_drop(DA(v)[i].p); |
| 303 | mp_drop(DA(v)[i].n); |
| 304 | mp_drop(DA(v)[i].r); |
| 305 | } |
| 306 | DA_DESTROY(v); |
| 307 | } |
| 308 | |
| 309 | |
| 310 | /*----- Primitive polynomials ---------------------------------------------*/ |
| 311 | |
| 312 | static int primitivep(fact_v *f, mp *p) |
| 313 | { |
| 314 | gfreduce r; |
| 315 | unsigned i; |
| 316 | mp *x = MP_NEW; |
| 317 | int rc = 1; |
| 318 | |
| 319 | if (!gf_irreduciblep(p)) |
| 320 | return (0); |
| 321 | #ifdef DEBUG |
| 322 | MP_PRINTX("p", p); |
| 323 | #endif |
| 324 | gfreduce_create(&r, p); |
| 325 | for (i = 0; i < DA_LEN(f); i++) { |
| 326 | x = gfreduce_exp(&r, x, MP_TWO, DA(f)[i].r); |
| 327 | #ifdef DEBUG |
| 328 | MP_PRINT(" r", DA(f)[i].r); |
| 329 | MP_PRINTX(" x", x); |
| 330 | #endif |
| 331 | if (MP_EQ(x, MP_ONE)) { |
| 332 | rc = 0; |
| 333 | break; |
| 334 | } |
| 335 | } |
| 336 | gfreduce_destroy(&r); |
| 337 | MP_DROP(x); |
| 338 | return (rc); |
| 339 | } |
| 340 | |
| 341 | static int dofip(fact_v *f, unsigned m, mp **p, unsigned n, unsigned x) |
| 342 | { |
| 343 | unsigned i; |
| 344 | |
| 345 | for (i = n + 1; i < x; i++) { |
| 346 | *p = mp_setbit(*p, *p, i); |
| 347 | if (n ? dofip(f, m, p, n - 1, i) : primitivep(f, *p)) |
| 348 | return (1); |
| 349 | *p = mp_clearbit(*p, *p, i); |
| 350 | } |
| 351 | return (0); |
| 352 | } |
| 353 | |
| 354 | static mp *fip(fact_v *f, unsigned m) |
| 355 | { |
| 356 | mp *p = MP_ONE; |
| 357 | unsigned n; |
| 358 | |
| 359 | p = mp_setbit(p, p, m); |
| 360 | n = 0; |
| 361 | while (!dofip(f, m, &p, n, m)) |
| 362 | n += 2; |
| 363 | return (p); |
| 364 | } |
| 365 | |
| 366 | static void findprim(unsigned m) |
| 367 | { |
| 368 | fact_v f = DA_INIT; |
| 369 | mp *q = mp_lsl(MP_NEW, MP_ONE, m); |
| 370 | mp *p; |
| 371 | unsigned i; |
| 372 | |
| 373 | q = mp_sub(q, q, MP_ONE); |
| 374 | factor(q, &f); |
| 375 | |
| 376 | #ifdef DEBUG |
| 377 | { |
| 378 | size_t i; |
| 379 | for (i = 0; i < DA_LEN(&f); i++) { |
| 380 | mp_writefile(DA(&f)[i].p, stdout, 10); |
| 381 | printf("^%u = ", DA(&f)[i].e); |
| 382 | mp_writefile(DA(&f)[i].n, stdout, 10); |
| 383 | fputs(" (", stdout); |
| 384 | mp_writefile(DA(&f)[i].r, stdout, 10); |
| 385 | fputs(")\n", stdout); |
| 386 | } |
| 387 | } |
| 388 | #endif |
| 389 | |
| 390 | p = fip(&f, m); |
| 391 | MP_PRINTX("p", p); |
| 392 | for (i = m + 1; i--; ) { |
| 393 | if (mp_testbit(p, i)) |
| 394 | printf("%u ", i); |
| 395 | } |
| 396 | putchar('\n'); |
| 397 | mp_drop(p); |
| 398 | mp_drop(q); |
| 399 | freefactors(&f); |
| 400 | } |
| 401 | |
| 402 | int main(int argc, char *argv[]) |
| 403 | { |
| 404 | ego(argv[0]); |
| 405 | rng = fibrand_create(0); |
| 406 | if (argc != 2) { |
| 407 | fprintf(stderr, "Usage: %s M\n", QUIS); |
| 408 | exit(1); |
| 409 | } |
| 410 | findprim(atoi(argv[1])); |
| 411 | return (0); |
| 412 | } |