optimized some parts of the series approximation

Author: jackyarndley

output.png

@@ -53,13 +53,13 @@ fn main() {
     println!("Zoom: {}", zoom);
 
     let mut renderer = FractalRenderer::new(
-        2000,
-        2000,
+        1000,
+        1000,
         zoom,
         iterations.parse::<usize>().unwrap(),
         center_re,
         center_im,
-        0.001,
+        0.01,
         false,
         32
     );

src/bin/main.rs

@@ -22,13 +22,13 @@ impl Perturbation {
 
                             if z_norm < reference.data[pixel.iteration - reference.start_iteration].z_tolerance {
                                 pixel.glitched = true;
-                                pixel.delta_current.reduce();
                                 break;
                             }
 
                             if z_norm > 1e16 {
                                 pixel.escaped = true;
-                                pixel.delta_current.reduce();
+                                pixel.delta_current.mantissa = pixel.delta_current.to_float();
+                                pixel.delta_current.exponent = 0;
                                 break;
                             }
                         }

src/math/perturbation.rs

@@ -61,9 +61,9 @@ impl SeriesApproximation {
 
         // Calculate the new coefficents
         for k in 2..=self.order {
-            let mut sum = ComplexExtended::new2(0.0, 0.0, 0);
+            let mut sum = self.coefficients[0] * self.coefficients[k];
 
-            for j in 0..=((k - 1) / 2) {
+            for j in 1..=((k - 1) / 2) {
                 sum += self.coefficients[j] * self.coefficients[k - j];
             }
             sum *= 2.0;
@@ -80,21 +80,19 @@ impl SeriesApproximation {
         // this section takes about half of the time
         for i in 0..self.original_probes.len() {
             // step the probe points using perturbation
-            self.current_probes[i] = self.coefficients[0] * self.current_probes[i] * 2.0 + self.current_probes[i] * self.current_probes[i] + self.original_probes[i];
+            self.current_probes[i] = self.current_probes[i] * (self.coefficients[0] * 2.0 + self.current_probes[i]);
+            self.current_probes[i] += self.original_probes[i];
+
+            // TODO this probably does not need to be done every iteration
             self.current_probes[i].reduce();
 
             // get the new approximations
-            let mut series_probe = ComplexExtended::new2(0.0, 0.0, 0);
-            let mut derivative_probe = ComplexExtended::new2(0.0, 0.0, 0);
+            let mut series_probe = self.next_coefficients[1] * self.approximation_probes[i][0];
+            let mut derivative_probe = self.next_coefficients[1] * self.approximation_probes_derivative[i][0];
 
-            for k in 1..=self.order {
+            for k in 2..=self.order {
                 series_probe += self.next_coefficients[k] * self.approximation_probes[i][k - 1];
-                derivative_probe += self.next_coefficients[k] * self.approximation_probes_derivative[i][k - 1] * k as f64;
-
-                if k % 16 == 0{
-                    series_probe.reduce();
-                    derivative_probe.reduce();
-                }
+                derivative_probe += self.next_coefficients[k] * self.approximation_probes_derivative[i][k - 1];
             };
 
             let relative_error = (self.current_probes[i] - series_probe).norm();
@@ -138,24 +136,23 @@ impl SeriesApproximation {
         self.original_probes.push(delta_probe);
         self.current_probes.push(delta_probe);
 
-        let mut delta_probe_n = Vec::with_capacity(self.order + 1);
-        let mut delta_probe_n_derivative = Vec::with_capacity(self.order + 1);
+        let mut delta_probe_n = delta_probe;
+
+        let mut delta_probe_n_2 = Vec::with_capacity(self.order + 1);
+        let mut delta_probe_n_derivative_2 = Vec::with_capacity(self.order + 1);
 
         // The first element will be 1, in order for the derivative to be calculated
-        delta_probe_n.push(delta_probe);
-        delta_probe_n_derivative.push(ComplexExtended::new2(1.0, 0.0, 0));
+        delta_probe_n_2.push(delta_probe_n);
+        delta_probe_n_derivative_2.push(ComplexExtended::new2(1.0, 0.0, 0));
 
         for i in 1..=self.order {
-            let mut delta_probe1 = delta_probe_n[i - 1] * delta_probe;
-            let mut delta_probe2 = delta_probe_n_derivative[i - 1] * delta_probe;
-            delta_probe1.reduce();
-            delta_probe2.reduce();
-            delta_probe_n.push(delta_probe1);
-            delta_probe_n_derivative.push(delta_probe2);
+            delta_probe_n_derivative_2.push((delta_probe_n * (i + 1) as f64));
+            delta_probe_n *= delta_probe;
+            delta_probe_n_2.push(delta_probe_n);
         }
 
-        self.approximation_probes.push(delta_probe_n);
-        self.approximation_probes_derivative.push(delta_probe_n_derivative);
+        self.approximation_probes.push(delta_probe_n_2);
+        self.approximation_probes_derivative.push(delta_probe_n_derivative_2);
     }
 
     // Get the current reference, and the current number of iterations done
@@ -181,31 +178,43 @@ impl SeriesApproximation {
     }
 
     pub fn evaluate(&self, point_delta: ComplexExtended) -> ComplexExtended {
-        let mut approximation = self.coefficients[1] * point_delta;
-        let mut original_point_n = point_delta * point_delta;
+        // 1907 ms packing opus 4K
+        // Horner's rule
+        let mut approximation = self.coefficients[self.order];
 
-        for k in 2..=self.order {
-            approximation += self.coefficients[k] * original_point_n;
-            original_point_n *= point_delta;
-
-            if k % 16 == 0 {
-                approximation.reduce();
-                original_point_n.reduce();
-            }
-        };
+        for k in (1..=(self.order - 1)).rev() {
+            approximation *= point_delta;
+            approximation += self.coefficients[k];
+        }
 
+        approximation *= point_delta;
+        approximation.reduce();
         approximation
     }
 
-    // pub fn evaluate_derivative(&self, point_delta: ComplexFixed<f64>) -> ComplexFixed<f64> {
+    // pub fn evaluate_derivative(&self, point_delta: ComplexExtended) -> FloatExtended {
     //     let mut original_point_derivative_n = ComplexExtended::new(1.0, 0, 0.0, 0);
     //     let mut approximation_derivative = ComplexExtended::new(0.0, 0, 0.0, 0);
-    //
+    
     //     for k in 1..=self.order {
     //         approximation_derivative += k as f64 * self.coefficients[k] * original_point_derivative_n;
     //         original_point_derivative_n *= ComplexExtended::new(point_delta.re, 0, point_delta.im, 0);
     //     };
-    //
+    
     //     approximation_derivative.to_float()
+
+
+    //     let mut approximation_derivative = self.coefficients[self.order] * self.order as f64;
+
+    //     for k in (1..=(self.order - 1)).rev() {
+    //         approximation *= point_delta;
+    //         approximation += self.coefficients[k] * k as f64;
+    //     }
+
+    //     approximation *= point_delta;
+    //     approximation.reduce();
+    //     approximation
+
+
     // }
 }

src/math/series_approximation.rs

@@ -69,6 +69,8 @@ impl FractalRenderer {
 
         println!("Rendering...");
 
+        // Series approximation currently has some overskipping issues
+        // this can be resolved by root finding and adding new probe points
         let mut series_approximation = SeriesApproximation::new(
             self.center_location.clone(),
             self.approximation_order,
@@ -120,9 +122,11 @@ impl FractalRenderer {
         println!("{:<14}{:>6} ms", "Iteration", time.elapsed().as_millis());
 
         let time = Instant::now();
-        Colouring::Iteration.run(&pixel_data, &mut self.image, self.maximum_iteration, delta_pixel, &reference);
+        Colouring::IterationSmooth.run(&pixel_data, &mut self.image, self.maximum_iteration, delta_pixel, &reference);
         println!("{:<14}{:>6} ms", "Coloring", time.elapsed().as_millis());
 
+        let time = Instant::now();
+
         // Remove all non-glitched points from the remaining points
         pixel_data.retain(|packet| {
             packet.glitched
@@ -157,7 +161,7 @@ impl FractalRenderer {
 
             Perturbation::iterate(&mut pixel_data, &r, r.current_iteration);
 
-            Colouring::Iteration.run(&pixel_data, &mut self.image, self.maximum_iteration, delta_pixel, &r);
+            Colouring::IterationSmooth.run(&pixel_data, &mut self.image, self.maximum_iteration, delta_pixel, &r);
 
             // Remove all non-glitched points from the remaining points
             pixel_data.retain(|packet| {

src/renderer.rs

@@ -5,7 +5,8 @@ use crate::math::Reference;
 pub enum Colouring {
     Iteration,
     IterationSquareRoot,
-    Histogram,
+    IterationSmooth,
+    IterationHistogram,
     // Distance
 }
 
@@ -56,21 +57,18 @@ impl Colouring {
             colours.push((red, green, blue))
         }
 
+        let escape_radius_ln = 1e16f64.ln();
+
         match self {
             Colouring::Iteration => {
-                // No smooth colouring at the moment
-
                 for pixel in pixel_data {
                     let (r, g, b) = if pixel.glitched && image.display_glitches {
                         (255, 0, 0)
                     } else if pixel.iteration >= maximum_iteration {
                         (0, 0, 0)
                     } else {
-                        let z_norm = (reference.data[pixel.iteration - reference.start_iteration].z_fixed + pixel.delta_current.to_float()).norm();
-                        let smooth = 1.0 - (z_norm.ln()).log2();
-
                         // 0.1656
-                        let hue = (7.0f64 * (pixel.iteration as f64 + smooth) + 8192.0) as usize % 8192;
+                        let hue = (250 * pixel.iteration) % 8192;
 
                         let colour = colours[hue];
 
@@ -98,7 +96,28 @@ impl Colouring {
                     image.plot(pixel.image_x, pixel.image_y, r, g, b);
                 }
             },
-            Colouring::Histogram => {
+            Colouring::IterationSmooth => {
+                for pixel in pixel_data {
+                    let (r, g, b) = if pixel.glitched && image.display_glitches {
+                        (255, 0, 0)
+                    } else if pixel.iteration >= maximum_iteration {
+                        (0, 0, 0)
+                    } else {
+                        let z_norm = (reference.data[pixel.iteration - reference.start_iteration].z_fixed + pixel.delta_current.mantissa).norm_sqr();
+                        let smooth = 1.0 - (z_norm.ln() / escape_radius_ln).log2();
+
+                        // 0.1656
+                        let hue = (250.0f64 * (pixel.iteration as f64 + smooth)) as usize % 8192;
+
+                        let colour = colours[hue];
+
+                        (colour.0 as u8, colour.1 as u8, colour.2 as u8)
+                    };
+
+                    image.plot(pixel.image_x, pixel.image_y, r, g, b);
+                }
+            },
+            Colouring::IterationHistogram => {
                 let mut iteration_counts = vec![0usize; maximum_iteration + 2];
 
                 for pixel in pixel_data {

src/util/colouring.rs

@@ -87,7 +87,6 @@ impl AddAssign for ComplexExtended {
                 self.mantissa += other.mantissa;
             }
         }
-        // self.reduce();
     }
 }
 
@@ -103,7 +102,6 @@ impl SubAssign for ComplexExtended {
             self.exponent = other.exponent;
             self.mantissa -= other.mantissa;
         }
-        self.reduce();
     }
 }
 
@@ -112,7 +110,6 @@ impl MulAssign<ComplexExtended> for ComplexExtended {
     fn mul_assign(&mut self, other: Self) {
         self.mantissa *= other.mantissa;
         self.exponent += other.exponent;
-        // self.reduce();
     }
 }