Merge pull request #14 from jackyarndley/development

added in optimized series approximation

Author: jackyarndley

.gitignore

@@ -1,4 +1,7 @@
 /target
 **/*.rs.bk
 .idea
-*.kfr
+*.kfr
+*.svg
+*.data
+*.old

Cargo.lock

@@ -13,6 +13,7 @@ rug = "1.9.0"
 itertools = "^0.9.0"
 
 #Additional commands that can improve performance (maybe by around 5-10%)
-#[profile.release]
+[profile.release]
 #lto = "fat"
 #codegen-units = 1
+#debug = true

Cargo.toml

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

output.png

@@ -1,4 +1,4 @@
-use crate::util::PixelData;
+use crate::util::{PixelData, FloatExp};
 
 use rayon::prelude::*;
 use crate::math::reference::Reference;
@@ -7,12 +7,12 @@ pub struct Perturbation {}
 
 impl Perturbation {
     pub fn iterate(pixel_data: &mut Vec<PixelData>, reference: &Reference, reference_current_iteration: usize) {
-        pixel_data.par_chunks_mut(4)
+        pixel_data.par_chunks_mut(64)
             .for_each(|pixel_data| {
                 for pixel in pixel_data {
                     let mut scaled_iterations = 0;
-                    let mut scaled_scale_factor_1 = 2.0f64.powi(pixel.delta_current.exponent);
-                    let mut scaled_delta_reference = 2.0f64.powi(pixel.delta_reference.exponent - pixel.delta_current.exponent) * pixel.delta_reference.mantissa;
+                    let mut scaled_scale_factor_1 = 1.0f64.ldexp(pixel.delta_current.exponent);
+                    let mut scaled_delta_reference = 1.0f64.ldexp(pixel.delta_reference.exponent - pixel.delta_current.exponent) * pixel.delta_reference.mantissa;
 
                     while pixel.iteration < reference_current_iteration {
                         let delta_current_float = scaled_scale_factor_1 * pixel.delta_current.mantissa;
@@ -43,8 +43,8 @@ impl Perturbation {
                                 // reset the scaled counter
                                 pixel.delta_current.reduce();
 
-                                scaled_scale_factor_1 = 2.0f64.powi(pixel.delta_current.exponent);
-                                scaled_delta_reference = 2.0f64.powi(pixel.delta_reference.exponent - pixel.delta_current.exponent) * pixel.delta_reference.mantissa;
+                                scaled_scale_factor_1 = 1.0f64.ldexp(pixel.delta_current.exponent);
+                                scaled_delta_reference = 1.0f64.ldexp(pixel.delta_reference.exponent - pixel.delta_current.exponent) * pixel.delta_reference.mantissa;
 
                                 scaled_iterations = 0;
                             },
@@ -61,8 +61,8 @@ impl Perturbation {
                                 if scaled_iterations > 250 {
                                     pixel.delta_current.reduce();
 
-                                    scaled_scale_factor_1 = 2.0f64.powi(pixel.delta_current.exponent);
-                                    scaled_delta_reference = 2.0f64.powi(pixel.delta_reference.exponent - pixel.delta_current.exponent) * pixel.delta_reference.mantissa;
+                                    scaled_scale_factor_1 = 1.0f64.ldexp(pixel.delta_current.exponent);
+                                    scaled_delta_reference = 1.0f64.ldexp(pixel.delta_reference.exponent - pixel.delta_current.exponent) * pixel.delta_reference.mantissa;
 
                                     scaled_iterations = 0;
                                 }

output2.png

@@ -46,88 +46,88 @@ impl SeriesApproximation {
         }
     }
 
-    pub fn step(&mut self) -> bool {
-        self.z.square_mut();
-        self.z += &self.c;
-
-        let z_fixed = to_fixed(&self.z);
-
-        if z_fixed.re.abs() < 1e-300 && z_fixed.im.abs() < 1e-300 {
-            println!("found slow at: {}", self.current_iteration);
-        }
+    pub fn run(&mut self) {
+        self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
+        self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
+        self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
+        self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
 
-        self.next_coefficients[0] = to_extended(&self.z);
-        self.next_coefficients[1] = self.coefficients[0] * self.coefficients[1] * 2.0 + ComplexExtended::new2(1.0, 0.0, 0);
+        let add_value = ComplexExtended::new2(1.0, 0.0, 0);
 
-        // Calculate the new coefficents
-        for k in 2..=self.order {
-            let mut sum = self.coefficients[0] * self.coefficients[k];
+        // Can be changed later into a better loop - this function could also return some more information
+        while self.current_iteration < self.maximum_iteration {
+            self.z.square_mut();
+            self.z += &self.c;
 
-            for j in 1..=((k - 1) / 2) {
-                sum += self.coefficients[j] * self.coefficients[k - j];
-            }
-            sum *= 2.0;
+            let z_fixed = to_fixed(&self.z);
 
-            // If even, we include the mid term as well
-            if k % 2 == 0 {
-                sum += self.coefficients[k / 2] * self.coefficients[k / 2];
+            if z_fixed.re.abs() < 1e-300 && z_fixed.im.abs() < 1e-300 {
+                println!("found slow at: {}", self.current_iteration);
             }
 
-            sum.reduce();
-            self.next_coefficients[k] = sum;
-        }
-
-        // 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.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 = self.next_coefficients[1] * self.approximation_probes[i][0];
-            let mut derivative_probe = self.next_coefficients[1] * self.approximation_probes_derivative[i][0];
+            self.next_coefficients[0] = to_extended(&self.z);
+            self.next_coefficients[1] = self.coefficients[0] * self.coefficients[1] * 2.0 + add_value;
+            self.next_coefficients[0].reduce();
+            self.next_coefficients[1].reduce();
 
+            // Calculate the new coefficents
             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];
-            };
+                let mut sum = self.coefficients[0] * self.coefficients[k];
 
-            let relative_error = (self.current_probes[i] - series_probe).norm_square();
-            let mut derivative = derivative_probe.norm_square();
+                for j in 1..=((k - 1) / 2) {
+                    sum += self.coefficients[j] * self.coefficients[k - j];
+                }
+                sum *= 2.0;
 
-            // 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);
-            }
+                // If even, we include the mid term as well
+                if k % 2 == 0 {
+                    sum += self.coefficients[k / 2] * self.coefficients[k / 2];
+                }
 
-            // Check that the error over the derivative is less than the pixel spacing
-            if relative_error / derivative > self.delta_pixel_square {
-                self.z -= &self.c;
-                self.z.sqrt_mut();
-                return false;
+                sum.reduce();
+                self.next_coefficients[k] = sum;
             }
-        }
 
-        // If the approximation is valid, continue
-        self.current_iteration += 1;
-
-        // Swap the two coefficients buffers
-        swap(&mut self.coefficients, &mut self.next_coefficients);
-        self.current_iteration < self.maximum_iteration
-    }
+            // 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.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 = 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 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];
+                };
+
+                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.mantissa = 1.0;
+                    derivative.exponent = 0;
+                }
+
+                // Check that the error over the derivative is less than the pixel spacing
+                if relative_error / derivative > self.delta_pixel_square {
+                    self.z -= &self.c;
+                    self.z.sqrt_mut();
+                    return;
+                }
+            }
 
-    pub fn run(&mut self) {
-        self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
-        self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
-        self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
-        self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
+            // If the approximation is valid, continue
+            self.current_iteration += 1;
 
-        // Can be changed later into a better loop - this function could also return some more information
-        while self.step() {
-            continue;
+            // Swap the two coefficients buffers
+            swap(&mut self.coefficients, &mut self.next_coefficients);
         }
     }
 
@@ -136,23 +136,23 @@ impl SeriesApproximation {
         self.original_probes.push(delta_probe);
         self.current_probes.push(delta_probe);
 
-        let mut delta_probe_n = delta_probe;
+        let mut current_value = 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);
+        let mut delta_n = Vec::with_capacity(self.order + 1);
+        let mut delta_derivative_n = Vec::with_capacity(self.order + 1);
 
         // The first element will be 1, in order for the derivative to be calculated
-        delta_probe_n_2.push(delta_probe_n);
-        delta_probe_n_derivative_2.push(ComplexExtended::new2(1.0, 0.0, 0));
+        delta_n.push(current_value);
+        delta_derivative_n.push(ComplexExtended::new2(1.0, 0.0, 0));
 
         for i in 1..=self.order {
-            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);
+            delta_derivative_n.push(current_value * (i + 1) as f64);
+            current_value *= delta_probe;
+            delta_n.push(current_value);
         }
 
-        self.approximation_probes.push(delta_probe_n_2);
-        self.approximation_probes_derivative.push(delta_probe_n_derivative_2);
+        self.approximation_probes.push(delta_n);
+        self.approximation_probes_derivative.push(delta_derivative_n);
     }
 
     // Get the current reference, and the current number of iterations done

src/bin/main.rs

@@ -2,11 +2,10 @@ use crate::util::{image::Image, ComplexFixed, ComplexArbitrary, PixelData, compl
 use crate::math::{SeriesApproximation, Perturbation};
 
 use std::time::Instant;
-use std::cmp::max;
+use std::cmp::{min, max};
 use std::f64::consts::LOG2_10;
 
 use rand::seq::SliceRandom;
-use itertools::Itertools;
 use rayon::prelude::*;
 
 pub struct FractalRenderer {
@@ -46,14 +45,21 @@ impl FractalRenderer {
             precision as u32,
             ComplexArbitrary::parse("(".to_owned() + center_real + "," + center_imag + ")").expect("Location is not valid!"));
 
+        let auto_approximation = if approximation_order == 0 {
+            let auto = (((image_width * image_height) as f64).log(1e6).powf(6.619) * 16.0f64) as usize;
+            min(max(auto, 3), 64)
+        } else {
+            approximation_order
+        };
+
         FractalRenderer {
             image_width,
             image_height,
             aspect,
             zoom,
             center_location,
             maximum_iteration,
-            approximation_order,
+            approximation_order: auto_approximation,
             glitch_tolerance,
             image: Image::new(image_width, image_height, display_glitches)
         }
@@ -95,27 +101,27 @@ impl FractalRenderer {
 
         let time = Instant::now();
 
-        let indices = (0..self.image_width).cartesian_product(0..self.image_height);
-
-        let mut pixel_data = indices.into_iter().par_bridge()
-                            .map(|(i, j)| {
-                                let element = ComplexFixed::new(i as f64 * delta_pixel + delta_top_left.re, j as f64 * delta_pixel + delta_top_left.im);
-                                let point_delta = ComplexExtended::new(element, -self.zoom.exponent);
-                                let new_delta = series_approximation.evaluate(point_delta);
-
-                                PixelData {
-                                    image_x: i,
-                                    image_y: j,
-                                    iteration: reference.start_iteration,
-                                    delta_centre: point_delta,
-                                    delta_reference: point_delta,
-                                    delta_start: new_delta,
-                                    delta_current: new_delta,
-                                    derivative_current: ComplexFixed::new(1.0, 0.0),
-                                    glitched: false,
-                                    escaped: false
-                                }
-                            }).collect::<Vec<PixelData>>();
+        let mut pixel_data = (0..(self.image_width * self.image_height)).into_par_iter()
+            .map(|index| {
+                let i = index % self.image_width;
+                let j = index / self.image_width;
+                let element = ComplexFixed::new(i as f64 * delta_pixel + delta_top_left.re, j as f64 * delta_pixel + delta_top_left.im);
+                let point_delta = ComplexExtended::new(element, -self.zoom.exponent);
+                let new_delta = series_approximation.evaluate(point_delta);
+
+                PixelData {
+                    image_x: i,
+                    image_y: j,
+                    iteration: reference.start_iteration,
+                    delta_centre: point_delta,
+                    delta_reference: point_delta,
+                    delta_start: new_delta,
+                    delta_current: new_delta,
+                    derivative_current: ComplexFixed::new(1.0, 0.0),
+                    glitched: false,
+                    escaped: false
+                }
+            }).collect::<Vec<PixelData>>();
 
         println!("{:<14}{:>6} ms", "Packing", time.elapsed().as_millis());
 
@@ -148,18 +154,15 @@ impl FractalRenderer {
 
             // this can be made faster, without having to do the series approximation again
             // this is done by storing more data in pixeldata2
-            pixel_data.chunks_mut(1)
-                        .for_each(|pixel_data| {
-                            for data in pixel_data {
-                                data.iteration = reference.start_iteration;
-                                data.glitched = false;
-                                data.escaped = false;
-                                data.delta_current = data.delta_start - delta_z;
-                                data.delta_reference = data.delta_centre - reference_wrt_sa;
-                                // might not need the evaluate here as if we store it separately, there is no need
-                                data.derivative_current = ComplexFixed::new(1.0, 0.0);
-                            }
-                        });
+            pixel_data.par_iter_mut()
+                .for_each(|pixel| {
+                    pixel.iteration = reference.start_iteration;
+                    pixel.glitched = false;
+                    pixel.delta_current = pixel.delta_start - delta_z;
+                    pixel.delta_reference = pixel.delta_centre - reference_wrt_sa;
+                        // might not need the evaluate here as if we store it separately, there is no need
+                        // data.derivative_current = ComplexFixed::new(1.0, 0.0);
+                });
 
             Perturbation::iterate(&mut pixel_data, &r, r.current_iteration);