<p>(Or really just poor conditioning in general.)</p>
<p>Recently, there has been some interest in correctly-rounded floating-point libraries (<a href="https://github.com/rutgers-apl/The-RLIBM-Project" rel="noreferrer">RLIBM</a> and <a href="https://core-math.gitlabpages.inria.fr/" rel="noreferrer">CORE-MATH</a>), which I think is great, but there is one thing that gives me pause.  Consider <span class="math-container">$\sin(2^{28})$</span> in single-precision.  The period of <span class="math-container">$\sin$</span> is <span class="math-container">$2\pi$</span>, which is appreciably smaller than an ulp at <span class="math-container">$2^{28}$</span> (16), so the result is unlikely to be meaningful.  Actually, there are two issues:</p>
<ol>
<li><p>The implementation of <span class="math-container">$\sin$</span> must go to a lot of extra effort to implement range reduction of such large arguments.</p>
</li>
<li><p>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.</p>
</li>
</ol>
<p>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?</p>
<p>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?</p>

<floating-point>
How to handle periodic functions of very large floats
CC BY-SA 4.0
2023-08-15T23:22:34.157Z
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.
kaya3
CC BY-SA 4.0
2023-08-15T23:48:08.233Z
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.
Moonchild
CC BY-SA 4.0
2023-08-16T13:49:36.067Z
"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](http://bbc.nvg.org/doc/BBCUserGuide-1.00.pdf)).
psmears
CC BY-SA 4.0
2023-08-16T15:36:38.873Z
_How to handle?_ The answer probably differs for numerical experts and the "general public".
Pablo H
<p>Some points to add in addition to already good answers.</p>
<p><strong>Degrees</strong></p>
<p>Should the angle use degree units, reduce the argument to a narrower range.  Common &quot;mod&quot; floating point functions are <em>exact</em>.</p>
<pre><code>// d in degrees
d = fmod(d, 360);
y = sin(d * (2*M_PI/360));
</code></pre>
<p><a href="//stackoverflow.com/q/31502120">Further reductions</a> are possible with <code>remquo()</code>.</p>
<p><strong>Further study</strong></p>
<p>The complexity for really good range reduction to <span class="math-container">$[-\pi, \pi]$</span> 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.</p>
<p><a href="https://redirect.cs.umbc.edu/%7Ephatak/645/supl/Ng-ArgReduction.pdf" rel="nofollow noreferrer" title="ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit"><em>Argument Reduction for Huge Arguments: Good to the Last Bit</em></a> is a useful read to form a high precision result for all <span class="math-container">$x$</span>.</p>
<p>It effectively discusses how to do a high-precision <span class="math-container">$\mathrm{mod}\,\pi$</span> (and how much extra precision is needed for 64-bit FP types).</p>
<p>Note: it identifies <span class="math-container">$x = 6{,}381{,}956{,}970{,}095{,}103 \cdot 2^{797}$</span> as a worst case.</p>
<hr />
<p><strong>What should the cut-off point be?</strong></p>
<p>Early trig function implementations introduced <code>TLOSS</code> (total loss of precision) and <code>PLOSS</code> (partial loss of precision). (Search these 2 keywords for many examples online. e.g. <a href="https://sites.ualberta.ca/dept/chemeng/AIX-43/share/man/info/C/a_doc_lib/libs/basetrf2/sin.htm" rel="nofollow noreferrer">AIX</a>.) They make sense if the angle is considered to represent a value ±1 <a href="https://en.wikipedia.org/wiki/Unit_in_the_last_place" rel="nofollow noreferrer">ULP</a>.  This model follows OP's line of thinking.  For <a href="https://en.wikipedia.org/wiki/Single-precision_floating-point_format" rel="nofollow noreferrer">binary32</a>, <code>TLOSS</code> tends to get raised about <span class="math-container">$|x| &gt; 2^{24}$</span> and <code>PLOSS</code> tends to get raised at <em>maybe</em> about <span class="math-container">$|x| &gt; 2^{12}$</span>.</p>
<p>More robust modern trig functions do not use the model of <span class="math-container">$x$</span> having error (&quot;the result is unlikely to be meaningful&quot;) and provide the best <code>sine(x)</code> as if <span class="math-container">$x$</span> is exact.  Thus no need for <code>PLOSS</code> nor <code>TLOSS</code>.</p>
<p>I recommend no cut-off.</p>
<p>I'll end with this rationale from <a href="https://redirect.cs.umbc.edu/%7Ephatak/645/supl/Ng-ArgReduction.pdf" rel="nofollow noreferrer" title="ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit">Good to the Last Bit</a>:</p>
<blockquote>
<p>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:</p>
<ol>
<li>It places an unnecessary burden on the user.</li>
<li>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.</li>
<li>It is a vendors responsibility to provide answers that are as correct as possible.</li>
</ol>
</blockquote>

