0

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.

1 Answers1

0

The precision for infinitely large numbers for floats are not available yet without using a lib. Therefore you should look into the decimal lib, you can even set the precision. Eg.

import decimal
decimal.getcontext().prec = 100
def pi():
    pi = decimal.Decimal(0)
    for k in range(350):
        pi += (decimal.Decimal(4)/(decimal.Decimal(8)*decimal.Decimal(k+1))...)

If you manage to force all the numbers to be integers, you don't need to worry about it

Adriano Martins
  • 1,788
  • 1
  • 19
  • 23
  • Hi Adriano, I am going to check the mpmath package, but I was trying to avoid these "brute force" solutions. In fact I was trying to learn a bit of numerical computing thinking. This is can be an art sometimes. – Mauricio Reis Mar 21 '18 at 02:52
  • Those are not _bruteforce_ alternatives, those are the only ways to work with a high precision on python. Also, I would recommend you to take a look on [NumPy](http://www.numpy.org/) - it is by far the best one out there. – Adriano Martins Mar 21 '18 at 18:12