fp64 Mercator tiles

// ═══ typeshade example — fp64 Web-Mercator tile pyramid ═══
//
// The actual tile-engine computation, at the zoom levels where f32 dies: a
// normalized Web-Mercator coordinate (Seoul, x ≈ 0.85272) is multiplied by
// 2^z and split with floor/fract into (tile index, in-tile uv). One f32 ulp
// of 0.85 is ~6e-8 — at tile-z 21 that is 27 PIXELS of a 640-px half, so the
// plain-f32 left half renders the city as blocky smears while the f64 right
// half keeps crisp tile borders through z 23. The scale 2^z is built by
// EXACT repeated doubling (a power of two is exact in f32; pow(2, z) need not
// be), so the multiply itself never eats the df64 budget. Tile indices at
// z ≤ 23 fit f32 integers (< 2^24) — only the FRACTION ever needed f64.
//
// 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,
  floor,
  fract,
  sin,
  dot,
  min,
  mix,
  step,
  smoothstep,
  toF32,
  toF64,
  f32T,
  u32,
  vec2fT,
  vec2f64T,
  Loop,
  Var,
  Let,
  uniformStruct,
  splitF64,
} from '../src/index.js';
import { VsOut, vs } from './_fullscreen.js';
import type { ShaderExample } from './_shared.js';

// Seoul in normalized Web-Mercator ([0,1)² — x = (lon+180)/360).
const SEOUL_X = 0.8527166666666667;
const SEOUL_Y = 0.3872522862452674;
const TILES_ACROSS = 3; // view width in tiles, per half

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,
    tile_z: f32T, // integer tile-pyramid level (slider, step 1)
    fp64: f32T, // toggle: 1 = split-screen f32 | f64 (canonical), 0 = all-f32
  },
);

