# How to handle periodic functions of very large floats

(Or really just poor conditioning in general.)

Recently, there has been some interest in correctly-rounded floating-point libraries (RLIBM and CORE-MATH), which I think is great, but there is one thing that gives me pause. Consider $$\sin(2^{28})$$ in single-precision. The period of $$\sin$$ is $$2\pi$$, which is appreciably smaller than an ulp at $$2^{28}$$ (16), so the result is unlikely to be meaningful. Actually, there are two issues:

1. The implementation of $$\sin$$ must go to a lot of extra effort to implement range reduction of such large arguments.

2. If the argument represents the result of a measurement or a computation—rounded even once—then the result will almost certainly be completely useless to the code that requested it; there is almost certainly a programming error.

When cases such as these can be trivially detected (as is the case for the sine of a large argument), does it make sense to signal an exception instead of returning a likely-meaningless result? Is there prior art on this?

What should the cutoff point be? Presumably, it should depend on the rounding interval of the argument, but that rounding interval will depend on the rounding mode; should that be taken into account?

• Arguably, one shouldn't use a floating-point representation at all for the input of a periodic function like sin; a fixed-point representation which wraps at the correct period would ensure not only equal precision across the whole domain, but also would not waste bits on redundant representations (e.g. floats can distinguish between PI, 3 * PI and 5 * PI, but there is no need for that if you're going to take the sin). A floating-point representation for the argument of sin might be useful if you specifically only want to compute it e.g. near 0 rather than across the whole domain. Commented Aug 15, 2023 at 23:22
• That sidesteps the problem, but does not solve it. The argument will likely be the result of some other floating-point computation, so now you effectively force the user to convert and range-reduce themselves, which is inconvenient and irregular. It also doesn't answer the more general question of what happens if you are able to detect pathologically bad conditioning in non-periodic cases. Commented Aug 15, 2023 at 23:48
• "Is there any prior art on this?" - One example was in BBC BASIC, which raised the error Accuracy lost in this situation (see eg page 462 of the manual). Commented Aug 16, 2023 at 13:49
• How to handle? The answer probably differs for numerical experts and the "general public". Commented Aug 16, 2023 at 15:36

I think this is in the category of "probable mistakes" ─ situations where you have evidence that the code is doing something the programmer didn't intend, but it's still possible they did intend it.

For example, it's plausible that a number theorist might actually want to compute the value of something like $$\sum_{k=0}^{28} \sin 2^k$$; sure, it looks strange to take the sin in radians of exact integers, but it is the kind of thing that number theorists are known to do. In that case, there is no "loss of precision" due to the scale of the argument: $$2^{28}$$ is the precise value that the programmer wants the sin of, and it is exactly representable as a single-precision float.

So I would say that if you do want to do something different when there is a potential loss of precision, then either it should be only a warning (and there should be a way to suppress it), or if you do raise an error at runtime then you should offer a way to explicitly do the calculation without that check (e.g. a separate sin_unchecked function).

For comparison, consider the following C or Java code, which is almost certainly a mistake:

while(do_thing());
do_other_thing();


The semicolon after the while condition is a statement that does nothing, which means the loop body does nothing and do_other_thing() is not part of the loop. This is a fairly common mistake in real code. But it's also not absolutely certain that it's a mistake; it's plausible to write while loops with empty bodies when the condition is the thing you want to loop, and it's plausible to have incorrect indentation.

This is what I mean by the category of "probable mistakes". Making it a warning is an option (and leaving the warning to a linter is another option; some real languages do). Making it an error is also an option ─ for example, you could require braces around the body of a while statement, so that an empty ; statement is not a valid body.

But because it could be intentional to have a loop with an empty body, the language should still have a way for the programmer to say "actually, I do want the loop body to be an empty statement". In this example, that's easily satisfied by allowing {} to be the empty statement instead of ;. But it's something you should consider, if you're going to make something an error.

