Real part of complex expMinusOne is bad when e^{-x} ~ cos y
- Dominant language
- Swift
- Stars
- 1.9k
- Forks
- 181
- Avg merge
- 2d 5h
- Merged PRs (30d)
- 6
Description
For example, in binary32:
x = 1.09590280055999755859375
y = 1.230000019073486328125
The real part of e^{x + iy} - 1, computed in high precision with Sollya and rounded to binary32, is
-3.4550861727211668039672076702117919921875e-8,
but the expression `fma(expm1(x),cos(y), -versin(y))` as used at https://github.com/apple/swift-numerics/blob/c26e60c6ba6e545cf8d3048a3cee554058f82db1/Sources/ComplexModule/ElementaryFunctions.swift#L117 gives (rounding each subexpression to binary32 in Sollya)
-6.9330411633927724324166774749755859375e-8,
in which every digit is wrong.
The problem is catastrophic cancellation when e^x cos y ~ 1. Factoring it to compute the addition with FMA still trips over catastrophic cancellation because (e^x - 1) cos y ~ versin y if e^x cos y ~ 1. Same goes for the other obvious factorings as e^x - 1 - e^x versin y, or e^{x + log cos y} - 1 = e^{x + log(1 - versin y)} - 1, although they do worse when x >>> 1 because the low resolution of floating-point values around y ~ acos(e^{-x}) gives no more opportunity for cancellation of versin y against e^x cos y once y = fl(pi/2).
I don't see any clear way to compute the error in either term short of computing the transcendental functions in extended precision. I guess that might be easy for binary32, but in binary64 you'd presumably have to do the transcendental functions in double-Dekker arithmetic. Maybe there's a clever algebraic trick to avoid the subtraction altogether, but I haven't thought of one.
---
Nit: There appears to be a typo in the comment at: https://github.com/apple/swift-numerics/blob/c26e60c6ba6e545cf8d3048a3cee554058f82db1/Sources/ComplexModule/ElementaryFunctions.swift#L86
It reads
```
exp(x) cosMinuxOne(y) + expMinusOne(y)
```
but it should read
```
exp(x) cosMinusOne(y) + expMinusOne(x)
```
Contributor guide
Research direction
Read Sources/ComplexModule/ElementaryFunctions.swift around lines 86 and 117, then reproduce the binary32 example from the issue. Determine how the real part of complex expMinusOne should avoid the reported cancellation and verify the corrected result against the high-precision value; also correct the noted comment typo.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- swift
- Domain
- tooling
- Issue type
- Bug
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Stale
- Clarity
- Mostly clear
- Newbie friendliness
- 35/100