patternMinor
Floating-point modular multiplication algorithm
Viewed 0 times
multiplicationpointfloatingmodularalgorithm
Problem
Is there a well-known algorithm for modular multiplication of floating-point numbers?
I would like to multiply some large angle in single precision (6-7 significant digits) and wrap it back to 360 degrees, without losing too many significant digits.
My current approach is to split the numbers
Does a more efficient algorithm exist?
Edit: Here's an example expression:
where
Unfortunately, double-precision math is not available.
I would like to multiply some large angle in single precision (6-7 significant digits) and wrap it back to 360 degrees, without losing too many significant digits.
float r = fmodf(a * b, 360);My current approach is to split the numbers
a and b into integer and decimal parts, get the remainder (a_int * b_int) % 360, and add it to the smaller decimal products.Does a more efficient algorithm exist?
Edit: Here's an example expression:
float L = fmodf(13.17639647 * d, 360);where
d is a decimal number between -36525 and 36525 (it could theoretically be a wider range).Unfortunately, double-precision math is not available.
Solution
A completely general solution for this is difficult to achieve if a wider floating-point format with at least twice the number of significand bits is not available. Normally, on common system platforms using C or C++ one could just use
A high-performance implementation is possible in the absence of
It is further assumed that the product $a \cdot b$ is restricted to $[-360 \cdot2^{22}, 360 \cdot 2^{22}]$, approximately $[-1.5\cdot 10^{9}, 1.5\cdot 10^{9}]$, such that $\lfloor \frac{a \cdot b}{360}\rfloor$ is exactly representable in a
With the help of FMA and the reciprocal $\frac{1}{360}$ one can quickly compute $\mathrm{nint}(\frac{a \cdot b}{360})$, where $\mathrm{nint}(\cdot)$ is the nearest integer function, and implement Kahan's accurate difference-of-products computation $ab-cd$. The remainder is computed as $a\cdot b - \mathrm{nint}(\frac{a \cdot b}{360}) \cdot 360$. To match
The final
`/*
Compute ab-cd with error
float rem = (float) fmod ((double)a * (double)b, 360.0); with float mapped to IEEE-754 binary32 type, double mapped to IEEE-754 binary64 types and therefore both the multiplication and fmod() being exact, with only a single rounding error incurred when rounding the result back to float.A high-performance implementation is possible in the absence of
double support if a few parameters and restrictions can be enforced. This answer assumes that a C or C++ program is executed on a system platform that supports IEEE-754 arithmetic and with the binary32 data type exposed as float, and that the default rounding mode in effect is "to nearest or even". It is further assumed that the fused multiply-add (FMA) operation is provided in hardware, and exposed (for the binary32 type) via the standard math function fmaf(). This latter assumption holds for most current systems based on the ARM, x86, and Power processor architectures, as well as commonly used GPUs.It is further assumed that the product $a \cdot b$ is restricted to $[-360 \cdot2^{22}, 360 \cdot 2^{22}]$, approximately $[-1.5\cdot 10^{9}, 1.5\cdot 10^{9}]$, such that $\lfloor \frac{a \cdot b}{360}\rfloor$ is exactly representable in a
float number. With this restriction in place, the desired result can be computed with an error of at most 2 ulp as follows.With the help of FMA and the reciprocal $\frac{1}{360}$ one can quickly compute $\mathrm{nint}(\frac{a \cdot b}{360})$, where $\mathrm{nint}(\cdot)$ is the nearest integer function, and implement Kahan's accurate difference-of-products computation $ab-cd$. The remainder is computed as $a\cdot b - \mathrm{nint}(\frac{a \cdot b}{360}) \cdot 360$. To match
fmod() semantics, the remainder must have the same sign as the dividend $ab$. This can be achieved by conditionally adding / subtracting $360$ if the condition is not met.The final
float result may round to exactly $360$. If this is problematic, an additional check should be added to the algorithm below to map this to zero.`/*
Compute ab-cd with error
Context
StackExchange Computer Science Q#156897, answer score: 3
Revisions (0)
No revisions yet.