• I definitely agree it could be useful to take the sine of 2^28, but it just doesn't make sense to me to use single floats for that, even though they can represent 2^28 exactly. What if you want to repeat the experiment with a base of 3 rather than 2?—then you're SOL. And the time taken for range reduction will increasingly dominate, so that you don't get much (if anything) in the way of a performance improvement. If you're doing such experiments, I can't think of a reason not to use extended-precision numbers. Commented Aug 16, 2023 at 19:06
• @Moonchild Yes, it wouldn't make much sense to use single-precision floats for that, I agree. But then again it doesn't make much sense to use single-precision floats for most things, except when you need a lot of them (e.g. in graphics applications where you are doing computations per pixel per 60th of a second) but in that case you usually don't want the whole thing to fail just because one calculation is suspect, and you don't want to pay extra to check that the magnitudes of the numbers you're taking sin of are small enough to be reasonable. Commented Aug 16, 2023 at 19:11
• The point I'm making is that floats can be very imprecise at those scales, but the magnitude of the ULP doesn't necessarily correspond with imprecision, because sometimes a float value is the exact number you intend it to be. So it's only evidence (albeit fairly good evidence) that it's a mistake, not definite proof. Commented Aug 16, 2023 at 19:14
• @kaya3: The expression x % TWO_PI would be precisely representable if TWO_PI was the name of a precisely representable value, but range reduction on sin(x) uses a modulus of 2π, which is not precisely representable, and thus x % π would generally not be precisely representable either. Commented Aug 16, 2023 at 21:49
• @kaya3: One thing that makes the notion of argument reduction for sin(x) particularly ironic, btw, is that the majority of code that would be asked to perform sin(x) for some float value is seeking to compute the sine of 2πy for some y, but is approximating that by asking for the sine of (26,353,590/4,194,304)y instead, and attempting to perform accurate range reduction mod 2π would make that a poorer approximation than would performing "less accurate" range reduction mod (26,353,590/4,194,304). Commented Aug 17, 2023 at 15:36

There's two major challenges that I see with this approach:

• The correct tolerance for "meaninglessness" is dependent on the use case
• Exceptional flow control is expensive from a user's perspective.

The tolerance issue is one I deal with often. There's no clear line between a good input and a bad input. I may want X significant digits in my output. How do I indicate that to the function like sin or cos? I've worked on systems where the differences between implementations of cosh play a major role. I've also worked on systems where all I really needed as a number, and if I get an arbitrary output for an absurd input, I'm fine with that.

Also remember that exceptions have to be handled to do any good. That puts a lot of cognitive load on the user's part. While I know signaling NaNs exist in IEEE-754, the major platforms I know of default to quiet NaNs because passing NaN around is deemed better than dealing with the exceptional behavior. I'd say this suggests it is "popular" to not raise exceptions for numerical issues which are worse than $$\sin(2^{38})$$. That's not to say its "right," but it does at least say that a lot of people would rather not have floating point exceptions thrown around. Use of exceptions here would likely do a lot to define your target user-base, and should probably be done in concert with similar numeric stability decisions in other aspects of the language.

• At the point when the size of the rounding interval exceeds the period of the function, I would argue that it's not giving you any information. (Before that point, I do agree that any specific choice is arbitrary.) Your points about handling exceptions seem to me to apply equally equally to bounds checking, yet I have heard no one say bounds checking is a bad idea (perhaps that sophisticated types should be used to preclude out-of-bounds accesses statically, but then you could do the same thing for float operations). Also, you already have to be prepared for the exception from sin(inf). Commented Aug 17, 2023 at 6:37
• "Also remember that exceptions have to be handled to do any good." Yes, but a sufficient exception handling can/should be done at some top level, as catch-all, telling the client that the request failed, with the failure reason (the exception instance) just being a minor detail. And personally, I'd prefer the fail-early behavior of floating point exceptions over continuing with useless NaN "calculations". Commented Aug 17, 2023 at 8:21
• @Moonchild I think a difference between this and bounds checking is that bounds checking occurs on a boundary that is very well agreed upon. There's a clear difference between last-item and one-past-end. In this case, there's a graceful degradation over the course of 64-bit's worth of dynamic range. Commented Aug 18, 2023 at 1:00

