I need to evaluate the finite sum of numbers which are increasing in absolute value, but are alternate. Problem is: the abolute values grow too fast and it starts accumulating numerical errors. These the functions definitions, one (Gj_full) straight to it and the other (Gj) recursively. fact_quo is a simple factorial function.
def fact_quo(n, m=1):
if (type(n) != int) or (type(m) != int):
raise TypeError("Arguments must be integers.")
if (n < 0) or (m < 0):
raise ValueError("n=" + str(n) + "\t m=" + str(m))
if (n == 0) or (n == 1) or (n == m):
return 1
if (n < m):
raise ValueError("n=" + str(n) + "\t m=" + str(m))
f = n
while n > (m+1):
n -= 1
f *= n
return f
def Gj_full(X, Y, Xl, Yl, j=0, coef=1):
if (X - Y + Xl + Yl) % 2 or X < Y or Y < j:
raise ValueError
u = (X - Y + Xl + Yl) // 2
v = coef * (2 ** (X - Y) * fact_quo(X, Y-j) * fact_quo(u+j, j) *
4 ** j * (-1) ** j)
w = 3 ** (u+j) * fact_quo(X-Y+j)
den2 = fact_quo(X) * fact_quo(Xl) * fact_quo(Yl)
z = (np.sqrt(fact_quo(X)) * np.sqrt(fact_quo(Y))
* np.sqrt(fact_quo(Xl)) * np.sqrt(fact_quo(Yl)))
return (v / (den2 * w)) * z
def Gj(X, Y, Xl, Yl, j=0, coef=1):
if (X - Y + Xl + Yl) % 2 or X < Y or Y < j:
raise ValueError
kX, kY, kXl, kYl, kj = X % 2, Y % 2, Xl % 2, Yl % 2, 0
g = coef * Gj_full(kX, kY, kXl, kYl, kj)
while kX < X:
u = (kX - kY + kXl + kYl) // 2
v = 4 * (u + kj + 1)
w = 3 * (kX + 2 - kY + kj) * (kX + 1 - kY + kj)
g *= (v / w) * np.sqrt(kX + 2) * np.sqrt(kX + 1)
kX += 2
while kXl < Xl:
u = (kX - kY + kXl + kYl) // 2
v = u + kj + 1
w = 3 * (kXl + 2) * (kXl + 1)
g *= (v / w) * np.sqrt(kXl + 2) * np.sqrt(kXl + 1)
kXl += 2
while kYl < Yl:
u = (kX - kY + kXl + kYl) // 2
v = u + kj + 1
w = 3 * (kYl + 2) * (kYl + 1)
g *= (v / w) * np.sqrt(kYl + 2) * np.sqrt(kYl + 1)
kYl += 2
while kY < Y:
u = (kX - kY + kXl + kYl) // 2
v = 3 * (kX - kY + kj) * (kX - kY - 1 + kj)
w = 4 * (kY + 2 - kj) * (kY + 1 - kj) * (u + kj)
g *= (v / w) * np.sqrt(kY + 2) * np.sqrt(kY + 1)
kY += 2
while kj < j:
u = (kX - kY + kXl + kYl) // 2
v = -4 * (kY - kj) * (u + kj + 1)
w = 3 * (kX - kY + kj + 1) * (kj + 1)
g *= (v / w)
kj += 1
return g
The (4/3) ** j and the factorials quicly increase the absolute value of the summing terms. The sum, however, are supposed to be smaller than 1. In fact, for X = Y and Xl = Yl = 0, the sum equals to (-1/3) ** X.
The precision for infinitely large numbers for floats are not available yet without using a lib. Therefore you should look into the
decimallib, you can even set the precision. Eg.If you manage to force all the numbers to be integers, you don't need to worry about it