SHARE
TWEET

Untitled

a guest Sep 19th, 2019 84 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import random
  2. from math import sqrt
  3. import numpy as np
  4.  
  5. def harmonic_series(stop_after=-1):
  6.     """
  7.  
  8.     :param stop_after:
  9.     :return:
  10.     """
  11.     sum = np.longdouble(0.0)
  12.     k = 1
  13.     while True:
  14.         sum += 1/k
  15.         yield sum
  16.         if k == stop_after:
  17.             break
  18.         k += 1
  19.  
  20. def alternating_harmonic_series(stop_after=-1):
  21.     """
  22.  
  23.     :param stop_after:
  24.     :return:
  25.     """
  26.     sum = np.longdouble(0.0)
  27.     k = 1
  28.     negative = False
  29.     while True:
  30.         if negative:
  31.             sum += (-1)/k
  32.         else:
  33.             sum += 1/k
  34.         yield sum
  35.         if k == stop_after:
  36.             break
  37.         k += 1
  38.         negative = not negative
  39.  
  40. def general_harmonic_series(a, b, alternating=False, stop_after=-1):
  41.     """
  42.  
  43.     :param a:
  44.     :param b:
  45.     :param alternating:
  46.     :param stop_after:
  47.     :return:
  48.     """
  49.     assert a != 0
  50.     assert b/a > 0.0
  51.     sum = np.longdouble(0.0)
  52.     n = 0
  53.     negative = False
  54.     while True:
  55.         if negative:
  56.             sum += (-1) / ((a*n) + b)
  57.         else:
  58.             sum += 1 / ((a*n) + b)
  59.         yield sum
  60.         if n == stop_after:
  61.             break
  62.         if alternating:
  63.             negative = not negative
  64.         n += 1
  65.  
  66. def leibniz_series(stop_after=-1):
  67.     """
  68.     Approaches pi/4
  69.     :param stop_after: if specified, stop after n-th element of series
  70.     :return: leibniz series generator object
  71.     """
  72.     sum = np.longdouble(0.0)
  73.     k = 0
  74.     negative = False
  75.     while True:
  76.         if negative:
  77.             sum += (-1) / ((2*k)+1)
  78.         else:
  79.             sum += 1 / ((2*k)+1)
  80.         yield sum
  81.         if k == stop_after:
  82.             break
  83.         k += 1
  84.         negative = not negative
  85.  
  86. def p_series(p, stop_after=-1):
  87.     """
  88.  
  89.         :param stop_after:
  90.         :return:
  91.     """
  92.     sum = np.longdouble(0.0)
  93.     k = 1
  94.     while True:
  95.         sum += 1 / (k ** p)
  96.         yield sum
  97.         if k == stop_after:
  98.             break
  99.         k += 1
  100.  
  101.  
  102. def random_harmonic_series(stop_after=-1):
  103.     """
  104.  
  105.     :param stop_after:
  106.     :return:
  107.     """
  108.     sum = np.longdouble(0.0)
  109.     k = 1
  110.     while True:
  111.         sum += random.choice([-1, 1]) / k
  112.         yield sum
  113.         if k == stop_after:
  114.             break
  115.         k += 1
  116.  
  117.  
  118. def euler_accelerator(series):
  119.     """
  120.     Lets a series converge faster
  121.     :param series:
  122.     :return:
  123.     """
  124.     s0 = series.__next__() # Sn-1
  125.     s1 = series.__next__() # Sn
  126.     s2 = series.__next__() # Sn+1
  127.     while True:
  128.         yield s2 - (sqrt(s2 - s1)) / (s0 -  2*s1  + s2)
  129.         s0, s1, s2 = s1, s2, series.__next__()
  130.  
  131.  
  132. def aitken_accelerator(series):
  133.     """
  134.     Lets a series converge even faster
  135.     :param series:
  136.     :return:
  137.     """
  138.     s0 = series.__next__()  # Sn-1
  139.     s1 = series.__next__()  # Sn
  140.     s2 = series.__next__()  # Sn+1
  141.     while True:
  142.         yield s2 - ( ((s2 - s1) ** 2) / ( s2 - (2*s1) + s0) )
  143.         s0, s1, s2 = s1, s2, series.__next__()
  144.  
  145. def multiplier(series,factor):
  146.     """
  147.     Multiplies every element of a generator by a factor
  148.     :param series: Generator object from where to obtain element
  149.     :param factor: factor to multiply each element with
  150.     :return: multiplied element generator object
  151.     """
  152.     for elem in series:
  153.         yield factor * elem
  154.  
  155.  
  156. def converger(series, precision_exponent=25):
  157.     """
  158.  
  159.     :param series:
  160.     :param precision_exponent:
  161.     :return:
  162.     """
  163.     precision = 10**(-precision_exponent)
  164.     print(f"Converging {series} to a delta precision of {precision}")
  165.     prev = np.inf
  166.     for idx, val in enumerate(series):
  167.         delta = abs(val - prev)
  168.         if delta < precision:
  169.             print(f"Reached target convergence precision at Step #{idx}")
  170.             return idx, val
  171.         prev = val
  172.  
  173.  
  174. if __name__ == '__main__':
  175.     to_converge = [ aitken_accelerator(leibniz_series())]
  176.     for result in map(converger, to_converge):
  177.         print(f"After #{result[0]} Steps, converged to a value of {result[1]*4}/4")
  178.         print(f"Difference to pi: {(result[1]*4) - np.pi}")
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top