import numpy as np

def gauss_quadrature_nodes_coeffs(n_nodes=2,a=-1.0,b=1.0):
  """
  Return the roots and coefficients for gauss quadrature.

  With n_nodes nodes we can integrate any polynomial of degree
  2*n_nodes-1 exactly.

  Taken from Numerical Analysis, Burden and Faires

  n_nodes: number of quadrature nodes. 
  a,b: left and right hand side bounds of the integral.
  """
  assert n_nodes >= 2
  assert n_nodes <= 16

  roots = [[0.5773502692,-0.5773502692], # 2 nodes
           [0.7745966692,0.0000000000,-0.7745966692], # 3 nodes
           [0.8611363116,0.3399810436,-0.3399810436,-0.8611363116],
           [0.9061798459,0.5384693101,0.0000000000,-0.5384693101,-0.9061798459],
           [-0.9324695142031520278123,-0.661209386466264513661,-0.2386191860831969086305,0.238619186083196908631,0.661209386466264513661,0.9324695142031520278123],
           [-0.949107912342758525,-0.74153118559939444,-0.405845151377397167,0.0,0.405845151377397167,0.74153118559939444,0.949107912342758525],
           [-0.960289856497536232,-0.79666647741362674,-0.525532409916328986,-0.183434642495649805,0.183434642495649805,0.525532409916328986,0.79666647741362674,0.960289856497536232],
            [-0.96816023950762609,-0.836031107326635794,-0.613371432700590397,-0.32425342340380893,0.0,0.32425342340380893,0.613371432700590397,0.836031107326635794,0.96816023950762609],
            [-0.97390652851717172,-0.86506336668898451,-0.679409568299024406,-0.433395394129247191,-0.14887433898163121,0.148874338981631211,0.433395394129247191,0.679409568299024406,0.865063366688984511,0.97390652851717172],
            [-0.978228658146056993,-0.887062599768095299,-0.730152005574049324,-0.519096129206811816,-0.26954315595234497,0.0,0.269543155952344972,0.519096129206811816,0.730152005574049324,0.887062599768095299,0.978228658146056993],
            [-0.981560634246719251,-0.904117256370474857,-0.769902674194304687,-0.587317954286617447,-0.367831498998180194,-0.125233408511468916,0.125233408511468916,0.367831498998180194,0.587317954286617447,0.769902674194304687,0.904117256370474857,0.981560634246719251],
            [-0.98418305471858815,-0.917598399222977965,-0.801578090733309913,-0.642349339440340221,-0.448492751036446853,-0.230458315955134794,0.0,0.23045831595513479,0.44849275103644685,0.642349339440340221,0.801578090733309913,0.917598399222977965,0.98418305471858815],
            [-0.986283808696812339,-0.928434883663573517,-0.827201315069764993,-0.68729290481168547,-0.515248636358154092,-0.31911236892788976,-0.108054948707343662,0.108054948707343662,0.31911236892788976,0.515248636358154092,0.68729290481168547,0.827201315069764993,0.928434883663573517,0.986283808696812339],
            [-0.987992518020485429,-0.9372733924007059,-0.848206583410427216,-0.724417731360170047,-0.570972172608538848,-0.39415134707756337,-0.201194093997434522,0.0,0.201194093997434522,0.39415134707756337,0.570972172608538848,0.724417731360170047,0.848206583410427216,0.937273392400705904,0.987992518020485429],
            [-0.989400934991649933,-0.944575023073232576,-0.865631202387831744,-0.75540440835500303,-0.61787624440264375,-0.458016777657227386,-0.281603550779258913,-0.09501250983763744,0.09501250983763744,0.281603550779258913,0.458016777657227386,0.617876244402643748,0.755404408355003034,0.86563120238783174,0.944575023073232576,0.989400934991649933]
           ]
  coeffs = [[1.0,1.0],
            [0.5555555556,0.8888888889,0.5555555556],
            [0.3478548451,0.6521451549,0.6521451549,0.3478548451],
            [0.2369268850,0.4786286705,0.5688888889,0.4786286705,0.2369268850],
            [0.1713244923791703450403,0.3607615730481386075698,0.4679139345726910473899,0.46791393457269104739,0.3607615730481386075698,0.1713244923791703450403],
            [0.129484966168869693,0.279705391489276668,0.381830050505118945,0.417959183673469388,0.38183005050511894,0.279705391489276668,0.129484966168869693],
            [0.101228536290376259,0.22238103445337447,0.313706645877887287,0.362683783378361983,0.36268378337836198,0.313706645877887287,0.222381034453374471,0.10122853629037626],
            [0.081274388361574412,0.1806481606948574,0.260610696402935462,0.31234707704000284,0.33023935500125976,0.31234707704000284,0.26061069640293546,0.180648160694857404,0.081274388361574412],
            [0.066671344308688138,0.14945134915058059,0.219086362515982044,0.26926671930999636,0.29552422471475287,0.29552422471475287,0.269266719309996355,0.21908636251598204,0.14945134915058059,0.066671344308688138],
            [0.055668567116173666,0.125580369464904625,0.186290210927734251,0.23319376459199048,0.26280454451024666,0.27292508677790063,0.262804544510246662,0.23319376459199048,0.18629021092773425,0.12558036946490462,0.055668567116173666],
            [0.047175336386511827,0.10693932599531843,0.160078328543346226,0.203167426723065922,0.23349253653835481,0.249147045813402785,0.24914704581340279,0.233492536538354809,0.203167426723065922,0.16007832854334623,0.106939325995318431,0.047175336386511827],
            [0.04048400476531588,0.092121499837728448,0.138873510219787239,0.178145980761945738,0.207816047536888502,0.22628318026289724,0.23255155323087391,0.22628318026289724,0.207816047536888502,0.17814598076194574,0.138873510219787239,0.09212149983772845,0.04048400476531588],
            [0.035119460331751863,0.08015808715976021,0.121518570687903185,0.157203167158193535,0.185538397477937814,0.205198463721295604,0.21526385346315779,0.21526385346315779,0.205198463721295604,0.18553839747793781,0.157203167158193535,0.121518570687903185,0.08015808715976021,0.035119460331751863],
            [0.0307532419961172684,0.0703660474881081247,0.107159220467171935,0.139570677926154315,0.166269205816993934,0.18616100001556221,0.198431485327111577,0.202578241925561273,0.19843148532711158,0.18616100001556221,0.166269205816993934,0.13957067792615431,0.10715922046717194,0.07036604748810812,0.030753241996117268],
            [0.027152459411754095,0.062253523938647893,0.095158511682492785,0.12462897125553387,0.149595988816576732,0.169156519395002538,0.18260341504492359,0.1894506104550685,0.189450610455068496,0.182603415044923589,0.169156519395002538,0.149595988816576732,0.12462897125553387,0.09515851168249278,0.062253523938647893,0.02715245941175409]
           ]

  # transform to other interval
  nodes = ((b-a)*np.array(roots[n_nodes-2]) + (b+a))/2.0
  coeffs = np.array(coeffs[n_nodes-2])*(b-a)/2.0

  return nodes,coeffs


