#!/usr/bin/env python3


# Truncation of Euler's Constant

def mainA():
    e_OEIS = 2.7182818284590452353
    N_max = 25
    N_range = np.arange(0, N_max)
    
    # a) Simple Relative Error Determination
    e_approx = np.empty(N_max)
    rel_err = np.empty(N_max)
    
    prev_e = 0
    for i, n in enumerate(N_range):
        e_approx[i] = prev_e + 1/np.math.factorial(n)
    
        rel_err[i] = abs( e_approx[i] - e_OEIS ) / e_OEIS
        prev_e = e_approx[i]
    

    ## Plot it all
    fig, ax = plt.subplots()
    ax.set_title("Relative Error for std float")
    ax.plot(N_range, rel_err)
    ax.set_yscale('log')
    ax.grid()
    ax.set_ylabel("Relative Error")
    ax.set_xlabel("N")
    
    plt.show()
    
def mainB():
    # b) Varying Floating Point Precision
    e_OEIS = 2.7182818284590452353
    N_max = 25
    N_range = np.arange(0, N_max)

    e_approx = np.zeros(N_max)
    rel_err = np.zeros(N_max)
    
    series_element_64 = np.empty(N_max, dtype=np.float64)
    series_element_32 = np.empty(N_max, dtype=np.float32)
    
    series64_diff = np.empty(N_max, dtype=np.float64)
    series32_diff = np.empty(N_max, dtype=np.float64)
    
    ## Determine the coefficients
    for i,n in enumerate(N_range):
        e_n = 1/np.math.factorial(n)
    
        series_element_64[i] = e_n # auto casts to float64
        series_element_32[i] = e_n # auto casts to float32
    
    ## Make a Cumulative Sum of the elements
    series64 = np.cumsum(series_element_64)
    series32 = np.cumsum(series_element_32)
    
    ## Determine the relative errror
    for i,n in enumerate(N_range):
        series64_diff[i] = abs( series64[i] - e_OEIS ) / e_OEIS
        series32_diff[i] = abs( series32[i] - e_OEIS ) / e_OEIS
    

    ## Plot it all
    fig, ax = plt.subplots()
    ax.set_title("Relative Error for float64 and float32")
    ax.plot(N_range, series64_diff, label="float64")
    ax.plot(N_range, series32_diff, label="float32")
    ax.set_yscale('log')
    ax.grid()
    ax.legend()
    ax.set_ylabel("Relative Error")
    ax.set_xlabel("N")
    
    plt.show()
    
    
def mainC():
    # c) Relative errors for different rounding accuracies
    round_accuracies = np.arange(1,6)

    e_OEIS = 2.7182818284590452353
    N_max = 25
    N_range = np.arange(0, N_max)
    
    # a) Simple Relative Error Determination
    e_approx = np.empty(N_max)
    
    rel_err = np.empty( (len(round_accuracies), N_max) )

    for i,d in enumerate(round_accuracies):
        prev_e = 0
        for j, n in enumerate(N_range):
            e_approx[i] = prev_e + round(1/np.math.factorial(n), d)
    
            rel_err[i,j] = abs( e_approx[i] - e_OEIS ) / e_OEIS
            prev_e = e_approx[i]
    

    ## Plot it all
    fig, ax = plt.subplots()
    ax.set_title("Relative Errors for std float with rounding at various digits d")
    for i,d in enumerate(round_accuracies):
        ax.plot(N_range, rel_err[i], label="d = {}".format(d))

    ax.set_yscale('log')
    ax.grid()
    ax.set_ylabel("Relative Error")
    ax.set_xlabel("N")
    ax.legend()
    
    plt.show()




if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt

    mainA()
    mainB()
    mainC()