How to handle periodic functions of very large floats

Question

(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?

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 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;
}

Answer 6

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:

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.

Answer 7

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.