If you claim IEEE-754 compliance, the code to handle $$\sin(2^{38})$$ has to exist, and memory is not so constrained that this is a huge problem these days. This implementation effort only has to be expended once, and there is an argument that programmers who try to evaluate $$\sin(2^{38})$$ deserve the performance they get.

Recent editions of the IEEE-754 standard also define functions such as sinPi, cosPi, tanPi, and their inverse versions (atan2Pi has an unfortunate name!), which are defined as, say:

$$\mathrm{sinPi}(x) = \sin(\pi x)$$

The theory is that range reduction to the range $$[-1,1]$$ is an extremely inexpensive operation, and is always exact.

• The main problem I'm concerned with is diagnosing erroneous code; range reduction is just an annoyance. Not super concerned with ieee-754 compliance for this issue. Hence, sinpi etc. isn't really helpful; it has a period of 2, so as an ulp of the argument gets close to 2, it suffers from the same precision loss. Commented Aug 16, 2023 at 5:19
• @Moonchild: Performing argument reduction on sin(16258053) would require subtracting about 16,258,052.99999988539857478731 from the operand. Computing sinpi(5175099) would require subtracting 5,175,098 from the operand. One of those computations would require substantially more precision than the other. Commented Aug 16, 2023 at 21:01
• Returning the fractional part of an IEEE-754 floating point number is an exact operation. It requires no rounding. Commented Aug 16, 2023 at 23:36
• "The theory is that range reduction to the range [−1,1] is an extremely inexpensive operation" --> excpept that it incurs rounding to go from X radians to [-1 .... 1] and that rounds then make sinPi(x) less meaningful for large x. Commented Aug 17, 2023 at 22:16
• Remember, IEEE-754 floating point numbers are a finite set, and over half of them are in the range $[-1,1]$. There is a sense in which this range is what floating point is optimised for, and functions like sinPi use that range more effectively than functions that work in radians. Commented Aug 18, 2023 at 0:52

IMHO, the underlying problem is that a floating-point number is only an approximation for some mathematical real number. From this point of view, it represents a range or interval of real numbers (with the range width being one ULP).

So, instead of ignoring that "interval" property, as our CPUs and libraries typically do, a more math-oriented approach would be to embrace it, and have functions like sin() accept and return a range value. Then e.g. sin(1e30) (with 1e30 denoting a range based on the IEEE double precision ULP) would give a range of [-1,1], thus telling the caller that the sine of the mathematical real number that got approximated to the double value can be anything.

Of course, that needs a whole new approach to math operations, so it might be too much for your project.

An approach better compatible with current libraries and CPUs might be to have an additional requiredPrecision argument to functions like sin(). Then

y = sin(x, 0.0001)


would either return a sine value with the guarantee that for any real number approximated by x, its mathematical sine value will be in the interval [y-0.0001, y+0.0001]. If the sin() function can't guarantee that, an exception would be appropriate. For the sin() function, a sufficient test can be to compare the ULP with the requiredPrecision argument, as the sine derivative is always in the range from -1 to 1.

This way, callers can express their precision requirements and get exceptions when the function can't fulfill that.

Without the callers specifying their requirements, any threshold you select will be arbitrary.

