fp64 Burning Ship

// ═══ @xgis/shader-dsl example — fp64 Burning Ship ═══
//
// The Burning Ship (z ← (|Re z| + i·|Im z|)² + c) is the fp64 family's
// showcase for `abs` ON THE f64 TYPE: the fold happens inside the
// extended-precision iteration, not after a narrowing — df64 abs is exact
// (negate both planes), so the fold costs no precision. The camera sits on
// the set's spike at c = −1.748 on the real axis (the axis itself never
// escapes — same real dynamics as the Mandelbrot needle), where the CPU
// check keeps 40–110 distinct escape bands per 48² window from a 1e-5 span
// all the way down to 1e-12. The plain-f32 left half collapses flat once the
// span drops under one ulp of 1.748 (~1e-7).
//
// The `_fp64` guard uniform is auto-injected by the lowering; the render
// harnesses bind it to 1.0f by probing the program for the Fp64Guard block.

import {
  fn,
  module,
  vec3,
  vec4,
  f32,
  f64,
  abs,
  pow,
  log2,
  max,
  mix,
  step,
  toF32,
  toF64,
  f32T,
  vec2fT,
  vec2f64T,
  If,
  Loop,
  Var,
  Let,
  u32,
  uniformStruct,
} from '../src/index.ts'
import { VsOut, vs } from './_fullscreen.ts'
import type { ShaderExample } from './_shared.ts'

// A spike point with escape structure at every depth (CPU-verified); y = 0
// keeps the never-escaping real axis mid-frame at every zoom.
const CENTER_X = -1.748
const CENTER_Y = 0
const ITER = 128

const U = uniformStruct(
  'Uniforms',
  { group: 0, binding: 0, as: 'u' },
  {
    center: vec2f64T, // one DF64Vec2 slot — host packs [hi.x, hi.y, lo.x, lo.y]
    resolution: vec2fT,
    zoom_exp: f32T, // view span = 10^-zoom_exp complex units
    fp64: f32T, // toggle: 1 = split-screen f32 | f64 (canonical), 0 = all-f32
  },
)

const fsShip = fn(
  'fs_ship',
  { vo: VsOut },
  (p) => {
    const span = Let(pow(f32(10.0), U.field.zoom_exp.neg()))
    const half = Let(p.vo.uv.x.mul(2.0))
    const sx = Let(half.sub(p.vo.uv.x.lt(0.5).select(0.0, 1.0)))
    const dx = Let(sx.sub(0.5).mul(span))
    const dy = Let(
      p.vo.uv.y.sub(0.5).mul(span).mul(U.field.resolution.y.div(U.field.resolution.x).mul(2.0)),
    )

    const it = Var(f32(0))
    const m2 = Var(f32(0)) // |z|² at escape (frozen once the guard fails)
    If(p.vo.uv.x.lt(0.5).or(U.field.fp64.lt(0.5)), () => {
      // f32 twin — SAME fold-and-square, center narrowed.
      const cx = Let(toF32(U.field.center.x).add(dx))
      const cy = Let(toF32(U.field.center.y).add(dy))
      const zx = Var(f32(0))
      const zy = Var(f32(0))
      Loop(
        u32(0),
        (j) => j.lt(u32(ITER)),
        () => {
          If(zx.mul(zx).add(zy.mul(zy)).le(16.0), () => {
            const nzx = Let(zx.mul(zx).sub(zy.mul(zy)).add(cx))
            zy.assign(abs(zx.mul(zy)).mul(2.0).add(cy))
            zx.assign(nzx)
            it.assign(it.add(1.0))
          })
        },
      )
      m2.assign(zx.mul(zx).add(zy.mul(zy)))
    }).else(() => {
      // f64 — |Re|, |Im| fold in extended precision (2|zx·zy| ≡ the |Im| fold
      // of the squared form: (|zx|+i|zy|)² has Im = 2|zx||zy| = 2|zx·zy|).
      const cx = Let(U.field.center.x.add(toF64(dx)))
      const cy = Let(U.field.center.y.add(toF64(dy)))
      const zx = Var(f64(0))
      const zy = Var(f64(0))
      Loop(
        u32(0),
        (j) => j.lt(u32(ITER)),
        () => {
          If(zx.mul(zx).add(zy.mul(zy)).le(16.0), () => {
            const nzx = Let(zx.mul(zx).sub(zy.mul(zy)).add(cx))
            zy.assign(abs(zx.mul(zy)).mul(2.0).add(cy))
            zx.assign(nzx)
            it.assign(it.add(1.0))
          })
        },
      )
      m2.assign(toF32(zx.mul(zx).add(zy.mul(zy))))
    })

    // Smooth escape time through an ember palette (dark hull → orange flame →
    // pale smoke); interior stays black. Same log₂ log₂ smoothing as
    // fp64-mandelbrot.ts.
    const sn = Let(it.sub(log2(max(log2(max(m2, 1.0001)), 0.0001))).add(1.0))
    const inside = Let(step(f32(ITER).sub(0.5), it))
    const s = Let(sn.div(ITER))
    const ease = Let(s.mul(s).mul(f32(3).sub(s.mul(2)))) // s²(3−2s) ember ramp
    const rgb = mix(
      mix(vec3(0.06, 0.02, 0.05), vec3(0.95, 0.45, 0.08), ease),
      vec3(1.0, 0.93, 0.75),
      s.mul(s),
    ).mul(f32(1).sub(inside))
    return vec4(rgb, f32(1))
  },
  { stage: 'fragment', retAttr: '@location(0)' },
)

