Payne-Hanek reduction modulo π/2 #
An argument x = m * 2 ^ e with m > 0 is reduced modulo π/2 from a window of the bits of
2/π. Let T be the cached lower bound of 2/π * 2 ^ N at a level of Constants with
N ≥ e + window + 6 bits, so that x * 2/π exceeds m * T * 2 ^ (e - N) by less than
2 m * 2 ^ (e - N). The bits of T of weight 2 ^ (N - e + 2) and above contribute multiples of
four, which change neither the quadrant nor the remainder, and the bits below the window contribute
less than m units of its last place. So the reduction multiplies m by the window bits of T
below that weight and keeps the product modulo 2 ^ window: one shift of T, one product and one
remainder, whatever the size of x * 2/π. The level is chosen from e: the window needs
b = max e 0 + window + 6 bits of 2/π, and the level read holds at most max 64 (2 b) bits,
computed once, on first use.
The result is a quadrant q and a fraction F with x * 2/π = 4 j + q + (F + η) / 2 ^ scale for
an integer j and some η with 0 ≤ η < error = 2 m (ArgumentReduction.Proof). The fraction has
magnitude at most half a quarter turn. When QuarterTurns.accepts bits holds, F has the sign of
F + η and relative error at most 2 ^ -bits. A rejected reduction is retried with a doubled
window: the error stays 2 m units while the units shrink, so some window is accepted for every
nonzero argument (ArgumentReduction.Termination). reduce bounds the number of doublings and
returns its first accepted attempt or else its last one, so it needs no search with a termination
proof. Callers check QuarterTurns.accepts on its result. fuelWithin? gives the most doublings
whose windows stay within a level of the cache.
QuarterTurns.argumentNumerator multiplies the fraction by the cached lower bound of π/2 to give
the reduced argument in radians, and rotate gives the sine and cosine of x from those of the
reduced argument. Binary sine and cosine read this reduction.
The method is that of M. H. Payne and R. N. Hanek, Radian reduction for trigonometric functions, ACM SIGNUM Newsletter 18(1), 1983, pages 19-24.
An argument x reduced modulo π/2: x * 2/π = 4 j + quadrant + (fraction + η) / 2 ^ scale for an
integer j and an unknown η with 0 ≤ η < error.
- quadrant : Fin 4
The quotient of
xbyπ/2, modulo four. - fraction : ℤ
The remainder in units of
2 ^ -scalequarter turns, below the unknown partη. - scale : ℕ
Fractional bits of
fraction. - error : ℕ
A strict upper bound of
η, in units of2 ^ -scalequarter turns.
Instances For
Instances For
Instances For
The cache level that reduceWindow reads: at least e + window + 6 and window + 6 bits.
Instances For
Reduce m * 2 ^ e modulo π/2 with window bits of 2/π.
With T the lower bound of 2/π * 2 ^ N at level windowLevel e window and
shift = N - e + 2 - window, piece holds the bits of T from shift to shift + window. The
product y = m * piece mod 2 ^ window is x * 2/π modulo four in units of 2 ^ -scale, with
scale = window - 2. Rounding y to the nearest unit 2 ^ scale gives the quadrant, and the
signed difference is the fraction. The bits of T below the window and the error of T make the
true remainder exceed the fraction by less than 2 m units.
Instances For
Whether a reduction keeps bits relative bits: error * 2 ^ (bits + 1) ≤ |fraction|. Then the
fraction has the sign of the true remainder and differs from it by at most 2 ^ -bits of it.
Instances For
The first window for a significand of p bits and bits relative bits. Its scale is
2 p + bits + 18. A remainder of f quarter turns from the nearest multiple of π/2 gives
|fraction| > |f| * 2 ^ scale - 2 m, and acceptance asks for 2 m * 2 ^ (bits + 1), so for
p ≤ bits + 1 the first window accepts every remainder of at least 2 ^ -(p + 16) quarter turns.
That covers the closest approaches to multiples of π/2 in the hard-case corpus of binary16,
bfloat16, binary32, binary64 and binary128.
Instances For
Window of attempt attempt: the start window doubled attempt times.
Instances For
Reduce m * 2 ^ e with the start window for a p-bit significand, doubling the window until
the reduction passes QuarterTurns.accepts bits or fuel doublings are used. Returns the first
accepted attempt, or the last attempt when every one is rejected, so callers that need bits
relative bits check QuarterTurns.accepts bits on the result. Every attempt reads at most level
windowLevel e (attemptWindow p bits fuel).
Instances For
Try attempt attempt and, while it is rejected, at most fuel attempts after it.
Instances For
The fuel of reduce within a resource limit: the largest count ≤ fuel whose window
attemptWindow p bits count reads a cache level of at most limit at exponent e, or none when
not even the first window does. The attempts before it have narrower windows, so they read at most
limit too.
Instances For
The reduced argument in radians as an integer over 2 ^ (scale + precision): the fraction times
the cached lower bound halfPiLower precision of π/2 * 2 ^ precision. It differs from the reduced
argument (fraction + η) / 2 ^ scale * π/2 by at most error / 2 ^ scale * π/2 plus
|fraction| / 2 ^ scale * 2 ^ (1 - precision).