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?

Gaskin

"round(cos(1000), 8) put it still returns a number written in scientific notation 1.2194074101485173e+225" It doesn't, actually. Scientific notation is just that... a notation. The result of round is a float, with a magnitude which is not decimal notation, or scientific notation, or any other. That comes as a result of how you print the value. – Zook

Look at the polynomial expansions of cosine. 1000/6.2... > 100, so your polynomial can't converge to the answer that far from zero. You need to take the argument modulo 2pi – Tinkle

You're worried about scientific notation, but entirely missing the fact that you have a number that's something times 10^255 – Tinkle

@Jared. I've updated my answer some more. – Tinkle

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.

Tinkle

When I add x %= 2.0 * math.pi to the beginning of function I get n -= ((x**c) / factorial(c)) OverflowError: int too large to convert to float – Gaskin

@Ironstone1_: you are taking way too many terms. factorial(172) is already larger than the largest float value. – Neoterize

Maybe Taylors series isn't the best method for calculating cos, I'll look for another. – Gaskin

@Ironstone1_. Most fast implementations use a bunch of hacks and something called Chebyshev Polynomials – Tinkle

@Ironstone1_. That being said, I've added a concrete implementation that's probably good enough for your immediate needs – Tinkle

Actually this happens with your original code when n > 85 also. – Crystallization

@YuriGinsburg. Well, (85*2))! ~= 7.26e306, while DBL_MAX ~= 1.80e308, so that's not surprising. Also, if you type "172 factorial" into Google, the result is "undefined". Heh. – Tinkle

@MadPhysicist Not sure what Google uses for computing factorial, but using Python I can easily get 100000! with no overflow. – Crystallization

@YuriGinsburg Google most likely uses IEEE754 double precision floats, while python supports infinite precision integers. The issue here isn't computing the factorial: it's computing 1/factorial as a float. – Tinkle

@YuriGinsburg. I've added a couple of examples to bypass the issue – Tinkle

@MadPhysicist This is good, the third function with the while loop, does appear to throw an error through if you enter cos(0) while abs(n / dn) > conv: ZeroDivisionError: float division by zero. – Gaskin

@MadPhysicist I understand the issue here ). But why Google uses double precision for factorial is still a mystery for me. – Crystallization

@YuriGinsburg. Because it was faster and easier to use floats for everything. Let's say that most people that are worried about factorials 86+ are not going to rely primarily on Google for that. – Tinkle

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.

Denominationalism

The function is not fine, and you can't magically turn 1.2*10^255 into a number between negative one and one. – Tinkle

Oops! I thought that read e-255 (hence, was a number quite close to zero). My mistake. – Denominationalism

Your answer makes a lot more sense now! – Tinkle

I'll leave my answer here, since I think the formatting information is still helpful, but yours is definitely the right answer to the question. – Denominationalism

I'll remove my downvote because I agree that the info is useful and I understand why you thought it was relevant now. – Tinkle