// `_fp64` guard lands at (group 0, binding 1) automatically.
const fp64BurningShipModule = module({
  funcs: [vs, fsShip],
  uses: [U, VsOut],
})

export const fp64BurningShip: ShaderExample = {
  id: 'fp64-burning-ship',
  title: 'fp64 Burning Ship',
  blurb:
    'The Burning Ship fractal (fold |Re z|, |Im z| before squaring) zoomed onto its needle spike — the fold runs as df64 abs INSIDE the extended-precision iteration, exact by construction. The plain-f32 left half collapses flat past a ~1e-7 span; the emulated-double right half keeps the flame filaments down to the df64 floor. Drag to pan, wheel to zoom, flip the fp64 toggle to collapse the right half in place.',
  category: 'generic',
  file: 'fp64-burning-ship.ts',
  module: fp64BurningShipModule,
  renderable: true,
  splitLabels: ['f32', 'f64 (emulated)'],
  controls: {
    center: {
      kind: 'pan2d',
      value: [CENTER_X, CENTER_Y],
      zoomExpField: 'zoom_exp',
      unitsPerWidth: 2,
    },
    resolution: { kind: 'resolution' },
    zoom_exp: {
      kind: 'slider',
      label: 'Zoom 10^-x',
      min: 0,
      max: 16,
      step: 0.05,
      value: 4,
      wheel: true,
    },
    fp64: { kind: 'toggle', label: 'fp64 emulation', value: true },
  },
}

struct Uniforms {
  center: DF64Vec2,
  resolution: vec2<f32>,
  zoom_exp: f32,
  fp64: f32,
}

struct VsOut {
  @builtin(position) pos: vec4<f32>,
  @location(0) uv: vec2<f32>,
}

struct DF64Vec2 {
  hi: vec2<f32>,
  lo: vec2<f32>,
}

@group(0) @binding(0) var<uniform> u: Uniforms;
@group(0) @binding(1) var _fp64: texture_2d<f32>;

@vertex
fn vs(@builtin(vertex_index) vi: u32) -> VsOut {
  let _cse0 = ((f32((vi & 1u)) * 4.0) - 1.0);
  let _cse1 = ((f32((vi >> 1u)) * 4.0) - 1.0);
  return VsOut(vec4<f32>(_cse0, _cse1, 0.0, 1.0), vec2<f32>(((_cse0 * 0.5) + 0.5), ((_cse1 * 0.5) + 0.5)));
}

