proof_engine/fractal/
newton.rs

1//! Newton fractal for arbitrary polynomials.
2
3/// A polynomial represented as coefficients [a₀, a₁, ..., aₙ] for a₀ + a₁x + ... + aₙxⁿ.
4#[derive(Debug, Clone)]
5pub struct Polynomial { pub coeffs: Vec<(f64, f64)> } // complex coefficients (re, im)
6
7impl Polynomial {
8    pub fn real(coeffs: &[f64]) -> Self {
9        Self { coeffs: coeffs.iter().map(|&c| (c, 0.0)).collect() }
10    }
11
12    pub fn eval(&self, z: (f64, f64)) -> (f64, f64) {
13        let mut result = (0.0, 0.0);
14        let mut zn = (1.0, 0.0);
15        for &c in &self.coeffs {
16            result = (result.0 + c.0 * zn.0 - c.1 * zn.1, result.1 + c.0 * zn.1 + c.1 * zn.0);
17            zn = (zn.0 * z.0 - zn.1 * z.1, zn.0 * z.1 + zn.1 * z.0);
18        }
19        result
20    }
21
22    pub fn derivative(&self) -> Polynomial {
23        if self.coeffs.len() <= 1 { return Polynomial { coeffs: vec![(0.0, 0.0)] }; }
24        let coeffs: Vec<(f64, f64)> = self.coeffs.iter().enumerate().skip(1)
25            .map(|(i, &(re, im))| (re * i as f64, im * i as f64)).collect();
26        Polynomial { coeffs }
27    }
28}
29
30/// Newton fractal renderer.
31pub struct NewtonFractal {
32    pub poly: Polynomial,
33    pub max_iter: u32,
34    pub tolerance: f64,
35}
36
37/// Result of one Newton pixel.
38#[derive(Debug, Clone, Copy)]
39pub struct NewtonPixel {
40    pub root_index: i32, // -1 if not converged
41    pub iterations: u32,
42    pub final_z: (f64, f64),
43}
44
45impl NewtonFractal {
46    pub fn new(poly: Polynomial, max_iter: u32) -> Self {
47        Self { poly, max_iter, tolerance: 1e-6 }
48    }
49
50    /// z³ - 1 = 0 (classic Newton fractal with 3 roots).
51    pub fn cubic_roots() -> Self {
52        Self::new(Polynomial::real(&[-1.0, 0.0, 0.0, 1.0]), 64)
53    }
54
55    pub fn compute_pixel(&self, z_re: f64, z_im: f64, roots: &[(f64, f64)]) -> NewtonPixel {
56        let dp = self.poly.derivative();
57        let mut z = (z_re, z_im);
58
59        for iter in 0..self.max_iter {
60            let fz = self.poly.eval(z);
61            let dfz = dp.eval(z);
62            let denom = dfz.0 * dfz.0 + dfz.1 * dfz.1;
63            if denom < 1e-20 { break; }
64            z.0 -= (fz.0 * dfz.0 + fz.1 * dfz.1) / denom;
65            z.1 -= (fz.1 * dfz.0 - fz.0 * dfz.1) / denom;
66
67            // Check convergence to known roots
68            for (ri, &root) in roots.iter().enumerate() {
69                let dr = z.0 - root.0;
70                let di = z.1 - root.1;
71                if dr * dr + di * di < self.tolerance * self.tolerance {
72                    return NewtonPixel { root_index: ri as i32, iterations: iter, final_z: z };
73                }
74            }
75        }
76        NewtonPixel { root_index: -1, iterations: self.max_iter, final_z: z }
77    }
78
79    /// Render the full Newton fractal.
80    pub fn render(&self, width: u32, height: u32, center: (f64, f64), zoom: f64, roots: &[(f64, f64)]) -> Vec<NewtonPixel> {
81        let scale = 2.0 / zoom;
82        let aspect = width as f64 / height as f64;
83        let mut pixels = Vec::with_capacity((width * height) as usize);
84        for py in 0..height { for px in 0..width {
85            let z_re = center.0 + (px as f64 / width as f64 - 0.5) * scale * aspect;
86            let z_im = center.1 + (py as f64 / height as f64 - 0.5) * scale;
87            pixels.push(self.compute_pixel(z_re, z_im, roots));
88        }}
89        pixels
90    }
91}
92
93#[cfg(test)]
94mod tests {
95    use super::*;
96    #[test]
97    fn cubic_roots_converges() {
98        let nf = NewtonFractal::cubic_roots();
99        let roots = vec![(1.0, 0.0), (-0.5, 0.866), (-0.5, -0.866)];
100        // Point near root 0
101        let p = nf.compute_pixel(0.9, 0.0, &roots);
102        assert_eq!(p.root_index, 0);
103    }
104    #[test]
105    fn polynomial_eval() {
106        let p = Polynomial::real(&[1.0, 0.0, 1.0]); // 1 + x²
107        let (re, im) = p.eval((2.0, 0.0));
108        assert!((re - 5.0).abs() < 1e-10);
109    }
110    #[test]
111    fn derivative_correct() {
112        let p = Polynomial::real(&[0.0, 0.0, 0.0, 1.0]); // x³
113        let dp = p.derivative(); // 3x²
114        let (re, _) = dp.eval((2.0, 0.0));
115        assert!((re - 12.0).abs() < 1e-10);
116    }
117}
proof_engine/fractal/newton.rs

proof_engine/fractal/
newton.rs