CC BY-SA 4.0
2023-08-17T22:08:06.907Z
2023-08-19T17:10:06.257Z
2023-08-19T17:10:06.257Z
chux
CC BY-SA 4.0
2023-12-20T20:06:32.260Z
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.
supercat
CC BY-SA 4.0
2023-12-21T17:19:17.717Z
@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.
chux
CC BY-SA 4.0
2023-12-21T18:10:51.750Z
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?
supercat
CC BY-SA 4.0
2023-12-21T20:50:38.773Z
@supercat 1) Sounds like a good question to post here or another site. Comments here too restrictive. 2) [Good to the Last Bit](http://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf) goes into that.  3) What is the point is which you would stop providing useful answers?  That choice could be quite arbitrary.
chux
CC BY-SA 4.0
2023-12-21T21:00:12.110Z
@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.
chux
CC BY-SA 4.0
2023-12-22T23:05:56.690Z
Even if `sinf` is never invoked with values larger than pi/2, an embedded program that calls that function would have no alternative but to include the entire `sinf` routine in ROM, using up code space that might be more usefully employed for other purposes.  Using `fmod` would do nothing to improve that situation, and might in fact make it worse if the implementation would have no other use for `fmod` [a general-purpose implementation would likely need to link `fmod`, but an implementation with specialized range reduction may--as its one consolation--not need to link that function).
supercat
CC BY-SA 4.0
2023-12-23T07:30:59.953Z
@supercat  I disagree that the code space concern is valid for in a tight code space situation as one can readily code small `my_sin()` with only a few ULP error in the primary range instead of a lengthy high quality one.  I look forward to seeing your posted answer to OP's question where you can well point out its benefits.
chux
<p>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.</p>

CC BY-SA 4.0
2023-08-17T22:51:40.490Z
2023-08-17T22:51:40.490Z
gnasher729
CC BY-SA 4.0
2023-08-18T16:13:17.967Z
I think this is what reduce_fast in my answer does?
qwr
CC BY-SA 4.0
2023-08-26T11:17:27.610Z
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.
gnasher729
<p>I think this is in the category of &quot;probable mistakes&quot; ─ 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.</p>
<p>For example, it's plausible that a number theorist might actually want to compute the value of something like <span class="math-container">$\sum_{k=0}^{28} \sin 2^k$</span>; 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 &quot;loss of precision&quot; due to the scale of the argument: <span class="math-container">$2^{28}$</span> is the precise value that the programmer wants the sin of, and it is exactly representable as a single-precision float.</p>
<p>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 <code>sin_unchecked</code> function).</p>
<hr />
<p>For comparison, consider the following C or Java code, which is almost certainly a mistake:</p>
<pre class="lang-c prettyprint-override"><code>while(do_thing());
    do_other_thing();
