surv/mandelbrot
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: Add arbitrary precision numbers at deep zooms

Browse Source
This commit is contained in:
surv
2026-09-25 15:30:05 +02:00
parent fbf0ef64af
commit bfd3d18f3d
11 changed files with 954 additions and 104 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/shaders/iterate_uniforms.wgsl
+4
View File
@@ -48,6 +48,10 @@ struct Uniforms {
// Kind-switch morph weight: each step is (1 - w)*f_kind + w*f_morph_from;
// 0 = no morph. Only read by MORPH pipelines (see mandelbrot.wgsl).
morph_w: f32,
// Deep views: binary exponent E of the view scale. `span` and `dc_offset`
// are uploaded multiplied by 2^-E so they stay in f32's range; 0 = not
// deep (plain f32 values). Only read by DEEP pipelines (mandelbrot.wgsl).
scale_exp: i32,
// Complex binomial coefficients C(complex_power, k) for k = 1..16, two per
// vec4 (k odd in .xy, k even in .zw), for the Complex Multibrot delta
// series. Precomputed on the CPU since they only depend on the power.
src/shaders/mandelbrot.wgsl
+587 -9
View File
@@ -11,6 +11,11 @@
// true value |y| drops below the delta |e|, or the reference runs out, we reset
// the reference index to 0 and carry the full value as the new delta (valid
// because X_0 = 0).
//
// Past ~1e30 zoom the deltas themselves leave f32's exponent range (smallest
// normal ~1.2e-38), so `DEEP` pipelines start each pixel in a rescaled form,
// e = w * 2^s with an f32 mantissa `w` and an i32 exponent `s`, and hand over
// to the plain f32 loop once the delta is big enough (see `iterate_sample`).
@group(0) @binding(0) var<uniform> u: Uniforms;
@group(0) @binding(1) var<storage, read> ref_orbit: array<vec2<f32>>;
@@ -20,6 +25,11 @@
// Only read by the adaptive-AA refine pass (`fs_refine`): the 1-sample-per-
// pixel data texture written by `fs_data`, which decides where to supersample.
@group(1) @binding(0) var coarse_tex: texture_2d<f32>;
// Per-point binary exponents of the reference orbit (`RefOrbit::exps`): the
// true X[m] is ref_orbit[m] * 2^ref_exp[m]. Non-zero only for points below
// f32's range, which only deep references contain; only `DEEP` pipelines read
// it (see `ref_at`).
@group(0) @binding(3) var<storage, read> ref_exp: array<i32>;
// Pipeline-overridable specialization constants, set per pipeline from the
// uniforms' `kind` / `is_julia` / `de_coloring` (see `PipelineKey` in
@@ -35,6 +45,10 @@ override DE: bool = false;
// second kind, `u.morph_from`. That one is a runtime value (it only lives for
// the length of the animation), so only MORPH pipelines pay for its branches.
override MORPH: bool = false;
// Deep view (`u.scale_exp != 0`): the per-pixel offset `dc`, the pixel size
// and `u.span` / `u.dc_offset` are all in units of 2^scale_exp, and each pixel
// starts in the rescaled deep phase (see `iterate_sample`).
override DEEP: bool = false;
struct VsOut {
@builtin(position) pos: vec4<f32>,
@@ -82,14 +96,18 @@ fn diffabs(c: f32, d: f32) -> f32 {
// tiny, but that only perturbs `s` by a relative f32 epsilon, and the result
// is `e * s`, so the delta keeps full relative precision.
fn multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: u32) -> vec2<f32> {
let y = z + e;
return cmul(e, multibrot_sum(z, z + e, p));
}
// The sum in `multibrot_delta`, sum_{k=0}^{p-1} y^k Z^{p-1-k} with y = Z+e.
fn multibrot_sum(z: vec2<f32>, y: vec2<f32>, p: u32) -> vec2<f32> {
var s = vec2<f32>(1.0, 0.0);
var zj = vec2<f32>(1.0, 0.0);
for (var j: u32 = 1u; j < p; j = j + 1u) {
zj = cmul(zj, z); // Z^j
s = cmul(s, y) + zj;
}
return cmul(e, s);
return s;
}
// Maximum number of terms in `complex_multibrot_delta`'s series (matches the
@@ -242,6 +260,324 @@ fn fprime(z: vec2<f32>) -> vec2<f32> {
return d;
}
// ---- Deep (rescaled) phase helpers -------------------------------------
//
// A deep value is an f32 mantissa times 2^exponent. All scaling is by exact
// powers of two, so it never rounds.
const LN2: f32 = 0.6931471805599453;
// Where `ldexp_sat` saturates. Only ever compared against small values
// (see `diffabs_scaled`), and far enough from f32's max that doubling it,
// or squaring a value of the size it's compared with, stays finite.
const LDEXP_SAT: f32 = 1.2676506e30; // 2^100
// x * 2^k for any k. WGSL's `ldexp` is only defined for exponents inside
// f32's range, so go through `frexp`: results below the smallest normal
// flush to 0, results above 2^100 saturate to ±2^100.
fn ldexp_sat(x: f32, k: i32) -> f32 {
if x == 0.0 {
return 0.0;
}
let f = frexp(x);
let ex = f.exp + k;
if ex > 100 {
return select(-LDEXP_SAT, LDEXP_SAT, x > 0.0);
}
if ex < -125 {
return 0.0;
}
return ldexp(f.fract, ex);
}
fn ldexp2_sat(v: vec2<f32>, k: i32) -> vec2<f32> {
return vec2<f32>(ldexp_sat(v.x, k), ldexp_sat(v.y, k));
}
// A complex number as mantissa * 2^e, with the mantissa's larger component in
// [0.5, 1) (or exactly 0, with e = 0).
struct Fe {
m: vec2<f32>,
e: i32,
};
fn fe_make(v: vec2<f32>, e: i32) -> Fe {
let a = max(abs(v.x), abs(v.y));
if a == 0.0 {
return Fe(vec2<f32>(0.0, 0.0), 0);
}
let k = frexp(a).exp;
return Fe(ldexp2_sat(v, -k), e + k);
}
// Reference point X[m] as f32 (points stored normalized flush to 0 here;
// only deep references have any, see `ref_exp`).
fn ref_at(m: u32) -> vec2<f32> {
let x = ref_orbit[m];
if DEEP {
let k = ref_exp[m];
if k != 0 {
return ldexp2_sat(x, k);
}
}
return x;
}
// X[m] with its full exponent range (deep phase only).
fn ref_fe(m: u32) -> Fe {
return fe_make(ref_orbit[m], ref_exp[m]);
}
// Complex log of m * 2^e (m != 0).
fn clog_fe(m: vec2<f32>, e: i32) -> vec2<f32> {
return vec2<f32>(0.5 * log(dot(m, m)) + f32(e) * LN2, atan2(m.y, m.x));
}
fn cexp(a: vec2<f32>) -> vec2<f32> {
return exp(a.x) * vec2<f32>(cos(a.y), sin(a.y));
}
// diffabs(c, 2^s * d) / 2^s = diffabs(c / 2^s, d): diffabs is positively
// homogeneous. Saturating c / 2^s is harmless: once it dwarfs |d| the result
// is just ±d.
fn diffabs_scaled(c: f32, d: f32, s: i32) -> f32 {
return diffabs(ldexp_sat(c, -s), d);
}
// Sum of two `Fe`s, at the larger one's exponent.
fn fe_add(a: Fe, b: Fe) -> Fe {
if a.m.x == 0.0 && a.m.y == 0.0 {
return b;
}
if b.m.x == 0.0 && b.m.y == 0.0 {
return a;
}
let e = max(a.e, b.e);
return fe_make(ldexp2_sat(a.m, a.e - e) + ldexp2_sat(b.m, b.e - e), e);
}
// One deep step's delta, (f(X + e) - f(X)) / 2^t: the mantissa `w` and its
// exponent `t`. `t` is the input scale s except next to the critical point,
// where the linear part of the step vanishes and the result is ~e^2, far
// below the input's scale.
struct DeepStep {
w: vec2<f32>,
t: i32,
};
// The per-kind formula of a deep step: (f(X + e) - f(X)) / 2^t for the
// delta e = w * 2^s, given X measured in units of 2^u (`x` = X / 2^u) and
// `sc` = 2^se, se = s - u, the delta's scale in those units. The result's
// scale is t = s + (p-1)·u for a degree-p kind, since every kind here but
// Lambda is p-homogeneous in (X, e) jointly (`diffabs_scaled` rescales the
// fold-point comparisons the same way). Usually u = 0 (x = X, sc = 2^s,
// t = s); see `deep_step_kind` for when it isn't. Lambda always gets u = 0.
// Every kind needs an arm here too.
fn advance_delta_scaled_kind(kind: u32, x: vec2<f32>, w: vec2<f32>, sc: f32, se: i32) -> vec2<f32> {
if kind == KIND_BURNING_SHIP {
let base = 2.0 * cmul(x, w) + sc * cmul(w, w);
let dp = x.x * w.y + x.y * w.x + sc * w.x * w.y;
return vec2<f32>(base.x, 2.0 * diffabs_scaled(x.x * x.y, dp, se));
} else if kind == KIND_TRICORN {
let cx = conj(x);
let cw = conj(w);
return 2.0 * cmul(cx, cw) + sc * cmul(cw, cw);
} else if kind == KIND_MULTIBROT {
return cmul(w, multibrot_sum(x, x + sc * w, clamp(u.power, 2u, 8u)));
} else if kind == KIND_CELTIC {
let sq = 2.0 * cmul(x, w) + sc * cmul(w, w);
return vec2<f32>(diffabs_scaled(x.x * x.x - x.y * x.y, sq.x, se), sq.y);
} else if kind == KIND_BUFFALO {
let sq = 2.0 * cmul(x, w) + sc * cmul(w, w);
return vec2<f32>(diffabs_scaled(x.x * x.x - x.y * x.y, sq.x, se),
-diffabs_scaled(2.0 * x.x * x.y, sq.y, se));
} else if kind == KIND_PERPENDICULAR {
let sq = 2.0 * cmul(x, w) + sc * cmul(w, w);
let da = diffabs_scaled(x.y, w.y, se); // (|Y + ey| - |Y|) / 2^s
let abs_yf = abs(x.y) + sc * da; // |Y + ey| * 2^(s-t)
return vec2<f32>(sq.x, -2.0 * (x.x * da + w.x * abs_yf));
} else if kind == KIND_LAMBDA {
let t = vec2<f32>(1.0 - 2.0 * x.x - sc * w.x, -2.0 * x.y - sc * w.y);
return cmul(u.lambda_l, cmul(w, t));
}
return 2.0 * cmul(x, w) + sc * cmul(w, w); // Mandelbrot (and Phoenix square part)
}
// Below this exponent a reference point counts as next to the critical point
// 0 (see `deep_step_kind`). Above it, the e^2 terms 2^s * w^2 can only flush
// to 0 when they are below 2^-50 of the linear ones.
const DEEP_X_NEAR_LOG2: i32 = -60;
// One deep step of `kind` for the delta e = w * 2^s from X (`x` as f32, `xf`
// at full range). Usually the input scale is kept (t = s). But when X is
// tiny (next to the critical point, e.g. at a minibrot's period), the linear
// term vanishes and the step's value is ~e^p, which would flush to 0 at
// scale s. There every z^p-like kind is p-homogeneous in (X, e) jointly, so
// both are measured in units of 2^u (u = the larger one's exponent) and the
// result lands at t = s + (p-1)·u (`ue` below).
fn deep_step_kind(kind: u32, x: vec2<f32>, xf: Fe, w: vec2<f32>, sc: f32, s: i32) -> DeepStep {
if kind == KIND_COMPLEX_MULTIBROT {
return complex_multibrot_step(xf, w, s);
}
let x_zero = xf.m.x == 0.0 && xf.m.y == 0.0;
// Lambda's critical point is 1/2 and its step has a constant linear
// term (λ·e), so it never needs this.
if kind == KIND_LAMBDA || (!x_zero && xf.e >= DEEP_X_NEAR_LOG2) {
return DeepStep(advance_delta_scaled_kind(kind, x, w, sc, s), s);
}
let kw = deep_log2(w, vec2<f32>(0.0, 0.0));
if kw == DEEP_ZERO {
return DeepStep(w, s); // e = 0: f(X) - f(X)
}
var ue = s + kw;
if !x_zero {
ue = max(ue, xf.e);
}
var deg = 2;
if kind == KIND_MULTIBROT {
deg = i32(clamp(u.power, 2u, 8u));
}
let xk = ldexp2_sat(xf.m, xf.e - ue);
let se = s - ue;
let dw = advance_delta_scaled_kind(kind, xk, w, ldexp_sat(1.0, se), se);
return DeepStep(dw, s + (deg - 1) * ue);
}
// Scaled `complex_multibrot_delta` with X as a full-range `Fe`. Same series
// as the f32 version, rewritten as X^(p-1) * w * sum_k C(p,k) r^(k-1)
// (r = e/X) so nothing is formed at the delta's true scale; X^(p-1)'s own
// exponent goes into the result's `t`. When |e/X| >= 0.5, X is itself tiny
// (|X| <= 2|e|), so both terms of the direct form are taken in log space at
// the larger one's scale. Branch-cut crossings leave the deep phase before
// stepping (`deep_cut_crossing`).
fn complex_multibrot_step(xf: Fe, w: vec2<f32>, s: i32) -> DeepStep {
let p = u.complex_power;
let wf = fe_make(w, s);
if wf.m.x == 0.0 && wf.m.y == 0.0 {
return DeepStep(vec2<f32>(0.0, 0.0), s);
}
if xf.m.x == 0.0 && xf.m.y == 0.0 {
// e^p.
let l = cmul(p, clog_fe(wf.m, wf.e));
let k = i32(floor(l.x / LN2));
return DeepStep(cexp(l - vec2<f32>(f32(k) * LN2, 0.0)), k);
}
let r = ldexp2_sat(cdiv(wf.m, xf.m), wf.e - xf.e);
if dot(r, r) < 0.25 {
var acc = cm_coef(1u);
var rk = r; // r^(k-1)
for (var k: u32 = 2u; k <= COMPLEX_MULTIBROT_TERMS; k = k + 1u) {
acc = acc + cmul(cm_coef(k), rk);
rk = cmul(rk, r);
if dot(rk, rk) < 1e-18 * dot(acc, acc) {
break;
}
}
// X^(p-1) = cexp(l) = cexp(l - k·ln2) * 2^k.
let l = cmul(p - vec2<f32>(1.0, 0.0), clog_fe(xf.m, xf.e));
let k = i32(floor(l.x / LN2));
let x_pm1 = cexp(l - vec2<f32>(f32(k) * LN2, 0.0));
return DeepStep(cmul(cmul(w, x_pm1), acc), s + k);
}
// (X + e)^p - X^p, both in log space (X + e may be exactly 0).
let xs = ldexp2_sat(xf.m, xf.e - s);
let yf = fe_make(xs + w, s);
let lb = cmul(p, clog_fe(xf.m, xf.e));
var k = i32(floor(lb.x / LN2));
var la = vec2<f32>(0.0, 0.0);
let y_zero = yf.m.x == 0.0 && yf.m.y == 0.0;
if !y_zero {
la = cmul(p, clog_fe(yf.m, yf.e));
k = max(k, i32(floor(la.x / LN2)));
}
let kl = vec2<f32>(f32(k) * LN2, 0.0);
var ya = vec2<f32>(0.0, 0.0);
if !y_zero {
ya = cexp(la - kl);
}
return DeepStep(ya - cexp(lb - kl), k);
}
// Deep-phase twin of `advance_delta` (same morph blend, at the larger of the
// two kinds' output scales).
fn advance_delta_scaled(x: vec2<f32>, xf: Fe, w: vec2<f32>, sc: f32, s: i32) -> DeepStep {
let a = deep_step_kind(KIND, x, xf, w, sc, s);
if MORPH {
let b = deep_step_kind(u.morph_from, x, xf, w, sc, s);
let t = max(a.t, b.t);
return DeepStep(mix(ldexp2_sat(a.w, a.t - t), ldexp2_sat(b.w, b.t - t), u.morph_w), t);
}
return a;
}
// Whether this step would take Complex Multibrot's X + e across the branch
// cut (see `complex_multibrot_delta`). The delta then jumps to the size of X,
// so the deep phase ends and the f32 loop takes the step.
fn deep_cut_crossing(xf: Fe, w: vec2<f32>, s: i32) -> bool {
let cm = KIND == KIND_COMPLEX_MULTIBROT || (MORPH && u.morph_from == KIND_COMPLEX_MULTIBROT);
if !cm {
return false;
}
let yn = xf.m + ldexp2_sat(w, s - xf.e); // (X + e) / 2^xe
return xf.m.x < 0.0 && ((xf.m.y < 0.0) != (yn.y < 0.0));
}
// Binary exponent of the largest component of a pair of complex mantissas
// (`frexp` convention: |x| < 2^k), or DEEP_ZERO when both are exactly 0.
const DEEP_ZERO: i32 = -100000;
fn deep_log2(a: vec2<f32>, b: vec2<f32>) -> i32 {
let m = max(max(abs(a.x), abs(a.y)), max(abs(b.x), abs(b.y)));
if m == 0.0 {
return DEEP_ZERO;
}
return frexp(m).exp;
}
// f'(y) at the full value y = X + w * 2^s, as an `Fe`. Plain `fprime` in f32
// unless y is below f32's comfortable range, which only happens next to the
// critical point 0 (after a rebase), where X is itself tiny: then y is
// formed in the 2^s-scaled domain and f' ~ p·y^(p-1) keeps its exponent.
fn deep_fprime(xf: Fe, w: vec2<f32>, s: i32, yt: vec2<f32>) -> Fe {
if MORPH || max(abs(yt.x), abs(yt.y)) >= DEEP_TINY {
return Fe(fprime(yt), 0);
}
let yf = fe_make(ldexp2_sat(xf.m, xf.e - s) + w, s);
if yf.m.x == 0.0 && yf.m.y == 0.0 {
return Fe(fprime(vec2<f32>(0.0, 0.0)), 0);
}
if KIND == KIND_MULTIBROT {
let p = clamp(u.power, 2u, 8u);
var ym = yf.m; // m^(p-1)
for (var k: u32 = 2u; k < p; k = k + 1u) {
ym = cmul(ym, yf.m);
}
return fe_make(f32(p) * ym, yf.e * i32(p - 1u));
} else if KIND == KIND_LAMBDA {
return Fe(fprime(vec2<f32>(0.0, 0.0)), 0); // λ(1 - 2y) ~ λ
} else if KIND == KIND_COMPLEX_MULTIBROT {
// p·y^(p-1) = p·exp(l), l = (p-1)·ln y; keep exp(Re l)'s exponent.
let p = u.complex_power;
let l = cmul(p - vec2<f32>(1.0, 0.0), clog_fe(yf.m, yf.e));
let k = i32(floor(l.x / LN2));
return fe_make(cmul(p, cexp(vec2<f32>(l.x - f32(k) * LN2, l.y))), k);
}
// z^2-like kinds: 2y (|f'| = |2y| for the abs variants too, see fprime).
return Fe(2.0 * yf.m, yf.e);
}
// The deep phase hands over to the f32 loop once the delta's magnitude
// reaches 2^DEEP_EXIT_LOG2: by then |e|^2 is still a normal f32 and the
// pixel offset dc (< 2^-99 on deep views) is below f32 rounding of e. The DE
// derivative only has to be a comfortably normal f32 (DEEP_EXIT_DZ_LOG2).
const DEEP_EXIT_LOG2: i32 = -48;
const DEEP_EXIT_DZ_LOG2: i32 = -100;
// Below this, a full value y goes through `deep_fprime`'s extended path.
const DEEP_TINY: f32 = 7.888609e-31; // 2^-100
// The mantissas are renormalized once their exponent drifts past ±this.
const DEEP_RENORM_LOG2: i32 = 16;
// A reference point can only matter for rebasing when it's within this many
// binades above the delta's scale (|w| < 2^DEEP_RENORM_LOG2).
const DEEP_NEAR_LOG2: i32 = 24;
// Weight of the Phoenix kind's p*z_{n-1} term in the current map: 1 for plain
// Phoenix, its morph share while switching to/from Phoenix, else 0.
fn phoenix_weight() -> f32 {
@@ -331,13 +667,14 @@ struct Sample {
// Perturbation iterate a single sample. `offset` is the per-pixel offset in
// complex units. For Mandelbrot it is the c-plane offset added every step (delta
// starts at 0); for Julia it is the z-plane offset that seeds the initial delta
// (c is fixed, so nothing is added per step).
// (c is fixed, so nothing is added per step). In `DEEP` pipelines both
// `offset` and `px` are in units of 2^u.scale_exp.
fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// Loop invariants, read once instead of on every iteration.
let max_iter = u.max_iter;
let bailout_sq = u.bailout_sq;
let ref_len = u.ref_len;
let z0 = ref_orbit[0]; // reference start (0 for Mandelbrot, center for Julia)
let z0 = ref_at(0u); // reference start (0 for Mandelbrot, center for Julia)
// Main cardioid / period-2 bulb bypass: those points never escape, so skip
// iterating them (they'd otherwise all burn the full max_iter). `offset` is
@@ -345,7 +682,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// orbit itself, since X_1 = X_0^2 + C_ref = C_ref. That's only f32-accurate,
// so skip the test once a pixel is smaller than that error (deep zoom),
// where it could misclassify pixels right at the boundary.
if KIND == KIND_MANDELBROT && !MORPH && !IS_JULIA && ref_len > 1u && px > 1e-6 {
if KIND == KIND_MANDELBROT && !MORPH && !IS_JULIA && !DEEP && ref_len > 1u && px > 1e-6 {
let c = ref_orbit[1] + offset;
let xq = c.x - 0.25;
let q = xq * xq + c.y * c.y;
@@ -361,6 +698,8 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// seeds the delta and nothing is added per step.
var step_add = offset;
var e = vec2<f32>(0.0, 0.0);
// Pixel size in complex units (`px` is pre-scaled in DEEP pipelines).
var px_t = px;
// Orbit derivative for distance estimation, pre-multiplied by the pixel
// size `px`. For the set plane it is px·d/dc (starts at 0, gains +px each
// step); for Julia it is px·d/dz0 (starts at px). The raw derivative grows
@@ -371,7 +710,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
if IS_JULIA {
step_add = vec2<f32>(0.0, 0.0);
e = offset;
dzs = vec2<f32>(px, 0.0);
dzs = vec2<f32>(px_t, 0.0);
}
// Previous-iterate state for the Phoenix two-term recurrence (delta of
// y_{n-1}, and its scaled derivative for DE). Both start at 0 (y_{-1} = 0).
@@ -403,8 +742,247 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
var cand_d2 = 3.0e38;
var mult2_cand = 1.0;
// Deep phase (DEEP pipelines only): iterate the delta as e = w * 2^sx,
// with `w` an f32 mantissa kept near 1 by renormalizing and `sx` an i32
// exponent, while |e| is too small for f32 (it starts at the pixel offset,
// ~2^scale_exp). The DE derivative is linear in the same way and carried
// as dzs = v * 2^sv, with its own exponent: near the critical point (after
// a rebase) f' is tiny, so dzs and e can drift far apart. Once both are
// big enough (or a rebase makes the delta large), the state is converted
// to plain f32 and the loop below carries on from the same n / m.
if DEEP {
let scale_e = u.scale_exp;
var sx = scale_e;
var sv = scale_e;
var sc = ldexp_sat(1.0, sx);
var w = vec2<f32>(0.0, 0.0);
var v = vec2<f32>(0.0, 0.0);
var w_prev = vec2<f32>(0.0, 0.0);
var v_prev = vec2<f32>(0.0, 0.0);
// Per-step additions (dc and px for the set plane), in units of 2^sx
// and 2^sv.
var d = offset;
var pd = px;
if IS_JULIA {
w = offset;
v = vec2<f32>(px, 0.0);
d = vec2<f32>(0.0, 0.0);
pd = 0.0;
}
let z0f = ref_fe(0u);
var xf = z0f;
var xt = z0;
var xf_old = xf;
var xt_old = xt;
var w_old = w;
var rebase_exit = false;
var cut_exit = false;
loop {
// Full value y = X + e as f32 (e flushes to 0 when negligible).
let yt = xt + ldexp2_sat(w, sx);
if dot(yt, yt) > bailout_sq {
escaped = true;
break;
}
if n >= max_iter {
break;
}
if deep_cut_crossing(xf, w, sx) {
cut_exit = true;
break;
}
if DE {
let fp = deep_fprime(xf, w, sx, yt);
if fp.e == 0 {
var v_new = cmul(fp.m, v);
if !IS_JULIA {
v_new.x = v_new.x + pd;
}
if phoenix_w > 0.0 {
v_new = v_new + phoenix_w * cmul(u.phoenix_p, v_prev);
v_prev = v;
}
v = v_new;
} else {
// f' is below f32's range (next to the critical point):
// sum the terms at the largest one's scale and move
// there, like the delta below.
var acc = fe_make(cmul(fp.m, v), fp.e + sv);
if !IS_JULIA {
acc = fe_add(acc, fe_make(vec2<f32>(px, 0.0), scale_e));
}
if phoenix_w > 0.0 {
acc = fe_add(acc, fe_make(phoenix_w * cmul(u.phoenix_p, v_prev), sv));
}
if acc.m.x == 0.0 && acc.m.y == 0.0 {
if phoenix_w > 0.0 {
v_prev = v;
}
v = acc.m;
} else {
if phoenix_w > 0.0 {
v_prev = ldexp2_sat(v, sv - acc.e);
}
v = acc.m;
sv = acc.e;
if !IS_JULIA {
pd = ldexp_sat(px, scale_e - sv);
}
}
}
}
w_old = w;
let st = advance_delta_scaled(xt, xf, w, sc, sx);
if st.t == sx {
w = st.w + d;
if phoenix_w > 0.0 {
w = w + phoenix_w * cmul(u.phoenix_p, w_prev);
w_prev = w_old;
}
} else {
// The step's value is at another scale (next to the
// critical point it is ~e^2, far below 2^sx): add dc and the
// Phoenix term at the largest addend's scale and move there.
var acc = fe_make(st.w, st.t);
if !IS_JULIA {
acc = fe_add(acc, fe_make(offset, scale_e));
}
if phoenix_w > 0.0 {
acc = fe_add(acc, fe_make(phoenix_w * cmul(u.phoenix_p, w_prev), sx));
}
if acc.m.x == 0.0 && acc.m.y == 0.0 {
w = acc.m;
if phoenix_w > 0.0 {
w_prev = w_old;
}
} else {
// w_old (the pre-step delta) is still needed at the new
// scale: it's the next previous delta, and a rebase reads it.
w_old = ldexp2_sat(w_old, sx - acc.e);
if phoenix_w > 0.0 {
w_prev = w_old;
}
w = acc.m;
sx = acc.e;
sc = ldexp_sat(1.0, sx);
if !IS_JULIA {
d = ldexp2_sat(offset, scale_e - sx);
}
}
}
m = m + 1u;
n = n + 1u;
if m >= ref_len {
escaped = true; // see the f32 loop's reference-exhausted case
break;
}
xf_old = xf;
xt_old = xt;
xf = ref_fe(m);
xt = ldexp2_sat(xf.m, xf.e);
// Rebase test |X + e| < |e|, in units of 2^sx. Only possible when
// |X| is within a few binades of |e| (|w| < 2^DEEP_RENORM_LOG2).
if xf.e - sx < DEEP_NEAR_LOG2 {
let q = ldexp2_sat(xf.m, xf.e - sx) + w;
if dot(q, q) < dot(w, w) {
// The new delta y - X[0] stays tiny only if X[0] is
// (always, for the set plane), and for Phoenix only if the
// previous full value y_{n-1} (its new previous delta) is.
let z0_small = (z0f.m.x == 0.0 && z0f.m.y == 0.0)
|| z0f.e - sx < DEEP_NEAR_LOG2;
let prev_small = phoenix_w == 0.0 || xf_old.e - sx < DEEP_NEAR_LOG2;
if !(z0_small && prev_small) {
rebase_exit = true;
break;
}
if phoenix_w > 0.0 {
w_prev = ldexp2_sat(xf_old.m, xf_old.e - sx) + w_old;
}
w = q - ldexp2_sat(z0f.m, z0f.e - sx);
xf = z0f;
xt = z0;
m = 0u;
}
}
// Leave once both the delta and the derivative fit in f32 (an
// exactly-zero one always does); otherwise renormalize the
// mantissas when they drift (exact: powers of two only).
let kw = deep_log2(w, w_prev);
let kv = deep_log2(v, v_prev);
let w_ok = kw == DEEP_ZERO || sx + kw > DEEP_EXIT_LOG2;
let v_ok = !DE || kv == DEEP_ZERO || sv + kv > DEEP_EXIT_DZ_LOG2;
if w_ok && v_ok {
break;
}
if kw != DEEP_ZERO && abs(kw) > DEEP_RENORM_LOG2 {
w = ldexp2_sat(w, -kw);
w_prev = ldexp2_sat(w_prev, -kw);
sx = sx + kw;
sc = ldexp_sat(1.0, sx);
if !IS_JULIA {
d = ldexp2_sat(offset, scale_e - sx);
}
}
if DE && kv != DEEP_ZERO && abs(kv) > DEEP_RENORM_LOG2 {
v = ldexp2_sat(v, -kv);
v_prev = ldexp2_sat(v_prev, -kv);
sv = sv + kv;
if !IS_JULIA {
pd = ldexp_sat(px, scale_e - sv);
}
}
}
// Hand over to the f32 loop. dc and px may flush to 0 here: they
// are below f32 rounding of the (now large enough) delta and
// derivative.
px_t = ldexp_sat(px, scale_e);
if !IS_JULIA {
step_add = ldexp2_sat(offset, scale_e);
}
e = ldexp2_sat(w, sx);
e_prev = ldexp2_sat(w_prev, sx);
dzs = ldexp2_sat(v, sv);
dzs_prev = ldexp2_sat(v_prev, sv);
xm = xt;
if rebase_exit {
// Rebase in plain f32: the new delta (y - X[0], and for Phoenix
// the previous full value) is no longer tiny.
if phoenix_w > 0.0 {
e_prev = xt_old + ldexp2_sat(w_old, sx);
}
e = (xt + e) - z0;
xm = z0;
m = 0u;
}
if cut_exit && e.y == 0.0 && w.y != 0.0 {
// Keep the side of the cut the pixel is on even if e.y flushed
// to 0: that decides the branch in the next (f32) step.
e.y = select(-1.17549435e-38, 1.17549435e-38, w.y > 0.0);
}
z = xm + e;
z2 = dot(z, z);
// Restart periodicity detection from here (check_at must stay ahead
// of n, or no window would ever close). Nothing is saved until the
// first window closes: `z` here is at an arbitrary phase, and the
// pixel still shadows the reference (exactly periodic when it's a
// minibrot nucleus), so comparing against it would flag exterior
// pixels as interior (see PERIOD_FIRST_CHECK). The sentinel is far
// outside the bailout radius, so no return can match it.
z_saved = vec2<f32>(1e18, 1e18);
z_cand = z;
while check_at <= n {
check_at = check_at * 2u;
}
}
loop {
if z2 > bailout_sq {
// (`escaped` may already be set by the deep phase.)
if escaped || z2 > bailout_sq {
escaped = true;
break;
}
@@ -428,7 +1006,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
if DE {
var dzs_new = cmul(fp, dzs);
if !IS_JULIA {
dzs_new.x = dzs_new.x + px;
dzs_new.x = dzs_new.x + px_t;
}
if phoenix_w > 0.0 {
dzs_new = dzs_new + phoenix_w * cmul(u.phoenix_p, dzs_prev);
@@ -457,7 +1035,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
escaped = true;
break;
}
xm = ref_orbit[m];
xm = ref_at(m);
z = xm + e;
z2 = dot(z, z);
if z2 < dot(e, e) {