@fragment
fn fs_ship(vo: VsOut) -> @location(0) vec4<f32> {
  let _licm0 = vec2<f32>(16.0, 0.0);
  let _licm1 = vec2<f32>(2.0, 0.0);
  let _cse0 = (vo.uv.x < 0.5);
  let _cse1 = vec2<f32>(u.center.hi.x, u.center.lo.x);
  let _cse2 = vec2<f32>(u.center.hi.y, u.center.lo.y);
  let _cse3 = vec2<f32>(0.0, 0.0);
  let _v0 = pow(10.0, (-u.zoom_exp));
  let _v1 = (vo.uv.x * 2.0);
  let _v2 = (_v1 - select(1.0, 0.0, _cse0));
  let _v3 = ((_v2 - 0.5) * _v0);
  let _v4 = (((vo.uv.y - 0.5) * _v0) * ((u.resolution.y / u.resolution.x) * 2.0));
  var _v5: f32 = 0.0;
  var _v6: f32 = 0.0;
  if ((_cse0 || (u.fp64 < 0.5))) {
    let _v7 = (df64_narrow(_cse1) + _v3);
    let _v8 = (df64_narrow(_cse2) + _v4);
    var _v9: f32 = 0.0;
    var _v10: f32 = 0.0;
    for (var _v11: u32 = 0u; (_v11 < 128u); _v11 = (_v11 + 1u)) {
      if ((((_v9 * _v9) + (_v10 * _v10)) <= 16.0)) {
        let _v12 = (((_v9 * _v9) - (_v10 * _v10)) + _v7);
        _v10 = ((abs((_v9 * _v10)) * 2.0) + _v8);
        _v9 = _v12;
        _v5 = (_v5 + 1.0);
      }
    }
    _v6 = ((_v9 * _v9) + (_v10 * _v10));
  } else {
    let _v13 = df64_add(_cse1, vec2<f32>(_v3, 0.0));
    let _v14 = df64_add(_cse2, vec2<f32>(_v4, 0.0));
    var _v15: vec2<f32> = _cse3;
    var _v16: vec2<f32> = _cse3;
    for (var _v17: u32 = 0u; (_v17 < 128u); _v17 = (_v17 + 1u)) {
      if (df64_le(df64_add(df64_mul(_v15, _v15), df64_mul(_v16, _v16)), _licm0)) {
        let _v18 = df64_add(df64_sub(df64_mul(_v15, _v15), df64_mul(_v16, _v16)), _v13);
        _v16 = df64_add(df64_mul(df64_abs(df64_mul(_v15, _v16)), _licm1), _v14);
        _v15 = _v18;
        _v5 = (_v5 + 1.0);
      }
    }
    _v6 = df64_narrow(df64_add(df64_mul(_v15, _v15), df64_mul(_v16, _v16)));
  }
  let _v19 = ((_v5 - log2(max(log2(max(_v6, 1.0001)), 0.0001))) + 1.0);
  let _v20 = step(127.5, _v5);
  let _v21 = (_v19 / 128.0);
  let _v22 = ((_v21 * _v21) * (3.0 - (_v21 * 2.0)));
  return vec4<f32>((mix(mix(vec3<f32>(0.06, 0.02, 0.05), vec3<f32>(0.95, 0.45, 0.08), _v22), vec3<f32>(1.0, 0.93, 0.75), (_v21 * _v21)) * (1.0 - _v20)), 1.0);
}

fn df64_twoSum(a: f32, b: f32) -> vec2<f32> {
  let _cse0 = textureLoad(_fp64, vec2<i32>(0, 0), 0).x;
  let _v0 = (a + b);
  let _v1 = (((_v0 * _cse0) - a) * _cse0);
  let _v2 = (((((a - ((_v0 - _v1) * _cse0)) * _cse0) * _cse0) * _cse0) + (b - _v1));
  return vec2<f32>(_v0, _v2);
}