</code></pre>
<p>The semicolon after the <code>while</code> condition is a statement that does nothing, which means the loop body does nothing and <code>do_other_thing()</code> is not part of the loop. This is a fairly common mistake in real code. But it's also not <em>absolutely</em> certain that it's a mistake; it's plausible to write <code>while</code> loops with empty bodies when the condition is the thing you want to loop, and it's plausible to have incorrect indentation.</p>
<p>This is what I mean by the category of &quot;probable mistakes&quot;. 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 <code>while</code> statement, so that an empty <code>;</code> statement is not a valid body.</p>
<p>But because it <em>could</em> be intentional to have a loop with an empty body, the language should still have a way for the programmer to say &quot;actually, I do want the loop body to be an empty statement&quot;. In this example, that's easily satisfied by allowing <code>{}</code> to be the empty statement instead of <code>;</code>. But it's something you should consider, if you're going to make something an error.</p>

CC BY-SA 4.0
2023-08-16T13:42:24.263Z
2023-08-16T19:27:13.140Z
2023-08-16T19:27:13.140Z
kaya3
CC BY-SA 4.0
2023-08-16T19:06:22.163Z
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.
Moonchild
CC BY-SA 4.0
2023-08-16T19:11:14.097Z
@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.
kaya3
CC BY-SA 4.0
2023-08-16T19:14:18.413Z
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.
kaya3
CC BY-SA 4.0
2023-08-16T19:18:01.967Z
You already have to pay for that extra check, because large floats are range-reduced differently.  And there is already the possibility of an exception, in the case when you pass in infinity.
Moonchild
CC BY-SA 4.0
2023-08-16T19:25:17.873Z
@Moonchild `x % TWO_PI` gives correctly-rounded results on floats regardless of their scale, and `Math.sin(Double.POSITIVE_INFINITY)` doesn't throw an exception in Java for example (it returns NaN).
kaya3
CC BY-SA 4.0
2023-08-16T19:30:52.957Z
In which language does x % TWO_PI work?  And how is that implemented?  You can see for example special code in rlibm [here](https://github.com/rutgers-apl/The-RLIBM-Project/blob/main/libm/sinf.c#L97) which needs 256 bits of 1/pi for large arguments, but much less for smaller arguments.  That's just for single floats (which not-so-coincidentally have an exponent range of ~256); double floats can get much larger.
Moonchild
CC BY-SA 4.0
2023-08-16T20:29:36.640Z
@Moonchild The remainder operation is one of the [required operations](https://en.wikipedia.org/wiki/IEEE_754#Required_operations) in the IEEE 754 spec. I wouldn't be able to tell you how it's implemented, though. I'm not sure I know of any languages which *don't* have it; in C it's the function `fmod` in `math.h`, in other languages (Java, Javascript, Python, Rust...) it's the `%` operator.
kaya3
CC BY-SA 4.0
2023-08-16T21:49:18.493Z
@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.
supercat
CC BY-SA 4.0
2023-08-16T21:56:07.283Z
@supercat That's a very good point; so it seems there is no extra cost for checking large numbers to throw an exception. But it's still doubtful that an exception is desirable, particularly for applications which do many such calculations and don't mind if a few of them are meaningless due to imprecision or domain errors.
kaya3
CC BY-SA 4.0
2023-08-16T22:05:45.607Z
@kaya3: When performing sinpi(x), the range reduction can be formed by adding or subtracting to/from x a value that is precisely representable if x is likewise.  If x is so big that the difference between representable values exceeds 2, the range reduction will yield precisely zero.  When computing sin(x), the amount to be added or subtracted to perform range reduction accurate to within 0.5ulp of the result would need precision far beyond that of x.
supercat
CC BY-SA 4.0
2023-08-17T15:36:50.667Z
@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).
supercat
CC BY-SA 4.0
2023-08-17T15:44:44.980Z
yeah. There's a reason `sinpi` is recommended by recent IEEE versions.
Oscar Smith
<p>There's two major challenges that I see with this approach:</p>
<ul>
<li>The correct tolerance for &quot;meaninglessness&quot; is dependent on the use case</li>
<li>Exceptional flow control is expensive from a user's perspective.</li>
</ul>
<p>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 <em>major</em> 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.</p>
<p>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 &quot;popular&quot; to not raise exceptions for numerical issues which are worse than <span class="math-container">$\sin(2^{38})$</span>.  That's not to say its &quot;right,&quot; 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.</p>

CC BY-SA 4.0
2023-08-17T02:35:29.310Z
2023-08-17T02:35:29.310Z
Cort Ammon
CC BY-SA 4.0
2023-08-18T01:00:54.627Z
@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.
Cort Ammon
CC BY-SA 4.0
2023-08-17T06:37:43.607Z
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).
Moonchild
CC BY-SA 4.0
2023-08-17T08:21:37.653Z
"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".
Ralf Kleberhoff
<p>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).</p>
<p>So, instead of ignoring that &quot;interval&quot; property, as our CPUs and libraries typically do, a more math-oriented approach would be to embrace it, and have functions like <code>sin()</code> accept and return a range value. Then e.g. <code>sin(1e30)</code> (with <code>1e30</code> denoting a range based on the IEEE double precision ULP) would give a range of <code>[-1,1]</code>, thus telling the caller that the sine of the mathematical real number that got approximated to the <code>double</code> value can be anything.</p>
<p>Of course, that needs a whole new approach to math operations, so it might be too much for your project.</p>
<p>An approach better compatible with current libraries and CPUs might be to have an additional <code>requiredPrecision</code> argument to functions like <code>sin()</code>. Then</p>
<pre><code>y = sin(x, 0.0001)
</code></pre>
<p>would either return a sine value with the guarantee that for any real number approximated by <code>x</code>, its mathematical sine value will be in the interval <code>[y-0.0001, y+0.0001]</code>. If the <code>sin()</code> function can't guarantee that, an exception would be appropriate. For the <code>sin()</code> function, a sufficient test can be to compare the ULP with the <code>requiredPrecision</code> argument, as the sine derivative is always in the range from -1 to 1.</p>
<p>This way, callers can express their precision requirements and get exceptions when the function can't fulfill that.</p>
<p>Without the callers specifying their requirements, any threshold you select will be arbitrary.</p>

