<s> from . functions import defun , defun_wrapped <EOL> @ defun_wrapped <EOL> def _erf_complex ( ctx , z ) : <EOL> z2 = ctx . square_exp_arg ( z , - <NUM_LIT:1> ) <EOL> v = ( <NUM_LIT:2> / ctx . sqrt ( ctx . pi ) ) * z * ctx . hyp1f1 ( ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , z2 ) <EOL> if not ctx . _re ( z ) : <EOL> v = ctx . _im ( v ) * ctx . j <EOL> return v <EOL> @ defun_wrapped <EOL> def _erfc_complex ( ctx , z ) : <EOL> if ctx . re ( z ) > <NUM_LIT:2> : <EOL> z2 = ctx . square_exp_arg ( z ) <EOL> nz2 = ctx . fneg ( z2 , exact = True ) <EOL> v = ctx . exp ( nz2 ) / ctx . sqrt ( ctx . pi ) * ctx . hyperu ( ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:1> , <NUM_LIT:2> ) , z2 ) <EOL> else : <EOL> v = <NUM_LIT:1> - ctx . _erf_complex ( z ) <EOL> if not ctx . _re ( z ) : <EOL> v = <NUM_LIT:1> + ctx . _im ( v ) * ctx . j <EOL> return v <EOL> @ defun <EOL> def erf ( ctx , z ) : <EOL> z = ctx . convert ( z ) <EOL> if ctx . _is_real_type ( z ) : <EOL> try : <EOL> return ctx . _erf ( z ) <EOL> except NotImplementedError : <EOL> pass <EOL> if ctx . _is_complex_type ( z ) and not z . imag : <EOL> try : <EOL> return type ( z ) ( ctx . _erf ( z . real ) ) <EOL> except NotImplementedError : <EOL> pass <EOL> return ctx . _erf_complex ( z ) <EOL> @ defun <EOL> def erfc ( ctx , z ) : <EOL> z = ctx . convert ( z ) <EOL> if ctx . _is_real_type ( z ) : <EOL> try : <EOL> return ctx . _erfc ( z ) <EOL> except NotImplementedError : <EOL> pass <EOL> if ctx . _is_complex_type ( z ) and not z . imag : <EOL> try : <EOL> return type ( z ) ( ctx . _erfc ( z . real ) ) <EOL> except NotImplementedError : <EOL> pass <EOL> return ctx . _erfc_complex ( z ) <EOL> @ defun <EOL> def square_exp_arg ( ctx , z , mult = <NUM_LIT:1> , reciprocal = False ) : <EOL> prec = ctx . prec * <NUM_LIT:4> + <NUM_LIT:20> <EOL> if reciprocal : <EOL> z2 = ctx . fmul ( z , z , prec = prec ) <EOL> z2 = ctx . fdiv ( ctx . one , z2 , prec = prec ) <EOL> else : <EOL> z2 = ctx . fmul ( z , z , prec = prec ) <EOL> if mult != <NUM_LIT:1> : <EOL> z2 = ctx . fmul ( z2 , mult , exact = True ) <EOL> return z2 <EOL> @ defun_wrapped <EOL> def erfi ( ctx , z ) : <EOL> if not z : <EOL> return z <EOL> z2 = ctx . square_exp_arg ( z ) <EOL> v = ( <NUM_LIT:2> / ctx . sqrt ( ctx . pi ) * z ) * ctx . hyp1f1 ( ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , z2 ) <EOL> if not ctx . _re ( z ) : <EOL> v = ctx . _im ( v ) * ctx . j <EOL> return v <EOL> @ defun_wrapped <EOL> def erfinv ( ctx , x ) : <EOL> xre = ctx . _re ( x ) <EOL> if ( xre != x ) or ( xre < - <NUM_LIT:1> ) or ( xre > <NUM_LIT:1> ) : <EOL> return ctx . bad_domain ( "<STR_LIT>" ) <EOL> x = xre <EOL> if not x : return x <EOL> if x == <NUM_LIT:1> : return ctx . inf <EOL> if x == - <NUM_LIT:1> : return ctx . ninf <EOL> if abs ( x ) < <NUM_LIT> : <EOL> a = <NUM_LIT> * x ** <NUM_LIT:3> + <NUM_LIT> * x <EOL> else : <EOL> u = ctx . ln ( <NUM_LIT:2> / ctx . pi / ( abs ( x ) - <NUM_LIT:1> ) ** <NUM_LIT:2> ) <EOL> a = ctx . sign ( x ) * ctx . sqrt ( u - ctx . ln ( u ) ) / ctx . sqrt ( <NUM_LIT:2> ) <EOL> ctx . prec += <NUM_LIT:10> <EOL> return ctx . findroot ( lambda t : ctx . erf ( t ) - x , a ) <EOL> @ defun_wrapped <EOL> def npdf ( ctx , x , mu = <NUM_LIT:0> , sigma = <NUM_LIT:1> ) : <EOL> sigma = ctx . convert ( sigma ) <EOL> return ctx . exp ( - ( x - mu ) ** <NUM_LIT:2> / ( <NUM_LIT:2> * sigma ** <NUM_LIT:2> ) ) / ( sigma * ctx . sqrt ( <NUM_LIT:2> * ctx . pi ) ) <EOL> @ defun_wrapped <EOL> def ncdf ( ctx , x , mu = <NUM_LIT:0> , sigma = <NUM_LIT:1> ) : <EOL> a = ( x - mu ) / ( sigma * ctx . sqrt ( <NUM_LIT:2> ) ) <EOL> if a < <NUM_LIT:0> : <EOL> return ctx . erfc ( - a ) / <NUM_LIT:2> <EOL> else : <EOL> return ( <NUM_LIT:1> + ctx . erf ( a ) ) / <NUM_LIT:2> <EOL> @ defun_wrapped <EOL> def betainc ( ctx , a , b , x1 = <NUM_LIT:0> , x2 = <NUM_LIT:1> , regularized = False ) : <EOL> if x1 == x2 : <EOL> v = <NUM_LIT:0> <EOL> elif not x1 : <EOL> if x1 == <NUM_LIT:0> and x2 == <NUM_LIT:1> : <EOL> v = ctx . beta ( a , b ) <EOL> else : <EOL> v = x2 ** a * ctx . hyp2f1 ( a , <NUM_LIT:1> - b , a + <NUM_LIT:1> , x2 ) / a <EOL> else : <EOL> m , d = ctx . nint_distance ( a ) <EOL> if m <= <NUM_LIT:0> : <EOL> if d < - ctx . prec : <EOL> h = + ctx . eps <EOL> ctx . prec *= <NUM_LIT:2> <EOL> a += h <EOL> elif d < - <NUM_LIT:4> : <EOL> ctx . prec -= d <EOL> s1 = x2 ** a * ctx . hyp2f1 ( a , <NUM_LIT:1> - b , a + <NUM_LIT:1> , x2 ) <EOL> s2 = x1 ** a * ctx . hyp2f1 ( a , <NUM_LIT:1> - b , a + <NUM_LIT:1> , x1 ) <EOL> v = ( s1 - s2 ) / a <EOL> if regularized : <EOL> v /= ctx . beta ( a , b ) <EOL> return v <EOL> @ defun <EOL> def gammainc ( ctx , z , a = <NUM_LIT:0> , b = None , regularized = False ) : <EOL> regularized = bool ( regularized ) <EOL> z = ctx . convert ( z ) <EOL> if a is None : <EOL> a = ctx . zero <EOL> lower_modified = False <EOL> else : <EOL> a = ctx . convert ( a ) <EOL> lower_modified = a != ctx . zero <EOL> if b is None : <EOL> b = ctx . inf <EOL> upper_modified = False <EOL> else : <EOL> b = ctx . convert ( b ) <EOL> upper_modified = b != ctx . inf <EOL> if not ( upper_modified or lower_modified ) : <EOL> if regularized : <EOL> if ctx . re ( z ) < <NUM_LIT:0> : <EOL> return ctx . inf <EOL> elif ctx . re ( z ) > <NUM_LIT:0> : <EOL> return ctx . one <EOL> else : <EOL> return ctx . nan <EOL> return ctx . gamma ( z ) <EOL> if a == b : <EOL> return ctx . zero <EOL> if ctx . re ( a ) > ctx . re ( b ) : <EOL> return - ctx . gammainc ( z , b , a , regularized ) <EOL> if upper_modified and lower_modified : <EOL> return + ctx . _gamma3 ( z , a , b , regularized ) <EOL> elif lower_modified : <EOL> return ctx . _upper_gamma ( z , a , regularized ) <EOL> elif upper_modified : <EOL> return ctx . _lower_gamma ( z , b , regularized ) <EOL> @ defun <EOL> def _lower_gamma ( ctx , z , b , regularized = False ) : <EOL> if ctx . isnpint ( z ) : <EOL> return type ( z ) ( ctx . inf ) <EOL> G = [ z ] * regularized <EOL> negb = ctx . fneg ( b , exact = True ) <EOL> def h ( z ) : <EOL> T1 = [ ctx . exp ( negb ) , b , z ] , [ <NUM_LIT:1> , z , - <NUM_LIT:1> ] , [ ] , G , [ <NUM_LIT:1> ] , [ <NUM_LIT:1> + z ] , b <EOL> return ( T1 , ) <EOL> return ctx . hypercomb ( h , [ z ] ) <EOL> @ defun <EOL> def _upper_gamma ( ctx , z , a , regularized = False ) : <EOL> if ctx . isint ( z ) : <EOL> try : <EOL> if regularized : <EOL> if ctx . isnpint ( z ) : <EOL> return type ( z ) ( ctx . zero ) <EOL> orig = ctx . prec <EOL> try : <EOL> ctx . prec += <NUM_LIT:10> <EOL> return ctx . _gamma_upper_int ( z , a ) / ctx . gamma ( z ) <EOL> finally : <EOL> ctx . prec = orig <EOL> else : <EOL> return ctx . _gamma_upper_int ( z , a ) <EOL> except NotImplementedError : <EOL> pass <EOL> nega = ctx . fneg ( a , exact = True ) <EOL> G = [ z ] * regularized <EOL> try : <EOL> def h ( z ) : <EOL> r = z - <NUM_LIT:1> <EOL> return [ ( [ ctx . exp ( nega ) , a ] , [ <NUM_LIT:1> , r ] , [ ] , G , [ <NUM_LIT:1> , - r ] , [ ] , <NUM_LIT:1> / nega ) ] <EOL> return ctx . hypercomb ( h , [ z ] , force_series = True ) <EOL> except ctx . NoConvergence : <EOL> def h ( z ) : <EOL> T1 = [ ] , [ <NUM_LIT:1> , z - <NUM_LIT:1> ] , [ z ] , G , [ ] , [ ] , <NUM_LIT:0> <EOL> T2 = [ - ctx . exp ( nega ) , a , z ] , [ <NUM_LIT:1> , z , - <NUM_LIT:1> ] , [ ] , G , [ <NUM_LIT:1> ] , [ <NUM_LIT:1> + z ] , a <EOL> return T1 , T2 <EOL> return ctx . hypercomb ( h , [ z ] ) <EOL> @ defun <EOL> def _gamma3 ( ctx , z , a , b , regularized = False ) : <EOL> pole = ctx . isnpint ( z ) <EOL> if regularized and pole : <EOL> return ctx . zero <EOL> try : <EOL> ctx . prec += <NUM_LIT:15> <EOL> T1 = ctx . gammainc ( z , a , regularized = regularized ) <EOL> T2 = ctx . gammainc ( z , b , regularized = regularized ) <EOL> R = T1 - T2 <EOL> if ctx . mag ( R ) - max ( ctx . mag ( T1 ) , ctx . mag ( T2 ) ) > - <NUM_LIT:10> : <EOL> return R <EOL> if not pole : <EOL> T1 = ctx . gammainc ( z , <NUM_LIT:0> , b , regularized = regularized ) <EOL> T2 = ctx . gammainc ( z , <NUM_LIT:0> , a , regularized = regularized ) <EOL> R = T1 - T2 <EOL> if <NUM_LIT:1> : <EOL> return R <EOL> finally : <EOL> ctx . prec -= <NUM_LIT:15> <EOL> raise NotImplementedError <EOL> @ defun_wrapped <EOL> def expint ( ctx , n , z ) : <EOL> if ctx . isint ( n ) and ctx . _is_real_type ( z ) : <EOL> try : <EOL> return ctx . _expint_int ( n , z ) <EOL> except NotImplementedError : <EOL> pass <EOL> if ctx . isnan ( n ) or ctx . isnan ( z ) : <EOL> return z * n <EOL> if z == ctx . inf : <EOL> return <NUM_LIT:1> / z <EOL> if z == <NUM_LIT:0> : <EOL> if ctx . re ( n ) <= <NUM_LIT:1> : <EOL> return type ( z ) ( ctx . inf ) <EOL> else : <EOL> return ctx . one / ( n - <NUM_LIT:1> ) <EOL> if n == <NUM_LIT:0> : <EOL> return ctx . exp ( - z ) / z <EOL> if n == - <NUM_LIT:1> : <EOL> return ctx . exp ( - z ) * ( z + <NUM_LIT:1> ) / z ** <NUM_LIT:2> <EOL> return z ** ( n - <NUM_LIT:1> ) * ctx . gammainc ( <NUM_LIT:1> - n , z ) <EOL> @ defun_wrapped <EOL> def li ( ctx , z , offset = False ) : <EOL> if offset : <EOL> if z == <NUM_LIT:2> : <EOL> return ctx . zero <EOL> return ctx . ei ( ctx . ln ( z ) ) - ctx . ei ( ctx . ln2 ) <EOL> if not z : <EOL> return z <EOL> if z == <NUM_LIT:1> : <EOL> return ctx . ninf <EOL> return ctx . ei ( ctx . ln ( z ) ) <EOL> @ defun <EOL> def ei ( ctx , z ) : <EOL> try : <EOL> return ctx . _ei ( z ) <EOL> except NotImplementedError : <EOL> return ctx . _ei_generic ( z ) <EOL> @ defun_wrapped <EOL> def _ei_generic ( ctx , z ) : <EOL> if z == ctx . inf : <EOL> return z <EOL> if z == ctx . ninf : <EOL> return ctx . zero <EOL> if ctx . mag ( z ) > <NUM_LIT:1> : <EOL> try : <EOL> r = ctx . one / z <EOL> v = ctx . exp ( z ) * ctx . hyper ( [ <NUM_LIT:1> , <NUM_LIT:1> ] , [ ] , r , <EOL> maxterms = ctx . prec , force_series = True ) / z <EOL> im = ctx . _im ( z ) <EOL> if im > <NUM_LIT:0> : <EOL> v += ctx . pi * ctx . j <EOL> if im < <NUM_LIT:0> : <EOL> v -= ctx . pi * ctx . j <EOL> return v <EOL> except ctx . NoConvergence : <EOL> pass <EOL> v = z * ctx . hyp2f2 ( <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:2> , z ) + ctx . euler <EOL> if ctx . _im ( z ) : <EOL> v += <NUM_LIT:0.5> * ( ctx . log ( z ) - ctx . log ( ctx . one / z ) ) <EOL> else : <EOL> v += ctx . log ( abs ( z ) ) <EOL> return v <EOL> @ defun <EOL> def e1 ( ctx , z ) : <EOL> try : <EOL> return ctx . _e1 ( z ) <EOL> except NotImplementedError : <EOL> return ctx . expint ( <NUM_LIT:1> , z ) <EOL> @ defun <EOL> def ci ( ctx , z ) : <EOL> try : <EOL> return ctx . _ci ( z ) <EOL> except NotImplementedError : <EOL> return ctx . _ci_generic ( z ) <EOL> @ defun_wrapped <EOL> def _ci_generic ( ctx , z ) : <EOL> if ctx . isinf ( z ) : <EOL> if z == ctx . inf : return ctx . zero <EOL> if z == ctx . ninf : return ctx . pi * <NUM_LIT> <EOL> jz = ctx . fmul ( ctx . j , z , exact = True ) <EOL> njz = ctx . fneg ( jz , exact = True ) <EOL> v = <NUM_LIT:0.5> * ( ctx . ei ( jz ) + ctx . ei ( njz ) ) <EOL> zreal = ctx . _re ( z ) <EOL> zimag = ctx . _im ( z ) <EOL> if zreal == <NUM_LIT:0> : <EOL> if zimag > <NUM_LIT:0> : v += ctx . pi * <NUM_LIT> <EOL> if zimag < <NUM_LIT:0> : v -= ctx . pi * <NUM_LIT> <EOL> if zreal < <NUM_LIT:0> : <EOL> if zimag >= <NUM_LIT:0> : v += ctx . pi * <NUM_LIT> <EOL> if zimag < <NUM_LIT:0> : v -= ctx . pi * <NUM_LIT> <EOL> if ctx . _is_real_type ( z ) and zreal > <NUM_LIT:0> : <EOL> v = ctx . _re ( v ) <EOL> return v <EOL> @ defun <EOL> def si ( ctx , z ) : <EOL> try : <EOL> return ctx . _si ( z ) <EOL> except NotImplementedError : <EOL> return ctx . _si_generic ( z ) <EOL> @ defun_wrapped <EOL> def _si_generic ( ctx , z ) : <EOL> if ctx . isinf ( z ) : <EOL> if z == ctx . inf : return <NUM_LIT:0.5> * ctx . pi <EOL> if z == ctx . ninf : return - <NUM_LIT:0.5> * ctx . pi <EOL> if ctx . mag ( z ) >= - <NUM_LIT:1> : <EOL> jz = ctx . fmul ( ctx . j , z , exact = True ) <EOL> njz = ctx . fneg ( jz , exact = True ) <EOL> v = ( - <NUM_LIT> ) * ( ctx . ei ( jz ) - ctx . ei ( njz ) ) <EOL> zreal = ctx . _re ( z ) <EOL> if zreal > <NUM_LIT:0> : <EOL> v -= <NUM_LIT:0.5> * ctx . pi <EOL> if zreal < <NUM_LIT:0> : <EOL> v += <NUM_LIT:0.5> * ctx . pi <EOL> if ctx . _is_real_type ( z ) : <EOL> v = ctx . _re ( v ) <EOL> return v <EOL> else : <EOL> return z * ctx . hyp1f2 ( ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , - <NUM_LIT> * z * z ) <EOL> @ defun_wrapped <EOL> def chi ( ctx , z ) : <EOL> nz = ctx . fneg ( z , exact = True ) <EOL> v = <NUM_LIT:0.5> * ( ctx . ei ( z ) + ctx . ei ( nz ) ) <EOL> zreal = ctx . _re ( z ) <EOL> zimag = ctx . _im ( z ) <EOL> if zimag > <NUM_LIT:0> : <EOL> v += ctx . pi * <NUM_LIT> <EOL> elif zimag < <NUM_LIT:0> : <EOL> v -= ctx . pi * <NUM_LIT> <EOL> elif zreal < <NUM_LIT:0> : <EOL> v += ctx . pi * <NUM_LIT> <EOL> return v <EOL> @ defun_wrapped <EOL> def shi ( ctx , z ) : <EOL> if ctx . mag ( z ) >= - <NUM_LIT:1> : <EOL> nz = ctx . fneg ( z , exact = True ) <EOL> v = <NUM_LIT:0.5> * ( ctx . ei ( z ) - ctx . ei ( nz ) ) <EOL> zimag = ctx . _im ( z ) <EOL> if zimag > <NUM_LIT:0> : v -= <NUM_LIT> * ctx . pi <EOL> if zimag < <NUM_LIT:0> : v += <NUM_LIT> * ctx . pi <EOL> return v <EOL> else : <EOL> return z * ctx . hyp1f2 ( ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , <NUM_LIT> * z * z ) <EOL> @ defun_wrapped <EOL> def fresnels ( ctx , z ) : <EOL> if z == ctx . inf : <EOL> return ctx . mpf ( <NUM_LIT:0.5> ) <EOL> if z == ctx . ninf : <EOL> return ctx . mpf ( - <NUM_LIT:0.5> ) <EOL> return ctx . pi * z ** <NUM_LIT:3> / <NUM_LIT:6> * ctx . hyp1f2 ( ( <NUM_LIT:3> , <NUM_LIT:4> ) , ( <NUM_LIT:3> , <NUM_LIT:2> ) , ( <NUM_LIT:7> , <NUM_LIT:4> ) , - ctx . pi ** <NUM_LIT:2> * z ** <NUM_LIT:4> / <NUM_LIT:16> ) <EOL> @ defun_wrapped <EOL> def fresnelc ( ctx , z ) : <EOL> if z == ctx . inf : <EOL> return ctx . mpf ( <NUM_LIT:0.5> ) <EOL> if z == ctx . ninf : <EOL> return ctx . mpf ( - <NUM_LIT:0.5> ) <EOL> return z * ctx . hyp1f2 ( ( <NUM_LIT:1> , <NUM_LIT:4> ) , ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:5> , <NUM_LIT:4> ) , - ctx . pi ** <NUM_LIT:2> * z ** <NUM_LIT:4> / <NUM_LIT:16> ) </s>
