fp64 (Extended-Precision Shader Math)

luma.gl emulates extended precision by representing one scalar as two f32 values, high + low. This double-single representation preserves small differences between large coordinates and supports arithmetic beyond ordinary f32 precision. It is not native IEEE 754 binary64 arithmetic: it offers up to roughly 48 significant bits, about 14 decimal digits, within essentially the exponent range of f32.

Choose a module

Need	Module or technique	Shader language
Add, subtract, multiply, divide, square root	fp64arithmetic	GLSL and WGSL
Exponential, logarithm, trigonometry, vector and matrix helpers	fp64, which depends on fp64arithmetic	GLSL vertex shaders; WGSL projection subset below
Subtract large binary64 coordinates before converting to a local value	Raw binary64 subtraction helpers in fp64arithmetic	WGSL
Render small offsets around a known origin	Subtract the origin on the CPU, then upload local f32 values	Any

Use the narrowest technique that satisfies the application's error budget. For example, a map can preserve accurate global positions while using ordinary f32 arithmetic after subtracting a nearby origin. Carrying double-single values through the entire shader is often unnecessary.

The fp64 module provides the full GLSL vertex library and the WGSL projection subset described below, including exponential, logarithm, and trigonometric functions. For GLSL fragment arithmetic, use fp64arithmetic, which supplies both vertex and fragment source.

Represent and upload values

A double-single scalar occupies a vec2 in GLSL or vec2f in WGSL:

value = high + low
high  = roundToFloat32(value)
low   = roundToFloat32(value - high)

Split the original JavaScript number before converting it to a Float32Array. Splitting a value that has already been rounded to f32 cannot recover the missing digits.

import {fp64ify, fp64LowPart} from '@luma.gl/shadertools';

const coordinate = 10_000_000.125;
const packedCoordinate = new Float32Array(2);
fp64ify(coordinate, packedCoordinate);
// packedCoordinate contains [10_000_000, 0.125].

const high = Math.fround(coordinate);
const low = fp64LowPart(coordinate);

The CPU helpers are also available as methods on fp64 and fp64arithmetic:

Helper	Result
fp64ify(value, output?, startIndex = 0)	Writes high and low terms at consecutive indices and returns the output array. Allocates an ordinary array when no output is supplied.
fp64LowPart(value)	Returns value - Math.fround(value); uploading it as f32 rounds the residual.
fp64ifyMatrix4(matrix)	Returns a 32-element Float32Array with adjacent high/low terms for each element, transposing the input with output[(row * 4 + column) * 2] split from matrix[column * 4 + row].

A two-dimensional double-single position needs four floats, for example [xHigh, xLow, yHigh, yLow]. A pair is one extended-precision scalar, not a spatial vector. Keep both terms until the final conversion or origin-relative result. Adding them with an ordinary shader + immediately rounds the result back to f32.

Use in a shader

For arithmetic only, register the smaller module:

import {fp64arithmetic} from '@luma.gl/shadertools';

const modules = [fp64arithmetic];

These snippets show shader expressions after the module has been assembled into the shader:

// Both terms must be supplied from the original coordinates.
vec2 position = vec2(10000000.0, 0.125);
vec2 origin = vec2(10000000.0, 0.0);
vec2 offset = sub_fp64(position, origin);
vec2 scaledOffset = mul_fp64(offset, vec2(2.0, 0.0));
float localPosition = scaledOffset.x + scaledOffset.y; // 0.25

let position = vec2f(10000000.0, 0.125);
let origin = vec2f(10000000.0, 0.0);
let offset = sub_fp64(position, origin);
let scaledOffset = mul_fp64(offset, vec2f(2.0, 0.0));
let localPosition = scaledOffset.x + scaledOffset.y;

To use the full GLSL vertex function library, register modules: [fp64]. It includes fp64arithmetic as a dependency. Engine module uniform management supplies the helper defaults; custom compute and direct pipeline integrations must ensure that their uniforms are bound and uploaded before execution. See uniform setup.