CC BY-SA 4.0
2023-08-16T07:11:54.413Z
2023-08-16T07:28:52.193Z
2023-08-16T07:28:52.193Z
Ralf Kleberhoff
CC BY-SA 4.0
2023-12-20T20:02:13.987Z
@qwr: What's unfortunate is that C's treatment of variadic functions discourages implementations from supporting sufficient precision in the manner advocated by Kahan.  Having programs perform sequences of operations in extended precision and rounding *at well defined points* could often be more accurate *and* more efficient both in hardware that's designed around that principle *and in implementations without hardware floating-point support*.  The concept wasn't simply adopted based upon the particular design of the 8087.
supercat
CC BY-SA 4.0
2023-08-16T07:28:48.857Z
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).
Moonchild
CC BY-SA 4.0
2023-08-16T07:33:49.697Z
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.
Moonchild
CC BY-SA 4.0
2023-08-16T15:07:32.963Z
@Moonchild: In what non-contrived scenarios would a function that is guaranteed to precisely compute sin(x) within 0.5ulp be more useful than one which is guaranteed to precisely compute sin(x') for some x' within 0.5ulp of x?
supercat
CC BY-SA 4.0
2023-08-16T19:02:05.317Z
I can't think of any (aside from reproducibility), but I don't see what the significance is.
Moonchild
CC BY-SA 4.0
2023-08-17T14:08:56.970Z
2^28 as a float could really represent 2^28 exactly, not a range based on ulp.
qwr
CC BY-SA 4.0
2023-08-17T14:22:07.617Z
@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.
Ralf Kleberhoff
CC BY-SA 4.0
2023-08-17T14:37:59.000Z
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
<p>For a real-world comparison, glibc's <a href="https://github.com/bminor/glibc/blob/master/sysdeps/ieee754/flt-32/s_sinf.c" rel="nofollow noreferrer"><code>sinf</code></a> does range reduce, with different cases for small and large examples.</p>
<p>That means <code>sinf(268435456.f) == sinf(4.5460540f)</code>, 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.</p>
<pre class="lang-c prettyprint-override"><code>/* 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.
*/
</code></pre>
<p><a href="https://github.com/bminor/glibc/blob/4290aed05135ae4c0272006442d147f2155e70d7/sysdeps/ieee754/flt-32/s_sincosf.h#L48" rel="nofollow noreferrer"><code>reduce_fast</code></a>:</p>
<pre class="lang-c prettyprint-override"><code>/* 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| &lt;= 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-&gt;hpi_inv;
  *np = converttoint (r);
  return x - roundtoint (r) * p-&gt;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-&gt;hpi_inv;
  int n = ((int32_t)r + 0x800000) &gt;&gt; 24;
  *np = n;
  return x - n * p-&gt;hpi;
#endif
}
</code></pre>
<p><a href="https://github.com/bminor/glibc/blob/4290aed05135ae4c0272006442d147f2155e70d7/sysdeps/ieee754/flt-32/s_sincosf.h#L75" rel="nofollow noreferrer"><code>reduce_large</code></a>: (I haven't benchmarked but this looks quite fast)</p>
<pre class="lang-c prettyprint-override"><code>/* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
   XI is a reinterpreted float and must be &gt;= 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-&gt;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 = &amp;__inv_pio4[(xi &gt;&gt; 26) &amp; 15];
  int shift = (xi &gt;&gt; 23) &amp; 7;
  uint64_t n, res0, res1, res2;

  xi = (xi &amp; 0xffffff) | 0x800000;
  xi &lt;&lt;= shift;

  res0 = xi * arr[0];
  res1 = (uint64_t)xi * arr[4];
  res2 = (uint64_t)xi * arr[8];
  res0 = (res2 &gt;&gt; 32) | (res0 &lt;&lt; 32);
  res0 += res1;

  n = (res0 + (1ULL &lt;&lt; 61)) &gt;&gt; 62;
  res0 -= n &lt;&lt; 62;
  double x = (int64_t)res0;
  *np = n;
  return x * pi63;
}
</code></pre>

CC BY-SA 4.0
2023-08-16T20:26:01.013Z
2023-08-16T21:07:27.777Z
2023-08-16T21:07:27.777Z
qwr
CC BY-SA 4.0
2023-08-18T14:53:36.783Z
`reduce_large()` is quite intriguing.  I have trouble finding the value of `__inv_pio4[]` in the link.  Could you please post it here?
chux
CC BY-SA 4.0
2023-08-18T16:03:28.907Z
@chux-ReinstateMonica https://github.com/bminor/glibc/blob/4290aed05135ae4c0272006442d147f2155e70d7/sysdeps/ieee754/flt-32/s_sincosf_data.c#L66
qwr
<p>If you claim IEEE-754 compliance, the code to handle <span class="math-container">$\sin(2^{38})$</span> 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 <span class="math-container">$\sin(2^{38})$</span> deserve the performance they get.</p>
<p>Recent editions of the IEEE-754 standard also define functions such as <code>sinPi</code>, <code>cosPi</code>, <code>tanPi</code>, and their inverse versions (<code>atan2Pi</code> has an unfortunate name!), which are defined as, say:</p>
<p><span class="math-container">$$\mathrm{sinPi}(x) = \sin(\pi x)$$</span></p>
<p>The theory is that range reduction to the range <span class="math-container">$[-1,1]$</span> is an extremely inexpensive operation, and is <em>always</em> exact.</p>

CC BY-SA 4.0
2023-08-16T03:34:04.260Z
2023-08-16T03:34:04.260Z
Pseudonym
CC BY-SA 4.0
2023-08-17T22:16:33.420Z
"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`.
chux
CC BY-SA 4.0
2023-08-18T00:43:23.670Z
@chux-ReinstateMonica Of course, it does, assuming that angles start in radians. That's sometimes the case and sometimes not. It's usually true in science, but it's not as common in engineering as you might think. I currently work in GIS, and we still use degrees. Those working in robotics and games often think in terms of "fractions of a turn". People working with Fourier series often have to evaluate things like $\int_{-L}^{L} f(x)\cos\left( \frac{n\pi x}{L} \right) \, dx$.
Pseudonym
CC BY-SA 4.0
2023-08-18T00:52:07.040Z
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.
Pseudonym
CC BY-SA 4.0
2023-12-20T19:58:55.220Z
@Pseudonym: Even in scientific applications, I find the claim dubious.  One might work with an angle that is a precise integer multiple of 1/n radians for some integer n, but unless n is a power of two, computing `sin(t * (1/n))` isn't apt to be meaningfully more precise than computing `sin2pi(t * (1/2πn))` .  In both cases, precision may be enhanced by performing argument reduction *before* the multiply.  The fact that code is computing the sine of an angle much larger than 2π would generally indicate that the author forfeited a chance to improve precision.
supercat
CC BY-SA 4.0
2023-12-20T21:42:00.797Z
@supercat If you were only ever evaluating sine and cosine, yes. But in scientific applications, you routinely work with functions, their analytic derivatives and their (sometimes partial) analytic integrals, so the factor of some multiple of $\pi$ always comes back.
Pseudonym
CC BY-SA 4.0
2023-12-21T02:45:43.233Z
@Pseudonym: If something has a period of e.g. 2π seconds, and a simulation is using a step size of 1 microsecond, then the period of the original signal, its first derivative, second derivative, etc. will all be 2000000π seconds, meaning one can compute phase angle as either `(step * (1/1000000.0)` radians, or `step * (1/ 2000000π)`  circles.
supercat
CC BY-SA 4.0
2023-12-21T07:07:22.447Z
@supercat I'm very confused by this, and perhaps we should take this to chat. My point about derivatives etc is merely that the derivative of $\cos x$ is $-\sin x$, but the derivative of $\mathrm{cosPi}\,x$ is $-\pi\,\mathrm{sinPi}\,x$. If all you're doing is calculating rotations or a pure periodic function, fine. But if you need to do anything that resembles analysis or closed-form calculus, working in terms of `cosPi` and `sinPi` doesn't necessarily fix anything.
Pseudonym
CC BY-SA 4.0
2023-08-16T05:19:16.663Z
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.
Moonchild
CC BY-SA 4.0
2023-08-16T21:01:03.383Z
@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.
supercat
CC BY-SA 4.0
2023-08-16T23:36:38.097Z
Returning the fractional part of an IEEE-754 floating point number is an exact operation. It requires no rounding.
Pseudonym
# How to handle periodic functions of very large floats
Posted by: username_0 on 2023-08-15
(Or really just poor conditioning in general.)

Recently, there has been some interest in correctly-rounded floating-point libraries ([RLIBM](https://github.com/rutgers-apl/The-RLIBM-Project) and [CORE-MATH](https://core-math.gitlabpages.inria.fr/)), 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?
**Tags:** floating-point
> **Comments:**
> **username_1**: 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.  
> **username_0**: 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.  
> **username_2**: "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](http://bbc.nvg.org/doc/BBCUserGuide-1.00.pdf)).  
> **username_3**: _How to handle?_ The answer probably differs for numerical experts and the "general public".  
---
## 7 answers
### Answer by username_1 on 2023-08-16. Score: 13
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.
> **Comments:**
> **username_0** (2023-08-16): 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.  
> **username_1** (2023-08-16): @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.  
> **username_1** (2023-08-16): 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.  
> **username_0** (2023-08-16): You already have to pay for that extra check, because large floats are range-reduced differently.  And there is already the possibility of an exception, in the case when you pass in infinity.  
> **username_1** (2023-08-16): @Moonchild `x % TWO_PI` gives correctly-rounded results on floats regardless of their scale, and `Math.sin(Double.POSITIVE_INFINITY)` doesn't throw an exception in Java for example (it returns NaN).  
> **username_0** (2023-08-16): In which language does x % TWO_PI work?  And how is that implemented?  You can see for example special code in rlibm [here](https://github.com/rutgers-apl/The-RLIBM-Project/blob/main/libm/sinf.c#L97) which needs 256 bits of 1/pi for large arguments, but much less for smaller arguments.  That's just for single floats (which not-so-coincidentally have an exponent range of ~256); double floats can get much larger.  
> **username_1** (2023-08-16): @Moonchild The remainder operation is one of the [required operations](https://en.wikipedia.org/wiki/IEEE_754#Required_operations) in the IEEE 754 spec. I wouldn't be able to tell you how it's implemented, though. I'm not sure I know of any languages which *don't* have it; in C it's the function `fmod` in `math.h`, in other languages (Java, Javascript, Python, Rust...) it's the `%` operator.  
> **username_5** (2023-08-16): @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.  
> **username_1** (2023-08-16): @supercat That's a very good point; so it seems there is no extra cost for checking large numbers to throw an exception. But it's still doubtful that an exception is desirable, particularly for applications which do many such calculations and don't mind if a few of them are meaningless due to imprecision or domain errors.  
> **username_5** (2023-08-16): @kaya3: When performing sinpi(x), the range reduction can be formed by adding or subtracting to/from x a value that is precisely representable if x is likewise.  If x is so big that the difference between representable values exceeds 2, the range reduction will yield precisely zero.  When computing sin(x), the amount to be added or subtracted to perform range reduction accurate to within 0.5ulp of the result would need precision far beyond that of x.  
> **username_5** (2023-08-17): @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).  
> **username_8** (2023-08-17): yeah. There's a reason `sinpi` is recommended by recent IEEE versions.  
---
### Answer by username_9 on 2023-08-17. Score: 12
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.
> **Comments:**
> **username_0** (2023-08-17): 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).  
> **username_10** (2023-08-17): "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".  
> **username_9** (2023-08-18): @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.  
---
### Answer by username_11 on 2023-08-16. Score: 8
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.
> **Comments:**
> **username_0** (2023-08-16): 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.  
> **username_5** (2023-08-16): @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.  
> **username_11** (2023-08-16): Returning the fractional part of an IEEE-754 floating point number is an exact operation. It requires no rounding.  
> **username_4** (2023-08-17): "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`.  
> **username_11** (2023-08-18): @chux-ReinstateMonica Of course, it does, assuming that angles start in radians. That's sometimes the case and sometimes not. It's usually true in science, but it's not as common in engineering as you might think. I currently work in GIS, and we still use degrees. Those working in robotics and games often think in terms of "fractions of a turn". People working with Fourier series often have to evaluate things like $\int_{-L}^{L} f(x)\cos\left( \frac{n\pi x}{L} \right) \, dx$.  
> **username_11** (2023-08-18): 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.  
> **username_5** (2023-12-20): @Pseudonym: Even in scientific applications, I find the claim dubious.  One might work with an angle that is a precise integer multiple of 1/n radians for some integer n, but unless n is a power of two, computing `sin(t * (1/n))` isn't apt to be meaningfully more precise than computing `sin2pi(t * (1/2πn))` .  In both cases, precision may be enhanced by performing argument reduction *before* the multiply.  The fact that code is computing the sine of an angle much larger than 2π would generally indicate that the author forfeited a chance to improve precision.  
> **username_11** (2023-12-20): @supercat If you were only ever evaluating sine and cosine, yes. But in scientific applications, you routinely work with functions, their analytic derivatives and their (sometimes partial) analytic integrals, so the factor of some multiple of $\pi$ always comes back.  
> **username_5** (2023-12-21): @Pseudonym: If something has a period of e.g. 2π seconds, and a simulation is using a step size of 1 microsecond, then the period of the original signal, its first derivative, second derivative, etc. will all be 2000000π seconds, meaning one can compute phase angle as either `(step * (1/1000000.0)` radians, or `step * (1/ 2000000π)`  circles.  
> **username_11** (2023-12-21): @supercat I'm very confused by this, and perhaps we should take this to chat. My point about derivatives etc is merely that the derivative of $\cos x$ is $-\sin x$, but the derivative of $\mathrm{cosPi}\,x$ is $-\pi\,\mathrm{sinPi}\,x$. If all you're doing is calculating rotations or a pure periodic function, fine. But if you need to do anything that resembles analysis or closed-form calculus, working in terms of `cosPi` and `sinPi` doesn't necessarily fix anything.  
---
### Answer by username_10 on 2023-08-16. Score: 6
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.
> **Comments:**
> **username_0** (2023-08-16): 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).  
> **username_0** (2023-08-16): 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.  
> **username_5** (2023-08-16): @Moonchild: In what non-contrived scenarios would a function that is guaranteed to precisely compute sin(x) within 0.5ulp be more useful than one which is guaranteed to precisely compute sin(x') for some x' within 0.5ulp of x?  
> **username_0** (2023-08-16): I can't think of any (aside from reproducibility), but I don't see what the significance is.  
> **username_7** (2023-08-17): 2^28 as a float could really represent 2^28 exactly, not a range based on ulp.  
> **username_10** (2023-08-17): @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.  
> **username_7** (2023-08-17): 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.  
> **username_5** (2023-12-20): @qwr: What's unfortunate is that C's treatment of variadic functions discourages implementations from supporting sufficient precision in the manner advocated by Kahan.  Having programs perform sequences of operations in extended precision and rounding *at well defined points* could often be more accurate *and* more efficient both in hardware that's designed around that principle *and in implementations without hardware floating-point support*.  The concept wasn't simply adopted based upon the particular design of the 8087.  
---
### Answer by username_7 on 2023-08-16. Score: 5
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.
*/

reduce_fast:
/* 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;
}
> **Comments:**
> **username_4** (2023-08-18): `reduce_large()` is quite intriguing.  I have trouble finding the value of `__inv_pio4[]` in the link.  Could you please post it here?  
> **username_7** (2023-08-18): @chux-ReinstateMonica https://github.com/bminor/glibc/blob/4290aed05135ae4c0272006442d147f2155e70d7/sysdeps/ieee754/flt-32/s_sincosf_data.c#L66  
---
### Answer by username_4 on 2023-08-17. Score: 3
Some points to add in addition to already good answers.
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:

It places an unnecessary burden on the user.
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.
It is a vendors responsibility to provide answers that are as correct as possible.
> **Comments:**
> **username_5** (2023-12-20): 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.  
> **username_4** (2023-12-21): @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.  
> **username_5** (2023-12-21): 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?  
> **username_4** (2023-12-21): @supercat 1) Sounds like a good question to post here or another site. Comments here too restrictive. 2) [Good to the Last Bit](http://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf) goes into that.  3) What is the point is which you would stop providing useful answers?  That choice could be quite arbitrary.  
> **username_4** (2023-12-21): @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.  
> **username_5** (2023-12-22): Even if `sinf` is never invoked with values larger than pi/2, an embedded program that calls that function would have no alternative but to include the entire `sinf` routine in ROM, using up code space that might be more usefully employed for other purposes.  Using `fmod` would do nothing to improve that situation, and might in fact make it worse if the implementation would have no other use for `fmod` [a general-purpose implementation would likely need to link `fmod`, but an implementation with specialized range reduction may--as its one consolation--not need to link that function).  
> **username_4** (2023-12-23): @supercat  I disagree that the code space concern is valid for in a tight code space situation as one can readily code small `my_sin()` with only a few ULP error in the primary range instead of a lengthy high quality one.  I look forward to seeing your posted answer to OP's question where you can well point out its benefits.  
---
### Answer by username_6 on 2023-08-17. Score: 1
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.
> **Comments:**
> **username_7** (2023-08-18): I think this is what reduce_fast in my answer does?  
> **username_6** (2023-08-26): 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.  
---
