<p>You could try the Taylor series expansion for <span class="math-container">$1/x$</span>:
<span class="math-container">$$
\frac{1}{x + a} = \frac{1}{x} - \frac{a}{x^2} + \frac{a^2}{x^3} - \frac{a^3}{x^4} + \cdots
$$</span>
when <span class="math-container">$a$</span> is small compared to <span class="math-container">$x$</span>. You can probably get away with just a handful of terms (maybe even just the first two).</p>

CC BY-SA 4.0
2022-09-29T22:43:46.993Z
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!
purefanatic
CC BY-SA 4.0
2022-09-29T23:04:11.493Z
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.
Nathan Reed
CC BY-SA 4.0
2022-09-29T23:07:52.497Z
Oh BTW, looks like there's a paper on this that might be of interest: https://hal.archives-ouvertes.fr/hal-00957379/document
Nathan Reed
CC BY-SA 4.0
2022-09-30T23:13:27.533Z
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.
purefanatic
CC BY-SA 4.0
2022-09-30T23:14:05.863Z
I'll read the paper tomorrow and let you know my findings, it looks perfect for my use case.
purefanatic
CC BY-SA 4.0
2022-09-27T23:55:39.590Z
2022-09-27T23:55:39.590Z
Nathan Reed
<p>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 &quot;coarse&quot; and &quot;fine&quot; parts:</p>
<pre><code>// buf is a Float32Array to be used as a vertex buffer, input is a Float64Array.
for (let i = 0; i &lt; input.length; ++i) {
    buf[i * 2] = input[i];
    buf[i * 2 + 1] = input[i] - buf[i * 2];
}
</code></pre>
<p>The y axis range that shall be visible is passed in the same manner:</p>
<pre><code>// 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);
</code></pre>
<p>Then, in the vertex shader I evaluate the y axis position by joining coarse- and fine-grained parts in the correct order:</p>
<pre><code>float y = ((n.x - range.x) + (n.y - range.y)) / range.z;
</code></pre>
<p>For points that are inside the given range, y is the correct position on a normalized y axis between 0 and 1.</p>
<p>This works really well. However, I also need to find the correct position on an <em>inversed</em> 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 <code>input</code> → <code>1 / input</code>, <code>rangeMin</code> → <code>1 / rangeMax</code> and <code>rangeMax</code> → <code>1 / rangeMin</code>. 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.</p>

<mathematics><algorithm><webgl><float>
Plot inverse of emulated double-precision floats
purefanatic
# Plot inverse of emulated double-precision floats

(Asked by: username_0 on 2022-09-26)

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 `input` → `1 / input`, `rangeMin` → `1 / rangeMax` and `rangeMax` → `1 / 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.

**Tags:** mathematics, algorithm, webgl, float


---

## 1 answers

### Answer by username_1 on 2022-09-27. Score: 0

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).

> **Comments:**
> **username_0**: 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!  
> **username_1**: 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.  
> **username_1**: Oh BTW, looks like there's a paper on this that might be of interest: https://hal.archives-ouvertes.fr/hal-00957379/document  
> **username_0**: 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.  
> **username_0**: I'll read the paper tomorrow and let you know my findings, it looks perfect for my use case.  

---

