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?

  • The package you need is Decimal, and see the related [question](https://stackoverflow.com/questions/28081091/what-is-the-largest-number-the-decimal-class-can-handle) – Kate Melnykova Sep 18 '20 at 22:24
  • If I just import it, is that enough? what else do I need to do? –  Sep 18 '20 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: https://docs.python.org/3.8/library/decimal.html – Kate Melnykova Sep 18 '20 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. – Kate Melnykova Sep 18 '20 at 23:12
  • Besides the standard libary's `decimal` module, there's also the third-party module [`mpmath`](http://mpmath.org/) which you can (also) get and install via [pypi](https://pypi.org/project/mpmath/) that you might find easier to use. – martineau Sep 18 '20 at 23:59
  • Greats, thanks so much everyone, let me try the suggestions. –  Sep 19 '20 at 02:53
  • @martineau I tried mpmath, it sucks. Apparently their decimal places stop at 50. Maybe its just that I don't know how to use it, but when I tried I was let down that the result was wrong. Then I realized the precision was very poor, pi at 50 decimal places. Not sure how to specify a better precision since these guys design these packages for their own use, not somebody else's, and it only makes sense to them. –  Sep 19 '20 at 03:23
  • Decimal doesn't have the zeta function, its formulas are very limited, just the very basic ones, really just arithmetic. –  Sep 19 '20 at 03:36
  • I probably will use BigFloat, am trying to install this package now and it contains a lot of number theory formulas that I need. –  Sep 19 '20 at 03:43
  • Your experience with mpmath surprises me because its precision is supposed to be arbtrary — for example under **Features** on its web page says "Almost any calculation can be performed just as well at 10-digit or 1000-digit precision" — so it sounds like you might not have been using it properly. – martineau Sep 19 '20 at 05:56
  • Yes, it's possible, go ahead and give it a try if you can. The above formula requires a lot of precision, so if you get it right you should see a value ~4.95 for x=12. I'm guessing that either the package has a limitation, or that I don't know how to use. As I said, at times only the creators know, after all it's their baby and I have no idea what was on their mind when they wrote those things. That's why I came here to ask. –  Sep 19 '20 at 06:40
  • Even language at times can be imprecise, so I need them or someone to show me how to do it. But I guess that package might not really be the guy, it seems BigFloat and Gmpy2 are the guys (they use GNU libraries GMP and MPFR, for precise floating point arithmetic with correct rounding), but unfortunately I can't install them on my computer, Murphy's law. –  Sep 19 '20 at 06:40
  • @kate-melnykova Decimal doesn't work cause it doesn't have the zeta function pre-built, and it doesn't have Bernoulli numbers either. –  Sep 19 '20 at 06:46
  • See the code with decimal below. Decimal does not support zeta-function, but math package (and, I guess, scipy) support it. – Kate Melnykova Sep 19 '20 at 15:35
  • Thanks, I will try your solution shortly. Since I've spent a long time researching how to install GMP and MPFR for use with BigFloat, I thought I might as well finish it. Good to learn Python too at the same time. –  Sep 19 '20 at 22:35
  • @kate-melnykova Your code is giving me inf as answer, I'm puzzled. Should I increase the precision? I tested it in Colab. My bad, your sum is not converging. Now since you changed my original formula I can't tell the source of the problem. I know it works in Mathematica, I think my formula I posted above should be correct as well. –  Sep 19 '20 at 22:42

1 Answers1

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
Kate Melnykova
  • 1,863
  • 1
  • 5
  • 17