Numerical implementation of the formula for the partial sum of a geometric series when ratio equals one. See https://math.stackexchange.com/questions/4288069/formula-for-the-partial-sum-of-a-geometric-series-when-ratio-equals-one

comparison_partial_sum_formulas.py

# Comparison of three numerical implementation of the formula for the partial sum of a geometric series
# when ratio equals (or is close to) one.
# See https://math.stackexchange.com/questions/4288069/formula-for-the-partial-sum-of-a-geometric-series-when-ratio-equals-one

from math import expm1, log1p

def SN_def(a,N):
    return sum(a**n for n in range(N))

def SN_form(a,N):
    return (1-a**N)/(1-a)

def SN_form2(eps,N):
    return -expm1(N*log1p(-eps))/eps

def SN_approx(eps,N):
    return N*(1-(N-1)/2*eps)
 
# Compare the implementations:

N = 100
eps_list = [10**(-i) for i in [1, 3, 5, 7, 9, 11, 13, 15, 17, 19]]
SN_fmt = lambda x: '{:20.15f}'.format(x)
txt_fmt = lambda s: '{:>20s}'.format(s)
sep = ', '

print(' '*7, end=sep)
print(txt_fmt('Definition with sum'), end=sep)
print(txt_fmt('Formula with expm1'), end=sep)
print(txt_fmt('Approx for small ε'), end=sep)
print(txt_fmt('Classical formula'))

for eps in eps_list:
    a = 1-eps
    print('ε={:.0e}'.format(eps), end=sep)
    print(SN_fmt(SN_def(a,N)), end=sep)
    print(SN_fmt(SN_form2(eps,N)), end=sep)
    print(SN_fmt(SN_approx(eps,N)), end=sep)
    try:
        print(SN_fmt(SN_form(a,N)))
    except ZeroDivisionError:
        print(txt_fmt('zero div. error'))