• A float is not its rounding interval, and although explicit interval arithmetic is definitely useful in some cases, it can also overestimate error (see kahan 'how futile are mindless...'), as well as harming performance unacceptably (factor of at least 2, if not more). Commented Aug 16, 2023 at 7:28
• I think that, once you reach the point where the size of the rounding interval is greater than the period of the function you approximate, there's a case to be made that stopping there is not 'arbitrary'. But I agree that any particular choice of stopping point prior to that one will probably be arbitrary. Commented Aug 16, 2023 at 7:33
• 2^28 as a float could really represent 2^28 exactly, not a range based on ulp.
– qwr
Commented Aug 17, 2023 at 14:08
• @qwr Of course it can, but it can as well be the result of a computation that mathematically would have produced something different, hopefully within +/- 0.5ulp, thus being a less-than-perfect representation of that desired number. From looking at the float number, you cannot tell. Commented Aug 17, 2023 at 14:22
• Well interval arithmetic is well-studied. I don't think it belongs in the C math library for many usecases with reasons I read in a presentation by Kahan, that often working with sufficient precision is enough and interval arithmetic can greatly overestimate error.
– qwr
Commented Aug 17, 2023 at 14:37

For a real-world comparison, glibc's sinf does range reduce, with different cases for small and large examples.

That means sinf(268435456.f) == sinf(4.5460540f), where the first number is 2^28 exactly and the second is 2^28 mod 2pi calculated with as much precision as a float can store.

/* Fast sinf implementation.  Worst-case ULP is 0.5607, maximum relative
error is 0.5303 * 2^-23.  A single-step range reduction is used for
small values.  Large inputs have their range reduced using fast integer
arithmetic.
*/

/* Fast range reduction using single multiply-subtract.  Return the modulo of
X as a value between -PI/4 and PI/4 and store the quadrant in NP.
The values for PI/2 and 2/PI are accessed via P.  Since PI/2 as a double
is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
the result is accurate for |X| <= 120.0.  */
static inline double
reduce_fast (double x, const sincos_t *p, int *np)
{
double r;
#if TOINT_INTRINSICS
/* Use fast round and lround instructions when available.  */
r = x * p->hpi_inv;
*np = converttoint (r);
return x - roundtoint (r) * p->hpi;
#else
/* Use scaled float to int conversion with explicit rounding.
hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
This avoids inaccuracies introduced by truncating negative values.  */
r = x * p->hpi_inv;
int n = ((int32_t)r + 0x800000) >> 24;
*np = n;
return x - n * p->hpi;
#endif
}


