Division by Zero error in calculating series

Question

I am trying to compute a series, and I am running into an issue that I don't know why is occurring.

"RuntimeWarning: divide by zero encountered in double_scalars"

When I checked the code, it didn't seem to have any singularities, so I am confused. Here is the code currently(log stands for natural logarithm)(edit: extending code if that helps):

from numpy import pi, log


#Create functions to calculate the sums

def phi(z: int):
    k = 0
    phi = 0

    #Loop through 1000 times to try to approximate the series value as if it went to infinity

    while k <= 100:
        phi += ((1/(k+1)) - (1/(k+(2*z))))
        k += 1

    return phi

def psi(z: int):
    psi = 0
    k = 1

    while k <= 101:
        psi += ((log(k))/( k**(2*z)))
        k += 1
    
    return psi

def sig(z: int):
    sig = 0
    k = 1

    while k <= 101:
        sig += ((log(k))**2)/(k^(2*z))
        k += 1

    return sig

def beta(z: int):
    beta = 0
    k = 1

    while k <= 101:
        beta += (1/(((2*z)+k)^2))
        k += 1
    
    return beta


#Create the formula to approximate the value. For higher accuracy, either calculate more derivatives of Bernoulli numbers or increase the boundry of k.

def Bern(z :int):
    #Define Euler–Mascheroni constant

    c = 0.577215664901532860606512

    #Begin computations (only approximation)

    B = (pi/6) * (phi(1) - c - 2 * log(2 * pi) - 1) - z * ((pi/6) * ((phi(1)- c - (2 * log(2 * pi)) - 1) * (phi(1) - c) + beta(1) - 2 * psi(1)) - 2 * (psi(1) * (phi(1) - c) + sig(1) + 2 * psi(1) * log(2 * pi)))

    #output

    return B

A = int(input("Choose any value: "))
print("The answer is", Bern(A + 1))

Any help would be much appreciated.

Answer 1

are you sure you need a ^ bitwise exclusive or operator instead of **? I've tried to run your code with input parameter z = 1. And on a second iteration the result of k^(2*z) was equal to 0, so where is from zero division error come from (2^2*1 = 0).