| 1 | #include <math.h> |
| 2 | |
| 3 | #include "ec.h" |
| 4 | #include "fibrand.h" |
| 5 | #include "field.h" |
| 6 | #include "mpbarrett.h" |
| 7 | #include "mpint.h" |
| 8 | #include "mprand.h" |
| 9 | #include "pgen.h" |
| 10 | #include "primetab.h" |
| 11 | |
| 12 | #include "factor.h" |
| 13 | |
| 14 | static grand *rng = 0; |
| 15 | |
| 16 | DA_DECL(mp_v, mp *); |
| 17 | |
| 18 | typedef struct fbe { |
| 19 | unsigned long p; |
| 20 | unsigned e; |
| 21 | mp *n; |
| 22 | } fbe; |
| 23 | DA_DECL(fbe_v, fbe); |
| 24 | |
| 25 | static void mpv_free(mp_v *v) |
| 26 | { |
| 27 | size_t i; |
| 28 | |
| 29 | for (i = 0; i < DA_LEN(v); i++) |
| 30 | mp_drop(DA(v)[i]); |
| 31 | DA_DESTROY(v); |
| 32 | } |
| 33 | |
| 34 | static mp *smallfactors(mp *x, mp_v *v) |
| 35 | { |
| 36 | mp f; |
| 37 | mpw fw; |
| 38 | mp *y = MP_NEW, *z = MP_NEW; |
| 39 | unsigned i; |
| 40 | |
| 41 | MP_COPY(x); |
| 42 | mp_build(&f, &fw, &fw + 1); |
| 43 | #ifdef DEBUG |
| 44 | MP_PRINT("target", x); |
| 45 | #endif |
| 46 | for (i = 0; i < NPRIME; i++) { |
| 47 | again: |
| 48 | fw = primetab[i]; |
| 49 | mp_div(&y, &z, x, &f); |
| 50 | if (MP_ZEROP(z)) { |
| 51 | #ifdef DEBUG |
| 52 | MP_PRINT("factor", &f); |
| 53 | #endif |
| 54 | MP_DROP(x); |
| 55 | x = y; |
| 56 | y = MP_NEW; |
| 57 | #ifdef DEBUG |
| 58 | MP_PRINT("remainder", x); |
| 59 | #endif |
| 60 | DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i])); |
| 61 | goto again; |
| 62 | } |
| 63 | } |
| 64 | mp_drop(y); |
| 65 | mp_drop(z); |
| 66 | return (x); |
| 67 | } |
| 68 | |
| 69 | #ifdef POLLARDRHO |
| 70 | static mp *pollardrho(mp *x) |
| 71 | { |
| 72 | unsigned i; |
| 73 | mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW; |
| 74 | mpbarrett mb; |
| 75 | |
| 76 | mpbarrett_create(&mb, x); |
| 77 | y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y); |
| 78 | for (i = 0; i < 1024; i++) { |
| 79 | y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y); |
| 80 | z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z); |
| 81 | z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z); |
| 82 | t = mp_sub(t, y, z); |
| 83 | mp_gcd(&g, 0, 0, x, t); |
| 84 | if (!MP_EQ(g, MP_ONE)) |
| 85 | goto done; |
| 86 | } |
| 87 | MP_DROP(g); |
| 88 | g = 0; |
| 89 | |
| 90 | done: |
| 91 | mpbarrett_destroy(&mb); |
| 92 | MP_DROP(y); |
| 93 | MP_DROP(z); |
| 94 | MP_DROP(t); |
| 95 | #ifdef DEBUG |
| 96 | MP_PRINT("rho factor", g); |
| 97 | #endif |
| 98 | return (g); |
| 99 | } |
| 100 | #endif |
| 101 | |
| 102 | static void addfbe(fbe_v *v, mp *B, unsigned long p) |
| 103 | { |
| 104 | fbe e; |
| 105 | unsigned w; |
| 106 | mp *pp = mp_fromulong(MP_NEW, p); |
| 107 | mp *qq = MP_NEW, *rr = MP_NEW; |
| 108 | mp *t; |
| 109 | |
| 110 | w = 0; |
| 111 | qq = MP_COPY(pp); |
| 112 | for (;;) { |
| 113 | w++; |
| 114 | rr = mp_mul(rr, qq, pp); |
| 115 | if (MP_CMP(rr, >=, B)) break; |
| 116 | t = rr; rr = qq; qq = t; |
| 117 | } |
| 118 | e.p = p; |
| 119 | e.e = w; |
| 120 | e.n = qq; |
| 121 | MP_DROP(rr); |
| 122 | MP_DROP(pp); |
| 123 | DA_PUSH(v, e); |
| 124 | } |
| 125 | |
| 126 | static void makefactorbase(mp *x, fbe_v *v) |
| 127 | { |
| 128 | mp *B; |
| 129 | unsigned nb; |
| 130 | unsigned long max; |
| 131 | unsigned long i; |
| 132 | size_t j; |
| 133 | |
| 134 | nb = sqrt(mp_bits(x)) * 1.5; |
| 135 | B = mp_lsl(MP_NEW, MP_ONE, nb); |
| 136 | if (nb < 32) |
| 137 | max = 1 << nb; |
| 138 | else |
| 139 | max = 0xfffffff0; |
| 140 | |
| 141 | addfbe(v, B, 2); |
| 142 | for (i = 3; i < max; i++) { |
| 143 | for (j = 0; j < DA_LEN(v); j++) { |
| 144 | if (i % DA(v)[j].p == 0) |
| 145 | goto composite; |
| 146 | } |
| 147 | addfbe(v, B, i); |
| 148 | composite:; |
| 149 | } |
| 150 | |
| 151 | #ifdef DEBUG |
| 152 | for (j = 0; j < DA_LEN(v); j++) { |
| 153 | printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e); |
| 154 | mp_writefile(DA(v)[j].n, stdout, 10); |
| 155 | putchar('\n'); |
| 156 | } |
| 157 | #endif |
| 158 | } |
| 159 | |
| 160 | static void freefactorbase(fbe_v *v) |
| 161 | { |
| 162 | size_t i; |
| 163 | |
| 164 | for (i = 0; i < DA_LEN(v); i++) |
| 165 | mp_drop(DA(v)[i].n); |
| 166 | DA_DESTROY(v); |
| 167 | } |
| 168 | |
| 169 | static mp *ecm(mp *x) |
| 170 | { |
| 171 | fbe_v fb = DA_INIT; |
| 172 | field *f = field_prime(x); |
| 173 | ec_curve *c; |
| 174 | ec p; |
| 175 | mp *a = MP_NEW, *b = MP_NEW; |
| 176 | mp *g = MP_NEW; |
| 177 | unsigned i; |
| 178 | size_t j; |
| 179 | |
| 180 | makefactorbase(x, &fb); |
| 181 | for (;;) { |
| 182 | a = mprand_range(a, x, rng, 0); |
| 183 | b = mprand_range(b, x, rng, 0); |
| 184 | c = ec_primeproj(f, a, b); |
| 185 | for (i = 0; i < 1024; i++) { |
| 186 | EC_CREATE(&p); |
| 187 | p.x = mprand_range(p.x, x, rng, 0); |
| 188 | p.y = mprand_range(p.y, x, rng, 0); |
| 189 | p.z = MP_ONE; |
| 190 | for (j = 0; j < DA_LEN(&fb); j++) |
| 191 | ec_imul(c, &p, &p, DA(&fb)[j].n); |
| 192 | if (EC_ATINF(&p)) |
| 193 | continue; |
| 194 | mp_gcd(&g, 0, 0, x, p.z); |
| 195 | if (!MP_EQ(g, MP_ONE)) |
| 196 | goto done; |
| 197 | EC_DESTROY(&p); |
| 198 | } |
| 199 | ec_destroycurve(c); |
| 200 | } |
| 201 | |
| 202 | done: |
| 203 | MP_DROP(a); MP_DROP(b); |
| 204 | EC_DESTROY(&p); |
| 205 | ec_destroycurve(c); |
| 206 | F_DESTROY(f); |
| 207 | freefactorbase(&fb); |
| 208 | #ifdef DEBUG |
| 209 | MP_PRINT("ecm factor", g); |
| 210 | #endif |
| 211 | return (g); |
| 212 | } |
| 213 | |
| 214 | static void dofactor(mp *x, mp_v *v) |
| 215 | { |
| 216 | mp *f; |
| 217 | |
| 218 | x = smallfactors(x, v); |
| 219 | if (MP_EQ(x, MP_ONE)) return; |
| 220 | |
| 221 | #ifdef POLLARDRHO |
| 222 | for (;;) { |
| 223 | if (pgen_primep(x, rng)) |
| 224 | goto done; |
| 225 | f = pollardrho(x); |
| 226 | if (!f) break; |
| 227 | dofactor(f, v); |
| 228 | mp_div(&x, 0, x, f); |
| 229 | MP_DROP(f); |
| 230 | } |
| 231 | #endif |
| 232 | |
| 233 | for (;;) { |
| 234 | if (pgen_primep(x, rng)) |
| 235 | goto done; |
| 236 | f = ecm(x); |
| 237 | if (!f) break; |
| 238 | dofactor(f, v); |
| 239 | mp_div(&x, 0, x, f); |
| 240 | MP_DROP(f); |
| 241 | } |
| 242 | |
| 243 | done: |
| 244 | DA_PUSH(v, x); |
| 245 | } |
| 246 | |
| 247 | static int cmpmp(const void *a, const void *b) |
| 248 | { |
| 249 | mp *const *aa = a, *const *bb = b; |
| 250 | return (mp_cmp(*aa, *bb)); |
| 251 | } |
| 252 | |
| 253 | void factor(mp *x, fact_v *v) |
| 254 | { |
| 255 | fact f; |
| 256 | mp_v fv = DA_INIT; |
| 257 | mp *e = MP_NEW; |
| 258 | size_t i; |
| 259 | |
| 260 | if (!rng) rng = fibrand_create(0); |
| 261 | dofactor(x, &fv); |
| 262 | qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp); |
| 263 | for (i = 0; i < DA_LEN(&fv); ) { |
| 264 | f.p = MP_COPY(DA(&fv)[i]); |
| 265 | for (f.e = 1, i++; |
| 266 | i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]); |
| 267 | f.e++, i++) ; |
| 268 | e = mp_fromuint(e, f.e); |
| 269 | f.n = mp_exp(MP_NEW, f.p, e); |
| 270 | f.r = MP_NEW; |
| 271 | mp_div(&f.r, 0, x, f.p); |
| 272 | DA_PUSH(v, f); |
| 273 | } |
| 274 | mp_drop(e); |
| 275 | mpv_free(&fv); |
| 276 | } |
| 277 | |
| 278 | void freefactors(fact_v *v) |
| 279 | { |
| 280 | size_t i; |
| 281 | |
| 282 | for (i = 0; i < DA_LEN(v); i++) { |
| 283 | mp_drop(DA(v)[i].p); |
| 284 | mp_drop(DA(v)[i].n); |
| 285 | mp_drop(DA(v)[i].r); |
| 286 | } |
| 287 | DA_DESTROY(v); |
| 288 | } |