feanor_math/algorithms/
ec_factor.rs

1use std::sync::atomic::AtomicU64;
2
3use super::int_factor::is_prime_power;
4use crate::algorithms::eea::signed_gcd;
5use crate::algorithms::sqr_mul;
6use crate::computation::*;
7use crate::divisibility::*;
8use crate::homomorphism::Homomorphism;
9use crate::integer::*;
10use crate::ordered::OrderedRingStore;
11use crate::pid::PrincipalIdealRingStore;
12use crate::primitive_int::StaticRing;
13use crate::ring::*;
14use crate::rings::finite::*;
15use crate::rings::zn::*;
16use crate::seq::VectorFn;
17use crate::{MAX_PROBABILISTIC_REPETITIONS, algorithms};
18
19type Point<R> = (El<R>, El<R>, El<R>);
20
21fn square<R>(Zn: &R, x: &El<R>) -> El<R>
22where
23    R: RingStore,
24{
25    let mut result: <<R as RingStore>::Type as RingBase>::Element = Zn.clone_el(x);
26    Zn.square(&mut result);
27    return result;
28}
29
30#[allow(unused)]
31fn point_eq<R>(Zn: &R, P: &Point<R>, Q: &Point<R>) -> bool
32where
33    R: RingStore,
34    R::Type: ZnRing,
35{
36    let factor_quo = if !Zn.is_zero(&Q.0) {
37        if Zn.is_zero(&P.0) {
38            return false;
39        }
40        (&P.0, &Q.0)
41    } else if !Zn.is_zero(&Q.1) {
42        if Zn.is_zero(&P.1) {
43            return false;
44        }
45        (&P.1, &Q.1)
46    } else {
47        assert!(!Zn.is_zero(&Q.2));
48        if Zn.is_zero(&P.2) {
49            return false;
50        }
51        (&P.2, &Q.2)
52    };
53    if !Zn.is_unit(factor_quo.1) {
54        let factor_of_n = signed_gcd(
55            Zn.integer_ring().clone_el(Zn.modulus()),
56            Zn.smallest_positive_lift(Zn.clone_el(factor_quo.1)),
57            Zn.integer_ring(),
58        );
59        let Zn_new = zn_big::Zn::new(
60            BigIntRing::RING,
61            int_cast(
62                Zn.integer_ring().checked_div(Zn.modulus(), &factor_of_n).unwrap(),
63                BigIntRing::RING,
64                Zn.integer_ring(),
65            ),
66        );
67        let red_map = ZnReductionMap::new(Zn, &Zn_new).unwrap();
68        if (Zn_new.is_zero(&red_map.map_ref(&Q.0))
69            && Zn_new.is_zero(&red_map.map_ref(&Q.1))
70            && Zn_new.is_zero(&red_map.map_ref(&Q.2)))
71            || (Zn_new.is_zero(&red_map.map_ref(&P.0))
72                && Zn_new.is_zero(&red_map.map_ref(&P.1))
73                && Zn_new.is_zero(&red_map.map_ref(&P.2)))
74        {
75            if (Zn_new.is_zero(&red_map.map_ref(&P.0))
76                && Zn_new.is_zero(&red_map.map_ref(&P.1))
77                && Zn_new.is_zero(&red_map.map_ref(&P.2)))
78                != (Zn_new.is_zero(&red_map.map_ref(&Q.0))
79                    && Zn_new.is_zero(&red_map.map_ref(&Q.1))
80                    && Zn_new.is_zero(&red_map.map_ref(&Q.2)))
81            {
82                return false;
83            }
84        } else if !point_eq(
85            &Zn_new,
86            &(red_map.map_ref(&P.0), red_map.map_ref(&P.1), red_map.map_ref(&P.2)),
87            &(red_map.map_ref(&Q.0), red_map.map_ref(&Q.1), red_map.map_ref(&Q.2)),
88        ) {
89            return false;
90        }
91
92        let Zn_new = zn_big::Zn::new(
93            BigIntRing::RING,
94            int_cast(factor_of_n, BigIntRing::RING, Zn.integer_ring()),
95        );
96        let red_map = ZnReductionMap::new(Zn, &Zn_new).unwrap();
97        if (Zn_new.is_zero(&red_map.map_ref(&Q.0))
98            && Zn_new.is_zero(&red_map.map_ref(&Q.1))
99            && Zn_new.is_zero(&red_map.map_ref(&Q.2)))
100            || (Zn_new.is_zero(&red_map.map_ref(&P.0))
101                && Zn_new.is_zero(&red_map.map_ref(&P.1))
102                && Zn_new.is_zero(&red_map.map_ref(&P.2)))
103        {
104            if (Zn_new.is_zero(&red_map.map_ref(&P.0))
105                && Zn_new.is_zero(&red_map.map_ref(&P.1))
106                && Zn_new.is_zero(&red_map.map_ref(&P.2)))
107                != (Zn_new.is_zero(&red_map.map_ref(&Q.0))
108                    && Zn_new.is_zero(&red_map.map_ref(&Q.1))
109                    && Zn_new.is_zero(&red_map.map_ref(&Q.2)))
110            {
111                return false;
112            }
113        } else if !point_eq(
114            &Zn_new,
115            &(red_map.map_ref(&P.0), red_map.map_ref(&P.1), red_map.map_ref(&P.2)),
116            &(red_map.map_ref(&Q.0), red_map.map_ref(&Q.1), red_map.map_ref(&Q.2)),
117        ) {
118            return false;
119        }
120        return true;
121    }
122    let factor = Zn.checked_div(factor_quo.0, factor_quo.1).unwrap();
123    if !Zn.is_unit(&factor) {
124        return false;
125    }
126    return Zn.eq_el(&P.0, &Zn.mul_ref(&factor, &Q.0))
127        && Zn.eq_el(&P.1, &Zn.mul_ref(&factor, &Q.1))
128        && Zn.eq_el(&P.2, &Zn.mul_ref(&factor, &Q.2));
129}
130
131#[inline(never)]
132fn edcurve_add<R>(Zn: &R, d: &El<R>, P: Point<R>, Q: &Point<R>) -> Point<R>
133where
134    R: RingStore,
135    R::Type: ZnRing,
136{
137    let (Px, Py, Pz) = P;
138    let (Qx, Qy, Qz) = Q;
139
140    let PxQx = Zn.mul_ref(&Px, Qx);
141    let PyQy = Zn.mul_ref(&Py, Qy);
142    let PzQz = Zn.mul_ref_snd(Pz, Qz);
143
144    let PzQz_sqr = square(Zn, &PzQz);
145    let dPxPyQxQy = Zn.mul_ref_snd(Zn.mul_ref(&PxQx, &PyQy), d);
146
147    let u1 = Zn.add_ref(&PzQz_sqr, &dPxPyQxQy);
148    let u2 = Zn.sub(PzQz_sqr, dPxPyQxQy);
149
150    let result = (
151        Zn.mul_ref_fst(
152            &PzQz,
153            Zn.mul_ref_snd(Zn.add(Zn.mul_ref_snd(Px, Qy), Zn.mul_ref_snd(Py, Qx)), &u2),
154        ),
155        Zn.mul(PzQz, Zn.mul_ref_snd(Zn.sub(PyQy, PxQx), &u1)),
156        Zn.mul(u1, u2),
157    );
158    debug_assert!(is_on_curve(Zn, d, &result));
159    return result;
160}
161
162#[inline(never)]
163fn edcurve_double<R>(Zn: &R, d: &El<R>, P: &Point<R>) -> Point<R>
164where
165    R: RingStore,
166    R::Type: ZnRing,
167{
168    let (Px, Py, Pz) = P;
169
170    let PxPy = Zn.mul_ref(Px, Py);
171    let Px_sqr = square(Zn, Px);
172    let Py_sqr = square(Zn, Py);
173    let Pz_sqr = square(Zn, Pz);
174    let Pz_pow4 = square(Zn, &Pz_sqr);
175    let d_PxPy_sqr = Zn.mul_ref_snd(Zn.mul_ref(&Px_sqr, &Py_sqr), d);
176
177    let u1 = Zn.add_ref(&Pz_pow4, &d_PxPy_sqr);
178    let u2 = Zn.sub(Pz_pow4, d_PxPy_sqr);
179
180    let result = (
181        Zn.mul_ref_fst(&Pz_sqr, Zn.mul_ref_snd(Zn.add_ref(&PxPy, &PxPy), &u2)),
182        Zn.mul_ref_fst(&Pz_sqr, Zn.mul_ref_snd(Zn.sub(Py_sqr, Px_sqr), &u1)),
183        Zn.mul(u1, u2),
184    );
185    debug_assert!(is_on_curve(Zn, d, &result));
186    return result;
187}
188
189fn ec_pow<R>(base: &Point<R>, d: &El<R>, power: &El<BigIntRing>, Zn: &R) -> Point<R>
190where
191    R: RingStore,
192    R::Type: ZnRing,
193{
194    let copy_point = |(x, y, z): &Point<R>| (Zn.clone_el(x), Zn.clone_el(y), Zn.clone_el(z));
195    let ZZ = BigIntRing::RING;
196
197    sqr_mul::generic_pow_shortest_chain_table(
198        copy_point(base),
199        power,
200        ZZ,
201        |P| Ok(edcurve_double(Zn, d, P)),
202        |P, Q| Ok(edcurve_add(Zn, d, copy_point(Q), P)),
203        |P| copy_point(P),
204        (Zn.zero(), Zn.one(), Zn.one()),
205    )
206    .unwrap_or_else(|x| x)
207}
208
209fn is_on_curve<R>(Zn: &R, d: &El<R>, P: &Point<R>) -> bool
210where
211    R: RingStore,
212    R::Type: ZnRing,
213{
214    let (x, y, z) = &P;
215    let x_sqr = square(Zn, x);
216    let y_sqr = square(Zn, y);
217    let z_sqr = square(Zn, z);
218    Zn.eq_el(
219        &Zn.mul_ref_snd(Zn.add_ref(&x_sqr, &y_sqr), &z_sqr),
220        &Zn.add(Zn.mul_ref(&z_sqr, &z_sqr), Zn.mul_ref_fst(d, Zn.mul(x_sqr, y_sqr))),
221    )
222}
223
224const POW_COST_CONSTANT: f64 = 0.1;
225
226/// returns `(ln_B, ln_attempts)`
227fn optimize_parameters(ln_p: f64, ln_n: f64) -> (f64, f64) {
228    let pow_cost_constant = POW_COST_CONSTANT;
229    let ln_cost_per_attempt = |ln_B: f64| ln_B + ln_B.ln() + pow_cost_constant * ln_n.ln();
230    let ln_cost_per_attempt_diff = |ln_B: f64| 1.0 + 1.0 / ln_B;
231    let ln_attempts = |ln_B: f64| {
232        let u = ln_p / ln_B;
233        u * (1.0 + 2.0f64.ln()) * u.ln() - u
234    };
235    let ln_attempts_diff = |ln_B: f64| {
236        let u = ln_p / ln_B;
237        let u_diff = -ln_p / (ln_B * ln_B);
238        u_diff * (1.0 + 2.0f64.ln()) * u.ln() + u * (1.0 + 2.0f64.ln()) * u_diff / u - u_diff
239    };
240    let f = |ln_B: f64| ln_cost_per_attempt(ln_B) - ln_attempts(ln_B);
241    let f_diff = |ln_B: f64| ln_cost_per_attempt_diff(ln_B) - ln_attempts_diff(ln_B);
242
243    let mut ln_B = (ln_p * ln_p.ln()).sqrt();
244    for _ in 0..10 {
245        ln_B = ln_B - f(ln_B) / f_diff(ln_B);
246    }
247    return (ln_B, ln_attempts(ln_B));
248}
249
250/// Optimizes the parameters to find a factor of size roughly `p`; `p` should be at most sqrt(n)
251fn lenstra_ec_factor_base<R, F, Controller>(
252    Zn: R,
253    log2_p: usize,
254    mut rng: F,
255    controller: Controller,
256) -> Result<Option<El<<R::Type as ZnRing>::IntegerRing>>, Controller::Abort>
257where
258    R: RingStore + Copy + Send + Sync,
259    El<R>: Send,
260    R::Type: ZnRing + DivisibilityRing,
261    F: FnMut() -> u64 + Send,
262    Controller: ComputationController,
263{
264    controller.run_computation(
265        format_args!(
266            "ec_factor(log2(n)={}, log2(p)={})",
267            Zn.integer_ring().abs_log2_ceil(Zn.modulus()).unwrap(),
268            log2_p
269        ),
270        |controller| {
271            let ZZ = BigIntRing::RING;
272            assert!(ZZ.is_leq(&ZZ.power_of_two(log2_p * 2), &Zn.size(&ZZ).unwrap()));
273            let log2_n = ZZ.abs_log2_ceil(&Zn.size(&ZZ).unwrap()).unwrap();
274            let ln_p = log2_p as f64 * 2.0f64.ln();
275            let (ln_B, ln_attempts) = optimize_parameters(ln_p, log2_n as f64 * 2.0f64.ln());
276            // after this many random curves, we expect to have found a factor with high
277            // probability, unless there is no factor of size about `log2_size`
278            let attempts = ln_attempts.exp() as usize;
279            log_progress!(controller, "(attempts={})", attempts);
280
281            let log2_B = ln_B / 2.0f64.ln();
282            assert!(log2_B <= i128::MAX as f64);
283
284            let primes =
285                algorithms::erathostenes::enumerate_primes(&StaticRing::<i128>::RING, &(1i128 << (log2_B as u64)));
286            let power_factorization = primes
287                .iter()
288                .map(|p| {
289                    (
290                        *p,
291                        log2_B.ceil() as usize / StaticRing::<i128>::RING.abs_log2_ceil(p).unwrap(),
292                    )
293                })
294                .collect::<Vec<_>>();
295            let power = ZZ.prod(
296                power_factorization
297                    .iter()
298                    .map(|(p, e)| ZZ.pow(ZZ.coerce(&StaticRing::<i128>::RING, *p), *e)),
299            );
300            let power_ref = &power;
301
302            let computation = ShortCircuitingComputation::new();
303
304            let base_rng_value = rng();
305            let rng_seed = AtomicU64::new(1);
306            let rng_seed_ref = &rng_seed;
307
308            computation
309                .handle(controller.clone())
310                .join_many((0..attempts).map_fn(move |_| {
311                    move |handle: ShortCircuitingComputationHandle<_, _>| {
312                        let mut rng = oorandom::Rand64::new(
313                            ((rng_seed_ref.fetch_add(1, std::sync::atomic::Ordering::Relaxed) as u128) << 64)
314                                | base_rng_value as u128,
315                        );
316                        let (x, y) = (
317                            Zn.random_element(|| rng.rand_u64()),
318                            Zn.random_element(|| rng.rand_u64()),
319                        );
320                        let (x_sqr, y_sqr) = (square(&Zn, &x), square(&Zn, &y));
321                        if let Some(d) =
322                            Zn.checked_div(&Zn.sub(Zn.add_ref(&x_sqr, &y_sqr), Zn.one()), &Zn.mul(x_sqr, y_sqr))
323                        {
324                            let P = (x, y, Zn.one());
325                            debug_assert!(is_on_curve(&Zn, &d, &P));
326                            let result = ec_pow(&P, &d, power_ref, &Zn).0;
327                            if !Zn.is_unit(&result) && !Zn.is_zero(&result) {
328                                return Ok(Some(result));
329                            }
330                        }
331                        log_progress!(handle, ".");
332                        checkpoint!(handle);
333                        return Ok(None);
334                    }
335                }));
336
337            if let Some(result) = computation.finish()? {
338                return Ok(Some(
339                    Zn.integer_ring()
340                        .ideal_gen(&Zn.smallest_positive_lift(result), Zn.modulus()),
341                ));
342            } else {
343                log_progress!(controller, "(no_factor)");
344                return Ok(None);
345            }
346        },
347    )
348}
349
350/// Given `Z/nZ`, tries to find a factor of `n` of size at most `2^min_factor_bound_log2`.
351/// If such a factor exists, the function is likely to successfully find it (with probability
352/// `1 - c^repetitions`, for a constant `0 < c < 1`). Otherwise, it is likely to return `None`.
353///
354/// Note that the returned value can be any nontrivial factor of `n`, and does not have to
355/// be bounded by `2^min_factor_bound_log2`. The function is more likely to find small factors,
356/// but can, in rare cases, find other factors as well.
357///
358/// # Explanation of logging output
359///
360/// If the passed computation controller accepts logging, it will receive the following symbols:
361///  - `c(m)` means that `m` random curves will be tried using the current parameters
362///  - `.` means an elliptic curve was tried and did not yield a factor of `n`
363#[stability::unstable(feature = "enable")]
364pub fn lenstra_ec_factor_small<R, Controller>(
365    Zn: R,
366    min_factor_bound_log2: usize,
367    repetitions: usize,
368    controller: Controller,
369) -> Result<Option<El<<R::Type as ZnRing>::IntegerRing>>, Controller::Abort>
370where
371    R: ZnRingStore + DivisibilityRingStore + Copy + Send + Sync,
372    El<R>: Send,
373    R::Type: ZnRing + DivisibilityRing,
374    Controller: ComputationController,
375{
376    assert!(!algorithms::miller_rabin::is_prime_base(Zn, 10));
377    assert!(is_prime_power(Zn.integer_ring(), Zn.modulus()).is_none());
378    let mut rng = oorandom::Rand64::new(Zn.integer_ring().default_hash(Zn.modulus()) as u128);
379
380    for log2_size in (16..min_factor_bound_log2).step_by(8) {
381        if let Some(factor) = lenstra_ec_factor_base(Zn, log2_size, || rng.rand_u64(), controller.clone())? {
382            return Ok(Some(factor));
383        }
384    }
385    for _ in 0..repetitions {
386        if let Some(factor) = lenstra_ec_factor_base(Zn, min_factor_bound_log2, || rng.rand_u64(), controller.clone())?
387        {
388            return Ok(Some(factor));
389        }
390    }
391    return Ok(None);
392}
393
394#[stability::unstable(feature = "enable")]
395pub fn lenstra_ec_factor<R, Controller>(
396    Zn: R,
397    controller: Controller,
398) -> Result<El<<R::Type as ZnRing>::IntegerRing>, Controller::Abort>
399where
400    R: ZnRingStore + DivisibilityRingStore + Copy + Send + Sync,
401    El<R>: Send,
402    R::Type: ZnRing + DivisibilityRing,
403    Controller: ComputationController,
404{
405    assert!(!algorithms::miller_rabin::is_prime_base(Zn, 10));
406    assert!(is_prime_power(Zn.integer_ring(), Zn.modulus()).is_none());
407    let ZZ = BigIntRing::RING;
408    let log2_N = ZZ.abs_log2_floor(&Zn.size(&ZZ).unwrap()).unwrap();
409    let mut rng = oorandom::Rand64::new(Zn.integer_ring().default_hash(Zn.modulus()) as u128);
410
411    // we first try to find smaller factors
412    for log2_size in (16..(log2_N / 2)).step_by(8) {
413        if let Some(factor) = lenstra_ec_factor_base(Zn, log2_size, || rng.rand_u64(), controller.clone())? {
414            return Ok(factor);
415        }
416    }
417    // this is now the general case
418    for _ in 0..MAX_PROBABILISTIC_REPETITIONS {
419        if let Some(factor) = lenstra_ec_factor_base(Zn, log2_N / 2, || rng.rand_u64(), controller.clone())? {
420            return Ok(factor);
421        }
422    }
423    unreachable!()
424}
425
426#[cfg(test)]
427use std::time::Instant;
428
429#[cfg(test)]
430use test::Bencher;
431
432#[cfg(test)]
433use crate::rings::rust_bigint::*;
434#[cfg(test)]
435use crate::rings::zn::zn_64::Zn;
436
437#[test]
438fn test_ec_factor() {
439    let n = 65537 * 65539;
440    let actual = lenstra_ec_factor(&Zn::new(n as u64), TEST_LOG_PROGRESS).unwrap_or_else(no_error);
441    assert!(actual != 1 && actual != n && n % actual == 0);
442}
443
444#[bench]
445fn bench_ec_factor_mersenne_number_58(bencher: &mut Bencher) {
446    let bits = 58;
447    let n = ((1i64 << bits) + 1) / 5;
448    let ring = Zn::new(n as u64);
449
450    bencher.iter(|| {
451        let p = lenstra_ec_factor(&ring, TEST_LOG_PROGRESS).unwrap_or_else(no_error);
452        assert!(n > 0 && n != 1 && n != p);
453        assert!(n % p == 0);
454    });
455}
456
457#[test]
458#[ignore]
459fn test_ec_factor_large() {
460    let ZZbig = BigIntRing::RING;
461    #[cfg(not(feature = "parallel"))]
462    let controller = TEST_LOG_PROGRESS;
463    #[cfg(feature = "parallel")]
464    let controller = RunMultithreadedLogProgress;
465
466    let n: i128 = 1073741827 * 71316922984999;
467
468    let begin = Instant::now();
469    let p = StaticRing::<i128>::RING.coerce(
470        &ZZbig,
471        lenstra_ec_factor(
472            &zn_big::Zn::new(&ZZbig, ZZbig.coerce(&StaticRing::<i128>::RING, n)),
473            controller.clone(),
474        )
475        .unwrap_or_else(no_error),
476    );
477    let end = Instant::now();
478    println!("Done in {} ms", (end - begin).as_millis());
479    assert!(p == 1073741827 || p == 71316922984999);
480
481    let n: i128 = 1152921504606847009 * 2305843009213693967;
482
483    let begin = Instant::now();
484    let p = StaticRing::<i128>::RING.coerce(
485        &ZZbig,
486        lenstra_ec_factor(
487            &zn_big::Zn::new(&ZZbig, ZZbig.coerce(&StaticRing::<i128>::RING, n)),
488            controller,
489        )
490        .unwrap_or_else(no_error),
491    );
492    let end = Instant::now();
493    println!("Done in {} ms", (end - begin).as_millis());
494    assert!(p == 1152921504606847009 || p == 2305843009213693967);
495}
496
497#[test]
498#[ignore]
499fn test_compute_partial_factorization() {
500    let ZZbig = BigIntRing::RING;
501    let n = int_cast(
502        RustBigintRing::RING.get_ring().parse("5164499756173817179311838344006023748659411585658447025661318713081295244033682389259290706560275662871806343945494986751", 10).unwrap(),
503        ZZbig,
504        RustBigintRing::RING
505    );
506
507    let Zn = zn_big::Zn::new(ZZbig, ZZbig.clone_el(&n));
508    let begin = Instant::now();
509    let factor = lenstra_ec_factor_small(&Zn, 50, 1, TEST_LOG_PROGRESS)
510        .unwrap_or_else(no_error)
511        .unwrap();
512    let end = Instant::now();
513    println!("Done in {} ms", (end - begin).as_millis());
514    ZZbig.println(&factor);
515    assert!(!ZZbig.is_one(&factor));
516    assert!(!ZZbig.eq_el(&factor, &n));
517    assert!(ZZbig.divides(&n, &factor));
518}
feanor_math/algorithms/ec_factor.rs

feanor_math/algorithms/
ec_factor.rs
