At the end of my post on how not to calculate cosine I said that the first step in calculating cosine, particularly cosine of a large number, would be to do range reduction. This post will present a simple range reduction method by Cody and Waite that is adequate for moderately large arguments.

If you want to compute the sine or cosine of an angle x you could start by reducing x mod 2π since that would not change the result. However, accurately reducing a number mod 2π is not trivial; that’s why range reduction is an area of algorithm development.

Range reduction mod π/2

Even better would be to reduce x mod π/2. Reducing to a smaller range means that power series method, and other methods such as rational approximation, will be more efficient.

So suppose you can find an integer k such that

x − k π/2 = y

where 0 ≤ y ≤ π/2. Then sin(x) is ±sin(y) or ±cos(y), depending on k mod 4 equals 0, 1, 2, or 3.

from math import *

def reduced_sin(x, k):
     match k % 4:
        case 0: return sin(x)
        case 1: return cos(x)
        case 2: return -sin(x)
        case 3: return -cos(x)

Naive range reduction

Now let’s set x = 500. Then k = 318 because that’s the multiple of π/2 we need to subtract to bring x into range, and the sine of x should be the negative of the sine of the reduced value y because 318 = 2 mod 4.

The following code computes sin(x) with naive range reduction

def naive_sin(x):
    k = floor(x / (pi/2))
    y = x % (pi/2)
    return reduced_sin(y, k)

and when x = 500 the error is on the order of 1.7 × 10⁻¹⁴.

Better range reduction

The value of k above is fine, but we’d like to calculate y more accurately. The following code is much better.

def Cody_Waite_sin(x):
    C1 = 1686629713 / 2**30
    C2 = 4701928774853425 / 2**86

    k = floor(x / (pi/2))
    y = (x - k*C1) - k*C2
    return reduced_sin(y, k)

This will compute sin(500) to full machine precision. What kind of magic is this?

The trick is that the exact value of C1 + C2 equals π/2 to more precision than is possible in a single float [1]. You can confirm, with bc or some other extended precision software, that the difference between C1 + C2 and π/2 is roughly 2⁻⁸⁸, while the limit of float precision is 2⁻⁵².

If we compute

y = x - k*(C1 + C2)

then we’re doing the same calculation as naive_sin and will get the same error. But if we compute

y = (x - k*C1) - k*C2

we will get a more accurate result, provided x isn’t too large.

You can use the following code to play around and see how large x can be before errors start to creep in. For small enough x, like 500, the Cody and Waite sine returns full precision. For larger x it’s better than naive sine but does not return full precision. And for large enough x it completely breaks down.

def compare(x):
    y0 = naive_sin(x) 
    y1 = Cody_Waite_sin(x)
    y2 = sin(x)
    print("Naive error:     ", y2 - y0)
    print("Cody Waite error:", y2 - y1)

Now this may seem circular since we’re using math.sin as our gold standard. However, this function is calling the sine function on your CPU, which is using sophisticated range reduction to compute its result accurately down to the last bit, assuming you run the code on a computer that’s less than 40 years old.

The Cody and Waite algorithm is inadequate for large x, but it’s a good place to begin studying range reduction. It shows there are clever ways of squeezing out more precision than seems possible.

[1] The numerator n₁ of C1 is ⌊2³⁰ π/2⌋. The numerator n₂ of C2 is the solution to

2⁸⁶⁻³⁰ n₁ + n₂ = ⌊2⁸⁶ π/2⌋.

In a footnote to the previous post, I said that Python’s math library can calculate the logarithm of extremely large numbers but not the cosine. This post will expand on that comment.

In this post I’ll use n = 200! as my example rather than 1000! because this value of N is larger than the largest representable floating point number but small enough to be more convenient to work with.

Suppose someone calculates 200! for you:

78865786736479050355236321393218506229513597768717326329474253324435\
94499634033429203042840119846239041772121389196388302576427902426371\
05061926624952829931113462857270763317237396988943922445621451664240\
25403329186413122742829485327752424240757390324032125740557956866022\
60319041703240623517008587961789222227896237038973747200000000000000\
00000000000000000000000000000000000

You could now calculate log(n) using

n = 7.886578673647905 × 10³⁷⁴

and so

log(n) = log(7.886578673647905 × 10³⁷⁴)
= log(7.886578673647905) + 374 log(10) = 863.2319871924055.

The key thing that makes this possible is that the least significant digits of n only affect the least significant digits of log(n). In the calculation above I kept the first 16 digits of n. Python couldn’t make use of any more digits, and had no need of any more digits, in order to produce the logarithm to machine precision.

Cosine doesn’t work that way. The cosine of n depends on the remainder when n is divided by 2π, and that remainder depends on every single digit of n. I’ll illustrate that below.

Using bc -l and setting the scale to 400, I can calculated n then calculate

cos(n + 10ⁱ)

for i running from 0 to 374, tweaking each digit one at a time. (Except when a digit is a 9 and the addition results in a carry.)

    n = 1
    for (i = 1; i <= 200; i++) n *= i
    scale = 400
    for (i = 1; i <= 374; i++) {
        x = c(n+10^i)
        scale = 16
        print x/1, "\n"
        scale = 400
    }

Here’s what a plot of the results look like.

The value of cos(n) is about −0.985, but the values above are all over the map. We can look at the range by projecting all the points over to the left edge then rotating a quarter turn:

The remarkable thing about this image is that there are a few gaps, i.e. a few values the cosine does not take on.

Here’s a more sophisticated way to look at it. The sequence 10ⁱ mod 2π is dense in [0, 2π], and so by going far enough out in the sequence, we can find a value that shifts the phase of n by any desired amount within any given tolerance.

Every digit in n matters, and changing any digit can change the value of cosine to be essentially any value. You cannot calculate the cosine of an enormous number without using some kind of extended precision arithmetic. There are clever range reduction algorithms that minimize the amount of extended arithmetic necessary, but extended arithmetic cannot be completely eliminated.