Problem 66 "Diophantine equation"

Consider quadratic Diophantine equations of the form:

x2 – Dy2 = 1

For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:

32 – 2×22 = 1
22 – 3×12 = 1
92 – 5×42 = 1
52 – 6×22 = 1
82 – 7×32 = 1

Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.

問 66 「ディオファントス方程式」

次の形式の, 2次のディオファントス方程式を考えよう:

x2 – Dy2 = 1

たとえば D=13 のとき, xを最小にする解は x is 6492 – 13×1802 = 1.である.

D が平方数(square)のとき, 正整数のなかに解は存在しないと考えられる.

D = {2, 3, 5, 6, 7} に対して x を最小にする解は次のようになる:

32 – 2×22 = 1
22 – 3×12 = 1
92 – 5×42 = 1
52 – 6×22 = 1
82 – 7×32 = 1

したがって, D ≤ 7 に対して x を最小にする解を考えると, D=5 のとき x は最大である.

D ≤ 1000 に対する x を最小にする解で, x が最大になるような D の値を見つけよ.


Non primitive Pythagorean triples and parabolas

Continued fractions and their periods

x^2 - Dy^2 = 1 ⬄ √D ≈ x/y 

mod fraction {
    trait MixedFraction {
        fn integer_part(&self) -> u32;
    }

    pub struct QuadraticIrrational {
        rational_part: u32,
        irrational_part: u32,
        downscaling_ratio: u32,
        pub integer_part: u32,
        pub iteration: u32,
        integers: Vec<u32>,
    }

    impl MixedFraction for QuadraticIrrational {
        fn integer_part(&self) -> u32 {
            let mixed_fraction 
                = self.rational_part as f64 
                + (self.irrational_part as f64).sqrt();
            (mixed_fraction / self.downscaling_ratio as f64) as u32
        }
    }

    impl QuadraticIrrational {
        pub fn new(square_free_integer: u32) -> Self {
            assert!(
                ((square_free_integer as f64).sqrt() as u32).pow(2)
                != square_free_integer
            );
            let mut q = Self {
                rational_part: 0,
                irrational_part: square_free_integer,
                downscaling_ratio: 1,
                integer_part: u32::MAX,
                iteration: 0,
                integers: vec![],
            };
            q.integer_part = q.integer_part();
            q.integers.push(q.integer_part);
            q
        }
        pub fn next_iteration(&mut self) {
            self.iteration += 1;
            self.rational_part
                = self.downscaling_ratio 
                * self.integer_part 
                - self.rational_part;
            self.downscaling_ratio
                = (self.irrational_part - self.rational_part.pow(2))
                / self.downscaling_ratio;
            self.integer_part = self.integer_part();
            self.integers.push(self.integer_part);
        }
        pub fn last_numerator_approx(&self) -> f64 {
            let mut result = self.integers[0] as f64;
            let mut term_before_last = 1f64;
            for a in (&self.integers[1..self.integers.len() - 1])
                .iter()
                .map(|&a| a as f64)
            {
                let t = result.clone();
                result *= a;
                result += term_before_last;
                term_before_last = t;
            }
            result
        }
        pub fn last_denominator_approx(&self) -> f64 {
            let mut result = 1f64;
            let mut term_before_last = 0f64;
            for a in (&self.integers[1..self.integers.len() - 1])
                .iter()
                .map(|&a| a as f64)
            {
                let t = result.clone();
                result *= a;
                result += term_before_last;
                term_before_last = t;
            }
            result
        }
    }
}

fn main() {
    let mut max_numerator_aprox = 0f64;
    let mut max_d = 0u32;
    for n in 2..=1000 {
        if ((n as f64).sqrt() as u32).pow(2) == n {
            continue;
        }
        let mut surd = fraction::QuadraticIrrational::new(n);
        let integer_part_orig = surd.integer_part;
        while {
            surd.next_iteration();
            surd.integer_part != integer_part_orig * 2 || surd.iteration % 2 != 0
        } {}
        let last_numerator_approx = surd.last_numerator_approx();
        if last_numerator_approx > max_numerator_aprox {
            max_numerator_aprox = last_numerator_approx;
            max_d = n;
        }
    }

    println!("{}", max_d);
    assert_eq!(max_d, 661);
}