fp64 relative-to-center

// ═══ typeshade example — fp64 relative-to-center rendering ═══
//
// THE technique every planet-scale engine ships (RTC / relative-to-eye,
// Cesium & co.): a survey marker sits at world (10⁸+3.7, 5×10⁷+2.3), and the
// only value the shader ever needs is delta = marker − eye — a SMALL number
// that f32 carries perfectly, PROVIDED the subtraction happens in extended
// precision FIRST. The f64 right half subtracts in df64 and narrows the
// result: a crisp reticle. The f32 left half narrows first and subtracts
// after — both operands quantize to ulp(10⁸) = 8 world units, the delta
// comes out in 8-unit steps, and the reticle snaps to a coarse grid, sits
// visibly OFF the true position, and squares off (drag to feel it jump).
// This is the same subtract-then-narrow discipline as the map package's
// extended-precision tile origins — demoted to a single reticle.
//
// 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,
  vec2,
  vec3,
  vec4,
  f32,
  pow,
  fract,
  abs,
  min,
  max,
  mix,
  length,
  smoothstep,
  exp,
  toF32,
  toF64,
  f32T,
  vec2fT,
  vec2f64T,
  Let,
  uniformStruct,
} from '../src/index.js';
import { VsOut, vs } from './_fullscreen.js';
import type { ShaderExample } from './_shared.js';

const U = uniformStruct(
  'Uniforms',
  { group: 0, binding: 0, as: 'u' },
  {
    center: vec2f64T, // camera/eye — one DF64Vec2 slot [hi.x, hi.y, lo.x, lo.y]
    mark: vec2f64T, // marker world position (eye + a few fractional units)
    resolution: vec2fT,
    zoom_exp: f32T, // view span = 10^-zoom_exp world units (negative = zoom out)
    fp64: f32T, // toggle: 1 = split-screen f32 | f64 (canonical), 0 = all-f32
  },
);

const fsRtc = fn(
  'fs_rtc',
  { 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 isF32 = Let(p.vo.uv.x.lt(0.5).or(U.field.fp64.lt(0.5)));
    // f64 path — subtract FIRST (df64, exact by Sterbenz-style cancellation),
    // narrow the small result after. The pixel's world position is
    // eye + offset; delta = (eye + offset) − marker.
    const ex64 = Let(toF32(U.field.center.x.add(toF64(dx)).sub(U.field.mark.x)));
    const ey64 = Let(toF32(U.field.center.y.add(toF64(dy)).sub(U.field.mark.y)));
    // f32 twin — narrow FIRST, subtract after: at 10⁸ both operands live on the
    // 8-unit ulp grid, so the delta quantizes to 8-unit steps.
    const ex32 = Let(toF32(U.field.center.x).add(dx).sub(toF32(U.field.mark.x)));
    const ey32 = Let(toF32(U.field.center.y).add(dy).sub(toF32(U.field.mark.y)));
    const ex = Let(isF32.select(ex32, ex64));
    const ey = Let(isF32.select(ey32, ey64));

    // Radar-style reticle in world units, everything scaled by the span so the
    // picture is zoom-invariant: rings every span/8, a crosshair, a hot dot.
    const rw = Let(span.mul(0.125)); // ring spacing
    const r = Let(length(vec2(ex, ey)));
    const pixw = Let(span.div(U.field.resolution.x.mul(0.5)));
    const tri = Let(
      abs(fract(r.div(rw)).sub(0.5))
        .neg()
        .add(0.5)
        .mul(rw),
    ); // dist to ring
    const ring = Let(f32(1).sub(smoothstep(f32(0), pixw.mul(1.6).add(1e-9), tri)));
    const cross = Let(
      f32(1).sub(smoothstep(f32(0), pixw.mul(1.4).add(1e-9), min(abs(ex), abs(ey)))),
    );
    const dotGlow = Let(exp(r.div(pixw.mul(6.0).add(1e-9)).neg()));
    const vignette = Let(max(f32(0), f32(1).sub(r.div(span.mul(0.75)))));

    const bg = Let(mix(vec3(0.01, 0.04, 0.02), vec3(0.02, 0.09, 0.045), vignette));
    const rgb = bg
      .add(vec3(0.1, 0.75, 0.3).mul(ring.mul(0.8)))
      .add(vec3(0.12, 0.9, 0.4).mul(cross.mul(0.55)))
      .add(vec3(1.0, 0.45, 0.25).mul(dotGlow));
    return vec4(rgb, f32(1));
  },
  { stage: 'fragment', retAttr: '@location(0)' },
);

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

