基于GPU.js实现Mandelbrot集的高精度浮点数方案咨询

Absolutely! Combining multiple standard floating-point numbers to build a higher-precision "composite float" is a perfect fit for your needs—you don’t need full arbitrary-precision libraries, just a targeted boost that works within GPU.js’s constraints. Here’s a practical, GPU-friendly approach tailored to Mandelbrot sets:

Core Idea: Split Values into High/Low Precision Components #

Instead of using a single Number (64-bit double-precision in JS), represent each value (real/imaginary parts of z and c) as two 64-bit floats:

• value_high: The main, large-magnitude component of the number
• value_low: A small offset that captures the precision lost by value_high

For example, if you have a number like 0.1234567890123456789, value_high would hold 0.12345678901234568 (the closest 64-bit float), and value_low would hold the difference between the true value and value_high to make up the missing digits.

Implement High-Precision Operations for Mandelbrot #

Mandelbrot only relies on addition and squaring (since z = z² + c is the core iteration). Here’s how to implement these with your composite floats:

1. High-Precision Addition #

Adding two composite values (a_h, a_l) and (b_h, b_l):

• First add the high components: sum_h = a_h + b_h
• Calculate the error introduced by this addition to compute the low component:

  const temp = sum_h - a_h;
  const sum_l = (a_h - (sum_h - temp)) + (b_h - temp) + a_l + b_l;
  

• Normalize the result: if sum_l is large enough to affect sum_h, carry over the excess to keep sum_l as a small offset relative to sum_h.

2. High-Precision Squaring (Critical for Mandelbrot) #

Squaring a composite value (x_h, x_l) (used for both real and imaginary parts of z):

• Compute the three parts of the product:

  const xh_sq = x_h * x_h;
  const cross = x_h * x_l * 2;
  const xl_sq = x_l * x_l;
  

• Combine these into a new composite value:
  • Start with new_h = xh_sq + cross
  • Then compute new_l = (xh_sq - (new_h - cross)) + xl_sq
  • Normalize to keep new_l as a small offset from new_h

For complex squaring (z² = (zr² - zi²) + (2*zr*zi)i), apply the above operations to the real and imaginary components separately, using high-precision subtraction (just addition with a negative value) for the real part.

GPU.js Implementation Tips #

Since GPU.js works exclusively with numerical operations (no strings or external libraries), you can pass composite values as arrays or separate float parameters to your kernel:

Example Kernel Snippet #

Here’s a simplified GPU.js kernel that handles high-precision complex squaring and addition for Mandelbrot:

const gpu = new GPU();

const mandelbrotKernel = gpu.createKernel(function(crH, crL, ciH, ciL, maxIterations) {
  const x = this.thread.x;
  const y = this.thread.y;
  
  // Initialize z with c (high/low components)
  let zrH = crH[x][y];
  let zrL = crL[x][y];
  let ziH = ciH[x][y];
  let ziL = ciL[x][y];
  
  let iterations = 0;
  while (iterations < maxIterations) {
    // Compute zr² - zi² (real part of z²)
    // First calculate zr²
    const zrHSq = zrH * zrH;
    const zrCross = zrH * zrL * 2;
    const zrLSq = zrL * zrL;
    let zrSqH = zrHSq + zrCross;
    let zrSqL = (zrHSq - (zrSqH - zrCross)) + zrLSq;
    
    // Calculate zi²
    const ziHSq = ziH * ziH;
    const ziCross = ziH * ziL * 2;
    const ziLSq = ziL * ziL;
    let ziSqH = ziHSq + ziCross;
    let ziSqL = (ziHSq - (ziSqH - ziCross)) + ziLSq;
    
    // High-precision subtraction: zr² - zi²
    let newZrH = zrSqH - ziSqH;
    const tempSub = newZrH - zrSqH;
    let newZrL = (zrSqH - (newZrH - tempSub)) + (ziSqH - tempSub) * -1 + zrSqL - ziSqL;
    
    // Compute 2*zr*zi (imaginary part of z²)
    // First multiply zr and zi
    const zrZiHH = zrH * ziH;
    const zrZiHL = zrH * ziL;
    const zrZiLH = zrL * ziH;
    const zrZiLL = zrL * ziL;
    let prodH = zrZiHH + zrZiHL + zrZiLH;
    let prodL = (zrZiHH - (prodH - zrZiHL - zrZiLH)) + zrZiLL;
    // Multiply by 2
    prodH *= 2;
    prodL *= 2;
    
    // Add c to z² (new z = z² + c)
    // Update real part
    const sumZrH = newZrH + crH[x][y];
    const tempZrSum = sumZrH - newZrH;
    zrL = (newZrH - (sumZrH - tempZrSum)) + (crH[x][y] - tempZrSum) + newZrL + crL[x][y];
    zrH = sumZrH;
    
    // Update imaginary part
    const sumZiH = prodH + ciH[x][y];
    const tempZiSum = sumZiH - prodH;
    ziL = (prodH - (sumZiH - tempZiSum)) + (ciH[x][y] - tempZiSum) + prodL + ciL[x][y];
    ziH = sumZiH;
    
    // Check for escape (use high component for fast early exit)
    const magSqH = zrH * zrH + ziH * ziH;
    if (magSqH > 4.0) break;
    
    iterations++;
  }
  
  return iterations / maxIterations;
}).setOutput([1024, 1024]);

Key Optimizations #

• Early Escape Check: Use only the high-precision component to check if the magnitude exceeds 4—this avoids unnecessary low-precision calculations and speeds up the kernel.
• Normalization: After each operation, ensure the low component is always a small offset relative to the high component to prevent precision overlap.
• Test on CPU First: Validate your high-precision operations on the CPU first (using the same composite float logic) before porting to GPU.js—this makes debugging much easier.

Why This Works #

Using two 64-bit floats gives you roughly 128 bits of effective precision, which is more than enough to eliminate pixelization when zooming deep into the Mandelbrot set. This approach avoids the need for string-based precision libraries and plays perfectly with GPU.js’s numerical-only execution model.

内容的提问来源于stack exchange，提问作者Benedikt Vilji Magnússon