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