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 = nfor i in range(n-1, 0, -1):c *= ireturn cdef cos(x, i=100):c = 2n = 0for i in range(i):if i % 2 == 0:n += ((x**c) / factorial(c))else:n -= ((x**c) / factorial(c))c += 2return 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 mathdef cos(x, i=30):x %= 2 * math.pic = 2n = 0f = 2for i in range(i):if i % 2 == 0:n += x**c / felse:n -= x**c / fc += 2f *= 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 mathdef cos(x, i=30):x %= 2 * math.pin = 0dn = x**2 / 2for c in range(2, 2 * i + 2, 2):n += dndn *= -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 mathdef cos(x, conv=2**-53):x %= 2 * math.pic = 2n = 1.0dn = -x**2 / 2.0while abs(n / dn) > conv:n += dnc += 2dn *= -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.