feanor_math/algorithms/
ec_factor.rs

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

feanor_math/algorithms/
ec_factor.rs
