[Python-ideas] Fwd: Trigonometry in degrees
Tim Peters
tim.peters at gmail.com
Thu Jun 14 14:44:34 EDT 2018
I should note that numeric code "that works" is often much subtler than it
appears at first glance. So, for educational purposes, I'll point out some
of what _wasn't_ said about this crucial function:
[Tim]
> import mpmath
> from math import fmod
> # Return (n, x) such that:
> # 1. d degrees is equivalent to x + n*(pi/2) radians.
> # 2. x is an mpmath float in [-pi/4, pi/4].
> # 3. n is an integer in range(4).
> # There is one potential rounding error, when mpmath.radians() is
> # used to convert a number of degrees between -45 and 45. This is
> # done using the current mpmath precision.
> def treduce(d):
> d = fmod(d, 360.0)
> n = round(d / 90.0)
> assert -4 <= n <= 4
> d -= n * 90.0
> assert -45.0 <= d <= 45.0
> return n & 3, mpmath.radians(d)
>
> How do we know there is at most one rounding error in that? No, it's not
obvious.
That `fmod` is exact is guaranteed by the relevant standards, but most
people who write a libm don't get it right at first. There is no "cheap"
way to implement it correctly. It requires getting the effect of doing
exact integer division on, potentially, multi-thousand bit integers.
Assuming x > y > 0, a correct implementation of fmod(x, y) ends up in a
loop that goes around a number of times roughly equal to log2(x/y),
simulating one-bit-at-a-time long division For example, here's glibc's
implementation:
https://github.com/bminor/glibc/blob/master/sysdeps/ieee754/dbl-64/e_fmod.c
Don't expect that to be easy to follow either ;-)
Then how do we know that `d -= n * 90.0" is exact? That's not obvious
either. It follows from the "Sterbenz lemma", one version of which: if x
and y are non-zero floats of the same sign within a factor of 2 of each
other,
1/2 <= x/y <= 2 (mathematically)
then x-y is exactly representable as a float too. This is true regardless
of the floating-point base, or of rounding mode in use. In IEEE-754, it
doesn't even need weasel words to exempt underflowing cases.
That lemma needs to be applied by cases, for each of the possible values of
(the integer) `n`. It gets closest to failing for |n| = 1. For example,
if d is a tiny bit larger than 45, n is 1, and then d/90 is
(mathematically) very close to 1/2.
Which is another thing that needs to be shown: "if d is a tiny bit larger
than 45, n is 1". Why? It's certainly true if we were using infinite
precision, but we're not. The smallest representable double > 45 is 45 +
2**-47:
>>> d = 45 + 2**-47
>>> d
45.00000000000001
>>> _.hex()
'0x1.6800000000001p+5'
Then d/90.0 (the argument to round()) is, with infinite precision,
(45 + 2**-47)/90 =
0.5 + 2**-47/90
1 ULP with respect to 0.5 is 2**-53, so that in turn is equal to
0.5 + 2**-53/(90/64) =
0.5 + (64/90)*2**-53 =
0.5 + 0.71111111111... * 2**-53
Because the tail (0.711...) is greater than 0.5 ULP, it rounds up under
nearest-even rounding, to
0.5 + 2**-53
>>> d / 90
0.5000000000000001
>>> 0.5 + 2**-53
0.5000000000000001
and so Python's round() rounds it up to 1:
>>> round(_)
1
Note that it would _not_ be true if truncating "rounding" were done, so
round-nearest is a hidden assumption in the code.
Similar analysis needs to be done at values near the boundaries around all
possible values of `n`.
That `assert -45.0 <= d <= 45.0` can't fall then follows from all of that.
In all, a proof that the code is correct is much longer than the code
itself. That's typical. Alas, it's also typical that math library sources
rarely point out the subtleties.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-ideas/attachments/20180614/c65669ee/attachment-0001.html>
More information about the Python-ideas
mailing list