// luma.gl
// SPDX-License-Identifier: MIT
// Copyright (c) vis.gl contributors
import { glsl } from "../../../lib/glsl-utils/highlight.js";
export const fp64arithmeticShader = `\
uniform float ONE;
vec2 split(float a) {
const float SPLIT = 4097.0;
float t = a * SPLIT;
#if defined(LUMA_FP64_CODE_ELIMINATION_WORKAROUND)
float a_hi = t * ONE - (t - a);
float a_lo = a * ONE - a_hi;
#else
float a_hi = t - (t - a);
float a_lo = a - a_hi;
#endif
return vec2(a_hi, a_lo);
}
vec2 split2(vec2 a) {
vec2 b = split(a.x);
b.y += a.y;
return b;
}
vec2 quickTwoSum(float a, float b) {
#if defined(LUMA_FP64_CODE_ELIMINATION_WORKAROUND)
float sum = (a + b) * ONE;
float err = b - (sum - a) * ONE;
#else
float sum = a + b;
float err = b - (sum - a);
#endif
return vec2(sum, err);
}
vec2 twoSum(float a, float b) {
float s = (a + b);
#if defined(LUMA_FP64_CODE_ELIMINATION_WORKAROUND)
float v = (s * ONE - a) * ONE;
float err = (a - (s - v) * ONE) * ONE * ONE * ONE + (b - v);
#else
float v = s - a;
float err = (a - (s - v)) + (b - v);
#endif
return vec2(s, err);
}
vec2 twoSub(float a, float b) {
float s = (a - b);
#if defined(LUMA_FP64_CODE_ELIMINATION_WORKAROUND)
float v = (s * ONE - a) * ONE;
float err = (a - (s - v) * ONE) * ONE * ONE * ONE - (b + v);
#else
float v = s - a;
float err = (a - (s - v)) - (b + v);
#endif
return vec2(s, err);
}
vec2 twoSqr(float a) {
float prod = a * a;
vec2 a_fp64 = split(a);
#if defined(LUMA_FP64_CODE_ELIMINATION_WORKAROUND)
float err = ((a_fp64.x * a_fp64.x - prod) * ONE + 2.0 * a_fp64.x *
a_fp64.y * ONE * ONE) + a_fp64.y * a_fp64.y * ONE * ONE * ONE;
#else
float err = ((a_fp64.x * a_fp64.x - prod) + 2.0 * a_fp64.x * a_fp64.y) + a_fp64.y * a_fp64.y;
#endif
return vec2(prod, err);
}
vec2 twoProd(float a, float b) {
float prod = a * b;
vec2 a_fp64 = split(a);
vec2 b_fp64 = split(b);
float err = ((a_fp64.x * b_fp64.x - prod) + a_fp64.x * b_fp64.y +
a_fp64.y * b_fp64.x) + a_fp64.y * b_fp64.y;
return vec2(prod, err);
}
vec2 sum_fp64(vec2 a, vec2 b) {
vec2 s, t;
s = twoSum(a.x, b.x);
t = twoSum(a.y, b.y);
s.y += t.x;
s = quickTwoSum(s.x, s.y);
s.y += t.y;
s = quickTwoSum(s.x, s.y);
return s;
}
vec2 sub_fp64(vec2 a, vec2 b) {
vec2 s, t;
s = twoSub(a.x, b.x);
t = twoSub(a.y, b.y);
s.y += t.x;
s = quickTwoSum(s.x, s.y);
s.y += t.y;
s = quickTwoSum(s.x, s.y);
return s;
}
vec2 mul_fp64(vec2 a, vec2 b) {
vec2 prod = twoProd(a.x, b.x);
prod.y += a.x * b.y;
#if defined(LUMA_FP64_HIGH_BITS_OVERFLOW_WORKAROUND)
prod = split2(prod);
#endif
prod = quickTwoSum(prod.x, prod.y);
prod.y += a.y * b.x;
#if defined(LUMA_FP64_HIGH_BITS_OVERFLOW_WORKAROUND)
prod = split2(prod);
#endif
prod = quickTwoSum(prod.x, prod.y);
return prod;
}
vec2 div_fp64(vec2 a, vec2 b) {
float xn = 1.0 / b.x;
#if defined(LUMA_FP64_HIGH_BITS_OVERFLOW_WORKAROUND)
vec2 yn = mul_fp64(a, vec2(xn, 0));
#else
vec2 yn = a * xn;
#endif
float diff = (sub_fp64(a, mul_fp64(b, yn))).x;
vec2 prod = twoProd(xn, diff);
return sum_fp64(yn, prod);
}
vec2 sqrt_fp64(vec2 a) {
if (a.x == 0.0 && a.y == 0.0) return vec2(0.0, 0.0);
if (a.x < 0.0) return vec2(0.0 / 0.0, 0.0 / 0.0);
float x = 1.0 / sqrt(a.x);
float yn = a.x * x;
#if defined(LUMA_FP64_CODE_ELIMINATION_WORKAROUND)
vec2 yn_sqr = twoSqr(yn) * ONE;
#else
vec2 yn_sqr = twoSqr(yn);
#endif
float diff = sub_fp64(a, yn_sqr).x;
vec2 prod = twoProd(x * 0.5, diff);
#if defined(LUMA_FP64_HIGH_BITS_OVERFLOW_WORKAROUND)
return sum_fp64(split(yn), prod);
#else
return sum_fp64(vec2(yn, 0.0), prod);
#endif
}
`;