pell_equation/
lib.rs

1#![doc = include_str!("../README.md")]
2#![warn(clippy::pedantic)]
3#![allow(clippy::many_single_char_names)]
4
5mod continued_fraction;
6
7use continued_fraction::ContinuedFractionIterator;
8use rug::{integer::IntegerExt64, Complete};
9type Z = rug::Integer;
10
11#[non_exhaustive]
12#[derive(Debug, thiserror::Error)]
13/// Error type
14pub enum Error {
15    /// Input value is out of domain
16    #[error("Input value is out of domain")]
17    OutOfDomain,
18}
19
20/// Fundamental solution of `x^2 - d*y^2 = ±1`
21#[derive(Debug, Clone, PartialEq, Eq)]
22#[non_exhaustive]
23pub enum Solution {
24    /// Fundamental solution of `x^2 - d*y^2 = -1`
25    Negative(Z, Z),
26    /// Fundamental solution of `x^2 - d*y^2 = 1`
27    Positive(Z, Z),
28    /// Not exsist nontirivial solution
29    NotExist,
30    /// Result digits reached limit
31    ReachedLimit,
32}
33
34macro_rules! check_limit {
35    ($target:expr, $limit:ident) => {
36        if let Some(limit) = $limit {
37            if $target.significant_bits_64() > limit {
38                return Solution::ReachedLimit;
39            }
40        }
41    };
42}
43macro_rules! check_limit2 {
44    ($target1:expr, $target2:expr, $limit:ident) => {
45        if let Some(limit) = $limit {
46            if $target1.significant_bits_64() > limit || $target2.significant_bits_64() > limit {
47                return Solution::ReachedLimit;
48            }
49        }
50    };
51}
52macro_rules! gen_x_neg {
53    ($y:ident, $d:ident, $limit:ident) => {{
54        let x2 = Z::from($d * $y.square_ref().complete() - 1);
55        debug_assert!(x2.is_perfect_square());
56        let x = x2.sqrt();
57        check_limit2!(x, $y, $limit);
58        Solution::Negative(x, $y)
59    }};
60}
61macro_rules! gen_x_pos {
62    ($y:ident, $d:ident, $limit:ident) => {{
63        let x2 = Z::from($d * $y.square_ref().complete() + 1);
64        debug_assert!(x2.is_perfect_square());
65        let x = x2.sqrt();
66        check_limit2!(x, $y, $limit);
67        Solution::Positive(x, $y)
68    }};
69}
70
71/// Calculate continued fraction of √d
72///
73/// Calculate [simple continued fraction](https://en.wikipedia.org/wiki/Simple_continued_fraction)
74/// of √d.\
75/// ex : √2 = [1; 2, 2, 2, ...]
76/// ```
77/// use rug::Integer;
78/// let v = pell_equation::continued_fraction_of_sqrt(Integer::from(2)).unwrap();
79/// assert_eq!(v, vec![Integer::from(1), Integer::from(2)]);
80/// ```
81/// # Errors
82/// If d is negative returns `Err(Error::OutOfDomain)`.  
83pub fn continued_fraction_of_sqrt(d: Z) -> Result<Vec<Z>, Error> {
84    let iter = ContinuedFractionIterator::new(d)?;
85    Ok(iter.collect())
86}
87
88fn next_pq_small(d: i64, sd: i64, p: i64, q: i64) -> (bool, i64, i64, i64) {
89    let int1 = (sd + p) / q;
90    let p1 = int1 * q - p;
91    let p2 = p1 + q;
92    let norm1 = (i128::from(d) - i128::from(p1) * i128::from(p1)).abs();
93    let norm2 = (i128::from(d) - i128::from(p2) * i128::from(p2)).abs();
94    let epsiron = norm1 < norm2;
95    if epsiron {
96        debug_assert_eq!(norm1 % i128::from(q), 0);
97        let q1 = i64::try_from(norm1 / i128::from(q)).unwrap();
98        (epsiron, int1, p1, q1)
99    } else {
100        debug_assert_eq!(norm2 % i128::from(q), 0);
101        let q2 = i64::try_from(norm2 / i128::from(q)).unwrap();
102        (epsiron, int1 + 1, p2, q2)
103    }
104}
105
106fn solve_pell_aux_small(d: i64, limit: Option<u64>) -> Solution {
107    // http://www.numbertheory.org/php/nscf_pell_short.html
108    // http://www.numbertheory.org/PDFS/nscf_midpoint.pdf
109    let sd = d.isqrt();
110    let mut q_old = None;
111    let mut p = 0i64;
112    let mut q = 1i64;
113    let mut b_old = Z::ONE.clone();
114    let mut b_cur = Z::ZERO;
115    let mut e_old = true;
116    loop {
117        let (e_cur, int, p_new, q_new) = next_pq_small(d, sd, p, q);
118        let b_new = if e_old {
119            (int * &b_cur + &b_old).complete()
120        } else {
121            (int * &b_cur - &b_old).complete()
122        };
123        check_limit!(b_new, limit);
124
125        if p == p_new {
126            let t = if e_old { b_old } else { -b_old };
127            let b = b_cur * (b_new + t);
128            return gen_x_pos!(b, d, limit);
129        } else if q == q_new {
130            return if e_cur {
131                let b = b_new.square_ref() + b_cur.square();
132                gen_x_neg!(b, d, limit)
133            } else {
134                let b = b_new.square_ref() - b_cur.square();
135                gen_x_pos!(b, d, limit)
136            };
137        } else if let Some(q_old) = q_old {
138            if q + q_old == p && !e_old {
139                let b = Z::from(2 * b_cur.square() - &b_new * &b_old);
140                return gen_x_neg!(b, d, limit);
141            }
142        }
143
144        q_old = if q % 2 == 0 { Some(q / 2) } else { None };
145        p = p_new;
146        q = q_new;
147        b_old = b_cur;
148        b_cur = b_new;
149        e_old = e_cur;
150    }
151}
152
153fn next_pq_large(d: &Z, sd: &Z, p: &Z, q: &Z) -> (bool, Z, Z, Z) {
154    let int1 = (sd + p).complete() / q;
155    let p1 = (&int1 * q - p).complete();
156    let p2 = (&p1 + q).complete();
157    let norm1 = (d - p1.square_ref()).complete().abs();
158    let norm2 = (d - p2.square_ref()).complete().abs();
159    let epsiron = norm1 < norm2;
160    if epsiron {
161        debug_assert!(norm1.is_divisible(q));
162        let q1 = norm1.div_exact(q);
163        (epsiron, int1, p1, q1)
164    } else {
165        debug_assert!(norm2.is_divisible(q));
166        let q2 = norm2.div_exact(q);
167        (epsiron, int1 + 1, p2, q2)
168    }
169}
170
171fn solve_pell_aux_large(d: &Z, limit: Option<u64>) -> Solution {
172    // http://www.numbertheory.org/php/nscf_pell_short.html
173    // http://www.numbertheory.org/PDFS/nscf_midpoint.pdf
174    let sd = d.sqrt_ref().complete();
175    let mut q_old = None;
176    let mut p = Z::ZERO;
177    let mut q = Z::ONE.clone();
178    let mut b_old = Z::ONE.clone();
179    let mut b_cur = Z::ZERO;
180    let mut e_old = true;
181    loop {
182        let (e_cur, int, p_new, q_new) = next_pq_large(d, &sd, &p, &q);
183        let b_new = if e_old {
184            (&int * &b_cur + &b_old).complete()
185        } else {
186            (&int * &b_cur - &b_old).complete()
187        };
188        check_limit!(b_new, limit);
189
190        if p == p_new {
191            let t = if e_old { b_old } else { -b_old };
192            let b = b_cur * (b_new + t);
193            return gen_x_pos!(b, d, limit);
194        } else if q == q_new {
195            return if e_cur {
196                let b = b_new.square_ref() + b_cur.square();
197                gen_x_neg!(b, d, limit)
198            } else {
199                let b = b_new.square_ref() - b_cur.square();
200                gen_x_pos!(b, d, limit)
201            };
202        } else if let Some(q_old) = q_old {
203            if &q + q_old == p && !e_old {
204                let b = Z::from(2 * b_cur.square() - &b_new * &b_old);
205                return gen_x_neg!(b, d, limit);
206            }
207        }
208
209        q_old = if q.is_divisible_u(2) {
210            Some(q.div_exact_u(2))
211        } else {
212            None
213        };
214        p = p_new;
215        q = q_new;
216        b_old = b_cur;
217        b_cur = b_new;
218        e_old = e_cur;
219    }
220}
221
222fn solve_pell_aux(d: &Z, limit: Option<u64>) -> Solution {
223    if d.is_negative() || d.is_perfect_square() {
224        return Solution::NotExist;
225    }
226    if let Some(d) = d.to_i64() {
227        solve_pell_aux_small(d, limit)
228    } else {
229        solve_pell_aux_large(d, limit)
230    }
231}
232
233/// Calculate fundamental solution of `x^2 - d*y^2 = ±1`
234///
235/// If d is negative or perfect square returns `NotExist`.\
236/// If `x^2 - d*y^2 = -1` has nontrivial solution, returns its fundamental solution.\
237/// Otherwise returns fundamental solution of `x^2 - d*y^2 = 1`.
238/// ```
239/// use rug::Integer;
240/// let v = pell_equation::solve_pell(&Integer::from(2), None);
241/// assert_eq!(v, pell_equation::Solution::Negative(Integer::from(1), Integer::from(1)));
242/// let w = pell_equation::solve_pell(&Integer::from(3), None);
243/// assert_eq!(w, pell_equation::Solution::Positive(Integer::from(2), Integer::from(1)));
244/// let v = pell_equation::solve_pell(&Integer::from(109), Some(23));
245/// assert_eq!(v, pell_equation::Solution::ReachedLimit);
246/// let w = pell_equation::solve_pell(&Integer::from(109), Some(24));
247/// assert_eq!(w, pell_equation::Solution::Negative(Integer::from(8890182), Integer::from(851525)));
248/// ```
249#[must_use]
250pub fn solve_pell(d: &Z, limit: Option<u64>) -> Solution {
251    solve_pell_aux(d, limit)
252}
253
254/// Calculate fundamental solution of `x^2 - d*y^2 = -1`
255///
256/// If `x^2 - d*y^2 = -1` has nontrivial solution, returns its fundamental solution.\
257/// Otherwise returns `NotExist`.
258/// ```
259/// use rug::Integer;
260/// use pell_equation::Solution;
261/// let v = pell_equation::solve_pell_negative(&Integer::from(2), None);
262/// assert_eq!(v, Solution::Negative(Integer::from(1), Integer::from(1)));
263/// let w = pell_equation::solve_pell_negative(&Integer::from(3), None);
264/// assert_eq!(w, Solution::NotExist);
265/// let v = pell_equation::solve_pell_negative(&Integer::from(109), Some(23));
266/// assert_eq!(v, Solution::ReachedLimit);
267/// let w = pell_equation::solve_pell_negative(&Integer::from(109), Some(24));
268/// assert_eq!(w, Solution::Negative(Integer::from(8890182), Integer::from(851525)));
269/// ```
270#[must_use]
271pub fn solve_pell_negative(d: &Z, limit: Option<u64>) -> Solution {
272    match solve_pell_aux(d, limit) {
273        Solution::Positive(_, _) => Solution::NotExist,
274        x => x,
275    }
276}
277
278/// Calculate fundamental solution of `x^2 - d*y^2 = 1`
279///
280/// If `x^2 - d*y^2 = 1` has nontrivial solution, returns its fundamental solution.\
281/// Otherwise returns `NotExist`.
282/// ```
283/// use rug::Integer;
284/// use pell_equation::Solution;
285/// let v = pell_equation::solve_pell_positive(&Integer::from(2), None);
286/// assert_eq!(v, Solution::Positive(Integer::from(3), Integer::from(2)));
287/// let w = pell_equation::solve_pell_positive(&Integer::from(3), None);
288/// assert_eq!(w, Solution::Positive(Integer::from(2), Integer::from(1)));
289/// let v = pell_equation::solve_pell_positive(&Integer::from(109), Some(47));
290/// assert_eq!(v, Solution::ReachedLimit);
291/// let w = pell_equation::solve_pell_positive(&Integer::from(109), Some(48));
292/// assert_eq!(w, Solution::Positive(
293///     Integer::from(158070671986249i64),
294///     Integer::from(15140424455100i64))
295/// );
296/// ```
297#[must_use]
298pub fn solve_pell_positive(d: &Z, limit: Option<u64>) -> Solution {
299    match solve_pell(d, limit) {
300        Solution::Negative(x, y) => {
301            let y2 = Z::from(2 * (&x * &y).complete());
302            let x2 = x.square() + y.square() * d;
303            check_limit2!(x2, y2, limit);
304            Solution::Positive(x2, y2)
305        }
306        x => x,
307    }
308}
309
310#[cfg(test)]
311mod tests {
312    use super::*;
313    fn to_z(v: &[i32]) -> Vec<Z> {
314        v.iter().map(|x| Z::from(*x)).collect()
315    }
316    #[test]
317    fn test_continued_fraction_of_sqrt() {
318        // https://planetmath.org/tableofcontinuedfractionsofsqrtnfor1n102
319        let continued_fractions = vec![
320            (2, vec![1, 2]),
321            (3, vec![1, 1, 2]),
322            (5, vec![2, 4]),
323            (6, vec![2, 2, 4]),
324            (7, vec![2, 1, 1, 1, 4]),
325            (8, vec![2, 1, 4]),
326            (10, vec![3, 6]),
327            (11, vec![3, 3, 6]),
328            (12, vec![3, 2, 6]),
329            (13, vec![3, 1, 1, 1, 1, 6]),
330            (14, vec![3, 1, 2, 1, 6]),
331            (15, vec![3, 1, 6]),
332            (17, vec![4, 8]),
333            (18, vec![4, 4, 8]),
334            (19, vec![4, 2, 1, 3, 1, 2, 8]),
335            (20, vec![4, 2, 8]),
336            (21, vec![4, 1, 1, 2, 1, 1, 8]),
337            (22, vec![4, 1, 2, 4, 2, 1, 8]),
338            (23, vec![4, 1, 3, 1, 8]),
339            (24, vec![4, 1, 8]),
340            (26, vec![5, 10]),
341            (27, vec![5, 5, 10]),
342            (28, vec![5, 3, 2, 3, 10]),
343            (29, vec![5, 2, 1, 1, 2, 10]),
344            (30, vec![5, 2, 10]),
345            (31, vec![5, 1, 1, 3, 5, 3, 1, 1, 10]),
346            (32, vec![5, 1, 1, 1, 10]),
347            (33, vec![5, 1, 2, 1, 10]),
348            (34, vec![5, 1, 4, 1, 10]),
349            (35, vec![5, 1, 10]),
350            (37, vec![6, 12]),
351            (38, vec![6, 6, 12]),
352            (39, vec![6, 4, 12]),
353            (40, vec![6, 3, 12]),
354            (41, vec![6, 2, 2, 12]),
355            (42, vec![6, 2, 12]),
356            (43, vec![6, 1, 1, 3, 1, 5, 1, 3, 1, 1, 12]),
357            (44, vec![6, 1, 1, 1, 2, 1, 1, 1, 12]),
358            (45, vec![6, 1, 2, 2, 2, 1, 12]),
359            (46, vec![6, 1, 3, 1, 1, 2, 6, 2, 1, 1, 3, 1, 12]),
360            (47, vec![6, 1, 5, 1, 12]),
361            (48, vec![6, 1, 12]),
362            (50, vec![7, 14]),
363            (51, vec![7, 7, 14]),
364            (52, vec![7, 4, 1, 2, 1, 4, 14]),
365            (53, vec![7, 3, 1, 1, 3, 14]),
366            (54, vec![7, 2, 1, 6, 1, 2, 14]),
367            (55, vec![7, 2, 2, 2, 14]),
368            (56, vec![7, 2, 14]),
369            (57, vec![7, 1, 1, 4, 1, 1, 14]),
370            (58, vec![7, 1, 1, 1, 1, 1, 1, 14]),
371            (59, vec![7, 1, 2, 7, 2, 1, 14]),
372            (60, vec![7, 1, 2, 1, 14]),
373            (61, vec![7, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14]),
374            (62, vec![7, 1, 6, 1, 14]),
375            (63, vec![7, 1, 14]),
376            (65, vec![8, 16]),
377            (66, vec![8, 8, 16]),
378            (67, vec![8, 5, 2, 1, 1, 7, 1, 1, 2, 5, 16]),
379            (68, vec![8, 4, 16]),
380            (69, vec![8, 3, 3, 1, 4, 1, 3, 3, 16]),
381            (70, vec![8, 2, 1, 2, 1, 2, 16]),
382            (71, vec![8, 2, 2, 1, 7, 1, 2, 2, 16]),
383            (72, vec![8, 2, 16]),
384            (73, vec![8, 1, 1, 5, 5, 1, 1, 16]),
385            (74, vec![8, 1, 1, 1, 1, 16]),
386            (75, vec![8, 1, 1, 1, 16]),
387            (76, vec![8, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16]),
388            (77, vec![8, 1, 3, 2, 3, 1, 16]),
389            (78, vec![8, 1, 4, 1, 16]),
390            (79, vec![8, 1, 7, 1, 16]),
391            (80, vec![8, 1, 16]),
392            (82, vec![9, 18]),
393            (83, vec![9, 9, 18]),
394            (84, vec![9, 6, 18]),
395            (85, vec![9, 4, 1, 1, 4, 18]),
396            (86, vec![9, 3, 1, 1, 1, 8, 1, 1, 1, 3, 18]),
397            (87, vec![9, 3, 18]),
398            (88, vec![9, 2, 1, 1, 1, 2, 18]),
399            (89, vec![9, 2, 3, 3, 2, 18]),
400            (90, vec![9, 2, 18]),
401            (91, vec![9, 1, 1, 5, 1, 5, 1, 1, 18]),
402            (92, vec![9, 1, 1, 2, 4, 2, 1, 1, 18]),
403            (93, vec![9, 1, 1, 1, 4, 6, 4, 1, 1, 1, 18]),
404            (94, vec![9, 1, 2, 3, 1, 1, 5, 1, 8, 1, 5, 1, 1, 3, 2, 1, 18]),
405            (95, vec![9, 1, 2, 1, 18]),
406            (96, vec![9, 1, 3, 1, 18]),
407            (97, vec![9, 1, 5, 1, 1, 1, 1, 1, 1, 5, 1, 18]),
408            (98, vec![9, 1, 8, 1, 18]),
409            (99, vec![9, 1, 18]),
410            (101, vec![10, 20]),
411        ];
412        for (d, w) in continued_fractions {
413            let d = Z::from(d);
414            let v = continued_fraction_of_sqrt(d).unwrap();
415            assert_eq!(v, to_z(&w));
416        }
417    }
418    #[test]
419    fn test_continued_fraction_of_sqrt338() {
420        let v = continued_fraction_of_sqrt(Z::from(338)).unwrap();
421        assert_eq!(v, to_z(&[18, 2, 1, 1, 2, 36]));
422    }
423    #[test]
424    fn test_solve_pell() {
425        let v = solve_pell(&Z::from(653), None);
426        assert_eq!(
427            v,
428            Solution::Negative(Z::from(2_291_286_382u64), Z::from(89_664_965))
429        );
430    }
431    #[test]
432    fn test_solve_pell2() {
433        let v = solve_pell(&Z::from(641), None);
434        assert_eq!(
435            v,
436            Solution::Negative(Z::from(36_120_833_468u64), Z::from(1_426_687_145))
437        );
438    }
439    #[test]
440    fn test_solve_pell3() {
441        let Solution::Negative(x, y) = solve_pell(&Z::from(1021), None) else {
442            panic!("not negative")
443        };
444        assert_eq!(
445            x,
446            Z::from_str_radix("315217280372584882515030", 10).unwrap()
447        );
448        assert_eq!(y, Z::from_str_radix("9865001296666956406909", 10).unwrap());
449    }
450    fn test_pell_aux_aux(d: i32, expected: &Solution) {
451        let solution = solve_pell_aux(&Z::from(d), None);
452        assert_eq!(&solution, expected);
453    }
454    #[test]
455    #[allow(clippy::too_many_lines)]
456    fn test_pell_aux() {
457        use Solution::{Negative, NotExist, Positive};
458        let list = [
459            (1, NotExist),
460            (2, Negative(Z::from(1), Z::from(1))),
461            (3, Positive(Z::from(2), Z::from(1))),
462            (4, NotExist),
463            (5, Negative(Z::from(2), Z::from(1))),
464            (6, Positive(Z::from(5), Z::from(2))),
465            (7, Positive(Z::from(8), Z::from(3))),
466            (8, Positive(Z::from(3), Z::from(1))),
467            (9, NotExist),
468            (10, Negative(Z::from(3), Z::from(1))),
469            (11, Positive(Z::from(10), Z::from(3))),
470            (12, Positive(Z::from(7), Z::from(2))),
471            (13, Negative(Z::from(18), Z::from(5))),
472            (14, Positive(Z::from(15), Z::from(4))),
473            (15, Positive(Z::from(4), Z::from(1))),
474            (16, NotExist),
475            (17, Negative(Z::from(4), Z::from(1))),
476            (18, Positive(Z::from(17), Z::from(4))),
477            (19, Positive(Z::from(170), Z::from(39))),
478            (20, Positive(Z::from(9), Z::from(2))),
479            (21, Positive(Z::from(55), Z::from(12))),
480            (22, Positive(Z::from(197), Z::from(42))),
481            (23, Positive(Z::from(24), Z::from(5))),
482            (24, Positive(Z::from(5), Z::from(1))),
483            (25, NotExist),
484            (26, Negative(Z::from(5), Z::from(1))),
485            (27, Positive(Z::from(26), Z::from(5))),
486            (28, Positive(Z::from(127), Z::from(24))),
487            (29, Negative(Z::from(70), Z::from(13))),
488            (30, Positive(Z::from(11), Z::from(2))),
489            (31, Positive(Z::from(1520), Z::from(273))),
490            (32, Positive(Z::from(17), Z::from(3))),
491            (33, Positive(Z::from(23), Z::from(4))),
492            (34, Positive(Z::from(35), Z::from(6))),
493            (35, Positive(Z::from(6), Z::from(1))),
494            (36, NotExist),
495            (37, Negative(Z::from(6), Z::from(1))),
496            (38, Positive(Z::from(37), Z::from(6))),
497            (39, Positive(Z::from(25), Z::from(4))),
498            (40, Positive(Z::from(19), Z::from(3))),
499            (41, Negative(Z::from(32), Z::from(5))),
500            (42, Positive(Z::from(13), Z::from(2))),
501            (43, Positive(Z::from(3482), Z::from(531))),
502            (44, Positive(Z::from(199), Z::from(30))),
503            (45, Positive(Z::from(161), Z::from(24))),
504            (46, Positive(Z::from(24335), Z::from(3588))),
505            (47, Positive(Z::from(48), Z::from(7))),
506            (48, Positive(Z::from(7), Z::from(1))),
507            (49, NotExist),
508            (50, Negative(Z::from(7), Z::from(1))),
509            (51, Positive(Z::from(50), Z::from(7))),
510            (52, Positive(Z::from(649), Z::from(90))),
511            (53, Negative(Z::from(182), Z::from(25))),
512            (54, Positive(Z::from(485), Z::from(66))),
513            (55, Positive(Z::from(89), Z::from(12))),
514            (56, Positive(Z::from(15), Z::from(2))),
515            (57, Positive(Z::from(151), Z::from(20))),
516            (58, Negative(Z::from(99), Z::from(13))),
517            (59, Positive(Z::from(530), Z::from(69))),
518            (60, Positive(Z::from(31), Z::from(4))),
519            (61, Negative(Z::from(29718), Z::from(3805))),
520            (62, Positive(Z::from(63), Z::from(8))),
521            (63, Positive(Z::from(8), Z::from(1))),
522            (64, NotExist),
523            (65, Negative(Z::from(8), Z::from(1))),
524            (66, Positive(Z::from(65), Z::from(8))),
525            (67, Positive(Z::from(48842), Z::from(5967))),
526            (68, Positive(Z::from(33), Z::from(4))),
527            (69, Positive(Z::from(7775), Z::from(936))),
528            (70, Positive(Z::from(251), Z::from(30))),
529            (71, Positive(Z::from(3480), Z::from(413))),
530            (72, Positive(Z::from(17), Z::from(2))),
531            (73, Negative(Z::from(1068), Z::from(125))),
532            (74, Negative(Z::from(43), Z::from(5))),
533            (75, Positive(Z::from(26), Z::from(3))),
534            (76, Positive(Z::from(57799), Z::from(6630))),
535            (77, Positive(Z::from(351), Z::from(40))),
536            (78, Positive(Z::from(53), Z::from(6))),
537            (79, Positive(Z::from(80), Z::from(9))),
538            (80, Positive(Z::from(9), Z::from(1))),
539            (81, NotExist),
540            (82, Negative(Z::from(9), Z::from(1))),
541            (83, Positive(Z::from(82), Z::from(9))),
542            (84, Positive(Z::from(55), Z::from(6))),
543            (85, Negative(Z::from(378), Z::from(41))),
544            (86, Positive(Z::from(10405), Z::from(1122))),
545            (87, Positive(Z::from(28), Z::from(3))),
546            (88, Positive(Z::from(197), Z::from(21))),
547            (89, Negative(Z::from(500), Z::from(53))),
548            (90, Positive(Z::from(19), Z::from(2))),
549            (91, Positive(Z::from(1574), Z::from(165))),
550            (92, Positive(Z::from(1151), Z::from(120))),
551            (93, Positive(Z::from(12151), Z::from(1260))),
552            (94, Positive(Z::from(2_143_295), Z::from(221_064))),
553            (95, Positive(Z::from(39), Z::from(4))),
554            (96, Positive(Z::from(49), Z::from(5))),
555            (97, Negative(Z::from(5604), Z::from(569))),
556            (98, Positive(Z::from(99), Z::from(10))),
557            (99, Positive(Z::from(10), Z::from(1))),
558            (100, NotExist),
559            (101, Negative(Z::from(10), Z::from(1))),
560            (102, Positive(Z::from(101), Z::from(10))),
561            (103, Positive(Z::from(227_528), Z::from(22419))),
562            (104, Positive(Z::from(51), Z::from(5))),
563            (105, Positive(Z::from(41), Z::from(4))),
564            (106, Negative(Z::from(4005), Z::from(389))),
565            (107, Positive(Z::from(962), Z::from(93))),
566            (108, Positive(Z::from(1351), Z::from(130))),
567            (109, Negative(Z::from(8_890_182), Z::from(851_525))),
568            (110, Positive(Z::from(21), Z::from(2))),
569            (111, Positive(Z::from(295), Z::from(28))),
570            (112, Positive(Z::from(127), Z::from(12))),
571            (113, Negative(Z::from(776), Z::from(73))),
572            (114, Positive(Z::from(1025), Z::from(96))),
573            (115, Positive(Z::from(1126), Z::from(105))),
574        ];
575        for (d, expected) in list {
576            test_pell_aux_aux(d, &expected);
577        }
578    }
579    fn solve_pell_rcf(d: &Z, limit: Option<u64>) -> Solution {
580        if d.is_negative() || d.is_perfect_square() {
581            return Solution::NotExist;
582        }
583        let (negative, b) = {
584            let iter = ContinuedFractionIterator::new_cut_tail(d.clone())
585                .unwrap()
586                .skip(1);
587            let mut b_old = Z::ZERO;
588            let mut b_now = Z::ONE.clone();
589            let mut negative = true;
590            for ai in iter {
591                b_old += &ai * &b_now;
592                check_limit!(b_old, limit);
593                std::mem::swap(&mut b_old, &mut b_now);
594                negative ^= true;
595            }
596            (negative, b_now)
597        };
598        if negative {
599            gen_x_neg!(b, d, limit)
600        } else {
601            gen_x_pos!(b, d, limit)
602        }
603    }
604    #[test]
605    fn test_pell_aux_2() {
606        for d in 2..30000 {
607            let d = Z::from(d);
608            let s2 = solve_pell_aux(&d, None);
609            let s1 = solve_pell_rcf(&d, None);
610            assert_eq!(s1, s2);
611        }
612    }
613}
pell_equation/lib.rs

pell_equation/
lib.rs
