Precisión en GPU frente a CPU: por qué el mismo cálculo da otro resultado

Construyes una matriz en JavaScript, la subes a un uniform y el shader la usa: parece la misma matriz, y no lo es. JavaScript calcula todo en doble precisión y el shader trabaja con la versión truncada a treinta y dos bits, así que el Float32Array que hay en medio no es un contenedor, es una frontera con aduana. Y al otro lado de esa frontera la aritmética ni siquiera es determinista entre fabricantes: la especificación de WGSL le concede libertades al compilador que hacen que dos GPU distintas den resultados distintos, y las dos estén cumpliendo la norma.

f32 es todo lo que hay

En WGSL el tipo de coma flotante es f32. Es el tipo por defecto de cualquier literal decimal, es lo que hay dentro de vec3f y de mat4x4f, y es lo único con lo que puedes contar. No existe f64 en WGSL. No hay doble precisión en la GPU en WebGPU, ni con una feature, ni con una extensión, ni por acuerdo con el fabricante. Si tu algoritmo necesita más de veinticuatro bits de significando, tienes que construirlo tú a partir de varios f32 o reformular el problema para no necesitarlo.

El otro tipo de coma flotante es f16, y solo aparece si el dispositivo anuncia la feature shader-f16 y tú la pides al crear el dispositivo. Da diez bits de mantisa, once de significando, y llega hasta 65504. Es media precisión: sirve para pesos de una red, para acumular color en un rango controlado o para reducir el ancho de banda de un buffer, y no sirve absolutamente para nada relacionado con posiciones.

const adapter = await navigator.gpu.requestAdapter();
const soportaF16 = adapter.features.has('shader-f16');
const device = await adapter.requestDevice({
  requiredFeatures: soportaF16 ? ['shader-f16'] : [],
});

Mientras tanto, en el otro lado, JavaScript no tiene ningún tipo de coma flotante que no sea f64. Un number es un doble de sesenta y cuatro bits con cincuenta y tres de significando. Cuando compones una matriz de modelo-vista-proyección en la CPU, cada multiplicación y cada suma se hacen con veintinueve bits más de los que tendrá el shader. Esa matriz es, a efectos prácticos, exacta.

La conversión ocurre en un punto y solo en uno: cuando escribes en un Float32Array. Ahí cada number se redondea al f32 más cercano, y a partir de ese momento el valor que tenías deja de existir.

const datos = new Float32Array(16);
datos.set(matrizEnDoble);              // aqui se pierde la precision, no antes
device.queue.writeBuffer(uniformBuffer, 0, datos);

Esta distinción tiene una consecuencia operativa enorme y muy poco aprovechada: todo cálculo que puedas mover a la CPU sale gratis en precisión. No es que la CPU sea más rápida, es que su aritmética por defecto tiene mil millones de veces menos error absoluto. Componer matrices, calcular diferencias de posiciones, acumular tiempo, resolver una intersección crítica: si se hace una vez por fotograma y no una vez por vértice, hazlo arriba.

Las libertades del compilador

Aquí es donde la intuición de la CPU deja de valer. En JavaScript, a * b + c significa exactamente una multiplicación redondeada seguida de una suma redondeada, y el resultado es reproducible en cualquier máquina del planeta. En WGSL no.

La especificación permite que el compilador contraiga una multiplicación y una suma en una única operación de multiplicación-suma fusionada, un FMA, que calcula el producto exacto y lo suma a c con un solo redondeo al final. El hardware lo hace en la misma latencia que una multiplicación sola, así que el compilador lo aplica siempre que puede. Y el resultado es distinto.

Toma dos valores exactamente representables en f32:

let a = 1.000244140625;   // 1 + 2^-12
let c = 1.00048828125;    // 1 + 2^-11
let d = a * a - c;        // el compilador decide como evaluarlo

Sin contracción: el producto exacto es 1 + 2^-11 + 2^-24. Redondeado a f32, el término 2^-24 cae exactamente en la mitad del escalón —el ULP en esa octava es 2^-23—, así que hay empate, se redondea al par, y el resultado es 1 + 2^-11. Restarle c da cero.

Con contracción: el producto no se redondea. La resta se hace sobre el valor exacto y da 2^-24, que es representable, así que el resultado es 5,96e-8.

Cero o 5,96e-8, según el driver. Las dos respuestas cumplen la especificación. Y no es un caso de laboratorio: a*b - c*d es el determinante de una matriz de dos por dos, el producto vectorial, el discriminante de una ecuación de segundo grado y el test de orientación de un triángulo, que son precisamente las expresiones donde un cambio de signo cambia el comportamiento del algoritmo, no solo el último bit del color.

La especificación tampoco exige que el compilador preserve el orden de evaluación de una cadena aritmética en todos los casos, así que sumar cuatro términos puede reasociarse. En coma flotante la suma no es asociativa, de modo que (x + y) + (z + w) y ((x + y) + z) + w no dan el mismo valor. Súmalo todo: el mismo shader, con las mismas entradas, puede dar resultados distintos en dos GPU distintas, y las dos son conformes. Deja de esperar bit-exactitud entre dispositivos y empieza a diseñar para que no la necesites.

Si necesitas la fusión de forma explícita y controlada, WGSL tiene fma(a, b, c), que garantiza el redondeo único. Lo que no tiene es la operación inversa: no hay manera portable de prohibirle al compilador que contraiga una expresión que tú preferirías en dos pasos.

