Пастки чисел з рухомою комою в коді симуляцій

1. Бітовий формат IEEE 754

JavaScript використовує подвійну точність (64 біти) для всіх операцій із числами. Типове значення в GLSL WebGL — одинарна точність (32 біти). Обидва відповідають IEEE 754:

Провідна 1. у мантисі є неявною (не зберігається) — це прихований біт. Денормалізовані числа (див. §5) є винятком: експонента = 0 означає, що провідний біт дорівнює 0, а не 1.

2. NaN, нескінченність та −0

IEEE 754 визначає кілька спеціальних значень, які поширюються крізь арифметику несподіваними способами:

// Перевірка на NaN — єдиний правильний спосіб
if (Number.isNaN(x)) { /* ... */ }   // правильно
if (x !== x) { /* ... */ }             // теж правильно (властивість NaN)
if (x === NaN) { /* ... */ }           // НЕПРАВИЛЬНО — завжди false

// Відрізнити -0 від +0
const isNegZero = (x) => x === 0 && (1/x) === -Infinity;

3. Точність і машинний епсилон

Машинний епсилон (ε) — це найменше число, для якого 1 + ε ≠ 1:

Похибки точності накопичуються. Після N додавань очікувана відносна похибка становить O(N·ε). У покроковій за часом симуляції з 60 fps × 3600 с = 216 000 кроків похибки float32 сягають ~2.5% — достатньо, щоб орбіти помітно дрейфували.

4. Катастрофічне скорочення

Коли віднімають два майже рівні великі числа, значущі цифри скорочуються, і результат має набагато меншу точність, ніж вхідні дані:

// Приклад: відстань між частинками, що майже перекриваються
let x1 = 1000000.0001, x2 = 1000000.0002;
let dx = x2 - x1;   // = 0.0001, але у float32 зберігається як ~0.000099182
                     // Похибка: ~0.8% лише від скорочення

Виправлення: Перецентруйте координати навколо цікавої області. Для планетарних симуляцій зберігайте позиції відносно камери або локального початку, а не в абсолютному світовому просторі:

// Зі світового простору у відносний до камери (трюк float-origin)
const camX = camera.position.x;
const relX = particle.x - camX;  // мале число, без скорочення

Та сама проблема виникає у квадратних формулах. Використовуйте числово стабільну форму:

// НЕСТАБІЛЬНО для |b| >> |c/b|
x = (-b + Math.sqrt(b*b - 4*a*c)) / (2*a);

// СТАБІЛЬНО: уникає скорочення, коли b і корінь дискримінанта майже скорочуються
const q = -0.5 * (b + Math.sign(b) * Math.sqrt(b*b - 4*a*c));
x1 = q / a;
x2 = c / q;

// Компенсоване сумування Кагана
function kahanSum(arr) {
  let sum = 0, c = 0;
  for (const x of arr) {
    const y = x - c;
    const t = sum + y;
    c = (t - sum) - y;
    sum = t;
  }
  return sum;
}

5. Денормалізовані числа та антипереповнення

Коли значення float32 падає нижче за ~1.2×10⁻³⁸, апаратне забезпечення активує денормалізований (субнормальний) режим: провідний біт змінюється з 1 на 0, дозволяючи поступове антипереповнення коштом точності. Денормалізована арифметика може бути у 50–100 разів повільнішою на CPU (і часто скидається до нуля на GPU).

// Приклад скидання денормалізованих — внесок маси частинки незначний,
// але все одно запускає вдвічі повільніший денормалізований шлях щокадру
particle.vx *= dampingFactor;   // якщо vx = 1e-39 (float32) → денормалізоване

// Виправлення: обнуліть крихітні швидкості
const MIN_V = 1e-30;
if (Math.abs(particle.vx) < MIN_V) particle.vx = 0;

6. Залежність від порядку та неасоціативність

Додавання з рухомою комою не асоціативне: (a + b) + c ≠ a + (b + c) у загальному випадку. Це означає:

• Суми паралельної редукції дають різні результати залежно від порядку потоків — важливо для відтворюваних повторів симуляцій.
• Зміна прапорців оптимізації компілятора (наприклад, -ffast-math) може дати інші результати через переупорядкування операцій з рухомою комою.
• Рушії JavaScript можуть внутрішньо використовувати 80-бітні регістри x87, даючи результати, які недетерміновано відрізняються від справжнього float64 (рідко, але реально на x86).

7. Порівняння значень з рухомою комою

Ніколи не порівнюйте числа з рухомою комою через ===, якщо ви не знаєте напевно, що обидва значення задані однаковим літералом або отримані з однакового обчислення без проміжного округлення:

// НЕПРАВИЛЬНО: майже ніколи не істинно через округлення
if (a + b === c) { ... }

// Абсолютний епсилон — добре для значень поблизу нуля (напр., швидкостей)
const nearZero = (x, eps = 1e-9) => Math.abs(x) < eps;

// Відносний епсилон — добре для значень, далеких від нуля (напр., енергії)
const nearEqual = (a, b, eps = 1e-9) =>
  Math.abs(a - b) <= eps * Math.max(Math.abs(a), Math.abs(b), 1);

// Порівняння в ULP — найсуворіше (потребує DataView)
function ulpDiff(a, b) {
  const buf = new ArrayBuffer(8);
  const view = new DataView(buf);
  view.setFloat64(0, a); const ia = view.getBigInt64(0);
  view.setFloat64(0, b); const ib = view.getBigInt64(0);
  return Number(ia > ib ? ia - ib : ib - ia);
}

8. Точність GPU (GLSL mediump / highp)

WebGL 1.0 вимагає GLSL mediump (float16 або float32 залежно від драйвера) для фрагментних шейдерів. WebGL 2.0 типово використовує float32 і для вершинних, і для фрагментних. Завжди оголошуйте точність явно:

// На початку кожного шейдера GLSL
precision highp float;   // рекомендовано для симуляцій
precision highp int;

// Уникайте неявного mediump (типове у фрагментних шейдерах WebGL 1)

Поширені пастки точності GPU:

• Позиції у світовому просторі у float32: На рельєфі планетарного масштабу координати гравця у (100 000, 50 000) втрачають ~1 см точності з float32. Використовуйте double на CPU і надсилайте на GPU позиції відносно камери.
• Час як uniform: float uTime втратить субмілісекундну точність після ~16 секунд. Передавайте як два float (старший + молодший) або періодично скидайте до 0.
• Тригонометричні функції: sin(x) у GLSL не визначена для |x| > 2π на деяких драйверах. Обгортайте через mod(x, 6.283185307).
• Цілочисельні атомарні операції у WebGPU: Обчислювальні шейдери WebGPU підтримують атомарні операції i32/u32, але не з рухомою комою. Використовуйте цілочисельні акумулятори, масштабовані множником.