fn df64_quickTwoSum(a: f32, b: f32) -> vec2<f32> {
  let _cse0 = textureLoad(_fp64, vec2<i32>(0, 0), 0).x;
  let _v0 = ((a + b) * _cse0);
  let _v1 = (b - ((_v0 - a) * _cse0));
  return vec2<f32>(_v0, _v1);
}

fn df64_split(a: f32) -> vec2<f32> {
  let _cse0 = textureLoad(_fp64, vec2<i32>(0, 0), 0).x;
  let _v0 = (a * (_cse0 * 4097.0));
  let _v1 = ((_v0 * _cse0) - (_v0 - a));
  let _v2 = ((a * _cse0) - _v1);
  return vec2<f32>(_v1, _v2);
}

fn df64_twoProd(a: f32, b: f32) -> vec2<f32> {
  let _v0 = (a * b);
  let _v1 = df64_split(a);
  let _v2 = df64_split(b);
  let _v3 = (((((_v1.x * _v2.x) - _v0) + (_v1.x * _v2.y)) + (_v1.y * _v2.x)) + (_v1.y * _v2.y));
  return vec2<f32>(_v0, _v3);
}

fn df64_add(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
  var _v0: vec2<f32> = df64_twoSum(a.x, b.x);
  let _v1 = df64_twoSum(a.y, b.y);
  _v0.y = (_v0.y + _v1.x);
  _v0 = df64_quickTwoSum(_v0.x, _v0.y);
  _v0.y = (_v0.y + _v1.y);
  _v0 = df64_quickTwoSum(_v0.x, _v0.y);
  return _v0;
}

fn df64_sub(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
  return df64_add(a, (-b));
}

fn df64_mul(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
  var _v0: vec2<f32> = df64_twoProd(a.x, b.x);
  _v0.y = (_v0.y + (a.x * b.y));
  _v0 = df64_quickTwoSum(_v0.x, _v0.y);
  _v0.y = (_v0.y + (a.y * b.x));
  return df64_quickTwoSum(_v0.x, _v0.y);
}

fn df64_le(a: vec2<f32>, b: vec2<f32>) -> bool {
  return ((a.x < b.x) || ((a.x == b.x) && (a.y <= b.y)));
}

fn df64_abs(a: vec2<f32>) -> vec2<f32> {
  return select(a, (-a), (a.x < 0.0));
}

fn df64_narrow(a: vec2<f32>) -> f32 {
  return (a.x + a.y);
}

#version 300 es
precision highp float;
precision highp int;

struct VsOut {
  vec4 pos;
  vec2 uv;
};
out vec2 uv;

VsOut vs_impl(uint vi) {
  float _cse0 = ((float((vi & 1u)) * 4.0) - 1.0);
  float _cse1 = ((float((vi >> 1u)) * 4.0) - 1.0);
  return VsOut(vec4(_cse0, _cse1, 0.0, 1.0), vec2(((_cse0 * 0.5) + 0.5), ((_cse1 * 0.5) + 0.5)));
}

void main() {
  VsOut _out = vs_impl(uint(gl_VertexID));
  gl_Position = _out.pos;
  uv = _out.uv;
}

#version 300 es
precision highp float;
precision highp int;

struct VsOut {
  vec4 pos;
  vec2 uv;
};

struct DF64Vec2 {
  vec2 hi;
  vec2 lo;
};
layout(std140) uniform Uniforms {
  DF64Vec2 center;
  vec2 resolution;
  float zoom_exp;
  float fp64;
} u;

