# Plot inverse of emulated double-precision floats

I need to accurately plot a line chart using WebGL. The numbers have a precision of around 33 bits - that's too many to fit into a single-precision float's mantissa. WebGL does not support the double-precision floats that arrive in my program so I split each number into "coarse" and "fine" parts:

// buf is a Float32Array to be used as a vertex buffer, input is a Float64Array.
for (let i = 0; i < input.length; ++i) {
buf[i * 2] = input[i];
buf[i * 2 + 1] = input[i] - buf[i * 2];
}


The y axis range that shall be visible is passed in the same manner:

// rangeMin and rangeMax are 64-bit floats that define the plotting boundaries.
// rangeMin is split up just like the input numbers.
// The value range is only needed in single precision.
gl.uniform3f(range,
Math.fround(rangeMin),
rangeMin - Math.fround(rangeMin),
rangeMax - rangeMin);


Then, in the vertex shader I evaluate the y axis position by joining coarse- and fine-grained parts in the correct order:

float y = ((n.x - range.x) + (n.y - range.y)) / range.z;


For points that are inside the given range, y is the correct position on a normalized y axis between 0 and 1.

This works really well. However, I also need to find the correct position on an inversed y axis, i.e. an axis that linearly spans from 1 / rangeMax to 1 / rangeMin. I could easily create a second vertex buffer and reuse the existing math, only substituting input1 / input, rangeMin1 / rangeMax and rangeMax1 / rangeMin. This would waste a lot of memory on the duplicated buffer. I feel like there should be a relatively simple way to approximate such kind of inverse in the vertex shader, maybe involving one or two steps of Newton's method, but I can't seem to figure it out.

• What is the point of keeping maximum accuracy for graphics ? Typical screen coordinates have an accuracy of 10 to 12 bits. Even printing on large sheets in high resolution will never require more than 20 bits.
– user1703
Jun 27, 2023 at 9:26

You could try the Taylor series expansion for $$1/x$$: $$\frac{1}{x + a} = \frac{1}{x} - \frac{a}{x^2} + \frac{a^2}{x^3} - \frac{a^3}{x^4} + \cdots$$ when $$a$$ is small compared to $$x$$. You can probably get away with just a handful of terms (maybe even just the first two).
• This seems promising although I still struggle to actually apply it. I see that I can approximate $\frac{1}{x+a}-\frac{1}{x}$ using this Taylor series. But I'm not sure how to apply that to my problem. Furthermore even the first term $-\frac{a}{x^2}$ won't have many significant digits because it's obviously numerically unstable, and later terms didn't contribute to the result in my tests - however it might still be sufficient. I only don't know what to do with this approximation! Sep 29, 2022 at 22:43
• It sounds like you're on the right track - use $1/x$ as the "coarse" part and $-a/x^2$ as the "fine" part for the reciprocal. I'm not sure what you mean about it being numerically unstable? If $|a| \ll |x|$, then $|a/x^2| \ll |1/x|$ correspondingly, so it's going to be about the same ratio between the coarse and fine parts. It will have about as many significant digits as $a$ has to begin with. Sep 29, 2022 at 23:04
• This is really helpful input, thanks. I compared taking only $1/x$ with adding $-a/x^2$ as fine part, which only improves precision by around .36 bits on average in my tests. But it doesn't do so always, depending on input. After much pondering I think the rounding error in the division $1/x$ is to blame, so I am unsure if this approach can work. With "numerically unstable" I meant that because $a$ is so small compared to $x^2$, $a/x^2$ might already be close to underflowing but that's "obviously" wrong, sorry. Sep 30, 2022 at 23:13