import random
from math import sqrt
import numpy as np

def harmonic_series(stop_after=-1):
    """

    :param stop_after:
    :return:
    """
    sum = np.longdouble(0.0)
    k = 1
    while True:
        sum += 1/k
        yield sum
        if k == stop_after:
            break
        k += 1

def alternating_harmonic_series(stop_after=-1):
    """

    :param stop_after:
    :return:
    """
    sum = np.longdouble(0.0)
    k = 1
    negative = False
    while True:
        if negative:
            sum += (-1)/k
        else:
            sum += 1/k
        yield sum
        if k == stop_after:
            break
        k += 1
        negative = not negative

def general_harmonic_series(a, b, alternating=False, stop_after=-1):
    """

    :param a:
    :param b:
    :param alternating:
    :param stop_after:
    :return:
    """
    assert a != 0
    assert b/a > 0.0
    sum = np.longdouble(0.0)
    n = 0
    negative = False
    while True:
        if negative:
            sum += (-1) / ((a*n) + b)
        else:
            sum += 1 / ((a*n) + b)
        yield sum
        if n == stop_after:
            break
        if alternating:
            negative = not negative
        n += 1

def leibniz_series(stop_after=-1):
    """
    Approaches pi/4
    :param stop_after: if specified, stop after n-th element of series
    :return: leibniz series generator object
    """
    sum = np.longdouble(0.0)
    k = 0
    negative = False
    while True:
        if negative:
            sum += (-1) / ((2*k)+1)
        else:
            sum += 1 / ((2*k)+1)
        yield sum
        if k == stop_after:
            break
        k += 1
        negative = not negative

def p_series(p, stop_after=-1):
    """

        :param stop_after:
        :return:
    """
    sum = np.longdouble(0.0)
    k = 1
    while True:
        sum += 1 / (k ** p)
        yield sum
        if k == stop_after:
            break
        k += 1


def random_harmonic_series(stop_after=-1):
    """

    :param stop_after:
    :return:
    """
    sum = np.longdouble(0.0)
    k = 1
    while True:
        sum += random.choice([-1, 1]) / k
        yield sum
        if k == stop_after:
            break
        k += 1


def euler_accelerator(series):
    """
    Lets a series converge faster
    :param series:
    :return:
    """
    s0 = series.__next__() # Sn-1
    s1 = series.__next__() # Sn
    s2 = series.__next__() # Sn+1
    while True:
        yield s2 - (sqrt(s2 - s1)) / (s0 -  2*s1  + s2)
        s0, s1, s2 = s1, s2, series.__next__()


def aitken_accelerator(series):
    """
    Lets a series converge even faster
    :param series:
    :return:
    """
    s0 = series.__next__()  # Sn-1
    s1 = series.__next__()  # Sn
    s2 = series.__next__()  # Sn+1
    while True:
        yield s2 - ( ((s2 - s1) ** 2) / ( s2 - (2*s1) + s0) )
        s0, s1, s2 = s1, s2, series.__next__()

def multiplier(series,factor):
    """
    Multiplies every element of a generator by a factor
    :param series: Generator object from where to obtain element
    :param factor: factor to multiply each element with
    :return: multiplied element generator object
    """
    for elem in series:
        yield factor * elem


def converger(series, precision_exponent=25):
    """

    :param series:
    :param precision_exponent:
    :return:
    """
    precision = 10**(-precision_exponent)
    print(f"Converging {series} to a delta precision of {precision}")
    prev = np.inf
    for idx, val in enumerate(series):
        delta = abs(val - prev)
        if delta < precision:
            print(f"Reached target convergence precision at Step #{idx}")
            return idx, val
        prev = val


if __name__ == '__main__':
    to_converge = [ aitken_accelerator(leibniz_series())]
    for result in map(converger, to_converge):
        print(f"After #{result[0]} Steps, converged to a value of {result[1]*4}/4")
        print(f"Difference to pi: {(result[1]*4) - np.pi}")