uniform sampler2D _fp64;
vec2 df64_twoSum(float a, float b);
vec2 df64_quickTwoSum(float a, float b);
vec2 df64_split(float a);
vec2 df64_twoProd(float a, float b);
vec2 df64_add(vec2 a, vec2 b);
vec2 df64_sub(vec2 a, vec2 b);
vec2 df64_mul(vec2 a, vec2 b);
bool df64_le(vec2 a, vec2 b);
vec2 df64_abs(vec2 a);
float df64_narrow(vec2 a);
vec2 df64_twoSum(float a, float b) {
  float _cse0 = texelFetch(_fp64, ivec2(0, 0), 0).x;
  float _v0 = (a + b);
  float _v1 = (((_v0 * _cse0) - a) * _cse0);
  float _v2 = (((((a - ((_v0 - _v1) * _cse0)) * _cse0) * _cse0) * _cse0) + (b - _v1));
  return vec2(_v0, _v2);
}

vec2 df64_quickTwoSum(float a, float b) {
  float _cse0 = texelFetch(_fp64, ivec2(0, 0), 0).x;
  float _v0 = ((a + b) * _cse0);
  float _v1 = (b - ((_v0 - a) * _cse0));
  return vec2(_v0, _v1);
}

vec2 df64_split(float a) {
  float _cse0 = texelFetch(_fp64, ivec2(0, 0), 0).x;
  float _v0 = (a * (_cse0 * 4097.0));
  float _v1 = ((_v0 * _cse0) - (_v0 - a));
  float _v2 = ((a * _cse0) - _v1);
  return vec2(_v1, _v2);
}

vec2 df64_twoProd(float a, float b) {
  float _v0 = (a * b);
  vec2 _v1 = df64_split(a);
  vec2 _v2 = df64_split(b);
  float _v3 = (((((_v1.x * _v2.x) - _v0) + (_v1.x * _v2.y)) + (_v1.y * _v2.x)) + (_v1.y * _v2.y));
  return vec2(_v0, _v3);
}

vec2 df64_add(vec2 a, vec2 b) {
  vec2 _v0 = df64_twoSum(a.x, b.x);
  vec2 _v1 = df64_twoSum(a.y, b.y);
  _v0.y = (_v0.y + _v1.x);
  _v0 = df64_quickTwoSum(_v0.x, _v0.y);
  _v0.y = (_v0.y + _v1.y);
  _v0 = df64_quickTwoSum(_v0.x, _v0.y);
  return _v0;
}

vec2 df64_sub(vec2 a, vec2 b) {
  return df64_add(a, (-b));
}

vec2 df64_mul(vec2 a, vec2 b) {
  vec2 _v0 = df64_twoProd(a.x, b.x);
  _v0.y = (_v0.y + (a.x * b.y));
  _v0 = df64_quickTwoSum(_v0.x, _v0.y);
  _v0.y = (_v0.y + (a.y * b.x));
  return df64_quickTwoSum(_v0.x, _v0.y);
}

bool df64_le(vec2 a, vec2 b) {
  return ((a.x < b.x) || ((a.x == b.x) && (a.y <= b.y)));
}

vec2 df64_abs(vec2 a) {
  return ((a.x < 0.0) ? (-a) : a);
}

float df64_narrow(vec2 a) {
  return (a.x + a.y);
}
in vec2 uv;
layout(location = 0) out vec4 _ret;