@invariant y el depth prepass

Hay un atributo en WGSL para el único caso donde la reproducibilidad es obligatoria: @invariant, que solo puede aplicarse a @builtin(position).

Lo que garantiza es preciso y conviene no exagerarlo: si dos shaders calculan la posición con la misma expresión a partir de las mismas entradas, el compilador está obligado a producir el mismo resultado bit a bit en los dos. No hace que dos fórmulas distintas coincidan. Lo que impide es que el optimizador aplique una contracción a FMA en un pipeline y no en el otro, o que reordene una suma aquí y no allá porque el contexto era distinto.

struct SalidaVS {
  @invariant @builtin(position) posicion: vec4f,
  @location(0) uv: vec2f,
}

@vertex
fn vs(@location(0) local: vec3f, @location(1) uv: vec2f) -> SalidaVS {
  var salida: SalidaVS;
  salida.posicion = camara.proyeccion * (camara.modeloVista * vec4f(local, 1.0));
  salida.uv = uv;
  return salida;
}

Ese atributo es exactamente lo que necesita un depth prepass. El patrón consiste en dibujar toda la geometría opaca en una primera pasada escribiendo solo profundidad, y repetir el dibujado con depthWriteEnabled: false y depthCompare: 'equal', de forma que el fragment shader caro se ejecute una vez por píxel visible y no una vez por fragmento dibujado.

const prepass = device.createRenderPipeline({
  layout: layoutCompartido,
  vertex: { module: moduloProfundidad, entryPoint: 'vs' },
  depthStencil: { format: 'depth32float', depthWriteEnabled: true, depthCompare: 'less' },
});

const color = device.createRenderPipeline({
  layout: layoutCompartido,
  vertex: { module: moduloColor, entryPoint: 'vs' },
  fragment: { module: moduloColor, entryPoint: 'fs', targets: [{ format }] },
  depthStencil: { format: 'depth32float', depthWriteEnabled: false, depthCompare: 'equal' },
});

equal no admite aproximaciones: si el prepass escribió 0,4372918 y la pasada de color calcula 0,43729183, la comparación falla y el píxel no se dibuja. Sin @invariant, el compilador tiene libertad para optimizar los dos módulos de forma distinta —el primero no tiene fragment shader, así que el análisis de vida de variables cambia, y con él las decisiones de contracción—, y el resultado son huecos negros o parpadeo. El fallo es además dependiente del driver: funciona en tu máquina y falla en la del cliente.

Dos condiciones acompañan al atributo y son igual de obligatorias. Los dos vertex shaders tienen que calcular la posición con la misma expresión, literalmente la misma; una versión «simplificada» del prepass que sea matemáticamente equivalente no sirve. Y los uniforms tienen que contener los mismos bits en las dos pasadas, lo cual significa no recomponer la matriz entre ellas.

Interpolación y funciones trascendentes

Entre el vertex shader y el fragment shader hay una etapa que también introduce error y que casi nunca se contabiliza. La interpolación con corrección de perspectiva no interpola el atributo: interpola el atributo dividido por w y también 1/w, ambos linealmente en espacio de pantalla, y en el fragmento divide uno por otro para recuperar el valor. Esa división es la que amplifica el error, y lo amplifica más cuanto mayor es w, es decir, más lejos. Un vértice a diez kilómetros contribuye a un fragmento con un w enorme, y el cociente de dos cantidades pequeñas y ruidosas produce normales sucias y coordenadas de textura que bailan.

WGSL te deja elegir el modo con @interpolate. El valor por defecto es perspective; linear interpola en espacio de pantalla sin dividir, que es lo correcto para cantidades que ya viven en pantalla; y flat no interpola nada, toma el valor del vértice de referencia y por tanto no introduce ningún error. Para índices, identificadores de material y cualquier cosa que deba ser exacta, flat no es una optimización, es una condición de corrección.

struct Interpolado {
  @builtin(position) posicion: vec4f,
  @location(0) normal: vec3f,                        // perspective, por defecto
  @location(1) @interpolate(flat) indiceMaterial: u32, // exacto, sin interpolar
}

Y queda el último frente. Las funciones trascendentes de WGSL —sin, cos, exp, pow, log, inverseSqrt— no están garantizadas al ULP. La especificación no exige el resultado correctamente redondeado, y cada implementación usa la aproximación de su hardware, que puede ser una tabla, una serie polinómica o una instrucción dedicada de precisión reducida. Dos GPU pueden devolver valores distintos para sin(x) con el mismo x, y la diferencia puede llegar a varios ULP.

De ahí sale la regla operativa que cierra la lección: nunca compares floats con == en un shader salvo que sepas exactamente por qué. Las excepciones legítimas existen y son pocas: comparar contra un valor que tú mismo escribiste sin operar sobre él, comparar contra cero un contador que solo se ha incrementado con enteros, o el depthCompare: 'equal' de arriba, que funciona precisamente porque @invariant te da la garantía que en cualquier otro sitio no tienes. Para todo lo demás, compara con tolerancia: abs(a - b) < epsilon, con un epsilon que escales con la magnitud de los operandos, porque un epsilon absoluto de 1e-6 es generosísimo cerca del uno e inservible cerca de mil.