4 ms·
That's a neat trick, "cancelling" out catastrophic cancellation. Of course it doesn't work anymore once x is smaller than machine epsilon, whereas expm1(x) / x
by cgadski 3y ago
That's a neat trick, "cancelling" out catastrophic cancellation. Of course it doesn't work anymore once x is smaller than machine epsilon, whereas expm1(x) / x will continue to work.
In general herbie is pretty good at suggesting the right functions / rearrangements. For this example, it finds expm1: https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e00f46723916.f54ecc7459e2826747d647a0d98ec9e41c611ae9/graph.html https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e...
- zokier 3y agoGotta love that Herbie also provides just `1` as alternative. "66.8% accurate, 305.0× speedup" But the expression you get if you disable "numerics" ruleset (which includes expm1) is pretty wild also: (if (<= (exp x) 4504410275303423/4503599627370496) (+ 1 (+ (* 1/24 (pow x 3)) (* x (+ 1/2 (* x 1/6))))) (/ (+ (exp x) -1) (log (exp x)))))
- aragilar 3y agoThat doesn't surprise me, as before floats are involved you need to transform the series into something you can compute (doing this by hand is a good way to understand what's going on underneath). https://developers.redhat.com/blog/2015/01/02/improving-math-performance-in-glibc https://developers.redhat.com/blog/2015/01/02/improving-math... gives a basic intro to how these functions are actually implemented, and we can if you compute exp(x) - 1 near 0 that 1 + (small stuff) - 1 is always going to be a problem (it'd be nice if we had a sufficiently smart compiler that could understand how to inline these things, but then it'd probably do surprising things as well).