vec4 fs_ship_impl(VsOut vo) {
  vec2 _licm0 = vec2(16.0, 0.0);
  vec2 _licm1 = vec2(2.0, 0.0);
  bool _cse0 = (vo.uv.x < 0.5);
  vec2 _cse1 = vec2(u.center.hi.x, u.center.lo.x);
  vec2 _cse2 = vec2(u.center.hi.y, u.center.lo.y);
  vec2 _cse3 = vec2(0.0, 0.0);
  float _v0 = pow(10.0, (-u.zoom_exp));
  float _v1 = (vo.uv.x * 2.0);
  float _v2 = (_v1 - (_cse0 ? 0.0 : 1.0));
  float _v3 = ((_v2 - 0.5) * _v0);
  float _v4 = (((vo.uv.y - 0.5) * _v0) * ((u.resolution.y / u.resolution.x) * 2.0));
  float _v5 = 0.0;
  float _v6 = 0.0;
  if ((_cse0 || (u.fp64 < 0.5))) {
    float _v7 = (df64_narrow(_cse1) + _v3);
    float _v8 = (df64_narrow(_cse2) + _v4);
    float _v9 = 0.0;
    float _v10 = 0.0;
    for (uint _v11 = 0u; (_v11 < 128u); _v11 = (_v11 + 1u)) {
      if ((((_v9 * _v9) + (_v10 * _v10)) <= 16.0)) {
        float _v12 = (((_v9 * _v9) - (_v10 * _v10)) + _v7);
        _v10 = ((abs((_v9 * _v10)) * 2.0) + _v8);
        _v9 = _v12;
        _v5 = (_v5 + 1.0);
      }
    }
    _v6 = ((_v9 * _v9) + (_v10 * _v10));
  } else {
    vec2 _v13 = df64_add(_cse1, vec2(_v3, 0.0));
    vec2 _v14 = df64_add(_cse2, vec2(_v4, 0.0));
    vec2 _v15 = _cse3;
    vec2 _v16 = _cse3;
    for (uint _v17 = 0u; (_v17 < 128u); _v17 = (_v17 + 1u)) {
      if (df64_le(df64_add(df64_mul(_v15, _v15), df64_mul(_v16, _v16)), _licm0)) {
        vec2 _v18 = df64_add(df64_sub(df64_mul(_v15, _v15), df64_mul(_v16, _v16)), _v13);
        _v16 = df64_add(df64_mul(df64_abs(df64_mul(_v15, _v16)), _licm1), _v14);
        _v15 = _v18;
        _v5 = (_v5 + 1.0);
      }
    }
    _v6 = df64_narrow(df64_add(df64_mul(_v15, _v15), df64_mul(_v16, _v16)));
  }
  float _v19 = ((_v5 - log2(max(log2(max(_v6, 1.0001)), 0.0001))) + 1.0);
  float _v20 = step(127.5, _v5);
  float _v21 = (_v19 / 128.0);
  float _v22 = ((_v21 * _v21) * (3.0 - (_v21 * 2.0)));
  return vec4((mix(mix(vec3(0.06, 0.02, 0.05), vec3(0.95, 0.45, 0.08), _v22), vec3(1.0, 0.93, 0.75), (_v21 * _v21)) * (1.0 - _v20)), 1.0);
}

void main() {
  VsOut vo;
  vo.pos = gl_FragCoord;
  vo.uv = uv;
  _ret = fs_ship_impl(vo);
}

{
  "bindGroups": [
    {
      "group": 0,
      "entries": [
        {
          "group": 0,
          "binding": 0,
          "name": "u",
          "space": "uniform",
          "resourceKind": "uniform-buffer",
          "structName": "Uniforms"
        }
      ]
    }
  ],
  "uniforms": [
    {
      "name": "Uniforms",
      "size": 32,
      "align": 16,
      "fields": [
        {
          "name": "center",
          "type": "vec2<f64>",
          "offset": 0,
          "align": 16,
          "size": 16
        },
        {
          "name": "resolution",
          "type": "vec2<f32>",
          "offset": 16,
          "align": 8,
          "size": 8
        },
        {
          "name": "zoom_exp",
          "type": "f32",
          "offset": 24,
          "align": 4,
          "size": 4
        },
        {
          "name": "fp64",
          "type": "f32",
          "offset": 28,
          "align": 4,
          "size": 4
        }
      ]
    }
  ],
  "storage": [],
  "entries": [
    {
      "name": "vs",
      "stage": "vertex",
      "inputs": [
        "u32"
      ],
      "output": "struct:VsOut"
    },
    {
      "name": "fs_ship",
      "stage": "fragment",
      "inputs": [
        "struct:VsOut"
      ],
      "output": "vec4<f32>"
    }
  ],
  "overrides": []
}