export const fp64Rtc: ShaderExample = {
  id: 'fp64-rtc',
  title: 'fp64 relative-to-center',
  blurb:
    'Why planet-scale engines render relative-to-eye: a reticle marks a survey point a few fractional units from the eye. Drag the DISTANCE slider out from the origin — near 10⁶ both halves put the reticle dead centre, but past ~10⁷·² the eye and marker land on the same f32 ulp grid, their difference quantizes, and the f32 left half snaps the reticle off-target in whole-ulp jumps. The f64 right half subtracts eye from marker in extended precision FIRST, then narrows the small delta — crisp to 10⁹. Wheel drives the distance; flip the fp64 toggle to compare.',
  category: 'cartographic',
  file: 'fp64-rtc.ts',
  module: fp64RtcModule,
  renderable: true,
  splitLabels: ['f32', 'f64 (emulated)'],
  controls: {
    // Eye and marker sweep together: eye = 10^mag, marker = eye + (3.7, 2.3).
    // At mag = 8 this is the classic (10⁸, 5·10⁷) view; the fractional (3.7, 2.3)
    // is exactly what f32 loses once one ulp grows past it.
    center: { kind: 'logmag2d', magField: 'mag', base: [1, 0.5], offset: [0, 0] },
    mark: { kind: 'logmag2d', magField: 'mag', base: [1, 0.5], offset: [3.7, 2.3] },
    resolution: { kind: 'resolution' },
    mag: {
      kind: 'slider',
      label: 'Distance from origin 10^x',
      min: 4,
      max: 9,
      step: 0.02,
      value: 6,
      wheel: true,
    },
    // View span across each half (world units).
    zoom_exp: { kind: 'slider', label: 'Zoom 10^-x', min: -2, max: 5, step: 0.05, value: -1.4 },
    fp64: { kind: 'toggle', label: 'fp64 emulation', value: true },
  },
};

struct Uniforms {
  @align(16) center: DF64Vec2,
  @align(16) mark: 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_rtc(vo: VsOut) -> @location(0) vec4<f32> {
  let _fp64_g = textureLoad(_fp64, vec2<i32>(0, 0), 0).x;
  let _v0 = pow(10.0, (-u.zoom_exp));
  let _v1 = (vo.uv.x * 2.0);
  let _cse0 = (vo.uv.x < 0.5);
  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));
  let _v5 = (_cse0 || (u.fp64 < 0.5));
  let _cse1 = vec2<f32>(u.center.hi.x, u.center.lo.x);
  let _cse2 = vec2<f32>(u.mark.hi.x, u.mark.lo.x);
  let _cse6 = bitcast<f32>(bitcast<u32>(0.0));
  let _cse3 = vec2<f32>(_cse6, _cse6);
  let _v6 = df64_narrow(df64_sub(df64_add(_cse1, vec2<f32>(_v3, 0.0), _fp64_g), df64_add(_cse2, _cse3, _fp64_g), _fp64_g));
  let _cse4 = vec2<f32>(u.center.hi.y, u.center.lo.y);
  let _cse5 = vec2<f32>(u.mark.hi.y, u.mark.lo.y);
  let _v7 = df64_narrow(df64_sub(df64_add(_cse4, vec2<f32>(_v4, 0.0), _fp64_g), df64_add(_cse5, _cse3, _fp64_g), _fp64_g));
  let _v8 = ((df64_narrow(_cse1) + _v3) - df64_narrow(_cse2));
  let _v9 = ((df64_narrow(_cse4) + _v4) - df64_narrow(_cse5));
  let _v10 = select(_v6, _v8, _v5);
  let _v11 = select(_v7, _v9, _v5);
  let _v12 = (_v0 * 0.125);
  let _v13 = length(vec2<f32>(_v10, _v11));
  let _v14 = (_v0 / (u.resolution.x * 0.5));
  let _v15 = (((-abs((fract((_v13 / _v12)) - 0.5))) + 0.5) * _v12);
  let _v16 = (1.0 - smoothstep(0.0, ((_v14 * 1.6) + 1e-9), _v15));
  let _v17 = (1.0 - smoothstep(0.0, ((_v14 * 1.4) + 1e-9), min(abs(_v10), abs(_v11))));
  let _v18 = exp((-(_v13 / ((_v14 * 6.0) + 1e-9))));
  let _v19 = max(0.0, (1.0 - (_v13 / (_v0 * 0.75))));
  let _v20 = mix(vec3<f32>(0.01, 0.04, 0.02), vec3<f32>(0.02, 0.09, 0.045), _v19);
  return vec4<f32>((((_v20 + (vec3<f32>(0.1, 0.75, 0.3) * (_v16 * 0.8))) + (vec3<f32>(0.12, 0.9, 0.4) * (_v17 * 0.55))) + (vec3<f32>(1.0, 0.45, 0.25) * _v18)), 1.0);
}

