# romberg.py 						
# Section: 5.3	
'''
  The Romberg algoritlm produces a triangular array of 
  numbers, which are numerical estimates of the definite integral
   \integral_a^b f(x) dx

  Compute an integral romberg array

  Use the Trapezoid rule to get the first column

  The second and successive columns in the Romberg array 
  are generated by the Richardson extrapolation formula
'''

import math

SZ = 5

def f(x):
  return 12*(math.sqrt(1-x*x) - math.sqrt(3)*x)


def romberg(f, a, b):
  print("Range:", a, "to", b)

  r = [0]*(SZ)
  for i in range(SZ):
    r[i] = [0]*SZ
  h = b - a
  r[0][0] = 0.5 * h *(f(a) + f(b))   

  # computes row-by-row 
  for i in range(1, SZ):
     h *= 0.5
     tot = 0
     k = 1
     while k <= math.pow(2,i)-1:
       tot += f(a + k*h)
       k += 2

     r[i][0] = 0.5 * r[i-1][0] + tot * h  
     # Trapezoid rule for first column

     for j in range(1, i+1):
       r[i][j] = r[i][j-1] + \
            (r[i][j-1] - r[i-1][j-1]) /(math.pow(4,j)-1) 
       # Richardson extrapolation for other columns

  # print triangular array
  for i in range(SZ):
    for j in range(i+1):
      print(f"{r[i][j]:13.10f}", end = '')
    print()
  print()


if __name__ == "__main__":
  print("Integrate f(x) = 12*(sqrt(1-x*x) - sqrt(3)*x")
  romberg(f, 0.0, 0.5)   # calculates pi
  
  print("Integrate 4*sqrt(1-x*x)")
  romberg((lambda x: 4*math.sqrt(1-x*x)), 0, 1)   
                          # calculates pi poorly
  print("Integrate exp(x)")
  romberg(math.exp, -1.0, 1.0) # e^1 - e^-1 = 2.3504
  
  print("Integrate 1/x")
  romberg((lambda x: 1/x), 1.0, 2.0)  # ln 2 = 0.69314718