Skip to main content

nekolib/math/
sqrt_fraction_.rs

1//! 平方根の連分数展開。
2
3use super::sqrt;
4
5use sqrt::Sqrt;
6
7/// 平方根の連分数展開。
8///
9/// $\\sqrt{n}$ の連分数展開を $[a\_0; a\_1, a\_2, \\dots]$ とする。
10/// このとき、$a\_\\bullet$ を生成するイテレータを返す。
11///
12/// # Examples
13/// ```
14/// use nekolib::math::sqrt_fraction;
15///
16/// let a: Vec<_> = sqrt_fraction(3).take(30).map(|x| x as f64).collect();
17/// let r = a.into_iter().rev().reduce(|a1, a0| a0 + 1.0 / a1).unwrap();
18///
19/// assert!((r - 3.0_f64.sqrt()).abs() < 1.0e-16);
20/// ```
21pub fn sqrt_fraction(n: i128) -> impl Iterator<Item = i128> {
22    let (r, next) = sqrt_fraction_fn(n);
23    std::iter::successors(Some((r, r, 1)), move |&(_, b, c)| Some(next(b, c)))
24        .map(|(a, ..)| a)
25}
26
27/// 平方根の連分数展開。
28///
29/// $\\sqrt{n}$ の連分数展開の $i$ step 目が次のように表されるとする。
30/// $$ \\sqrt{n} = a\_0 +
31/// \\frac{1}{\\dots\\,+\\frac{1}{a\_{i-1}+\\frac{\\sqrt{n}-b\_{i-1}}{c\_{i-1}}}}.
32/// $$
33///
34/// $a\_0$ と、関数 $f: (b\_i, c\_i)\\mapsto (a\_{i+1}, b\_{i+1}, c\_{i+1})$
35/// を返す。
36/// ただし、$(a\_0, b\_0, c\_0)
37/// = (\\lfloor\\sqrt{n}\\rfloor, \\lfloor\\sqrt{n}\\rfloor, 1)$ である。
38///
39/// 実際の連分数展開が欲しいときは [`sqrt_fraction`] を用いればよい。
40/// 周期検出などをしたいときは $(b\_\\bullet, c\_\\bullet)$ が必要になる。
41///
42/// # Derivation
43/// `todo!()`
44///
45/// # Examples
46/// ```
47/// use nekolib::math::sqrt_fraction_fn;
48///
49/// let (a0, next) = sqrt_fraction_fn(3);
50/// assert_eq!(a0, 1);
51/// let (a1, b1, c1) = next(a0, 1);
52/// let (a2, b2, c2) = next(b1, c1);
53/// let (a3, b3, c3) = next(b2, c2);
54///
55/// assert_eq!((a1, b1, c1), (a3, b3, c3));  // sqrt(3) has period 2
56/// assert_eq!([a0, a1, a2], [1, 1, 2]);  // sqrt(3) = [1; (1, 2)]
57/// ```
58pub fn sqrt_fraction_fn(
59    n: i128,
60) -> (i128, impl Fn(i128, i128) -> (i128, i128, i128)) {
61    let r = n.sqrt();
62    let next = move |b, c| {
63        assert_eq!((n - b * b) % c, 0);
64        let c_ = (n - b * b) / c;
65        let a_ = (r + b) / c_;
66        let b_ = a_ * c_ - b;
67        (a_, b_, c_)
68    };
69    (r, next)
70}