<s> import math <EOL> from . . libmp . backend import xrange <EOL> class QuadratureRule ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , ctx ) : <EOL> self . ctx = ctx <EOL> self . standard_cache = { } <EOL> self . transformed_cache = { } <EOL> self . interval_count = { } <EOL> def clear ( self ) : <EOL> """<STR_LIT>""" <EOL> self . standard_cache = { } <EOL> self . transformed_cache = { } <EOL> self . interval_count = { } <EOL> def calc_nodes ( self , degree , prec , verbose = False ) : <EOL> r"""<STR_LIT>""" <EOL> raise NotImplementedError <EOL> def get_nodes ( self , a , b , degree , prec , verbose = False ) : <EOL> """<STR_LIT>""" <EOL> key = ( a , b , degree , prec ) <EOL> if key in self . transformed_cache : <EOL> return self . transformed_cache [ key ] <EOL> orig = self . ctx . prec <EOL> try : <EOL> self . ctx . prec = prec + <NUM_LIT:20> <EOL> if ( degree , prec ) in self . standard_cache : <EOL> nodes = self . standard_cache [ degree , prec ] <EOL> else : <EOL> nodes = self . calc_nodes ( degree , prec , verbose ) <EOL> self . standard_cache [ degree , prec ] = nodes <EOL> nodes = self . transform_nodes ( nodes , a , b , verbose ) <EOL> if key in self . interval_count : <EOL> self . transformed_cache [ key ] = nodes <EOL> else : <EOL> self . interval_count [ key ] = True <EOL> finally : <EOL> self . ctx . prec = orig <EOL> return nodes <EOL> def transform_nodes ( self , nodes , a , b , verbose = False ) : <EOL> r"""<STR_LIT>""" <EOL> ctx = self . ctx <EOL> a = ctx . convert ( a ) <EOL> b = ctx . convert ( b ) <EOL> one = ctx . one <EOL> if ( a , b ) == ( - one , one ) : <EOL> return nodes <EOL> half = ctx . mpf ( <NUM_LIT:0.5> ) <EOL> new_nodes = [ ] <EOL> if ctx . isinf ( a ) or ctx . isinf ( b ) : <EOL> if ( a , b ) == ( ctx . ninf , ctx . inf ) : <EOL> p05 = - half <EOL> for x , w in nodes : <EOL> x2 = x * x <EOL> px1 = one - x2 <EOL> spx1 = px1 ** p05 <EOL> x = x * spx1 <EOL> w *= spx1 / px1 <EOL> new_nodes . append ( ( x , w ) ) <EOL> elif a == ctx . ninf : <EOL> b1 = b + <NUM_LIT:1> <EOL> for x , w in nodes : <EOL> u = <NUM_LIT:2> / ( x + one ) <EOL> x = b1 - u <EOL> w *= half * u ** <NUM_LIT:2> <EOL> new_nodes . append ( ( x , w ) ) <EOL> elif b == ctx . inf : <EOL> a1 = a - <NUM_LIT:1> <EOL> for x , w in nodes : <EOL> u = <NUM_LIT:2> / ( x + one ) <EOL> x = a1 + u <EOL> w *= half * u ** <NUM_LIT:2> <EOL> new_nodes . append ( ( x , w ) ) <EOL> elif a == ctx . inf or b == ctx . ninf : <EOL> return [ ( x , - w ) for ( x , w ) in self . transform_nodes ( nodes , b , a , verbose ) ] <EOL> else : <EOL> raise NotImplementedError <EOL> else : <EOL> C = ( b - a ) / <NUM_LIT:2> <EOL> D = ( b + a ) / <NUM_LIT:2> <EOL> for x , w in nodes : <EOL> new_nodes . append ( ( D + C * x , C * w ) ) <EOL> return new_nodes <EOL> def guess_degree ( self , prec ) : <EOL> """<STR_LIT>""" <EOL> g = int ( <NUM_LIT:4> + max ( <NUM_LIT:0> , self . ctx . log ( prec / <NUM_LIT> , <NUM_LIT:2> ) ) ) <EOL> g += <NUM_LIT:2> <EOL> return g <EOL> def estimate_error ( self , results , prec , epsilon ) : <EOL> r"""<STR_LIT>""" <EOL> if len ( results ) == <NUM_LIT:2> : <EOL> return abs ( results [ <NUM_LIT:0> ] - results [ <NUM_LIT:1> ] ) <EOL> try : <EOL> if results [ - <NUM_LIT:1> ] == results [ - <NUM_LIT:2> ] == results [ - <NUM_LIT:3> ] : <EOL> return self . ctx . zero <EOL> D1 = self . ctx . log ( abs ( results [ - <NUM_LIT:1> ] - results [ - <NUM_LIT:2> ] ) , <NUM_LIT:10> ) <EOL> D2 = self . ctx . log ( abs ( results [ - <NUM_LIT:1> ] - results [ - <NUM_LIT:3> ] ) , <NUM_LIT:10> ) <EOL> except ValueError : <EOL> return epsilon <EOL> D3 = - prec <EOL> D4 = min ( <NUM_LIT:0> , max ( D1 ** <NUM_LIT:2> / D2 , <NUM_LIT:2> * D1 , D3 ) ) <EOL> return self . ctx . mpf ( <NUM_LIT:10> ) ** int ( D4 ) <EOL> def summation ( self , f , points , prec , epsilon , max_degree , verbose = False ) : <EOL> """<STR_LIT>""" <EOL> ctx = self . ctx <EOL> I = err = ctx . zero <EOL> for i in xrange ( len ( points ) - <NUM_LIT:1> ) : <EOL> a , b = points [ i ] , points [ i + <NUM_LIT:1> ] <EOL> if a == b : <EOL> continue <EOL> if ( a , b ) == ( ctx . ninf , ctx . inf ) : <EOL> _f = f <EOL> f = lambda x : _f ( - x ) + _f ( x ) <EOL> a , b = ( ctx . zero , ctx . inf ) <EOL> results = [ ] <EOL> for degree in xrange ( <NUM_LIT:1> , max_degree + <NUM_LIT:1> ) : <EOL> nodes = self . get_nodes ( a , b , degree , prec , verbose ) <EOL> if verbose : <EOL> print ( "<STR_LIT>" % ( ctx . nstr ( a ) , ctx . nstr ( b ) , degree , max_degree ) ) <EOL> results . append ( self . sum_next ( f , nodes , degree , prec , results , verbose ) ) <EOL> if degree > <NUM_LIT:1> : <EOL> err = self . estimate_error ( results , prec , epsilon ) <EOL> if err <= epsilon : <EOL> break <EOL> if verbose : <EOL> print ( "<STR_LIT>" , ctx . nstr ( err ) ) <EOL> I += results [ - <NUM_LIT:1> ] <EOL> if err > epsilon : <EOL> if verbose : <EOL> print ( "<STR_LIT>" , ctx . nstr ( err ) ) <EOL> return I , err <EOL> def sum_next ( self , f , nodes , degree , prec , previous , verbose = False ) : <EOL> r"""<STR_LIT>""" <EOL> return self . ctx . fdot ( ( w , f ( x ) ) for ( x , w ) in nodes ) <EOL> class TanhSinh ( QuadratureRule ) : <EOL> r"""<STR_LIT>""" <EOL> def sum_next ( self , f , nodes , degree , prec , previous , verbose = False ) : <EOL> """<STR_LIT>""" <EOL> h = self . ctx . mpf ( <NUM_LIT:2> ) ** ( - degree ) <EOL> if previous : <EOL> S = previous [ - <NUM_LIT:1> ] / ( h * <NUM_LIT:2> ) <EOL> else : <EOL> S = self . ctx . zero <EOL> S += self . ctx . fdot ( ( w , f ( x ) ) for ( x , w ) in nodes ) <EOL> return h * S <EOL> def calc_nodes ( self , degree , prec , verbose = False ) : <EOL> r"""<STR_LIT>""" <EOL> ctx = self . ctx <EOL> nodes = [ ] <EOL> extra = <NUM_LIT:20> <EOL> ctx . prec += extra <EOL> tol = ctx . ldexp ( <NUM_LIT:1> , - prec - <NUM_LIT:10> ) <EOL> pi4 = ctx . pi / <NUM_LIT:4> <EOL> t0 = ctx . ldexp ( <NUM_LIT:1> , - degree ) <EOL> if degree == <NUM_LIT:1> : <EOL> nodes . append ( ( ctx . zero , ctx . pi / <NUM_LIT:2> ) ) <EOL> h = t0 <EOL> else : <EOL> h = t0 * <NUM_LIT:2> <EOL> expt0 = ctx . exp ( t0 ) <EOL> a = pi4 * expt0 <EOL> b = pi4 / expt0 <EOL> udelta = ctx . exp ( h ) <EOL> urdelta = <NUM_LIT:1> / udelta <EOL> for k in xrange ( <NUM_LIT:0> , <NUM_LIT:20> * <NUM_LIT:2> ** degree + <NUM_LIT:1> ) : <EOL> c = ctx . exp ( a - b ) <EOL> d = <NUM_LIT:1> / c <EOL> co = ( c + d ) / <NUM_LIT:2> <EOL> si = ( c - d ) / <NUM_LIT:2> <EOL> x = si / co <EOL> w = ( a + b ) / co ** <NUM_LIT:2> <EOL> diff = abs ( x - <NUM_LIT:1> ) <EOL> if diff <= tol : <EOL> break <EOL> nodes . append ( ( x , w ) ) <EOL> nodes . append ( ( - x , w ) ) <EOL> a *= udelta <EOL> b *= urdelta <EOL> if verbose and k % <NUM_LIT> == <NUM_LIT> : <EOL> print ( "<STR_LIT>" , ctx . nstr ( - ctx . log ( diff , <NUM_LIT:10> ) / prec ) ) <EOL> ctx . prec -= extra <EOL> return nodes <EOL> class GaussLegendre ( QuadratureRule ) : <EOL> """<STR_LIT>""" <EOL> def calc_nodes ( self , degree , prec , verbose = False ) : <EOL> """<STR_LIT>""" <EOL> ctx = self . ctx <EOL> epsilon = ctx . ldexp ( <NUM_LIT:1> , - prec - <NUM_LIT:8> ) <EOL> orig = ctx . prec <EOL> ctx . prec = int ( prec * <NUM_LIT> ) <EOL> if degree == <NUM_LIT:1> : <EOL> x = ctx . mpf ( <NUM_LIT:3> ) / <NUM_LIT:5> <EOL> w = ctx . mpf ( <NUM_LIT:5> ) / <NUM_LIT:9> <EOL> nodes = [ ( - x , w ) , ( ctx . zero , ctx . mpf ( <NUM_LIT:8> ) / <NUM_LIT:9> ) , ( x , w ) ] <EOL> ctx . prec = orig <EOL> return nodes <EOL> nodes = [ ] <EOL> n = <NUM_LIT:3> * <NUM_LIT:2> ** ( degree - <NUM_LIT:1> ) <EOL> upto = n // <NUM_LIT:2> + <NUM_LIT:1> <EOL> for j in xrange ( <NUM_LIT:1> , upto ) : <EOL> r = ctx . mpf ( math . cos ( math . pi * ( j - <NUM_LIT> ) / ( n + <NUM_LIT:0.5> ) ) ) <EOL> while <NUM_LIT:1> : <EOL> t1 , t2 = <NUM_LIT:1> , <NUM_LIT:0> <EOL> for j1 in xrange ( <NUM_LIT:1> , n + <NUM_LIT:1> ) : <EOL> t3 , t2 , t1 = t2 , t1 , ( ( <NUM_LIT:2> * j1 - <NUM_LIT:1> ) * r * t1 - ( j1 - <NUM_LIT:1> ) * t2 ) / j1 <EOL> t4 = n * ( r * t1 - t2 ) / ( r ** <NUM_LIT:2> - <NUM_LIT:1> ) <EOL> t5 = r <EOL> a = t1 / t4 <EOL> r = r - a <EOL> if abs ( a ) < epsilon : <EOL> break <EOL> x = r <EOL> w = <NUM_LIT:2> / ( ( <NUM_LIT:1> - r ** <NUM_LIT:2> ) * t4 ** <NUM_LIT:2> ) <EOL> if verbose and j % <NUM_LIT:30> == <NUM_LIT:15> : <EOL> print ( "<STR_LIT>" % ( j , upto ) ) <EOL> nodes . append ( ( x , w ) ) <EOL> nodes . append ( ( - x , w ) ) <EOL> ctx . prec = orig <EOL> return nodes <EOL> class QuadratureMethods ( object ) : <EOL> def __init__ ( ctx , * args , ** kwargs ) : <EOL> ctx . _gauss_legendre = GaussLegendre ( ctx ) <EOL> ctx . _tanh_sinh = TanhSinh ( ctx ) <EOL> def quad ( ctx , f , * points , ** kwargs ) : <EOL> r"""<STR_LIT>""" <EOL> rule = kwargs . get ( '<STR_LIT>' , '<STR_LIT>' ) <EOL> if type ( rule ) is str : <EOL> if rule == '<STR_LIT>' : <EOL> rule = ctx . _tanh_sinh <EOL> elif rule == '<STR_LIT>' : <EOL> rule = ctx . _gauss_legendre <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" % rule ) <EOL> else : <EOL> rule = rule ( ctx ) <EOL> verbose = kwargs . get ( '<STR_LIT>' ) <EOL> dim = len ( points ) <EOL> orig = prec = ctx . prec <EOL> epsilon = ctx . eps / <NUM_LIT:8> <EOL> m = kwargs . get ( '<STR_LIT>' ) or rule . guess_degree ( prec ) <EOL> points = [ ctx . _as_points ( p ) for p in points ] <EOL> try : <EOL> ctx . prec += <NUM_LIT:20> <EOL> if dim == <NUM_LIT:1> : <EOL> v , err = rule . summation ( f , points [ <NUM_LIT:0> ] , prec , epsilon , m , verbose ) <EOL> elif dim == <NUM_LIT:2> : <EOL> v , err = rule . summation ( lambda x : rule . summation ( lambda y : f ( x , y ) , points [ <NUM_LIT:1> ] , prec , epsilon , m ) [ <NUM_LIT:0> ] , <EOL> points [ <NUM_LIT:0> ] , prec , epsilon , m , verbose ) <EOL> elif dim == <NUM_LIT:3> : <EOL> v , err = rule . summation ( lambda x : rule . summation ( lambda y : rule . summation ( lambda z : f ( x , y , z ) , points [ <NUM_LIT:2> ] , prec , epsilon , m ) [ <NUM_LIT:0> ] , <EOL> points [ <NUM_LIT:1> ] , prec , epsilon , m ) [ <NUM_LIT:0> ] , <EOL> points [ <NUM_LIT:0> ] , prec , epsilon , m , verbose ) <EOL> else : <EOL> raise NotImplementedError ( "<STR_LIT>" ) <EOL> finally : <EOL> ctx . prec = orig <EOL> if kwargs . get ( "<STR_LIT:error>" ) : <EOL> return + v , err <EOL> return + v <EOL> def quadts ( ctx , * args , ** kwargs ) : <EOL> """<STR_LIT>""" <EOL> kwargs [ '<STR_LIT>' ] = '<STR_LIT>' <EOL> return ctx . quad ( * args , ** kwargs ) <EOL> def quadgl ( ctx , * args , ** kwargs ) : <EOL> """<STR_LIT>""" <EOL> kwargs [ '<STR_LIT>' ] = '<STR_LIT>' <EOL> return ctx . quad ( * args , ** kwargs ) <EOL> def quadosc ( ctx , f , interval , omega = None , period = None , zeros = None ) : <EOL> r"""<STR_LIT>""" <EOL> a , b = ctx . _as_points ( interval ) <EOL> a = ctx . convert ( a ) <EOL> b = ctx . convert ( b ) <EOL> if [ omega , period , zeros ] . count ( None ) != <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if a == ctx . ninf and b == ctx . inf : <EOL> s1 = ctx . quadosc ( f , [ a , <NUM_LIT:0> ] , omega = omega , zeros = zeros , period = period ) <EOL> s2 = ctx . quadosc ( f , [ <NUM_LIT:0> , b ] , omega = omega , zeros = zeros , period = period ) <EOL> return s1 + s2 <EOL> if a == ctx . ninf : <EOL> if zeros : <EOL> return ctx . quadosc ( lambda x : f ( - x ) , [ - b , - a ] , lambda n : zeros ( - n ) ) <EOL> else : <EOL> return ctx . quadosc ( lambda x : f ( - x ) , [ - b , - a ] , omega = omega , period = period ) <EOL> if b != ctx . inf : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if not zeros : <EOL> if omega : <EOL> period = <NUM_LIT:2> * ctx . pi / omega <EOL> zeros = lambda n : n * period / <NUM_LIT:2> <EOL> n = <NUM_LIT:1> <EOL> s = ctx . quadgl ( f , [ a , zeros ( n ) ] ) <EOL> def term ( k ) : <EOL> return ctx . quadgl ( f , [ zeros ( k ) , zeros ( k + <NUM_LIT:1> ) ] ) <EOL> s += ctx . nsum ( term , [ n , ctx . inf ] ) <EOL> return s <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> import doctest <EOL> doctest . testmod ( ) </s>