if __name__ == "__main__":

  f = lambda x : x**2 -2*x - 3
  fint = lambda x: x**3 / 3 - x**2 -3*x # integral
  a = -1.0
  b = 1.0
  true_int = fint(b) - fint(a)
  nodes,coeffs = gauss_quadrature_nodes_coeffs(n_nodes=2,a=a,b=b)
  fX = f(nodes)
  print(fX @ coeffs)
  print(true_int)
 
  f = lambda x : x**2 -2*x - 3
  fint = lambda x: x**3 / 3 - x**2 -3*x # integral
  a = -7.0
  b = 9.21
  true_int = fint(b) - fint(a)
  nodes,coeffs = gauss_quadrature_nodes_coeffs(n_nodes=2,a=a,b=b)
  fX = f(nodes)
  print(fX @ coeffs)
  print(true_int)


  f = lambda x : 13*x**7 -102.8*x**5 + x**2 -2*x - 14.
  fint = lambda x: 13*x**8 / 8 - 102.8*x**6/6 + x**3 / 3 - x**2 -14*x # integral
  a = -1.0
  b = 4.0
  true_int = fint(b) - fint(a)
  n_nodes = 4
  nodes,coeffs = gauss_quadrature_nodes_coeffs(n_nodes=n_nodes,a=a,b=b)
  fX = f(nodes)
  print(fX @ coeffs)
  print(true_int)


  f = lambda x : 13*x**12 -102.8*x**5 + x**2 -2*x - 14.
  fint = lambda x: 13*x**13 / 13 - 102.8*x**6/6 + x**3 / 3 - x**2 -14*x # integral
  a = -1.0
  b = 4.0
  true_int = fint(b) - fint(a)
  n_nodes = 7
  nodes,coeffs = gauss_quadrature_nodes_coeffs(n_nodes=n_nodes,a=a,b=b)
  fX = f(nodes)
  print(fX @ coeffs)
  print(true_int)

  f = lambda x : 13*x**19 -102.8*x**5 + x**2 -2*x - 14.
  fint = lambda x: 13*x**20 / 20 - 102.8*x**6/6 + x**3 / 3 - x**2 -14*x # integral
  a = -1.0
  b = 4.0
  true_int = fint(b) - fint(a)
  n_nodes = 10
  nodes,coeffs = gauss_quadrature_nodes_coeffs(n_nodes=n_nodes,a=a,b=b)
  fX = f(nodes)
  print(fX @ coeffs)
  print(true_int)