reduce_large: (I haven't benchmarked but this looks quite fast)

/* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
Reduction uses a table of 4/PI with 192 bits of precision.  A 32x96->128 bit
multiply computes the exact 2.62-bit fixed-point modulo.  Since the result
can have at most 29 leading zeros after the binary point, the double
precision result is accurate to 33 bits.  */
static inline double
reduce_large (uint32_t xi, int *np)
{
const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
int shift = (xi >> 23) & 7;
uint64_t n, res0, res1, res2;

xi = (xi & 0xffffff) | 0x800000;
xi <<= shift;

res0 = xi * arr[0];
res1 = (uint64_t)xi * arr[4];
res2 = (uint64_t)xi * arr[8];
res0 = (res2 >> 32) | (res0 << 32);
res0 += res1;

n = (res0 + (1ULL << 61)) >> 62;
res0 -= n << 62;
double x = (int64_t)res0;
*np = n;
return x * pi63;
}

• reduce_large() is quite intriguing. I have trouble finding the value of __inv_pio4[] in the link. Could you please post it here? Commented Aug 18, 2023 at 14:53
• @chux-ReinstateMonica github.com/bminor/glibc/blob/…
– qwr
Commented Aug 18, 2023 at 16:03

Degrees

Should the angle use degree units, reduce the argument to a narrower range. Common "mod" floating point functions are exact.

// d in degrees
d = fmod(d, 360);
y = sin(d * (2*M_PI/360));


Further reductions are possible with remquo().

Further study

The complexity for really good range reduction to $$[-\pi, \pi]$$ for large arguments is about as difficult as computing in-range trig functions. The is no need to have a cut-off other than a smaller code-footprint.

Argument Reduction for Huge Arguments: Good to the Last Bit is a useful read to form a high precision result for all $$x$$.

It effectively discusses how to do a high-precision $$\mathrm{mod}\,\pi$$ (and how much extra precision is needed for 64-bit FP types).

Note: it identifies $$x = 6{,}381{,}956{,}970{,}095{,}103 \cdot 2^{797}$$ as a worst case.

What should the cut-off point be?

Early trig function implementations introduced TLOSS (total loss of precision) and PLOSS (partial loss of precision). (Search these 2 keywords for many examples online. e.g. AIX.) They make sense if the angle is considered to represent a value ±1 ULP. This model follows OP's line of thinking. For binary32, TLOSS tends to get raised about $$|x| > 2^{24}$$ and PLOSS tends to get raised at maybe about $$|x| > 2^{12}$$.

More robust modern trig functions do not use the model of $$x$$ having error ("the result is unlikely to be meaningful") and provide the best sine(x) as if $$x$$ is exact. Thus no need for PLOSS nor TLOSS.

I recommend no cut-off.

I'll end with this rationale from Good to the Last Bit:

It is often argued that being concerned about large arguments is unnecessary, because sophisticated users simply know better than to compute with large angles. It is our contention that this position is suboptimal, because:

1. It places an unnecessary burden on the user.
2. The consequences of producing incorrect (inaccurate) answers may be catastrophic; many people assume that computers can do arithmetic very well. While numerical analysts know better, not all programmers are numerical analysts, nor should they be.
3. It is a vendors responsibility to provide answers that are as correct as possible.
• Perhaps the Standard could have recognized "general-prupose" trig functions and "super-precise" ones, and allow implementations to either make both symbols identify the same super-precise functions, or implement the general-purpose functions in a faster but less precise way, but bundle reference implementations of the more precise functions which woudln't need to worry about speed since they'd almost never be used anyway. Commented Dec 20, 2023 at 20:06
• @supercat <have ... general-prupose" trig functions and "super-precise" ones">, Hmmm. C can already do that with a call for the float version for lower precision and the double one for higher precision. The float version might be faster than the double one. Commented Dec 21, 2023 at 17:19
• A function that returns the properly rounded double-precision sine of any aribitrary-selected number that's within 0.5ulp of a double-precision argument could be faster and simpler than one which always returns the properly rounded single-precision sine of the exact value passed. Consider that there are precisely representable single-precision values like 16367173 times 2⁷³ whose sine could not be computed to even within +/-1ulp without performing computations using a value of pi accurate to well over 100 bits of precision. How much usefulness would such precision add? Commented Dec 21, 2023 at 18:10
• @supercat 1) Sounds like a good question to post here or another site. Comments here too restrictive. 2) Good to the Last Bit goes into that. 3) What is the point is which you would stop providing useful answers? That choice could be quite arbitrary. Commented Dec 21, 2023 at 20:50
• @supercat Curious why should a a sinf(float x) provide a poor answer for 16367173 times 2⁷³ ? Aside from the complexity of writing code, why should large x provide a poor answer, especially if primary range x results are fast? If code wants to avoid extra computation of a good sinf() for large x, it can do the fmod() itself and live with the weaker results as x` grows. Commented Dec 21, 2023 at 21:00

If you have a fused multiply-add available then you can do range reduction in a reasonably wide range: you usually reduce by a multiple of pi/2. So find the floating point number p2 closest to pi/2 and the floating point number p’ closest to pi/2 - p2. Let n = x / p2 rounded to the nearest integer, calculate x - n * p2 using a fused multiply add, subtract n * p’ from x. Repeat after x is small. You’ll need to be careful with the implementation but you’ll get perfect results up to say 2^53 and then precision goes down.

• I think this is what reduce_fast in my answer does?
– qwr
Commented Aug 18, 2023 at 16:13
• No. By using fused multiply-add you get about 106 bit of precision, so you get perfect results for x up to 2^53, not 120. You can't write x - a*b and assume you get a fused multiply-add or fused multiply-subtract instruction. Commented Aug 26, 2023 at 11:17