fn df64_twoSum(a: f32, b: f32, _fp64_g: f32) -> vec2<f32> {
  let _v0 = (a + b);
  let _v1 = (((_v0 * _fp64_g) - a) * _fp64_g);
  let _v2 = (((a - ((_v0 - _v1) * _fp64_g)) * _fp64_g) + (b - _v1));
  return vec2<f32>(_v0, _v2);
}

fn df64_quickTwoSum(a: f32, b: f32, _fp64_g: f32) -> vec2<f32> {
  let _v0 = ((a + b) * _fp64_g);
  let _v1 = (b - ((_v0 - a) * _fp64_g));
  return vec2<f32>(_v0, _v1);
}

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

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

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

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

out vec2 uv;

void main() {
  uint vi = uint(gl_VertexID);
  float _cse0 = ((float((vi & 1u)) * 4.0) - 1.0);
  float _cse1 = ((float((vi >> 1u)) * 4.0) - 1.0);
  gl_Position = vec4(_cse0, _cse1, 0.0, 1.0);
  uv = vec2(((_cse0 * 0.5) + 0.5), ((_cse1 * 0.5) + 0.5));
}

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

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

uniform highp sampler2D _fp64;
vec2 df64_twoSum(float a, float b, float _fp64_g) {
  float _v0 = (a + b);
  float _v1 = (((_v0 * _fp64_g) - a) * _fp64_g);
  float _v2 = (((a - ((_v0 - _v1) * _fp64_g)) * _fp64_g) + (b - _v1));
  return vec2(_v0, _v2);
}

vec2 df64_quickTwoSum(float a, float b, float _fp64_g) {
  float _v0 = ((a + b) * _fp64_g);
  float _v1 = (b - ((_v0 - a) * _fp64_g));
  return vec2(_v0, _v1);
}

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

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

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

void main() {
  float _fp64_g = texelFetch(_fp64, ivec2(0, 0), 0).x;
  float _v0 = pow(10.0, (-u.zoom_exp));
  float _v1 = (uv.x * 2.0);
  bool _cse0 = (uv.x < 0.5);
  float _v2 = (_v1 - (_cse0 ? 0.0 : 1.0));
  float _v3 = ((_v2 - 0.5) * _v0);
  float _v4 = (((uv.y - 0.5) * _v0) * ((u.resolution.y / u.resolution.x) * 2.0));
  bool _v5 = (_cse0 || (u.fp64 < 0.5));
  vec2 _cse1 = vec2(u.center.hi.x, u.center.lo.x);
  vec2 _cse2 = vec2(u.mark.hi.x, u.mark.lo.x);
  float _cse6 = uintBitsToFloat(floatBitsToUint(0.0));
  vec2 _cse3 = vec2(_cse6, _cse6);
  float _v6 = df64_narrow(df64_sub(df64_add(_cse1, vec2(_v3, 0.0), _fp64_g), df64_add(_cse2, _cse3, _fp64_g), _fp64_g));
  vec2 _cse4 = vec2(u.center.hi.y, u.center.lo.y);
  vec2 _cse5 = vec2(u.mark.hi.y, u.mark.lo.y);
  float _v7 = df64_narrow(df64_sub(df64_add(_cse4, vec2(_v4, 0.0), _fp64_g), df64_add(_cse5, _cse3, _fp64_g), _fp64_g));
  float _v8 = ((df64_narrow(_cse1) + _v3) - df64_narrow(_cse2));
  float _v9 = ((df64_narrow(_cse4) + _v4) - df64_narrow(_cse5));
  float _v10 = (_v5 ? _v8 : _v6);
  float _v11 = (_v5 ? _v9 : _v7);
  float _v12 = (_v0 * 0.125);
  float _v13 = length(vec2(_v10, _v11));
  float _v14 = (_v0 / (u.resolution.x * 0.5));
  float _v15 = (((-abs((fract((_v13 / _v12)) - 0.5))) + 0.5) * _v12);
  float _v16 = (1.0 - smoothstep(0.0, ((_v14 * 1.6) + 1e-9), _v15));
  float _v17 = (1.0 - smoothstep(0.0, ((_v14 * 1.4) + 1e-9), min(abs(_v10), abs(_v11))));
  float _v18 = exp((-(_v13 / ((_v14 * 6.0) + 1e-9))));
  float _v19 = max(0.0, (1.0 - (_v13 / (_v0 * 0.75))));
  vec3 _v20 = mix(vec3(0.01, 0.04, 0.02), vec3(0.02, 0.09, 0.045), _v19);
  _ret = vec4((((_v20 + (vec3(0.1, 0.75, 0.3) * (_v16 * 0.8))) + (vec3(0.12, 0.9, 0.4) * (_v17 * 0.55))) + (vec3(1.0, 0.45, 0.25) * _v18)), 1.0);
}