const fsTiles = fn(
  'fs_tiles',
  { vo: VsOut },
  (p) => {
    // 2^z by exact repeated doubling — every step is a power of two, exact in
    // f32 up to 2^127; pow(2.0, z) is implementation-defined precision.
    const scale = Var(f32(1));
    Loop(
      u32(0),
      (j) => toF32(j).lt(U.field.tile_z),
      () => {
        scale.assign(scale.mul(2.0));
      },
    );
    const span = Let(f32(TILES_ACROSS).div(scale)); // world units per half-width

    const halfUv = Let(p.vo.uv.x.mul(2.0));
    const sx = Let(halfUv.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 — the multiply by an exact power of two costs nothing; floor
    // and fract on the f64 type split tile index from in-tile position.
    const tx = Let(U.field.center.x.add(toF64(dx)).mul(scale));
    const ty = Let(U.field.center.y.add(toF64(dy)).mul(scale));
    const tix64 = Let(toF32(floor(tx))); // < 2^24 at z ≤ 23 — exact as f32
    const tiy64 = Let(toF32(floor(ty)));
    const lx64 = Let(toF32(fract(tx)));
    const ly64 = Let(toF32(fract(ty)));
    // Checker parity in f64: at z 23 the SUM of the two indices can pass 2^24,
    // so fold to 0/1 before narrowing.
    const par64 = Let(toF32(fract(floor(tx).add(floor(ty)).mul(0.5))).mul(2.0));
    // f32 twin — SAME formulas, coordinate narrowed first: the in-tile
    // fraction moves in ulp(0.85)·2^z jumps (27 px at z 21).
    const tx32 = Let(toF32(U.field.center.x).add(dx).mul(scale));
    const ty32 = Let(toF32(U.field.center.y).add(dy).mul(scale));
    const tix32 = Let(floor(tx32));
    const tiy32 = Let(floor(ty32));
    const lx32 = Let(fract(tx32));
    const ly32 = Let(fract(ty32));
    const par32 = Let(fract(tix32.add(tiy32).mul(0.5)).mul(2.0));

    const tix = Let(isF32.select(tix32, tix64));
    const tiy = Let(isF32.select(tiy32, tiy64));
    const lx = Let(isF32.select(lx32, lx64));
    const ly = Let(isF32.select(ly32, ly64));
    const par = Let(isF32.select(par32, par64));

    // Map-ish styling: a per-tile pastel wash (hash of the tile index), the
    // checker parity as a second tone, crisp tile borders, and an in-tile
    // gradient so the interior visibly flows (or, on the left, visibly jumps).
    const h = Let(fract(sin(dot(vec2(tix, tiy), vec2(127.1, 311.7))).mul(43758.5453)));
    const tint = Let(
      mix(vec3(0.78, 0.84, 0.88), vec3(0.9, 0.87, 0.78), h).mul(mix(f32(0.88), f32(1.0), par)),
    );
    const grad = Let(mix(f32(0.92), f32(1.05), lx.mul(0.6).add(ly.mul(0.4))));
    const edge = Let(min(min(lx, f32(1).sub(lx)), min(ly, f32(1).sub(ly))));
    const pixw = Let(f32(TILES_ACROSS).div(U.field.resolution.x.mul(0.5))); // tile units / px
    const border = Let(smoothstep(f32(0), pixw.mul(1.6).add(1e-9), edge));
    const inkline = Let(step(edge, pixw.mul(6.0))); // wider label-side rule
    const rgb = tint
      .mul(grad)
      .mul(mix(f32(0.45), f32(1.0), border))
      .sub(vec3(0.05, 0.04, 0.02).mul(inkline));
    return vec4(rgb, f32(1));
  },
  { stage: 'fragment', retAttr: '@location(0)' },
);

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

// DF64Vec2 std140 buffer order is PLANE-major: [hi.x, hi.y, lo.x, lo.y].
const [HX, LX] = splitF64(SEOUL_X);
const [HY, LY] = splitF64(SEOUL_Y);

export const fp64MercatorTiles: ShaderExample = {
  id: 'fp64-mercator-tiles',
  title: 'fp64 Mercator tiles',
  blurb:
    'The tile-engine formula itself: normalized Web-Mercator (Seoul) × 2^z, floor → tile index, fract → in-tile uv. One f32 ulp of x ≈ 0.85 is 27 px at tile-z 21, so the plain-f32 left half smears tiles into blocks while the f64 right half keeps crisp borders through z 23. The 2^z scale is built by exact repeated doubling — the multiply is free; only the fraction ever needed extended precision. Wheel through zoom levels and flip the fp64 toggle to compare.',
  category: 'cartographic',
  file: 'fp64-mercator-tiles.ts',
  module: fp64MercatorTilesModule,
  renderable: true,
  splitLabels: ['f32', 'f64 (emulated)'],
  controls: {
    center: { kind: 'const', value: [HX, HY, LX, LY] },
    resolution: { kind: 'resolution' },
    tile_z: {
      kind: 'slider',
      label: 'Tile zoom z',
      min: 12,
      max: 23,
      step: 1,
      value: 16,
      wheel: true,
    },
    fp64: { kind: 'toggle', label: 'fp64 emulation', value: true },
  },
};

struct Uniforms {
  @align(16) center: DF64Vec2,
  resolution: vec2<f32>,
  tile_z: 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_tiles(vo: VsOut) -> @location(0) vec4<f32> {
  let _fp64_g = textureLoad(_fp64, vec2<i32>(0, 0), 0).x;
  let _licm0 = u.tile_z;
  var _v0: f32 = 1.0;
  for (var _v1: u32 = 0u; (f32(_v1) < _licm0); _v1 = (_v1 + 1u)) {
    _v0 = (_v0 * 2.0);
  }
  let _v2 = (3.0 / _v0);
  let _v3 = (vo.uv.x * 2.0);
  let _cse0 = (vo.uv.x < 0.5);
  let _v4 = (_v3 - select(1.0, 0.0, _cse0));
  let _v5 = ((_v4 - 0.5) * _v2);
  let _v6 = (((vo.uv.y - 0.5) * _v2) * ((u.resolution.y / u.resolution.x) * 2.0));
  let _v7 = (_cse0 || (u.fp64 < 0.5));
  let _cse1 = vec2<f32>(u.center.hi.x, u.center.lo.x);
  let _gv0 = vec2<f32>(_v0, 0.0);
  let _v8 = df64_mul(df64_add(_cse1, vec2<f32>(_v5, 0.0), _fp64_g), _gv0, _fp64_g);
  let _cse2 = vec2<f32>(u.center.hi.y, u.center.lo.y);
  let _v9 = df64_mul(df64_add(_cse2, vec2<f32>(_v6, 0.0), _fp64_g), _gv0, _fp64_g);
  let _gv1 = df64_floor(_v8, _fp64_g);
  let _v10 = df64_narrow(_gv1);
  let _gv2 = df64_floor(_v9, _fp64_g);
  let _v11 = df64_narrow(_gv2);
  let _cse4 = bitcast<f32>(bitcast<u32>(0.0));
  let _cse3 = vec2<f32>(_cse4, _cse4);
  let _v12 = df64_narrow(df64_fract(df64_add(_v8, _cse3, _fp64_g), _fp64_g));
  let _v13 = df64_narrow(df64_fract(df64_add(_v9, _cse3, _fp64_g), _fp64_g));
  let _v14 = (df64_narrow(df64_fract((df64_add(_gv1, _gv2, _fp64_g) * 0.5), _fp64_g)) * 2.0);
  let _v15 = ((df64_narrow(_cse1) + _v5) * _v0);
  let _v16 = ((df64_narrow(_cse2) + _v6) * _v0);
  let _v17 = floor(_v15);
  let _v18 = floor(_v16);
  let _v19 = fract(_v15);
  let _v20 = fract(_v16);
  let _v21 = (fract(((_v17 + _v18) * 0.5)) * 2.0);
  let _v22 = select(_v10, _v17, _v7);
  let _v23 = select(_v11, _v18, _v7);
  let _v24 = select(_v12, _v19, _v7);
  let _v25 = select(_v13, _v20, _v7);
  let _v26 = select(_v14, _v21, _v7);
  let _v27 = fract((sin(dot(vec2<f32>(_v22, _v23), vec2<f32>(127.1, 311.7))) * 43758.5453));
  let _v28 = (mix(vec3<f32>(0.78, 0.84, 0.88), vec3<f32>(0.9, 0.87, 0.78), _v27) * mix(0.88, 1.0, _v26));
  let _v29 = mix(0.92, 1.05, ((_v24 * 0.6) + (_v25 * 0.4)));
  let _v30 = min(min(_v24, (1.0 - _v24)), min(_v25, (1.0 - _v25)));
  let _v31 = (3.0 / (u.resolution.x * 0.5));
  let _v32 = smoothstep(0.0, ((_v31 * 1.6) + 1e-9), _v30);
  let _v33 = step(_v30, (_v31 * 6.0));
  return vec4<f32>((((_v28 * _v29) * mix(0.45, 1.0, _v32)) - (vec3<f32>(0.05, 0.04, 0.02) * _v33)), 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_split(a: f32, _fp64_g: f32) -> vec2<f32> {
  let _v0 = (a * (_fp64_g * 4097.0));
  let _v1 = ((_v0 * _fp64_g) - (_v0 - a));
  let _v2 = ((a * _fp64_g) - _v1);
  return vec2<f32>(_v1, _v2);
}

fn df64_twoProd(a: f32, b: f32, _fp64_g: f32) -> vec2<f32> {
  let _v0 = (a * b);
  let _v1 = df64_split(a, _fp64_g);
  let _v2 = df64_split(b, _fp64_g);
  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>, _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_mul(a: vec2<f32>, b: vec2<f32>, _fp64_g: f32) -> vec2<f32> {
  var _v0: vec2<f32> = df64_twoProd(a.x, b.x, _fp64_g);
  _v0.y = (_v0.y + (a.x * b.y));
  _v0 = df64_quickTwoSum(_v0.x, _v0.y, _fp64_g);
  _v0.y = (_v0.y + (a.y * b.x));
  return df64_quickTwoSum(_v0.x, _v0.y, _fp64_g);
}

fn df64_floor(a: vec2<f32>, _fp64_g: f32) -> vec2<f32> {
  let _v0 = floor(a.x);
  return select(vec2<f32>(_v0, 0.0), df64_quickTwoSum(_v0, floor(a.y), _fp64_g), (_v0 == a.x));
}

fn df64_fract(a: vec2<f32>, _fp64_g: f32) -> vec2<f32> {
  return df64_sub(a, df64_floor(a, _fp64_g), _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;
  vec2 resolution;
  float tile_z;
  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_split(float a, float _fp64_g) {
  float _v0 = (a * (_fp64_g * 4097.0));
  float _v1 = ((_v0 * _fp64_g) - (_v0 - a));
  float _v2 = ((a * _fp64_g) - _v1);
  return vec2(_v1, _v2);
}

vec2 df64_twoProd(float a, float b, float _fp64_g) {
  float _v0 = (a * b);
  vec2 _v1 = df64_split(a, _fp64_g);
  vec2 _v2 = df64_split(b, _fp64_g);
  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, 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);
}

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

vec2 df64_floor(vec2 a, float _fp64_g) {
  float _v0 = floor(a.x);
  return ((_v0 == a.x) ? df64_quickTwoSum(_v0, floor(a.y), _fp64_g) : vec2(_v0, 0.0));
}

vec2 df64_fract(vec2 a, float _fp64_g) {
  return df64_sub(a, df64_floor(a, _fp64_g), _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 _licm0 = u.tile_z;
  float _v0 = 1.0;
  for (uint _v1 = 0u; (float(_v1) < _licm0); _v1 = (_v1 + 1u)) {
    _v0 = (_v0 * 2.0);
  }
  float _v2 = (3.0 / _v0);
  float _v3 = (uv.x * 2.0);
  bool _cse0 = (uv.x < 0.5);
  float _v4 = (_v3 - (_cse0 ? 0.0 : 1.0));
  float _v5 = ((_v4 - 0.5) * _v2);
  float _v6 = (((uv.y - 0.5) * _v2) * ((u.resolution.y / u.resolution.x) * 2.0));
  bool _v7 = (_cse0 || (u.fp64 < 0.5));
  vec2 _cse1 = vec2(u.center.hi.x, u.center.lo.x);
  vec2 _gv0 = vec2(_v0, 0.0);
  vec2 _v8 = df64_mul(df64_add(_cse1, vec2(_v5, 0.0), _fp64_g), _gv0, _fp64_g);
  vec2 _cse2 = vec2(u.center.hi.y, u.center.lo.y);
  vec2 _v9 = df64_mul(df64_add(_cse2, vec2(_v6, 0.0), _fp64_g), _gv0, _fp64_g);
  vec2 _gv1 = df64_floor(_v8, _fp64_g);
  float _v10 = df64_narrow(_gv1);
  vec2 _gv2 = df64_floor(_v9, _fp64_g);
  float _v11 = df64_narrow(_gv2);
  float _cse4 = uintBitsToFloat(floatBitsToUint(0.0));
  vec2 _cse3 = vec2(_cse4, _cse4);
  float _v12 = df64_narrow(df64_fract(df64_add(_v8, _cse3, _fp64_g), _fp64_g));
  float _v13 = df64_narrow(df64_fract(df64_add(_v9, _cse3, _fp64_g), _fp64_g));
  float _v14 = (df64_narrow(df64_fract((df64_add(_gv1, _gv2, _fp64_g) * 0.5), _fp64_g)) * 2.0);
  float _v15 = ((df64_narrow(_cse1) + _v5) * _v0);
  float _v16 = ((df64_narrow(_cse2) + _v6) * _v0);
  float _v17 = floor(_v15);
  float _v18 = floor(_v16);
  float _v19 = fract(_v15);
  float _v20 = fract(_v16);
  float _v21 = (fract(((_v17 + _v18) * 0.5)) * 2.0);
  float _v22 = (_v7 ? _v17 : _v10);
  float _v23 = (_v7 ? _v18 : _v11);
  float _v24 = (_v7 ? _v19 : _v12);
  float _v25 = (_v7 ? _v20 : _v13);
  float _v26 = (_v7 ? _v21 : _v14);
  float _v27 = fract((sin(dot(vec2(_v22, _v23), vec2(127.1, 311.7))) * 43758.5453));
  vec3 _v28 = (mix(vec3(0.78, 0.84, 0.88), vec3(0.9, 0.87, 0.78), _v27) * mix(0.88, 1.0, _v26));
  float _v29 = mix(0.92, 1.05, ((_v24 * 0.6) + (_v25 * 0.4)));
  float _v30 = min(min(_v24, (1.0 - _v24)), min(_v25, (1.0 - _v25)));
  float _v31 = (3.0 / (u.resolution.x * 0.5));
  float _v32 = smoothstep(0.0, ((_v31 * 1.6) + 1e-9), _v30);
  float _v33 = step(_v30, (_v31 * 6.0));
  _ret = vec4((((_v28 * _v29) * mix(0.45, 1.0, _v32)) - (vec3(0.05, 0.04, 0.02) * _v33)), 1.0);
}