Merge pull request #8 from jackyarndley/extended

Extended

Author: jackyarndley

.gitignore

@@ -1,3 +1,4 @@
 /target
 **/*.rs.bk
-.idea
+.idea
+*.kfr

Cargo.lock

@@ -12,6 +12,8 @@ pbr = "1.0.2"
 rayon = "1.3.0"
 packed_simd = "0.3.3"
 rug = "1.7.0"
+float-extended = {path = "../rust-float-extended/"}
+itertools = "0.9.0"
 
 #Additional commands that can improve performance (maybe by around 5-10%)
 #[profile.release]

Cargo.toml

@@ -1,74 +1,67 @@
 use std::time::Instant;
 use rust_fractal::renderer::FractalRenderer;
+use std::fs;
+use std::io::{stdout, stdin, Write};
 
 fn main() {
     println!("Mandelbrot Renderer");
-    // let center = (
-    //     "-1.99996619445037030418434688506350579675531241540724851511761922944801584242342684381376129778868913812287046406560949864353810575744772166485672496092803920095332",
-    //     "+0.00000000000000000000000000000000030013824367909383240724973039775924987346831190773335270174257280120474975614823581185647299288414075519224186504978181625478529");
-    // let zoom = 2.3620330788506154104770818136626E157;
 
-   // let center = (
-   //     "0.0",
-   //     "0.0");
-   // let zoom = 1.0;
+    // TODO on some images, we may need to check the imaginary part of the floatexp to make sure that is not too large.
 
-   //  let center = (
-   //     "-0.749999987350877481864108888020013802837969258626230419972587823828734338471228477079750588709551510361714463695461745528645748607681279674273355384334270208362211787387351792878073779449767292692440",
-   //     "0.001000038688236832013124581230049849132759425863378894883003211011278068229551274712347955044740933397589760194545872789087012331273586364914484522575986336846199522726507205442204060303594956029930");
-   // let zoom = 3.7E191;
+    let mut s = String::new();
+    print!("File to render: ");
+    let _ = stdout().flush();
+    stdin().read_line(&mut s).expect("User did not enter a correct string.");
 
-   // let center = (
-   //     "-1.689573325432279612075987864633746594591438093139394112928000260",
-   //     "0.000000000000000000000000000000000145514706909258179374258");
-   // let zoom = 1.2980742146337048E34;
+    s.pop();
 
-    // let center = (
-    //     "-1.999999999138270118722935763129859470012913240693218269085000",
-    //     "-0.000000000000000000322680215517275822769282166130504217892606");
-    // let zoom = 1.63e030;
+    if !s.ends_with(".kfr") {
+        s.push_str(".kfr")
+    };
 
-    // let center = (
-    //     "-1.76904090125288168240284276810365472222460651367280300612031519864776551390840630119504817227707689276670981344397551593371805167279428144061791279408701383958080000000000000000000000002e+00",
-    //     "-3.1605400272558568759745268636087334444350515177571908056976983793515506960262149704114708190147368048572976192049103272074188848768835600195188510326583157513e-03");
-    // let zoom = 7.5E139;
+    println!("Reading: {}", s);
 
-    // let center = (
-    //     "-1.7685304554715107439359975226423287998937254896566074423088827415719065562223946222490318219590860313620613389226996388",
-    //     "0.0059445386042946892550190097831583153039526746987365071641624831873576346016435037438597231970921541448697901394734283301");
-    // let zoom = 2.00000000000E113;
+    let data = fs::read_to_string(s).expect("Unable to read file.");
 
-    // let center = (
-    //     "-0.74977496303986055825690204249855929432518723234763348374372443797716597482713646451452176727030299366505937952743219986590574037612687173829307656886118099999999999999999999999",
-    //     "-0.06638008421017090687747407166586874456952046519386815946449195256714882710925246541317934635067121915196760463342031197755105049354085337445454326026704199999999999999999999999");
-    // let zoom = 4.70017300819E149;
+    let mut center_re = "-0.75";
+    let mut center_im  = "0.0";
+    let mut zoom = "1E0";
+    let mut iterations = "1000";
 
-    // let center = (
-    //     "-1.99996619445037030418434688506350579675531241540724851511761922944801584242342684381376129778868913812287046406560949864353810575744772166485672496092803920095332176654846785710146671838403830455177782981641506704068474778442",
-    //     "-0.000000000000000000000000000000000300138243679093832407249730397759249873468311907733352701742572801204749756148235811856472992884140755192241865049781816254785289455481579902975133989545509287024987599289178526901586954397865");
-    // let zoom = 8.80307862787E167;
+    for line in data.lines() {
+        let mut parts = line.split_whitespace();
+        // println!("{}", line);
+        let temp = parts.next().unwrap();
 
-    let center = (
-        "-1.479796270901760024425279934770411525645551054432599517909807632824286254403907594526888466099962805022975196472549771681831234491695559852583955204986197762293872769474806903995564259040667568599770094383157857518790853771783314763231302599999999999999999999999998",
-        "-0.001199443952813447746281973233374468444560314114132538362037569205657422216739564521471119107626453330996365067987088146663639996715939831819152248618042255824652268918299630897525386638029428706473919823922522752497780934312003352081931299999999999999999999999998");
-    let zoom = 3.27799999998E235;
-
-    // let center = (
-    //     "-0.74962449737876168207862620426836138684529527812364025754481571424558286479672801698750203976779135608148687747196595174858125388297577788573082753210469806739721377901297451172624762450940529814713048377873612297220313016947039287469999999999999999",
-    //     "0.03427010874046016945172951545749474868051534672439111853174854738370595772935930171842999778713222794215453436628998200591914208207632674780780978496784843807690510401829301309069542198413574392425885166367231192087416338497065495509999999999999999");
-    // let zoom = 7.58065474756E227;
+        match temp {
+            "Re:" => {
+                center_re = parts.next().unwrap();
+            },
+            "Im:" => {
+                center_im = parts.next().unwrap();
+            }
+            "Zoom:" => {
+                zoom = parts.next().unwrap();
+            },
+            "Iterations:" => {
+                iterations = parts.next().unwrap();
+            }
+            _ => {}
+        }
+    }
 
+    println!("Zoom: {}", zoom);
 
     let mut renderer = FractalRenderer::new(
         1000,
         1000,
         zoom,
-        100000,
-        center.0,
-        center.1,
+        iterations.parse::<usize>().unwrap(),
+        center_re,
+        center_im,
         0.001,
-        false,
-        64
+        true,
+        16
     );
 
     let time = Instant::now();

output.png

@@ -1,5 +1,4 @@
-// pub mod renderer;
+#![feature(extra_log_consts)]
 pub mod renderer;
 pub mod util;
-pub mod colouring;
 pub mod math;

src/bin/main.rs

@@ -1,5 +1,6 @@
-// pub mod perturbation_f64;
-// pub mod perturbation_f64_vectorised;
-pub mod perturbation;
-pub mod reference;
-pub mod series_approximation;
+pub mod perturbation_double;
+pub mod perturbation_extended;
+pub mod reference_double;
+pub mod reference_extended;
+pub mod series_approximation_double;
+pub mod series_approximation_extended;

src/lib.rs

@@ -1,19 +1,16 @@
-use crate::math::reference::Reference;
-use crate::util::PixelData;
-
+use crate::math::reference_double::ReferenceDouble;
+use crate::util::PixelDataDouble;
 use rayon::prelude::*;
 
-pub struct Perturbation {}
+pub struct PerturbationDouble {}
 
-impl Perturbation {
-    pub fn iterate(pixel_data: &mut Vec<PixelData>, reference: &Reference, maximum_iteration: usize) {
-        pixel_data.par_chunks_mut(1)
+impl PerturbationDouble {
+    pub fn iterate(pixel_data: &mut Vec<PixelDataDouble>, reference: &ReferenceDouble, reference_current_iteration: usize) {
+        pixel_data.par_chunks_mut(4)
             .for_each(|pixel_data| {
                 for packet in pixel_data {
-                    // TODO investigate setting up a floatexp-type thing which only rescales every 500 or so iterations
-
                     // normal
-                    while packet.iteration < maximum_iteration {
+                    while packet.iteration < reference_current_iteration {
                         // This uses the difference between the starting iteration of the reference - can be used to skip some
                         let z = packet.delta_current + reference.z_reference[packet.iteration - reference.start_iteration];
                         let z_norm = z.norm_sqr();
@@ -24,13 +21,13 @@ impl Perturbation {
                             break;
                         }
 
-                        if z_norm > 65536.0 {
+                        if z_norm > 1e16 {
                             packet.escaped = true;
                             packet.delta_current = z;
                             break;
                         }
 
-                        packet.derivative_current = 2.0 * z * packet.derivative_current + 1.0;
+                        // packet.derivative_current = 2.0 * z * packet.derivative_current + 1.0;
                         packet.delta_current = 2.0 * reference.z_reference[packet.iteration - reference.start_iteration] * packet.delta_current + packet.delta_current * packet.delta_current + packet.delta_reference;
                         packet.iteration += 1;
                     }

src/math/mod.rs

@@ -0,0 +1,114 @@
+use crate::util::PixelDataExtended;
+use float_extended::util::FloatExp;
+
+use rayon::prelude::*;
+use crate::math::reference_extended::ReferenceExtended;
+use std::cmp::max;
+
+pub struct PerturbationExtended {}
+
+impl PerturbationExtended {
+    pub fn iterate(pixel_data: &mut Vec<PixelDataExtended>, reference: &ReferenceExtended, reference_current_iteration: usize) {
+        pixel_data.par_chunks_mut(4)
+            .for_each(|pixel_data| {
+                for packet in pixel_data {
+                    let original_reference = (packet.delta_reference, packet.p_initial);
+
+                    // normalise so that the current d0 is same as current dn
+                    packet.delta_reference *= 1.0.ldexp(packet.p_initial - packet.p_current);
+                    packet.p_initial = packet.p_current;
+                    let mut temp = 1.0.ldexp(packet.p_current);
+
+                    // normal but with added exponent
+                    while packet.iteration < reference_current_iteration {
+                        // packet.derivative_current = 2.0 * z * packet.derivative_current + 1.0;
+                        let z_reference = reference.z_reference[packet.iteration - reference.start_iteration];
+
+                        if packet.p_current > -500 {
+                            let mut z = 0.5 * z_reference.0;
+
+                            if z_reference.1 < 0 {
+                                z *= 1.0.ldexp(z_reference.1);
+                            }
+
+                            z += packet.delta_current * temp;
+                            let z_norm = z.norm_sqr();
+
+                            if z_norm < reference.z_tolerance[packet.iteration - reference.start_iteration] {
+                                packet.glitched = true;
+                                packet.delta_current = z;
+                                break;
+                            }
+
+                            if z_norm > 1e16 {
+                                packet.escaped = true;
+                                packet.delta_current = z;
+                                break;
+                            }
+
+                            // We have already done this check, we only need to check if the current value is in floatexp form
+                            if z_reference.1 < 0 {
+                                // we rescale back to the initial delta level
+
+                                // this is the new p value we will be using
+                                let new_p = max(packet.p_current + z_reference.1, original_reference.1);
+
+                                packet.delta_reference = original_reference.0 * 1.0.ldexp(original_reference.1 - new_p);
+
+                                let temp9 = z_reference.0 * 1.0.ldexp(packet.p_current + z_reference.1 - new_p) + packet.delta_current * 1.0.ldexp(2 * packet.p_current - new_p);
+
+                                packet.delta_current *= temp9;
+                                packet.delta_current += packet.delta_reference;
+                                packet.p_current = new_p;
+                                temp = 1.0.ldexp(packet.p_current);
+                            } else {
+                                // for the delta squared bit, as we want it in terms of the current p level, the powers are 2^(2 * pinit) / 2^(pinit)
+                                let temp2 = packet.delta_current * temp + z_reference.0;
+                                packet.delta_current *= temp2;
+                                packet.delta_current += packet.delta_reference;
+                            }
+                        } else {
+                            // We have already done this check, we only need to check if the current value is in floatexp form
+                            if z_reference.1 < 0 {
+                                // we rescale back to the initial delta level
+
+                                // this is the new p value we will be using
+                                let new_p = max(packet.p_current + z_reference.1, original_reference.1);
+
+                                packet.delta_reference = original_reference.0 * 1.0.ldexp(original_reference.1 - new_p);
+
+                                let temp9 = z_reference.0 * 1.0.ldexp(packet.p_current + z_reference.1 - new_p) + packet.delta_current * 1.0.ldexp(2 * packet.p_current - new_p);
+
+                                packet.delta_current *= temp9;
+                                packet.delta_current += packet.delta_reference;
+                                packet.p_current = new_p;
+                                temp = 1.0.ldexp(packet.p_current);
+                            } else {
+                                packet.delta_current *= z_reference.0;
+                                packet.delta_current += packet.delta_reference;
+                            }
+                        }
+
+                        packet.iteration += 1;
+
+                        // 2x time with 100 vs 1
+                        if packet.iteration % 400 == 0 {
+                            // here we should do the rescale, based on the real value
+                            let (temp_mantissa, added_exponent) = packet.delta_current.re.frexp();
+                            packet.delta_current.re = temp_mantissa;
+
+                            packet.delta_current.im = packet.delta_current.im.ldexp(-added_exponent);
+                            packet.p_current += added_exponent;
+                            packet.delta_reference.re = packet.delta_reference.re.ldexp(-added_exponent);
+                            packet.delta_reference.im = packet.delta_reference.im.ldexp(-added_exponent);
+                            temp = 1.0.ldexp(packet.p_current);
+                            // packet.derivative_current /= 1.0.ldexp(added_exponent);
+                        }
+                    }
+                }
+            });
+    }
+}
+
+
+