Merge pull request #13 from jackyarndley/development

add in minor optimizations

Author: jackyarndley

output.png

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

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

@@ -23,7 +23,7 @@ impl Reference {
             maximum_iteration,
             z,
             c,
-            data: Vec::new()
+            data: Vec::with_capacity(1000)
         }
     }
 

src/math/reference.rs

@@ -8,7 +8,7 @@ use std::mem::swap;
 pub struct SeriesApproximation {
     pub current_iteration: usize,
     maximum_iteration: usize,
-    delta_pixel: FloatExtended,
+    delta_pixel_square: FloatExtended,
     z: ComplexArbitrary,
     c: ComplexArbitrary,
     pub order: usize,
@@ -22,7 +22,7 @@ pub struct SeriesApproximation {
 }
 
 impl SeriesApproximation {
-    pub fn new(c: ComplexArbitrary, order: usize, maximum_iteration: usize, delta_pixel: FloatExtended, delta_top_left: ComplexExtended) -> Self {
+    pub fn new(c: ComplexArbitrary, order: usize, maximum_iteration: usize, delta_pixel_square: FloatExtended, delta_top_left: ComplexExtended) -> Self {
         let mut coefficients = vec![ComplexExtended::new2(0.0, 0.0, 0); order as usize + 1];
 
         coefficients[0] = to_extended(&c);
@@ -32,7 +32,7 @@ impl SeriesApproximation {
         SeriesApproximation {
             current_iteration: 1,
             maximum_iteration,
-            delta_pixel,
+            delta_pixel_square,
             z: c.clone(),
             c,
             order,
@@ -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,33 +80,31 @@ 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();
-            let mut derivative = derivative_probe.norm();
+            let relative_error = (self.current_probes[i] - series_probe).norm_square();
+            let mut derivative = derivative_probe.norm_square();
 
             // Check to make sure that the derivative is greater than or equal to 1
             if derivative.to_float() < 1.0 {
                 derivative = FloatExtended::new(1.0, 0);
             }
 
             // Check that the error over the derivative is less than the pixel spacing
-            if relative_error / derivative > self.delta_pixel {
+            if relative_error / derivative > self.delta_pixel_square {
                 self.z -= &self.c;
                 self.z.sqrt_mut();
                 return false;
@@ -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,11 +69,15 @@ impl FractalRenderer {
 
         println!("Rendering...");
 
+        let delta_pixel_extended = FloatExtended::new(delta_pixel, -self.zoom.exponent);
+
+        // 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,
             self.maximum_iteration,
-            FloatExtended::new(delta_pixel, -self.zoom.exponent),
+            delta_pixel_extended * delta_pixel_extended,
             ComplexExtended::new(delta_top_left, -self.zoom.exponent),
         );
 
@@ -120,9 +124,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);
+        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 +163,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);
+            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

@@ -1,15 +1,17 @@
 use crate::util::image::Image;
 use crate::util::PixelData;
+use crate::math::Reference;
 
 pub enum Colouring {
     Iteration,
     IterationSquareRoot,
-    Histogram,
+    IterationSmooth,
+    IterationHistogram,
     // Distance
 }
 
 impl Colouring {
-    pub fn run(&self, pixel_data: &Vec<PixelData>, image: &mut Image, maximum_iteration: usize, _delta_pixel: f64) {
+    pub fn run(&self, pixel_data: &Vec<PixelData>, image: &mut Image, maximum_iteration: usize, _delta_pixel: f64, reference: &Reference) {
         // Palette is temporarily here
         let mut colours = Vec::new();
 
@@ -55,18 +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 {
                         // 0.1656
-                        let hue = (250.0f64 * pixel.iteration as f64) as usize % 8192;
+                        let hue = (250 * pixel.iteration) % 8192;
 
                         let colour = colours[hue];
 
@@ -94,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();
     }
 }