Approximating cos using the Taylor series

Question

I'm using the Taylors series to calculate the cos of a number, with small numbers the function returns accurate results for example cos(5) gives 0.28366218546322663. But with larger numbers it returns inaccurate results such as cos(1000) gives 1.2194074101485173e+225

def factorial(n):
    c = n
    for i in range(n-1, 0, -1):
        c *= i
    return c

def cos(x, i=100):
    c = 2
    n = 0
    for i in range(i):
        if i % 2 == 0:
            n += ((x**c) / factorial(c))
        else:
            n -= ((x**c) / factorial(c))
        c += 2
    return 1 - n

I tried using round(cos(1000), 8) put it still returns a number written in scientific notation 1.2194074101485173e+225 with the e+ part. math.cos(1000) gives 0.5623790762907029, how can I round my numbers so they are the same as the math.cos method?

Answer 1

A McLaurin series uses Euler's ideas to approximate the value of a function using appropriate polynomials. The polynomials obviously diverge from a function like cos(x) because they all go towards infinity at some point, while cos doesn't. An order 100 polynomial can approximate at most 50 periods of the function on each side of zero. Since 50 * 2pi << 1000, your polynomial can't approximate cos(1000).

To get even close to a reasonable solution, the order of your polynomial must be at least x / pi. You can try to compute a polynomial of order 300+, but you're very likely to run into some major numerical issues because of the finite precision of floats and the enormity of factorials.

Instead, use the periodicity of cos(x) and add the following as the first line of your function:

x %= 2.0 * math.pi

You'll also want to limit the order of your polynomial to avoid problems with factorials that are too large to fit in a float. Furthermore, you can, and should compute your factorials by incrementing prior results instead of starting from scratch at every iteration. Here is a concrete example:

import math

def cos(x, i=30):
    x %= 2 * math.pi
    c = 2
    n = 0
    f = 2
    for i in range(i):
        if i % 2 == 0:
            n += x**c / f
        else:
            n -= x**c / f
        c += 2
        f *= c * (c - 1)
    return 1 - n

>>> print(cos(5), math.cos(5))
0.28366218546322663 0.28366218546322625

>>> print(cos(1000), math.cos(1000))
0.5623790762906707 0.5623790762907029

>>> print(cos(1000, i=86))
...
OverflowError: int too large to convert to float

You can further get away from numerical bottlenecks by noticing that the incremental product is x**2 / (c * (c - 1)). This is something that will remain well bounded for much larger i than you can support with a direct factorial:

import math

def cos(x, i=30):
    x %= 2 * math.pi
    n = 0
    dn = x**2 / 2
    for c in range(2, 2 * i + 2, 2):
        n += dn
        dn *= -x**2 / ((c + 1) * (c + 2))
    return 1 - n

>>> print(cos(5), math.cos(5))
0.28366218546322675 0.28366218546322625
>>> print(cos(1000), math.cos(1000))
0.5623790762906709 0.5623790762907029
>>> print(cos(1000, i=86), math.cos(1000))
0.5623790762906709 0.5623790762907029
>>> print(cos(1000, i=1000), math.cos(1000))
0.5623790762906709 0.5623790762907029

Notice that past a certain point, no matter how many loops you do, the result doesn't change. This is because now dn converges to zero, as Euler intended.

You can use this information to improve your loop even further. Since floats have finite precision (53 bits in the mantissa, to be specific), you can stop iteration when |dn / n| < 2**-53:

import math

def cos(x, conv=2**-53):
    x %= 2 * math.pi
    c = 2
    n = 1.0
    dn = -x**2 / 2.0
    while abs(n / dn) > conv:
        n += dn
        c += 2
        dn *= -x**2 / (c * (c - 1))
    return n

>>> print(cos2(5), math.cos(5))
0.28366218546322675 0.28366218546322625
>>> print(cos(1000), math.cos(1000))
0.5623790762906709 0.5623790762907029
>>> print(cos(1000, 1e-6), math.cos(1000))
0.5623792855306163 0.5623790762907029
>>> print(cos2(1000, 1e-100), math.cos(1000))
0.5623790762906709 0.5623790762907029

The parameter conv is not just the bound on |dn/n|. Since the following terms switch sign, it is also an upper bound on the overall precision of the result.

Answer 2

The number returned is simply a number; it has no sense of notation until you print it. If you're looking to control how the value is printed, let's suppose you're printing like so

print(cos(1000))

Then we can use format strings to control output

print("{:f}".format(cos(1000)))

If you're on Python 3.6 or newer, we can even interpolate it directly into the string literal.

print(f"{cos(1000):f}")

You can read the above links to see more details about the format mini-language (the language is the same between the two features). For instance, if you want to print a specific number of decimal places, you can request that too. We can print exactly three decimal places as follows

print("{:.3f}".format(cos(1000)))
print(f"{cos(1000):.3f}")

However, as Mad Physicist pointed out, there are some more mathematical problems with your code as well, so I strongly urge you to read his answer as well.