0

I am not even sure if that is the best wording for what I'm trying to achieve, but I am trying to evaluate a power series in Python, and as I'd imagined, both because I'm new to Python and because power series evaluations are tricky when using floating-point numbers, my calculation is crashing.

My power series is an approximation for the simple prime-power counting function, J(x). In Mathematica I can evaluate this power series (I've done it up to prime 151, it's very time-consuming), but even in Mathematica, if I use floating-point in the formulas (instead of symbolic and then convert to decimal), the result is wrong. My plan is to create an executable that I can leave running on my job's remote desktop on Unix or Windows.

This is the code I created in Python (using scipy):

x = 12 M = 9 * x soma = 0 for i in range(1, M+1): termo1 = (-1)**i * x * (2*pi*x)**(2*i) / (2*i+1) for j in range(1, i+1): termo2 = (-1)**j * log(zeta(2*j))/((2*pi)**(2*j) * factorial(2*i+1-2*j)) soma = soma + termo1 * termo2 soma = -4 * soma

When we're new to something, the stress to figure out how to do even simple things can be tough (I'm a newbie so just saying use Decimal is greek to me). This is the error that I'm getting:

File "D:/iTunes/Python/PyCharm/Zeta.py", line 31, in <module> termo1 = (-1)**i * x * (2*pi*x)**(2*i) / (2*i+1) OverflowError: (34, 'Result too large')

How do I fix this? Or even better, is there a package that assumes that I need great precision in all the functions (zeta, log, factorial, etc.) and that spares me the trouble of having to figure things out myself?

Improve this question
16
  • The package you need is Decimal, and see the related question Commented Sep 18, 2020 at 22:24
  • If I just import it, is that enough? what else do I need to do? Commented Sep 18, 2020 at 23:06
  • 1
    No, you need to initialize your numbers as decimals directly in the code. See documentation on how to use the package: docs.python.org/3.8/library/decimal.html Commented Sep 18, 2020 at 23:08
  • You may also get a huge improvement in time if you optimize (-1)**i, (2*pi)**(2*j), factorials, etc. Roughly speaking, Python will run the Taylor series to compute each of these quantities for every iteration. However, it is easy to compute each of the quantities above if you knew their values for the previous iteration. Commented Sep 18, 2020 at 23:12
  • Besides the standard libary's decimal module, there's also the third-party module mpmath which you can (also) get and install via pypi that you might find easier to use. Commented Sep 18, 2020 at 23:59

1 Answer 1

Reset to default
0

This code works on Google Colab. Please update the precision as needed. Note that multiplication and division makes a few last digits to be inaccurate.

import math from math import pi, log, factorial from scipy.special import zeta from decimal import Decimal, getcontext getcontext().prec = 100 x = Decimal(12) M = int(9 * x) soma = Decimal(0) power_of_one = Decimal(-1) power_of_two_pi_x = (2 * Decimal(pi) * x) * (2 * Decimal(pi) * x) print(type(power_of_two_pi_x)) for i in range(1, M+1): termo1 = power_of_one * x * power_of_two_pi_x / (2*i+1) power_of_one *= -1 power_of_two_pi_x *= (2 * Decimal(pi) * x) * (2 * Decimal(pi) * x) power_of_one_j = -1 power_of_two_pi = Decimal(4 * pi * pi) factori = factorial(2 * i - 1) for j in range(1, i+1): termo2 = power_of_one_j * Decimal(log(zeta(2*j)))/power_of_two_pi * Decimal(factori) power_of_one_j *= -1 power_of_two_pi *= Decimal(4 * pi * pi) factori //= 2 * i + 1 - 2 * j if i > j: factori //= 2 * i - 2* j soma += termo1 * termo2 soma = -4 * soma
Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

I see, in line termo2=…, I didn't put factorial and power of 2pi in the denominator. Sorry.