From 7f4f31a09f5c684ce835966f1d451cd60d8d4dfd Mon Sep 17 00:00:00 2001 From: supersurviveur Date: Tue, 15 Sep 2026 19:28:49 +0200 Subject: [PATCH] feat: Add more fractals --- src/app.rs | 47 ++++++++++- src/fractal/reference.rs | 156 ++++++++++++++++++++++++++++++++++-- src/fractal/renderer.rs | 6 +- src/fractal/share.rs | 19 ++++- src/shaders/mandelbrot.wgsl | 51 +++++++++++- src/worker.rs | 4 + 6 files changed, 269 insertions(+), 14 deletions(-) diff --git a/src/app.rs b/src/app.rs index a1de971..9d5a003 100644 --- a/src/app.rs +++ b/src/app.rs @@ -42,6 +42,10 @@ const KINDS: &[(FractalKind, &str)] = &[ (FractalKind::BurningShip, "Burning Ship"), (FractalKind::Tricorn, "Tricorn"), (FractalKind::Multibrot, "Multibrot"), + (FractalKind::Celtic, "Celtic"), + (FractalKind::Perpendicular, "Perpendicular"), + (FractalKind::Buffalo, "Buffalo"), + (FractalKind::Phoenix, "Phoenix"), ]; /// UI label for a fractal kind. @@ -99,6 +103,7 @@ struct RequestKey { half_height: f64, julia: bool, julia_c: (f64, f64), + phoenix_p: (f64, f64), iter: u32, kind: FractalKind, power: u32, @@ -123,6 +128,8 @@ pub struct FractalApp { /// Exponent for the Multibrot kind. power: u32, julia_c: (f64, f64), + /// Distortion constant `p` for the Phoenix kind (`z^2 + c + p·z_{n-1}`). + phoenix_p: (f64, f64), max_iterations: u32, /// When set, `max_iterations` tracks the zoom depth automatically (so deep /// zooms stay sharp without hand-tuning); the manual slider takes over when @@ -228,6 +235,7 @@ impl FractalApp { kind: FractalKind::Mandelbrot, power: 3, julia_c: (-0.8, 0.156), + phoenix_p: (-0.5, 0.0), max_iterations: 512, auto_iterations: true, color_scale: 0.15, @@ -272,6 +280,10 @@ impl FractalApp { "burningship" | "burning_ship" | "ship" => FractalKind::BurningShip, "tricorn" | "mandelbar" => FractalKind::Tricorn, "multibrot" | "multi" => FractalKind::Multibrot, + "celtic" => FractalKind::Celtic, + "perpendicular" | "perp" => FractalKind::Perpendicular, + "buffalo" => FractalKind::Buffalo, + "phoenix" => FractalKind::Phoenix, _ => FractalKind::Mandelbrot, }; if let Ok(p) = std::env::var("MANDEL_POWER") @@ -377,6 +389,7 @@ impl FractalApp { half_height: self.view.half_height, iterations: self.max_iterations, julia_c: self.julia_c, + phoenix_p: self.phoenix_p, color_scale: self.color_scale, color_offset: self.color_offset, palette: self.palette, @@ -393,6 +406,7 @@ impl FractalApp { self.kind = s.kind; self.power = s.power.clamp(2, 8); self.julia_c = s.julia_c; + self.phoenix_p = s.phoenix_p; self.color_scale = s.color_scale; self.color_offset = s.color_offset; self.palette = (s.palette as usize).min(PALETTE_NAMES.len() - 1) as u32; @@ -437,6 +451,10 @@ impl FractalApp { FractalKind::BurningShip => (-0.5, -0.5, 1.3), FractalKind::Tricorn => (-0.25, 0.0, 1.6), FractalKind::Multibrot => (0.0, 0.0, 1.5), + FractalKind::Celtic => (-0.5, 0.0, 1.6), + FractalKind::Perpendicular => (-0.5, 0.0, 1.5), + FractalKind::Buffalo => (-0.5, -0.5, 1.5), + FractalKind::Phoenix => (0.0, 0.0, 1.6), }; ViewState::with_center(big_from_f64(cr, 53), big_from_f64(ci, 53), hh) } @@ -448,6 +466,7 @@ impl FractalApp { half_height: self.view.half_height, julia: matches!(self.mode, FractalMode::Julia), julia_c: self.julia_c, + phoenix_p: self.phoenix_p, iter: self.max_iterations, kind: self.kind, power: self.power, @@ -470,6 +489,7 @@ impl FractalApp { }; if key.julia != matches!(self.mode, FractalMode::Julia) || key.julia_c != self.julia_c + || key.phoenix_p != self.phoenix_p || key.iter != self.max_iterations || key.kind != self.kind || key.power != self.power @@ -519,6 +539,7 @@ impl FractalApp { precision, kind: key.kind, power: key.power, + phoenix_p: key.phoenix_p, }); self.pending = true; } @@ -536,6 +557,7 @@ impl FractalApp { precision, key.kind, key.power, + key.phoenix_p, ) } else { compute_set_reference( @@ -545,6 +567,7 @@ impl FractalApp { precision, key.kind, key.power, + key.phoenix_p, ) }; self.apply_reference( @@ -580,8 +603,9 @@ impl FractalApp { kind: self.kind.shader_id(), power: self.power, dc_offset: self.dc_offset(), + phoenix_p: [self.phoenix_p.0 as f32, self.phoenix_p.1 as f32], de_coloring: self.de_coloring as u32, - _pad: 0, + _pad: [0, 0, 0], } } @@ -794,6 +818,22 @@ impl FractalApp { if self.kind == FractalKind::Multibrot { ui.add(egui::Slider::new(&mut self.power, 2..=8).text("power")); } + if self.kind == FractalKind::Phoenix { + ui.horizontal(|ui| { + ui.label("p ="); + ui.add( + egui::DragValue::new(&mut self.phoenix_p.0) + .speed(0.001) + .range(-2.0..=2.0), + ); + ui.add( + egui::DragValue::new(&mut self.phoenix_p.1) + .speed(0.001) + .range(-2.0..=2.0), + ); + ui.label("i"); + }); + } if self.kind != prev_kind { self.view = Self::default_view_for(self.mode, self.kind); } @@ -868,8 +908,9 @@ impl FractalApp { ui.checkbox(&mut self.de_coloring, "Distance shading") .on_hover_text( "Shade by distance to the set boundary (from the orbit derivative) \ - for crisp filaments at deep zoom. Exact for Mandelbrot/Multibrot, \ - approximate for Burning Ship/Tricorn.", + for crisp filaments at deep zoom. Exact for the holomorphic kinds \ + (Mandelbrot/Multibrot/Phoenix), approximate for the abs-based kinds \ + (Burning Ship/Tricorn/Celtic/Perpendicular/Buffalo).", ); ui.separator(); diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 697af30..2b118f6 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -10,7 +10,7 @@ //! * Mandelbrot-set: `z0 = 0`, `c = view center` (the c-plane point per pixel). //! * Julia-set: `z0 = view center`, `c = fractal constant` (fixed per view). -use crate::view::Big; +use crate::view::{Big, big_from_f64}; /// The iteration formula. Must be kept in sync with `advance_delta` and the /// `KIND_*` constants in the shader. @@ -24,6 +24,14 @@ pub enum FractalKind { Tricorn, /// `z -> z^power + c` (power >= 2). Multibrot, + /// `z -> |Re(z^2)| + i·Im(z^2) + c` (abs on the real output of the square). + Celtic, + /// `z -> (x^2 - y^2) - 2·x·|y|·i + c` (abs on the imaginary input). + Perpendicular, + /// `z -> |Re(z^2)| - |Im(z^2)|·i + c` (abs on both outputs). + Buffalo, + /// `z -> z^2 + c + p·z_{n-1}` (two-term recurrence; `p` is `phoenix_p`). + Phoenix, } impl FractalKind { @@ -34,6 +42,10 @@ impl FractalKind { FractalKind::BurningShip => 1, FractalKind::Tricorn => 2, FractalKind::Multibrot => 3, + FractalKind::Celtic => 4, + FractalKind::Perpendicular => 5, + FractalKind::Buffalo => 6, + FractalKind::Phoenix => 7, } } } @@ -55,12 +67,19 @@ pub fn compute_reference( precision: usize, kind: FractalKind, power: u32, + phoenix_p: (f64, f64), ) -> Vec<[f32; 2]> { let cr = c_re.clone().with_precision(precision).value(); let ci = c_im.clone().with_precision(precision).value(); let mut zr = z0_re.clone().with_precision(precision).value(); let mut zi = z0_im.clone().with_precision(precision).value(); + // Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0). + let mut zr_prev = big_zero(precision); + let mut zi_prev = big_zero(precision); + // Phoenix distortion constant `p` (a small fixed complex number). + let pr = big_from_f64(phoenix_p.0, precision); + let pi = big_from_f64(phoenix_p.1, precision); let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1); @@ -97,8 +116,37 @@ pub fn compute_reference( let (pr, pi) = complex_pow(&zr, &zi, power.max(2), precision); (pr + &cr, pi + &ci) } + FractalKind::Celtic => { + // |Re(z^2)| + i·Im(z^2): abs the real output of the square. + let re = big_abs(&zr.sqr() - &zi.sqr()) + &cr; + let im = ((&zr * &zi) << 1) + &ci; + (re, im) + } + FractalKind::Perpendicular => { + // (x^2 - y^2) - 2·x·|y| i: abs the imaginary input. + let re = &zr.sqr() - &zi.sqr() + &cr; + let im = &ci - ((&zr * &big_abs(zi.clone())) << 1); + (re, im) + } + FractalKind::Buffalo => { + // |Re(z^2)| - |Im(z^2)| i: abs both outputs. + let re = big_abs(&zr.sqr() - &zi.sqr()) + &cr; + let im = &ci - big_abs((&zr * &zi) << 1); + (re, im) + } + FractalKind::Phoenix => { + // z^2 + c + p·z_{n-1}. + let re2 = &zr.sqr() - &zi.sqr(); + let im2 = (&zr * &zi) << 1; + let pzr = &pr * &zr_prev - &pi * &zi_prev; + let pzi = &pr * &zi_prev + &pi * &zr_prev; + (re2 + &cr + pzr, im2 + &ci + pzi) + } }; + // Shift the previous iterate (only the Phoenix arm reads it). + zr_prev = zr; + zi_prev = zi; zr = new_zr.with_precision(precision).value(); zi = new_zi.with_precision(precision).value(); } @@ -139,10 +187,11 @@ pub fn compute_set_reference( precision: usize, kind: FractalKind, power: u32, + phoenix_p: (f64, f64), ) -> Vec<[f32; 2]> { let zero = big_zero(precision); compute_reference( - &zero, &zero, center_re, center_im, max_iter, precision, kind, power, + &zero, &zero, center_re, center_im, max_iter, precision, kind, power, phoenix_p, ) } @@ -156,7 +205,7 @@ mod tests { fn reference_matches_naive_f64() { let cr = Big::try_from(-0.75_f64).unwrap(); let ci = Big::try_from(0.1_f64).unwrap(); - let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Mandelbrot, 2); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Mandelbrot, 2, (0.0, 0.0)); // Independent naive f64 orbit. let (c_re, c_im) = (-0.75_f64, 0.1_f64); @@ -180,7 +229,8 @@ mod tests { fn interior_orbit_runs_full_length() { let cr = Big::try_from(-0.2_f64).unwrap(); let ci = Big::try_from(0.0_f64).unwrap(); - let points = compute_set_reference(&cr, &ci, 500, 120, FractalKind::Mandelbrot, 2); + let points = + compute_set_reference(&cr, &ci, 500, 120, FractalKind::Mandelbrot, 2, (0.0, 0.0)); assert_eq!(points.len(), 501, "interior orbit should not escape"); } @@ -189,7 +239,8 @@ mod tests { fn burning_ship_reference_matches_naive_f64() { let cr = Big::try_from(-1.75_f64).unwrap(); let ci = Big::try_from(-0.03_f64).unwrap(); - let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::BurningShip, 2); + let points = + compute_set_reference(&cr, &ci, 60, 200, FractalKind::BurningShip, 2, (0.0, 0.0)); let (c_re, c_im) = (-1.75_f64, -0.03_f64); let (mut zr, mut zi) = (0.0_f64, 0.0_f64); @@ -209,7 +260,8 @@ mod tests { fn multibrot3_reference_matches_naive_f64() { let cr = Big::try_from(0.3_f64).unwrap(); let ci = Big::try_from(0.2_f64).unwrap(); - let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Multibrot, 3); + let points = + compute_set_reference(&cr, &ci, 60, 200, FractalKind::Multibrot, 3, (0.0, 0.0)); let (c_re, c_im) = (0.3_f64, 0.2_f64); let (mut zr, mut zi) = (0.0_f64, 0.0_f64); @@ -242,6 +294,7 @@ mod tests { 200, FractalKind::Mandelbrot, 2, + (0.0, 0.0), ); let (mut zr, mut zi) = (0.15_f64, -0.1_f64); @@ -256,4 +309,95 @@ mod tests { zi = nzi; } } + + /// Celtic reference matches a naive f64 iteration: real = |x^2 - y^2| + cr. + #[test] + fn celtic_reference_matches_naive_f64() { + let cr = Big::try_from(-0.6_f64).unwrap(); + let ci = Big::try_from(0.4_f64).unwrap(); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Celtic, 2, (0.0, 0.0)); + + let (c_re, c_im) = (-0.6_f64, 0.4_f64); + let (mut zr, mut zi) = (0.0_f64, 0.0_f64); + for point in &points { + let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs())); + assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}"); + assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}"); + let nzr = (zr * zr - zi * zi).abs() + c_re; + let nzi = 2.0 * zr * zi + c_im; + zr = nzr; + zi = nzi; + } + } + + /// Perpendicular reference matches a naive f64 iteration: + /// real = x^2 - y^2 + cr, imag = -2·x·|y| + ci. + #[test] + fn perpendicular_reference_matches_naive_f64() { + let cr = Big::try_from(-0.7_f64).unwrap(); + let ci = Big::try_from(-0.2_f64).unwrap(); + let points = + compute_set_reference(&cr, &ci, 60, 200, FractalKind::Perpendicular, 2, (0.0, 0.0)); + + let (c_re, c_im) = (-0.7_f64, -0.2_f64); + let (mut zr, mut zi) = (0.0_f64, 0.0_f64); + for point in &points { + let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs())); + assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}"); + assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}"); + let nzr = zr * zr - zi * zi + c_re; + let nzi = -2.0 * zr * zi.abs() + c_im; + zr = nzr; + zi = nzi; + } + } + + /// Buffalo reference matches a naive f64 iteration: + /// real = |x^2 - y^2| + cr, imag = -|2·x·y| + ci. + #[test] + fn buffalo_reference_matches_naive_f64() { + let cr = Big::try_from(-1.2_f64).unwrap(); + let ci = Big::try_from(-0.35_f64).unwrap(); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Buffalo, 2, (0.0, 0.0)); + + let (c_re, c_im) = (-1.2_f64, -0.35_f64); + let (mut zr, mut zi) = (0.0_f64, 0.0_f64); + for point in &points { + let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs())); + assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}"); + assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}"); + let nzr = (zr * zr - zi * zi).abs() + c_re; + let nzi = -(2.0 * zr * zi).abs() + c_im; + zr = nzr; + zi = nzi; + } + } + + /// Phoenix reference matches a naive f64 two-term iteration + /// `z_{n+1} = z_n^2 + c + p·z_{n-1}` (z_0 = 0, z_{-1} = 0). + #[test] + fn phoenix_reference_matches_naive_f64() { + let cr = Big::try_from(0.5667_f64).unwrap(); + let ci = Big::try_from(0.0_f64).unwrap(); + let p = (-0.5_f64, 0.0_f64); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Phoenix, 2, p); + + let (c_re, c_im) = (0.5667_f64, 0.0_f64); + let (mut zr, mut zi) = (0.0_f64, 0.0_f64); + let (mut pr, mut pi) = (0.0_f64, 0.0_f64); // previous iterate + for point in &points { + let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs())); + assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}"); + assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}"); + // p·z_{n-1} = (p.0 + i p.1)(pr + i pi). + let pzr = p.0 * pr - p.1 * pi; + let pzi = p.0 * pi + p.1 * pr; + let nzr = zr * zr - zi * zi + c_re + pzr; + let nzi = 2.0 * zr * zi + c_im + pzi; + pr = zr; + pi = zi; + zr = nzr; + zi = nzi; + } + } } diff --git a/src/fractal/renderer.rs b/src/fractal/renderer.rs index 31c4bb0..bcf7eb4 100644 --- a/src/fractal/renderer.rs +++ b/src/fractal/renderer.rs @@ -43,10 +43,14 @@ pub struct Uniforms { /// or reused reference (computed at a slightly different center) still maps /// correctly. Added to every pixel's per-pixel offset. pub dc_offset: [f32; 2], + /// Distortion constant `p` for the Phoenix map (`z^2 + c + p·z_{n-1}`); + /// ignored by other kinds. Kept next to `dc_offset` so both `vec2`s land on + /// 8-byte boundaries, matching the shader's layout. + pub phoenix_p: [f32; 2], /// 0 = escape-time coloring, 1 = distance-estimation shading. pub de_coloring: u32, /// Padding to a 16-byte multiple (uniform buffer requirement). - pub _pad: u32, + pub _pad: [u32; 3], } /// Offscreen texture the fractal is rendered into, plus the bind group used to diff --git a/src/fractal/share.rs b/src/fractal/share.rs index 9fe9da9..3bf8080 100644 --- a/src/fractal/share.rs +++ b/src/fractal/share.rs @@ -21,6 +21,8 @@ pub struct ShareState { pub half_height: f64, pub iterations: u32, pub julia_c: (f64, f64), + /// Distortion constant for the Phoenix kind (ignored by others). + pub phoenix_p: (f64, f64), pub color_scale: f32, pub color_offset: f32, /// Palette index (`palette_id` in the shader). @@ -36,6 +38,10 @@ impl ShareState { FractalKind::BurningShip => "burning", FractalKind::Multibrot => "multi", FractalKind::Tricorn => "tricorn", + FractalKind::Celtic => "celtic", + FractalKind::Perpendicular => "perp", + FractalKind::Buffalo => "buffalo", + FractalKind::Phoenix => "phoenix", })); s.push_str(&format!("&pw={}", self.power)); s.push_str(&format!( @@ -43,6 +49,7 @@ impl ShareState { self.center_re, self.center_im, self.half_height, self.iterations )); s.push_str(&format!("&jr={}&ji={}", self.julia_c.0, self.julia_c.1)); + s.push_str(&format!("&px={}&py={}", self.phoenix_p.0, self.phoenix_p.1)); s.push_str(&format!( "&cs={}&co={}&pal={}", self.color_scale, self.color_offset, self.palette @@ -68,6 +75,10 @@ impl ShareState { "multi" => FractalKind::Multibrot, "burning" => FractalKind::BurningShip, "tricorn" => FractalKind::Tricorn, + "celtic" => FractalKind::Celtic, + "perp" => FractalKind::Perpendicular, + "buffalo" => FractalKind::Buffalo, + "phoenix" => FractalKind::Phoenix, _ => FractalKind::Mandelbrot }) .unwrap_or(FractalKind::Mandelbrot), @@ -80,6 +91,10 @@ impl ShareState { map.get("jr").and_then(|s| s.parse().ok()).unwrap_or(-0.8), map.get("ji").and_then(|s| s.parse().ok()).unwrap_or(0.156), ), + phoenix_p: ( + map.get("px").and_then(|s| s.parse().ok()).unwrap_or(-0.5), + map.get("py").and_then(|s| s.parse().ok()).unwrap_or(0.0), + ), color_scale: map.get("cs").and_then(|s| s.parse().ok()).unwrap_or(0.02), color_offset: map.get("co").and_then(|s| s.parse().ok()).unwrap_or(0.0), palette: map.get("pal").and_then(|s| s.parse().ok()).unwrap_or(0), @@ -95,13 +110,14 @@ mod tests { fn round_trip() { let s = ShareState { julia: true, - kind: FractalKind::Multibrot, + kind: FractalKind::Phoenix, power: 5, center_re: "-0.743643887037158704752191506114774".into(), center_im: "0.131825904205311970493132056385139".into(), half_height: 1.5e-20, iterations: 4000, julia_c: (-0.123, 0.745), + phoenix_p: (-0.5, 0.1), color_scale: 0.02, color_offset: 0.25, palette: 3, @@ -115,6 +131,7 @@ mod tests { assert_eq!(d.half_height, s.half_height); assert_eq!(d.iterations, s.iterations); assert_eq!(d.julia_c, s.julia_c); + assert_eq!(d.phoenix_p, s.phoenix_p); assert_eq!(d.palette, s.palette); } diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index fbec524..a24cae2 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -22,11 +22,14 @@ struct Uniforms { is_julia: u32, palette_id: u32, aa_level: u32, - // Iteration formula: 0 Mandelbrot, 1 Burning Ship, 2 Tricorn, 3 Multibrot. + // Iteration formula (see the KIND_* constants below). kind: u32, // Exponent for the Multibrot kind. power: u32, dc_offset: vec2, + // Distortion constant p for the Phoenix map (z^2 + c + p*z_{n-1}); unused + // by other kinds. Placed by dc_offset so both vec2s stay 8-byte aligned. + phoenix_p: vec2, // 0 = escape-time coloring, 1 = distance-estimation shading. de_coloring: u32, }; @@ -35,6 +38,10 @@ const KIND_MANDELBROT: u32 = 0u; const KIND_BURNING_SHIP: u32 = 1u; const KIND_TRICORN: u32 = 2u; const KIND_MULTIBROT: u32 = 3u; +const KIND_CELTIC: u32 = 4u; +const KIND_PERPENDICULAR: u32 = 5u; +const KIND_BUFFALO: u32 = 6u; +const KIND_PHOENIX: u32 = 7u; @group(0) @binding(0) var u: Uniforms; @group(0) @binding(1) var ref_orbit: array>; @@ -128,8 +135,25 @@ fn advance_delta(z: vec2, e: vec2) -> vec2 { return 2.0 * cmul(cz, ce) + cmul(ce, ce); } else if (u.kind == KIND_MULTIBROT) { return multibrot_delta(z, e, clamp(u.power, 2u, 8u)); + } else if (u.kind == KIND_CELTIC) { + // z^2 delta split: sq.x = delta of Re(z^2), sq.y = delta of Im(z^2). + // Celtic abs the real output, so |Re(z^2)| delta = diffabs(Re(Z^2), sq.x). + let sq = 2.0 * cmul(z, e) + cmul(e, e); + return vec2(diffabs(z.x * z.x - z.y * z.y, sq.x), sq.y); + } else if (u.kind == KIND_BUFFALO) { + // Abs both outputs: real |Re(z^2)|, imag -|Im(z^2)| (Im(Z^2) = 2 X Y). + let sq = 2.0 * cmul(z, e) + cmul(e, e); + return vec2(diffabs(z.x * z.x - z.y * z.y, sq.x), + -diffabs(2.0 * z.x * z.y, sq.y)); + } else if (u.kind == KIND_PERPENDICULAR) { + // real x^2 - y^2 (ordinary square delta), imag -2 x |y|. + // d(-2 x |y|) = -2[ X·(|Y+ey|-|Y|) + ex·|Y+ey| ]; diffabs gives |Y+ey|-|Y|. + let sq = 2.0 * cmul(z, e) + cmul(e, e); + let da = diffabs(z.y, e.y); // |Y + ey| - |Y| + let abs_yf = abs(z.y) + da; // |Y + ey| + return vec2(sq.x, -2.0 * (z.x * da + e.x * abs_yf)); } - return 2.0 * cmul(z, e) + cmul(e, e); // Mandelbrot + return 2.0 * cmul(z, e) + cmul(e, e); // Mandelbrot (and Phoenix square part) } // Derivative f'(Z) of the iteration map at the full value Z, used to propagate @@ -183,6 +207,10 @@ fn shade(offset: vec2, px: f32) -> vec3 { // (starts at 0, gains +1 each step); for Julia it is d/dz0 (starts at 1). var dz = vec2(0.0, 0.0); var dz_seed = vec2(1.0, 0.0); + // Previous-iterate state for the Phoenix two-term recurrence (delta of + // y_{n-1}, and its derivative for DE). Both start at 0 (y_{-1} = 0). + var e_prev = vec2(0.0, 0.0); + var dz_prev = vec2(0.0, 0.0); if (u.is_julia != 0u) { step_add = vec2(0.0, 0.0); e = offset; @@ -210,12 +238,24 @@ fn shade(offset: vec2, px: f32) -> vec3 { // Propagate the derivative of the full orbit (unaffected by rebasing, // which only re-expresses the same value). Only when DE is enabled. + // Phoenix's two-term map adds p·dz_{n-1} and carries the previous dz. if (u.de_coloring != 0u) { - dz = cmul(fprime(z), dz) + dz_seed; + var dz_new = cmul(fprime(z), dz) + dz_seed; + if (u.kind == KIND_PHOENIX) { + dz_new = dz_new + cmul(u.phoenix_p, dz_prev); + dz_prev = dz; + } + dz = dz_new; } // Advance the delta by this fractal's formula (+ dc for the set plane). + // Phoenix additionally adds p·e_{n-1} and carries the previous delta. + let e_old = e; e = advance_delta(xm, e) + step_add; + if (u.kind == KIND_PHOENIX) { + e = e + cmul(u.phoenix_p, e_prev); + e_prev = e_old; + } m = m + 1u; n = n + 1u; @@ -231,6 +271,11 @@ fn shade(offset: vec2, px: f32) -> vec3 { if (dot(y, y) < dot(e, e)) { // Rebase to index 0: carry the full value as the new delta. Valid // because y_n = X[0] + (y_n - X[0]); for Mandelbrot X[0]=0. + // Phoenix: after rebasing the implied previous reference is Y[-1]=0, + // so the previous delta becomes the full previous value y_n (= z). + if (u.kind == KIND_PHOENIX) { + e_prev = z; + } e = y - z0; m = 0u; } diff --git a/src/worker.rs b/src/worker.rs index 5ff4327..790efab 100644 --- a/src/worker.rs +++ b/src/worker.rs @@ -22,6 +22,8 @@ pub struct RefRequest { pub precision: usize, pub kind: FractalKind, pub power: u32, + /// Distortion constant for the Phoenix map (ignored by other kinds). + pub phoenix_p: (f64, f64), } pub struct RefResult { @@ -101,6 +103,7 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> { req.precision, req.kind, req.power, + req.phoenix_p, ) } else { compute_set_reference( @@ -110,6 +113,7 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> { req.precision, req.kind, req.power, + req.phoenix_p, ) } }