Is there a way to find the error in values like `.log(2)` and `.pi / 3`?

[Nevin]

Suppose I have a value x of some type T: Real, which is close to ln(2), and I want to know how far off it is. Calculating the difference y = x - .log(2) causes catastrophic cancelation, with only a few bits of precision remaining. We could regain the lost precision by subtracting from y an accurate value for ln(2) - .log(2), but I don’t know how to find it.

For comparison, if we’re interested in π rather than ln(2), we can find an accurate value for π - .pi like so:

extension Real {
  static var truePiMinusSelfPi: Self {
    // sin(π - x) ≈ x for small x, so use x = π - .pi
    return .sin(.pi)
  }
}

Which lets us calculate distance from π as:

let y = (x - .pi) - .truePiMinusSelfPi

If T is a radix-2 type, then this trivially generalizes to π/2, π/4, etc. But if we care about π/3, or if we’re writing generic code that doesn’t know the radix of T, it is not so obvious how to proceed.

Is there a good way to find ln(2) - .log(2), or π/3 - (.pi / 3), or other similar values?

[scanon (Steve Canon)]

The most general approach here would be to evaluate ln(2) via series approximation in head-tail arithmetic, giving you about twice as many bits as the format provides, and then compute x - head - tail, yielding an answer with roughly full precision.

[Nevin]

Thanks, I was attempting to address issue #158 (about the real part of expMinusOne) by finding particular points along the critical curve near which we can evaluate to high precision.

It technically worked, but it only really helped extremely close to those particular points, so I’m not sure it’s a promising direction after all.

---

As a quick sketch, consider the point (ln(2), π/3), where we have (e^x, cos(y)) = (2, 1/2).

Take a nearby point (x, y) = (ln(2) + δ, π/3 + ε). Then e^x = 2e^δ and:

cos(y) = (1/2)cos(ε) - (√3/2)sin(ε)
= (1 - 2sin²(ε/2) - sin(ε)√3) / 2

So:

e^x cos(y) - 1 = e^δ (1 - 2sin²(ε/2) - sin(ε)√3) - 1
= (e^δ - 1) - e^δ (2sin²(ε/2) + sin(ε)√3)

Thus, if we have accurate values for δ = x - ln(2), and ε = y - π/3, then we can evaluate the above expression to get an accurate value for e^x cos(y) - 1.

[scanon (Steve Canon)]

I'm not sure that would resolve the problem anyway; in particular close to the curve e^δcos(π/3+ε) = 1/2 you're back in the same cancellation scenario. It's a thorny problem.

[Nevin]

@scanon

close to the curve e^δcos(π/3+ε) = 1/2 you're back in the same cancellation scenario

Right, but ε and δ are small, so the errors have been reduced from the order of .ulpOfOne to the order of .ulpOfOne.ulp.

---

To give a little more detail, the expression we’re evaluating is:

(e^δ - 1) - e^δ (2sin²(ε/2) + sin(ε)√3)

The first term (e^δ - 1) ≈ δ, and evaluating it as .expMinusOne(δ) has an error around δ.ulp.

The last term sin(ε)√3 ≈ ε√3, and evaluating it as .sin(ε) * .sqrt(3) has an error around ε.ulp.

The term 2sin²(ε/2) ≈ ε²/2 is much smaller and doesn’t change the magnitude or error estimate.

The factor e^δ ≈ 1, so multiplying by .exp(δ) doesn’t change them either.

Thus we have something of order δ minus something of order ε. There may indeed be catastrophic cancelation when we subtract them, but the error is then of order max(δ.ulp, ε.ulp).