Scalar functions

Each arithmetic argument and return value is a high/low pair. The spelling of addition is sum_fp64, not add_fp64.

Function	Operation	Availability
sum_fp64(a, b)	Add two pairs	GLSL and WGSL
sub_fp64(a, b)	Subtract two pairs	GLSL and WGSL
mul_fp64(a, b)	Multiply two pairs	GLSL and WGSL
div_fp64(a, b)	Divide two pairs	GLSL and WGSL
sqrt_fp64(a)	Square root	GLSL and WGSL
exp_fp64(a)	Exponential	GLSL vertex and WGSL
log_fp64(a)	Natural logarithm	GLSL vertex and WGSL
sin_fp64(a), cos_fp64(a), tan_fp64(a)	Trigonometry, with angles in radians	GLSL vertex and WGSL
radians_fp64(degrees)	Convert degrees to radians	GLSL vertex and WGSL
mix_fp64(a, b, fraction)	Interpolate pairs with an ordinary float fraction	GLSL vertex and WGSL

The GLSL library also provides vec2_*_fp64, vec3_*_fp64, and vec4_*_fp64 helpers, including lengths, distances, dot products, and mat4_vec4_mul_fp64. These operate on arrays of high/low pairs; they do not accept ordinary vec3 or mat4 values as interchangeable inputs. See the module source for individual signatures.

For WGSL normalization, classification, comparisons, and raw binary64 subtraction, see the fp64arithmetic reference.

WGSL projection functions

Import fp64 as a shader module to use exp_fp64, log_fp64, sin_fp64, cos_fp64, tan_fp64, radians_fp64, and mix_fp64. Each scalar is a vec2f(high, low). These functions use the selected arithmetic mode of their fp64arithmetic dependency. GLSL-only vector helpers not listed here are not yet ported.

WGSL uses return values instead of GLSL output parameters:

• vec4_fp64(a: vec4f) -> array<vec2f, 4> converts four ordinary floats.
• vec4_dot_fp64(a: array<vec2f, 4>, b: array<vec2f, 4>) -> vec2f computes a dot product.
• mat4_vec4_mul_fp64(matrix: array<vec2f, 16>, a: array<vec2f, 4>) -> array<vec2f, 4> uses the same row-major pair ordering as the GLSL helper: output row i uses matrix entries i * 4 through i * 4 + 3. Transpose column-major input before splitting.

The WGSL transcendental functions target geographic projection, including the Mercator expression log(tan(PI / 4 + latitude / 2)). Trigonometric range reduction is intended for angles within one turn; large-angle accuracy is not guaranteed, and magnitudes above 1,000,000 radians or non-finite inputs return NaN. Tangent is ill-conditioned near its poles. Logarithm accepts positive inputs; zero and negative inputs return NaN, matching the GLSL helper. Positive infinity returns infinity. These functions retain the binary32 exponent range and double-single precision limits described below; they do not provide binary64 semantics.

GPU regression tests compare classic, hybrid, and integer-controlled execution with JavaScript binary64 references, including low-limb differences near one, geographic boundaries, nearby latitudes, and matrix cancellation. Their tolerances are regression criteria, not a universal precision guarantee across inputs or devices. Classic precision tests exclude Apple/Metal, where the module defaults to integer arithmetic because floating-point residuals are not preserved. Boundary checks and compilation still cover classic mode there. The numerical suite requires a hardware adapter, matching the existing arithmetic tests: software adapters validate WGSL module compilation, while native compute-pipeline compilation of the emulated arithmetic can be prohibitively slow.

Precision and range

Property	Native f32	Double-single	IEEE binary64
Significant bits	24	Up to about 48	53
Approximate decimal precision	7 digits	14 digits	15–16 digits
Storage per scalar	4 bytes	8 bytes	8 bytes
Exponent range	Binary32	Essentially binary32	Binary64
Arithmetic	Native shader operations	Emulated operations on pairs	Not available as a portable WebGL/WGSL scalar

These figures describe representation capacity, not a universal error bound for every operation or sequence. Cancellation, intermediate overflow, underflow, compiler transformations, and conditioning of the algorithm all affect accuracy. Extended precision cannot repair an ill-conditioned formula or input that was rounded before upload.

The low term becomes subnormal near underflow: for result magnitudes from about 2^-102 down toward 2^-126, precision tapers from roughly 48 toward 24 bits. Classic floating-point paths can lose subnormal intermediates when hardware flushes them to zero, even if a later result would be representable. A double-single algorithm is not a general implementation of binary64 rounding, NaN propagation, infinities, or signed-zero rules.

Historical classic-WebGL measurements on a 2015 AMD Radeon R9 M370X found errors below one ulp for addition/subtraction, approximately 1.5 ulps for multiplication, 2 for division, 2.6 for square root and exponential, 11.6 for logarithm, and 5 for trigonometry. These are measurements of that implementation and device, not portable guarantees. Logarithm accuracy also depends on the native logarithm used as a starting estimate.

WGSL arithmetic modes and backend behavior

The data representation and function names stay the same across modes:

Mode	Approach	Main tradeoff
Classic	Floating-point transforms recover rounding residuals	Depends on the compiler preserving the required evaluation steps; subnormal intermediates remain hardware-sensitive.
Integer-controlled	Integer operations reconstruct critical arithmetic and rounding	More predictable rounding, but larger shaders and higher compilation/execution cost.
Hybrid	Integer reconstruction of high-part residuals with native low-term accumulation	A throughput/precision compromise; not the fully controlled correctness path.

For the public fp64arithmetic WGSL module, assembly automatically enables integer-controlled arithmetic on Apple WebGPU adapters. Explicit defines can override this default. Higher-level kernels can choose their own policy: do not infer a projection kernel's mode solely from the module default or from the name "precise". See mode selection and the projection API.

Shader compilation is part of the cost. A large shader with repeated calls to integer-heavy helpers can expand substantially when the backend inlines them. Measure pipeline creation separately from steady-state dispatch time. A fast compile does not establish numerical accuracy, and a successful software-adapter test does not establish hardware behavior.

Live example and benchmark

The Mandelbrot panes show where f32 loses detail during a deep zoom. On Apple WebGPU, the live FP64 pane uses hybrid arithmetic; the benchmark separately exposes the integer-controlled path. It measures arithmetic workloads and reports numerical error together with timings. Run it on the devices and browsers used by the application rather than treating one GPU's result as a universal ranking.

Optional interactive GPU benchmarkExplore floating-point precision.Compare Mandelbrot rendering and compute precision when you are ready to use your GPU.

Validation and troubleshooting

• Coordinates still jitter: check that splitting happens before conversion to f32, that high and low terms arrive in the expected order, and that every precision-sensitive operation uses the pair helpers. Prefer a nearby origin when possible.
• Results are invalid or lose low terms: check the module uniform binding and upload, especially in custom compute paths. ONE must be 1.0 and SPLIT must be 4097.0.
• Only one browser/backend fails: compare classic and integer-controlled modes against a CPU reference; compiler transformations and subnormal handling differ. Explicitly record the selected mode, GPU, backend, browser, and driver with results.
• Compilation takes seconds or minutes: time assembly and pipeline creation independently, reduce repeated heavy helper call sites, and evaluate a different arithmetic mode against the same numerical reference before adopting it.
• A value is stored as vec2u: establish whether it contains raw binary64 words or bitcast double-single limbs. Those layouts are not interchangeable.

Validate representative inputs and edge cases: large nearby coordinates, cancellation, different signs, extreme exponents, subnormal low terms, and non-finite values if the application accepts them. Specify whether the reference is exact arithmetic, a sequence of binary64 operations, or an absolute/relative tolerance. Report maximum error as well as typical error.

For numerical derivations and alternatives, see GPU Floating-Point Precision Techniques.