<s> """<STR_LIT>""" <EOL> import math <EOL> from bisect import bisect <EOL> from . backend import xrange <EOL> from . backend import MPZ , MPZ_ZERO , MPZ_ONE , MPZ_TWO , MPZ_FIVE , BACKEND <EOL> from . libmpf import ( <EOL> round_floor , round_ceiling , round_down , round_up , <EOL> round_nearest , round_fast , <EOL> ComplexResult , <EOL> bitcount , bctable , lshift , rshift , giant_steps , sqrt_fixed , <EOL> from_int , to_int , from_man_exp , to_fixed , to_float , from_float , <EOL> from_rational , normalize , <EOL> fzero , fone , fnone , fhalf , finf , fninf , fnan , <EOL> mpf_cmp , mpf_sign , mpf_abs , <EOL> mpf_pos , mpf_neg , mpf_add , mpf_sub , mpf_mul , mpf_div , mpf_shift , <EOL> mpf_rdiv_int , mpf_pow_int , mpf_sqrt , <EOL> reciprocal_rnd , negative_rnd , mpf_perturb , <EOL> isqrt_fast <EOL> ) <EOL> from . libintmath import ifib <EOL> if BACKEND == '<STR_LIT>' : <EOL> EXP_COSH_CUTOFF = <NUM_LIT> <EOL> else : <EOL> EXP_COSH_CUTOFF = <NUM_LIT> <EOL> EXP_SERIES_U_CUTOFF = <NUM_LIT> <EOL> if BACKEND == '<STR_LIT>' : <EOL> COS_SIN_CACHE_PREC = <NUM_LIT> <EOL> else : <EOL> COS_SIN_CACHE_PREC = <NUM_LIT:200> <EOL> COS_SIN_CACHE_STEP = <NUM_LIT:8> <EOL> cos_sin_cache = { } <EOL> MAX_LOG_INT_CACHE = <NUM_LIT> <EOL> log_int_cache = { } <EOL> LOG_TAYLOR_PREC = <NUM_LIT> <EOL> LOG_TAYLOR_SHIFT = <NUM_LIT:9> <EOL> log_taylor_cache = { } <EOL> LOG_AGM_MAG_PREC_RATIO = <NUM_LIT:20> <EOL> ATAN_TAYLOR_PREC = <NUM_LIT> <EOL> ATAN_TAYLOR_SHIFT = <NUM_LIT:7> <EOL> atan_taylor_cache = { } <EOL> cache_prec_steps = [ <NUM_LIT> , <NUM_LIT> ] <EOL> for k in xrange ( <NUM_LIT:1> , bitcount ( LOG_TAYLOR_PREC ) + <NUM_LIT:1> ) : <EOL> cache_prec_steps += [ min ( <NUM_LIT:2> ** k , LOG_TAYLOR_PREC ) + <NUM_LIT:20> ] * <NUM_LIT:2> ** ( k - <NUM_LIT:1> ) <EOL> def constant_memo ( f ) : <EOL> """<STR_LIT>""" <EOL> f . memo_prec = - <NUM_LIT:1> <EOL> f . memo_val = None <EOL> def g ( prec , ** kwargs ) : <EOL> memo_prec = f . memo_prec <EOL> if prec <= memo_prec : <EOL> return f . memo_val >> ( memo_prec - prec ) <EOL> newprec = int ( prec * <NUM_LIT> + <NUM_LIT:10> ) <EOL> f . memo_val = f ( newprec , ** kwargs ) <EOL> f . memo_prec = newprec <EOL> return f . memo_val >> ( newprec - prec ) <EOL> g . __name__ = f . __name__ <EOL> g . __doc__ = f . __doc__ <EOL> return g <EOL> def def_mpf_constant ( fixed ) : <EOL> """<STR_LIT>""" <EOL> def f ( prec , rnd = round_fast ) : <EOL> wp = prec + <NUM_LIT:20> <EOL> v = fixed ( wp ) <EOL> if rnd in ( round_up , round_ceiling ) : <EOL> v += <NUM_LIT:1> <EOL> return normalize ( <NUM_LIT:0> , v , - wp , bitcount ( v ) , prec , rnd ) <EOL> f . __doc__ = fixed . __doc__ <EOL> return f <EOL> def bsp_acot ( q , a , b , hyperbolic ) : <EOL> if b - a == <NUM_LIT:1> : <EOL> a1 = MPZ ( <NUM_LIT:2> * a + <NUM_LIT:3> ) <EOL> if hyperbolic or a & <NUM_LIT:1> : <EOL> return MPZ_ONE , a1 * q ** <NUM_LIT:2> , a1 <EOL> else : <EOL> return - MPZ_ONE , a1 * q ** <NUM_LIT:2> , a1 <EOL> m = ( a + b ) // <NUM_LIT:2> <EOL> p1 , q1 , r1 = bsp_acot ( q , a , m , hyperbolic ) <EOL> p2 , q2 , r2 = bsp_acot ( q , m , b , hyperbolic ) <EOL> return q2 * p1 + r1 * p2 , q1 * q2 , r1 * r2 <EOL> def acot_fixed ( a , prec , hyperbolic ) : <EOL> """<STR_LIT>""" <EOL> N = int ( <NUM_LIT> * prec / math . log ( a ) + <NUM_LIT:20> ) <EOL> p , q , r = bsp_acot ( a , <NUM_LIT:0> , N , hyperbolic ) <EOL> return ( ( p + q ) << prec ) // ( q * a ) <EOL> def machin ( coefs , prec , hyperbolic = False ) : <EOL> """<STR_LIT>""" <EOL> extraprec = <NUM_LIT:10> <EOL> s = MPZ_ZERO <EOL> for a , b in coefs : <EOL> s += MPZ ( a ) * acot_fixed ( MPZ ( b ) , prec + extraprec , hyperbolic ) <EOL> return ( s >> extraprec ) <EOL> @ constant_memo <EOL> def ln2_fixed ( prec ) : <EOL> """<STR_LIT>""" <EOL> return machin ( [ ( <NUM_LIT> , <NUM_LIT> ) , ( - <NUM_LIT:2> , <NUM_LIT> ) , ( <NUM_LIT:8> , <NUM_LIT> ) ] , prec , True ) <EOL> @ constant_memo <EOL> def ln10_fixed ( prec ) : <EOL> """<STR_LIT>""" <EOL> return machin ( [ ( <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT:20> , <NUM_LIT> ) ] , prec , True ) <EOL> """<STR_LIT>""" <EOL> CHUD_A = MPZ ( <NUM_LIT> ) <EOL> CHUD_B = MPZ ( <NUM_LIT> ) <EOL> CHUD_C = MPZ ( <NUM_LIT> ) <EOL> CHUD_D = MPZ ( <NUM_LIT:12> ) <EOL> def bs_chudnovsky ( a , b , level , verbose ) : <EOL> """<STR_LIT>""" <EOL> if b - a == <NUM_LIT:1> : <EOL> g = MPZ ( ( <NUM_LIT:6> * b - <NUM_LIT:5> ) * ( <NUM_LIT:2> * b - <NUM_LIT:1> ) * ( <NUM_LIT:6> * b - <NUM_LIT:1> ) ) <EOL> p = b ** <NUM_LIT:3> * CHUD_C ** <NUM_LIT:3> // <NUM_LIT> <EOL> q = ( - <NUM_LIT:1> ) ** b * g * ( CHUD_A + CHUD_B * b ) <EOL> else : <EOL> if verbose and level < <NUM_LIT:4> : <EOL> print ( "<STR_LIT>" , a , b ) <EOL> mid = ( a + b ) // <NUM_LIT:2> <EOL> g1 , p1 , q1 = bs_chudnovsky ( a , mid , level + <NUM_LIT:1> , verbose ) <EOL> g2 , p2 , q2 = bs_chudnovsky ( mid , b , level + <NUM_LIT:1> , verbose ) <EOL> p = p1 * p2 <EOL> g = g1 * g2 <EOL> q = q1 * p2 + q2 * g1 <EOL> return g , p , q <EOL> @ constant_memo <EOL> def pi_fixed ( prec , verbose = False , verbose_base = None ) : <EOL> """<STR_LIT>""" <EOL> N = int ( prec / <NUM_LIT> / <NUM_LIT> + <NUM_LIT:2> ) <EOL> if verbose : <EOL> print ( "<STR_LIT>" , N ) <EOL> g , p , q = bs_chudnovsky ( <NUM_LIT:0> , N , <NUM_LIT:0> , verbose ) <EOL> sqrtC = isqrt_fast ( CHUD_C << ( <NUM_LIT:2> * prec ) ) <EOL> v = p * CHUD_C * sqrtC // ( ( q + CHUD_A * p ) * CHUD_D ) <EOL> return v <EOL> def degree_fixed ( prec ) : <EOL> return pi_fixed ( prec ) // <NUM_LIT> <EOL> def bspe ( a , b ) : <EOL> """<STR_LIT>""" <EOL> if b - a == <NUM_LIT:1> : <EOL> return MPZ_ONE , MPZ ( b ) <EOL> m = ( a + b ) // <NUM_LIT:2> <EOL> p1 , q1 = bspe ( a , m ) <EOL> p2 , q2 = bspe ( m , b ) <EOL> return p1 * q2 + p2 , q1 * q2 <EOL> @ constant_memo <EOL> def e_fixed ( prec ) : <EOL> """<STR_LIT>""" <EOL> N = int ( <NUM_LIT> * prec / math . log ( prec ) + <NUM_LIT:20> ) <EOL> p , q = bspe ( <NUM_LIT:0> , N ) <EOL> return ( ( p + q ) << prec ) // q <EOL> @ constant_memo <EOL> def phi_fixed ( prec ) : <EOL> """<STR_LIT>""" <EOL> prec += <NUM_LIT:10> <EOL> a = isqrt_fast ( MPZ_FIVE << ( <NUM_LIT:2> * prec ) ) + ( MPZ_ONE << prec ) <EOL> return a >> <NUM_LIT:11> <EOL> mpf_phi = def_mpf_constant ( phi_fixed ) <EOL> mpf_pi = def_mpf_constant ( pi_fixed ) <EOL> mpf_e = def_mpf_constant ( e_fixed ) <EOL> mpf_degree = def_mpf_constant ( degree_fixed ) <EOL> mpf_ln2 = def_mpf_constant ( ln2_fixed ) <EOL> mpf_ln10 = def_mpf_constant ( ln10_fixed ) <EOL> @ constant_memo <EOL> def ln_sqrt2pi_fixed ( prec ) : <EOL> wp = prec + <NUM_LIT:10> <EOL> return to_fixed ( mpf_log ( mpf_shift ( mpf_pi ( wp ) , <NUM_LIT:1> ) , wp ) , prec - <NUM_LIT:1> ) <EOL> @ constant_memo <EOL> def sqrtpi_fixed ( prec ) : <EOL> return sqrt_fixed ( pi_fixed ( prec ) , prec ) <EOL> mpf_sqrtpi = def_mpf_constant ( sqrtpi_fixed ) <EOL> mpf_ln_sqrt2pi = def_mpf_constant ( ln_sqrt2pi_fixed ) <EOL> def mpf_pow ( s , t , prec , rnd = round_fast ) : <EOL> """<STR_LIT>""" <EOL> ssign , sman , sexp , sbc = s <EOL> tsign , tman , texp , tbc = t <EOL> if ssign and texp < <NUM_LIT:0> : <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> if texp >= <NUM_LIT:0> : <EOL> return mpf_pow_int ( s , ( - <NUM_LIT:1> ) ** tsign * ( tman << texp ) , prec , rnd ) <EOL> if texp == - <NUM_LIT:1> : <EOL> if tman == <NUM_LIT:1> : <EOL> if tsign : <EOL> return mpf_div ( fone , mpf_sqrt ( s , prec + <NUM_LIT:10> , <EOL> reciprocal_rnd [ rnd ] ) , prec , rnd ) <EOL> return mpf_sqrt ( s , prec , rnd ) <EOL> else : <EOL> if tsign : <EOL> return mpf_pow_int ( mpf_sqrt ( s , prec + <NUM_LIT:10> , <EOL> reciprocal_rnd [ rnd ] ) , - tman , prec , rnd ) <EOL> return mpf_pow_int ( mpf_sqrt ( s , prec + <NUM_LIT:10> , rnd ) , tman , prec , rnd ) <EOL> c = mpf_log ( s , prec + <NUM_LIT:10> , rnd ) <EOL> return mpf_exp ( mpf_mul ( t , c ) , prec , rnd ) <EOL> def int_pow_fixed ( y , n , prec ) : <EOL> """<STR_LIT>""" <EOL> if n == <NUM_LIT:2> : <EOL> return ( y * y ) , <NUM_LIT:0> <EOL> bc = bitcount ( y ) <EOL> exp = <NUM_LIT:0> <EOL> workprec = <NUM_LIT:2> * ( prec + <NUM_LIT:4> * bitcount ( n ) + <NUM_LIT:4> ) <EOL> _ , pm , pe , pbc = fone <EOL> while <NUM_LIT:1> : <EOL> if n & <NUM_LIT:1> : <EOL> pm = pm * y <EOL> pe = pe + exp <EOL> pbc += bc - <NUM_LIT:2> <EOL> pbc = pbc + bctable [ int ( pm >> pbc ) ] <EOL> if pbc > workprec : <EOL> pm = pm >> ( pbc - workprec ) <EOL> pe += pbc - workprec <EOL> pbc = workprec <EOL> n -= <NUM_LIT:1> <EOL> if not n : <EOL> break <EOL> y = y * y <EOL> exp = exp + exp <EOL> bc = bc + bc - <NUM_LIT:2> <EOL> bc = bc + bctable [ int ( y >> bc ) ] <EOL> if bc > workprec : <EOL> y = y >> ( bc - workprec ) <EOL> exp += bc - workprec <EOL> bc = workprec <EOL> n = n // <NUM_LIT:2> <EOL> return pm , pe <EOL> def nthroot_fixed ( y , n , prec , exp1 ) : <EOL> start = <NUM_LIT:50> <EOL> try : <EOL> y1 = rshift ( y , prec - n * start ) <EOL> r = MPZ ( int ( y1 ** ( <NUM_LIT:1.0> / n ) ) ) <EOL> except OverflowError : <EOL> y1 = from_int ( y1 , start ) <EOL> fn = from_int ( n ) <EOL> fn = mpf_rdiv_int ( <NUM_LIT:1> , fn , start ) <EOL> r = mpf_pow ( y1 , fn , start ) <EOL> r = to_int ( r ) <EOL> extra = <NUM_LIT:10> <EOL> extra1 = n <EOL> prevp = start <EOL> for p in giant_steps ( start , prec + extra ) : <EOL> pm , pe = int_pow_fixed ( r , n - <NUM_LIT:1> , prevp ) <EOL> r2 = rshift ( pm , ( n - <NUM_LIT:1> ) * prevp - p - pe - extra1 ) <EOL> B = lshift ( y , <NUM_LIT:2> * p - prec + extra1 ) // r2 <EOL> r = ( B + ( n - <NUM_LIT:1> ) * lshift ( r , p - prevp ) ) // n <EOL> prevp = p <EOL> return r <EOL> def mpf_nthroot ( s , n , prec , rnd = round_fast ) : <EOL> """<STR_LIT>""" <EOL> sign , man , exp , bc = s <EOL> if sign : <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> if not man : <EOL> if s == fnan : <EOL> return fnan <EOL> if s == fzero : <EOL> if n > <NUM_LIT:0> : <EOL> return fzero <EOL> if n == <NUM_LIT:0> : <EOL> return fone <EOL> return finf <EOL> if not n : <EOL> return fnan <EOL> if n < <NUM_LIT:0> : <EOL> return fzero <EOL> return finf <EOL> flag_inverse = False <EOL> if n < <NUM_LIT:2> : <EOL> if n == <NUM_LIT:0> : <EOL> return fone <EOL> if n == <NUM_LIT:1> : <EOL> return mpf_pos ( s , prec , rnd ) <EOL> if n == - <NUM_LIT:1> : <EOL> return mpf_div ( fone , s , prec , rnd ) <EOL> rnd = reciprocal_rnd [ rnd ] <EOL> flag_inverse = True <EOL> extra_inverse = <NUM_LIT:5> <EOL> prec += extra_inverse <EOL> n = - n <EOL> if n > <NUM_LIT:20> and ( n >= <NUM_LIT> or prec < int ( <NUM_LIT> + <NUM_LIT> * n ** <NUM_LIT> ) ) : <EOL> prec2 = prec + <NUM_LIT:10> <EOL> fn = from_int ( n ) <EOL> nth = mpf_rdiv_int ( <NUM_LIT:1> , fn , prec2 ) <EOL> r = mpf_pow ( s , nth , prec2 , rnd ) <EOL> s = normalize ( r [ <NUM_LIT:0> ] , r [ <NUM_LIT:1> ] , r [ <NUM_LIT:2> ] , r [ <NUM_LIT:3> ] , prec , rnd ) <EOL> if flag_inverse : <EOL> return mpf_div ( fone , s , prec - extra_inverse , rnd ) <EOL> else : <EOL> return s <EOL> prec2 = prec + <NUM_LIT:2> * n - ( prec % n ) <EOL> if n > <NUM_LIT:10> : <EOL> prec2 += prec2 // <NUM_LIT:10> <EOL> prec2 = prec2 - prec2 % n <EOL> shift = bc - prec2 <EOL> sign1 = <NUM_LIT:0> <EOL> es = exp + shift <EOL> if es < <NUM_LIT:0> : <EOL> sign1 = <NUM_LIT:1> <EOL> es = - es <EOL> if sign1 : <EOL> shift += es % n <EOL> else : <EOL> shift -= es % n <EOL> man = rshift ( man , shift ) <EOL> extra = <NUM_LIT:10> <EOL> exp1 = ( ( exp + shift - ( n - <NUM_LIT:1> ) * prec2 ) // n ) - extra <EOL> rnd_shift = <NUM_LIT:0> <EOL> if flag_inverse : <EOL> if rnd == '<STR_LIT:u>' or rnd == '<STR_LIT:c>' : <EOL> rnd_shift = <NUM_LIT:1> <EOL> else : <EOL> if rnd == '<STR_LIT:d>' or rnd == '<STR_LIT:f>' : <EOL> rnd_shift = <NUM_LIT:1> <EOL> man = nthroot_fixed ( man + rnd_shift , n , prec2 , exp1 ) <EOL> s = from_man_exp ( man , exp1 , prec , rnd ) <EOL> if flag_inverse : <EOL> return mpf_div ( fone , s , prec - extra_inverse , rnd ) <EOL> else : <EOL> return s <EOL> def mpf_cbrt ( s , prec , rnd = round_fast ) : <EOL> """<STR_LIT>""" <EOL> return mpf_nthroot ( s , <NUM_LIT:3> , prec , rnd ) <EOL> def log_int_fixed ( n , prec , ln2 = None ) : <EOL> """<STR_LIT>""" <EOL> if n in log_int_cache : <EOL> value , vprec = log_int_cache [ n ] <EOL> if vprec >= prec : <EOL> return value >> ( vprec - prec ) <EOL> wp = prec + <NUM_LIT:10> <EOL> if wp <= LOG_TAYLOR_SHIFT : <EOL> if ln2 is None : <EOL> ln2 = ln2_fixed ( wp ) <EOL> r = bitcount ( n ) <EOL> x = n << ( wp - r ) <EOL> v = log_taylor_cached ( x , wp ) + r * ln2 <EOL> else : <EOL> v = to_fixed ( mpf_log ( from_int ( n ) , wp + <NUM_LIT:5> ) , wp ) <EOL> if n < MAX_LOG_INT_CACHE : <EOL> log_int_cache [ n ] = ( v , wp ) <EOL> return v >> ( wp - prec ) <EOL> def agm_fixed ( a , b , prec ) : <EOL> """<STR_LIT>""" <EOL> i = <NUM_LIT:0> <EOL> while <NUM_LIT:1> : <EOL> anew = ( a + b ) >> <NUM_LIT:1> <EOL> if i > <NUM_LIT:4> and abs ( a - anew ) < <NUM_LIT:8> : <EOL> return a <EOL> b = isqrt_fast ( a * b ) <EOL> a = anew <EOL> i += <NUM_LIT:1> <EOL> return a <EOL> def log_agm ( x , prec ) : <EOL> """<STR_LIT>""" <EOL> x2 = ( x * x ) >> prec <EOL> s = a = b = x2 <EOL> while a : <EOL> b = ( b * x2 ) >> prec <EOL> a = ( a * b ) >> prec <EOL> s += a <EOL> s += ( MPZ_ONE << prec ) <EOL> s = ( s * s ) >> ( prec - <NUM_LIT:2> ) <EOL> s = ( s * isqrt_fast ( x << prec ) ) >> prec <EOL> t = a = b = x <EOL> while a : <EOL> b = ( b * x2 ) >> prec <EOL> a = ( a * b ) >> prec <EOL> t += a <EOL> t = ( MPZ_ONE << prec ) + ( t << <NUM_LIT:1> ) <EOL> t = ( t * t ) >> prec <EOL> p = agm_fixed ( s , t , prec ) <EOL> return ( pi_fixed ( prec ) << prec ) // p <EOL> def log_taylor ( x , prec , r = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> for i in xrange ( r ) : <EOL> x = isqrt_fast ( x << prec ) <EOL> one = MPZ_ONE << prec <EOL> v = ( ( x - one ) << prec ) // ( x + one ) <EOL> sign = v < <NUM_LIT:0> <EOL> if sign : <EOL> v = - v <EOL> v2 = ( v * v ) >> prec <EOL> v4 = ( v2 * v2 ) >> prec <EOL> s0 = v <EOL> s1 = v // <NUM_LIT:3> <EOL> v = ( v * v4 ) >> prec <EOL> k = <NUM_LIT:5> <EOL> while v : <EOL> s0 += v // k <EOL> k += <NUM_LIT:2> <EOL> s1 += v // k <EOL> v = ( v * v4 ) >> prec <EOL> k += <NUM_LIT:2> <EOL> s1 = ( s1 * v2 ) >> prec <EOL> s = ( s0 + s1 ) << ( <NUM_LIT:1> + r ) <EOL> if sign : <EOL> return - s <EOL> return s <EOL> def log_taylor_cached ( x , prec ) : <EOL> """<STR_LIT>""" <EOL> n = x >> ( prec - LOG_TAYLOR_SHIFT ) <EOL> cached_prec = cache_prec_steps [ prec ] <EOL> dprec = cached_prec - prec <EOL> if ( n , cached_prec ) in log_taylor_cache : <EOL> a , log_a = log_taylor_cache [ n , cached_prec ] <EOL> else : <EOL> a = n << ( cached_prec - LOG_TAYLOR_SHIFT ) <EOL> log_a = log_taylor ( a , cached_prec , <NUM_LIT:8> ) <EOL> log_taylor_cache [ n , cached_prec ] = ( a , log_a ) <EOL> a >>= dprec <EOL> log_a >>= dprec <EOL> u = ( ( x - a ) << prec ) // a <EOL> v = ( u << prec ) // ( ( MPZ_TWO << prec ) + u ) <EOL> v2 = ( v * v ) >> prec <EOL> v4 = ( v2 * v2 ) >> prec <EOL> s0 = v <EOL> s1 = v // <NUM_LIT:3> <EOL> v = ( v * v4 ) >> prec <EOL> k = <NUM_LIT:5> <EOL> while v : <EOL> s0 += v // k <EOL> k += <NUM_LIT:2> <EOL> s1 += v // k <EOL> v = ( v * v4 ) >> prec <EOL> k += <NUM_LIT:2> <EOL> s1 = ( s1 * v2 ) >> prec <EOL> s = ( s0 + s1 ) << <NUM_LIT:1> <EOL> return log_a + s <EOL> def mpf_log ( x , prec , rnd = round_fast ) : <EOL> """<STR_LIT>""" <EOL> sign , man , exp , bc = x <EOL> if not man : <EOL> if x == fzero : return fninf <EOL> if x == finf : return finf <EOL> if x == fnan : return fnan <EOL> if sign : <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> wp = prec + <NUM_LIT:20> <EOL> if man == <NUM_LIT:1> : <EOL> if not exp : <EOL> return fzero <EOL> return from_man_exp ( exp * ln2_fixed ( wp ) , - wp , prec , rnd ) <EOL> mag = exp + bc <EOL> abs_mag = abs ( mag ) <EOL> if abs_mag <= <NUM_LIT:1> : <EOL> tsign = <NUM_LIT:1> - abs_mag <EOL> if tsign : <EOL> tman = ( MPZ_ONE << bc ) - man <EOL> else : <EOL> tman = man - ( MPZ_ONE << ( bc - <NUM_LIT:1> ) ) <EOL> tbc = bitcount ( tman ) <EOL> cancellation = bc - tbc <EOL> if cancellation > wp : <EOL> t = normalize ( tsign , tman , abs_mag - bc , tbc , tbc , '<STR_LIT:n>' ) <EOL> return mpf_perturb ( t , tsign , prec , rnd ) <EOL> else : <EOL> wp += cancellation <EOL> if abs_mag > <NUM_LIT> : <EOL> if bitcount ( abs_mag ) > wp : <EOL> return from_man_exp ( exp * ln2_fixed ( wp ) , - wp , prec , rnd ) <EOL> if wp <= LOG_TAYLOR_PREC : <EOL> m = log_taylor_cached ( lshift ( man , wp - bc ) , wp ) <EOL> if mag : <EOL> m += mag * ln2_fixed ( wp ) <EOL> else : <EOL> optimal_mag = - wp // LOG_AGM_MAG_PREC_RATIO <EOL> n = optimal_mag - mag <EOL> x = mpf_shift ( x , n ) <EOL> wp += ( - optimal_mag ) <EOL> m = - log_agm ( to_fixed ( x , wp ) , wp ) <EOL> m -= n * ln2_fixed ( wp ) <EOL> return from_man_exp ( m , - wp , prec , rnd ) <EOL> def mpf_log_hypot ( a , b , prec , rnd ) : <EOL> """<STR_LIT>""" <EOL> if not b [ <NUM_LIT:1> ] : <EOL> a , b = b , a <EOL> if not a [ <NUM_LIT:1> ] : <EOL> if not b [ <NUM_LIT:1> ] : <EOL> if a == b == fzero : <EOL> return fninf <EOL> if fnan in ( a , b ) : <EOL> return fnan <EOL> return finf <EOL> if a == fzero : <EOL> return mpf_log ( mpf_abs ( b ) , prec , rnd ) <EOL> if a == fnan : <EOL> return fnan <EOL> return finf <EOL> a2 = mpf_mul ( a , a ) <EOL> b2 = mpf_mul ( b , b ) <EOL> extra = <NUM_LIT:20> <EOL> h2 = mpf_add ( a2 , b2 , prec + extra ) <EOL> cancelled = mpf_add ( h2 , fnone , <NUM_LIT:10> ) <EOL> mag_cancelled = cancelled [ <NUM_LIT:2> ] + cancelled [ <NUM_LIT:3> ] <EOL> if cancelled == fzero or mag_cancelled < - extra // <NUM_LIT:2> : <EOL> h2 = mpf_add ( a2 , b2 , prec + extra - min ( a2 [ <NUM_LIT:2> ] , b2 [ <NUM_LIT:2> ] ) ) <EOL> return mpf_shift ( mpf_log ( h2 , prec , rnd ) , - <NUM_LIT:1> ) <EOL> def atan_newton ( x , prec ) : <EOL> if prec >= <NUM_LIT:100> : <EOL> r = math . atan ( int ( ( x >> ( prec - <NUM_LIT> ) ) ) / <NUM_LIT> ** <NUM_LIT> ) <EOL> else : <EOL> r = math . atan ( int ( x ) / <NUM_LIT> ** prec ) <EOL> prevp = <NUM_LIT:50> <EOL> r = MPZ ( int ( r * <NUM_LIT> ** <NUM_LIT> ) >> ( <NUM_LIT> - prevp ) ) <EOL> extra_p = <NUM_LIT:50> <EOL> for wp in giant_steps ( prevp , prec ) : <EOL> wp += extra_p <EOL> r = r << ( wp - prevp ) <EOL> cos , sin = cos_sin_fixed ( r , wp ) <EOL> tan = ( sin << wp ) // cos <EOL> a = ( ( tan - rshift ( x , prec - wp ) ) << wp ) // ( ( MPZ_ONE << wp ) + ( ( tan ** <NUM_LIT:2> ) >> wp ) ) <EOL> r = r - a <EOL> prevp = wp <EOL> return rshift ( r , prevp - prec ) <EOL> def atan_taylor_get_cached ( n , prec ) : <EOL> prec2 = ( <NUM_LIT:1> << ( bitcount ( prec - <NUM_LIT:1> ) ) ) + <NUM_LIT:20> <EOL> dprec = prec2 - prec <EOL> if ( n , prec2 ) in atan_taylor_cache : <EOL> a , atan_a = atan_taylor_cache [ n , prec2 ] <EOL> else : <EOL> a = n << ( prec2 - ATAN_TAYLOR_SHIFT ) <EOL> atan_a = atan_newton ( a , prec2 ) <EOL> atan_taylor_cache [ n , prec2 ] = ( a , atan_a ) <EOL> return ( a >> dprec ) , ( atan_a >> dprec ) <EOL> def atan_taylor ( x , prec ) : <EOL> n = ( x >> ( prec - ATAN_TAYLOR_SHIFT ) ) <EOL> a , atan_a = atan_taylor_get_cached ( n , prec ) <EOL> d = x - a <EOL> s0 = v = ( d << prec ) // ( ( a ** <NUM_LIT:2> >> prec ) + ( a * d >> prec ) + ( MPZ_ONE << prec ) ) <EOL> v2 = ( v ** <NUM_LIT:2> >> prec ) <EOL> v4 = ( v2 * v2 ) >> prec <EOL> s1 = v // <NUM_LIT:3> <EOL> v = ( v * v4 ) >> prec <EOL> k = <NUM_LIT:5> <EOL> while v : <EOL> s0 += v // k <EOL> k += <NUM_LIT:2> <EOL> s1 += v // k <EOL> v = ( v * v4 ) >> prec <EOL> k += <NUM_LIT:2> <EOL> s1 = ( s1 * v2 ) >> prec <EOL> s = s0 - s1 <EOL> return atan_a + s <EOL> def atan_inf ( sign , prec , rnd ) : <EOL> if not sign : <EOL> return mpf_shift ( mpf_pi ( prec , rnd ) , - <NUM_LIT:1> ) <EOL> return mpf_neg ( mpf_shift ( mpf_pi ( prec , negative_rnd [ rnd ] ) , - <NUM_LIT:1> ) ) <EOL> def mpf_atan ( x , prec , rnd = round_fast ) : <EOL> sign , man , exp , bc = x <EOL> if not man : <EOL> if x == fzero : return fzero <EOL> if x == finf : return atan_inf ( <NUM_LIT:0> , prec , rnd ) <EOL> if x == fninf : return atan_inf ( <NUM_LIT:1> , prec , rnd ) <EOL> return fnan <EOL> mag = exp + bc <EOL> if mag > prec + <NUM_LIT:20> : <EOL> return atan_inf ( sign , prec , rnd ) <EOL> if - mag > prec + <NUM_LIT:20> : <EOL> return mpf_perturb ( x , <NUM_LIT:1> - sign , prec , rnd ) <EOL> wp = prec + <NUM_LIT:30> + abs ( mag ) <EOL> if mag >= <NUM_LIT:2> : <EOL> x = mpf_rdiv_int ( <NUM_LIT:1> , x , wp ) <EOL> reciprocal = True <EOL> else : <EOL> reciprocal = False <EOL> t = to_fixed ( x , wp ) <EOL> if sign : <EOL> t = - t <EOL> if wp < ATAN_TAYLOR_PREC : <EOL> a = atan_taylor ( t , wp ) <EOL> else : <EOL> a = atan_newton ( t , wp ) <EOL> if reciprocal : <EOL> a = ( ( pi_fixed ( wp ) >> <NUM_LIT:1> ) + <NUM_LIT:1> ) - a <EOL> if sign : <EOL> a = - a <EOL> return from_man_exp ( a , - wp , prec , rnd ) <EOL> def mpf_atan2 ( y , x , prec , rnd = round_fast ) : <EOL> xsign , xman , xexp , xbc = x <EOL> ysign , yman , yexp , ybc = y <EOL> if not yman : <EOL> if y == fzero and x != fnan : <EOL> if mpf_sign ( x ) >= <NUM_LIT:0> : <EOL> return fzero <EOL> return mpf_pi ( prec , rnd ) <EOL> if y in ( finf , fninf ) : <EOL> if x in ( finf , fninf ) : <EOL> return fnan <EOL> if y == finf : <EOL> return mpf_shift ( mpf_pi ( prec , rnd ) , - <NUM_LIT:1> ) <EOL> return mpf_neg ( mpf_shift ( mpf_pi ( prec , negative_rnd [ rnd ] ) , - <NUM_LIT:1> ) ) <EOL> return fnan <EOL> if ysign : <EOL> return mpf_neg ( mpf_atan2 ( mpf_neg ( y ) , x , prec , negative_rnd [ rnd ] ) ) <EOL> if not xman : <EOL> if x == fnan : <EOL> return fnan <EOL> if x == finf : <EOL> return fzero <EOL> if x == fninf : <EOL> return mpf_pi ( prec , rnd ) <EOL> if y == fzero : <EOL> return fzero <EOL> return mpf_shift ( mpf_pi ( prec , rnd ) , - <NUM_LIT:1> ) <EOL> tquo = mpf_atan ( mpf_div ( y , x , prec + <NUM_LIT:4> ) , prec + <NUM_LIT:4> ) <EOL> if xsign : <EOL> return mpf_add ( mpf_pi ( prec + <NUM_LIT:4> ) , tquo , prec , rnd ) <EOL> else : <EOL> return mpf_pos ( tquo , prec , rnd ) <EOL> def mpf_asin ( x , prec , rnd = round_fast ) : <EOL> sign , man , exp , bc = x <EOL> if bc + exp > <NUM_LIT:0> and x not in ( fone , fnone ) : <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> wp = prec + <NUM_LIT:15> <EOL> a = mpf_mul ( x , x ) <EOL> b = mpf_add ( fone , mpf_sqrt ( mpf_sub ( fone , a , wp ) , wp ) , wp ) <EOL> c = mpf_div ( x , b , wp ) <EOL> return mpf_shift ( mpf_atan ( c , prec , rnd ) , <NUM_LIT:1> ) <EOL> def mpf_acos ( x , prec , rnd = round_fast ) : <EOL> sign , man , exp , bc = x <EOL> if bc + exp > <NUM_LIT:0> : <EOL> if x not in ( fone , fnone ) : <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> if x == fnone : <EOL> return mpf_pi ( prec , rnd ) <EOL> wp = prec + <NUM_LIT:15> <EOL> a = mpf_mul ( x , x ) <EOL> b = mpf_sqrt ( mpf_sub ( fone , a , wp ) , wp ) <EOL> c = mpf_div ( b , mpf_add ( fone , x , wp ) , wp ) <EOL> return mpf_shift ( mpf_atan ( c , prec , rnd ) , <NUM_LIT:1> ) <EOL> def mpf_asinh ( x , prec , rnd = round_fast ) : <EOL> wp = prec + <NUM_LIT:20> <EOL> sign , man , exp , bc = x <EOL> mag = exp + bc <EOL> if mag < - <NUM_LIT:8> : <EOL> if mag < - wp : <EOL> return mpf_perturb ( x , <NUM_LIT:1> - sign , prec , rnd ) <EOL> wp += ( - mag ) <EOL> q = mpf_sqrt ( mpf_add ( mpf_mul ( x , x ) , fone , wp ) , wp ) <EOL> q = mpf_add ( mpf_abs ( x ) , q , wp ) <EOL> if sign : <EOL> return mpf_neg ( mpf_log ( q , prec , negative_rnd [ rnd ] ) ) <EOL> else : <EOL> return mpf_log ( q , prec , rnd ) <EOL> def mpf_acosh ( x , prec , rnd = round_fast ) : <EOL> wp = prec + <NUM_LIT:15> <EOL> if mpf_cmp ( x , fone ) == - <NUM_LIT:1> : <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> q = mpf_sqrt ( mpf_add ( mpf_mul ( x , x ) , fnone , wp ) , wp ) <EOL> return mpf_log ( mpf_add ( x , q , wp ) , prec , rnd ) <EOL> def mpf_atanh ( x , prec , rnd = round_fast ) : <EOL> sign , man , exp , bc = x <EOL> if ( not man ) and exp : <EOL> if x in ( fzero , fnan ) : <EOL> return x <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> mag = bc + exp <EOL> if mag > <NUM_LIT:0> : <EOL> if mag == <NUM_LIT:1> and man == <NUM_LIT:1> : <EOL> return [ finf , fninf ] [ sign ] <EOL> raise ComplexResult ( "<STR_LIT>" ) <EOL> wp = prec + <NUM_LIT:15> <EOL> if mag < - <NUM_LIT:8> : <EOL> if mag < - wp : <EOL> return mpf_perturb ( x , sign , prec , rnd ) <EOL> wp += ( - mag ) <EOL> a = mpf_add ( x , fone , wp ) <EOL> b = mpf_sub ( fone , x , wp ) <EOL> return mpf_shift ( mpf_log ( mpf_div ( a , b , wp ) , prec , rnd ) , - <NUM_LIT:1> ) <EOL> def mpf_fibonacci ( x , prec , rnd = round_fast ) : <EOL> sign , man , exp , bc = x <EOL> if not man : <EOL> if x == fninf : <EOL> return fnan <EOL> return x <EOL> size = abs ( exp + bc ) <EOL> if exp >= <NUM_LIT:0> : <EOL> if size < <NUM_LIT:10> or size <= bitcount ( prec ) : <EOL> return from_int ( ifib ( to_int ( x ) ) , prec , rnd ) <EOL> wp = prec + size + <NUM_LIT:20> <EOL> a = mpf_phi ( wp ) <EOL> b = mpf_add ( mpf_shift ( a , <NUM_LIT:1> ) , fnone , wp ) <EOL> u = mpf_pow ( a , x , wp ) <EOL> v = mpf_cos_pi ( x , wp ) <EOL> v = mpf_div ( v , u , wp ) <EOL> u = mpf_sub ( u , v , wp ) <EOL> u = mpf_div ( u , b , prec , rnd ) <EOL> return u <EOL> def exponential_series ( x , prec , type = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> if x < <NUM_LIT:0> : <EOL> x = - x <EOL> sign = <NUM_LIT:1> <EOL> else : <EOL> sign = <NUM_LIT:0> <EOL> r = int ( <NUM_LIT:0.5> * prec ** <NUM_LIT:0.5> ) <EOL> xmag = bitcount ( x ) - prec <EOL> r = max ( <NUM_LIT:0> , xmag + r ) <EOL> extra = <NUM_LIT:10> + <NUM_LIT:2> * max ( r , - xmag ) <EOL> wp = prec + extra <EOL> x <<= ( extra - r ) <EOL> one = MPZ_ONE << wp <EOL> alt = ( type == <NUM_LIT:2> ) <EOL> if prec < EXP_SERIES_U_CUTOFF : <EOL> x2 = a = ( x * x ) >> wp <EOL> x4 = ( x2 * x2 ) >> wp <EOL> s0 = s1 = MPZ_ZERO <EOL> k = <NUM_LIT:2> <EOL> while a : <EOL> a //= ( k - <NUM_LIT:1> ) * k ; s0 += a ; k += <NUM_LIT:2> <EOL> a //= ( k - <NUM_LIT:1> ) * k ; s1 += a ; k += <NUM_LIT:2> <EOL> a = ( a * x4 ) >> wp <EOL> s1 = ( x2 * s1 ) >> wp <EOL> if alt : <EOL> c = s1 - s0 + one <EOL> else : <EOL> c = s1 + s0 + one <EOL> else : <EOL> u = int ( <NUM_LIT> * prec ** <NUM_LIT> ) <EOL> x2 = a = ( x * x ) >> wp <EOL> xpowers = [ one , x2 ] <EOL> for i in xrange ( <NUM_LIT:1> , u ) : <EOL> xpowers . append ( ( xpowers [ - <NUM_LIT:1> ] * x2 ) >> wp ) <EOL> sums = [ MPZ_ZERO ] * u <EOL> k = <NUM_LIT:2> <EOL> while a : <EOL> for i in xrange ( u ) : <EOL> a //= ( k - <NUM_LIT:1> ) * k <EOL> if alt and k & <NUM_LIT:2> : sums [ i ] -= a <EOL> else : sums [ i ] += a <EOL> k += <NUM_LIT:2> <EOL> a = ( a * xpowers [ - <NUM_LIT:1> ] ) >> wp <EOL> for i in xrange ( <NUM_LIT:1> , u ) : <EOL> sums [ i ] = ( sums [ i ] * xpowers [ i ] ) >> wp <EOL> c = sum ( sums ) + one <EOL> if type == <NUM_LIT:0> : <EOL> s = isqrt_fast ( c * c - ( one << wp ) ) <EOL> if sign : <EOL> v = c - s <EOL> else : <EOL> v = c + s <EOL> for i in xrange ( r ) : <EOL> v = ( v * v ) >> wp <EOL> return v >> extra <EOL> else : <EOL> pshift = wp - <NUM_LIT:1> <EOL> for i in xrange ( r ) : <EOL> c = ( ( c * c ) >> pshift ) - one <EOL> s = isqrt_fast ( abs ( ( one << wp ) - c * c ) ) <EOL> if sign : <EOL> s = - s <EOL> return ( c >> extra ) , ( s >> extra ) <EOL> def exp_basecase ( x , prec ) : <EOL> """<STR_LIT>""" <EOL> if prec > EXP_COSH_CUTOFF : <EOL> return exponential_series ( x , prec , <NUM_LIT:0> ) <EOL> r = int ( prec ** <NUM_LIT:0.5> ) <EOL> prec += r <EOL> s0 = s1 = ( MPZ_ONE << prec ) <EOL> k = <NUM_LIT:2> <EOL> a = x2 = ( x * x ) >> prec <EOL> while a : <EOL> a //= k ; s0 += a ; k += <NUM_LIT:1> <EOL> a //= k ; s1 += a ; k += <NUM_LIT:1> <EOL> a = ( a * x2 ) >> prec <EOL> s1 = ( s1 * x ) >> prec <EOL> s = s0 + s1 <EOL> u = r <EOL> while r : <EOL> s = ( s * s ) >> prec <EOL> r -= <NUM_LIT:1> <EOL> return s >> u <EOL> def exp_expneg_basecase ( x , prec ) : <EOL> """<STR_LIT>""" <EOL> if prec > EXP_COSH_CUTOFF : <EOL> cosh , sinh = exponential_series ( x , prec , <NUM_LIT:1> ) <EOL> return cosh + sinh , cosh - sinh <EOL> a = exp_basecase ( x , prec ) <EOL> b = ( MPZ_ONE << ( prec + prec ) ) // a <EOL> return a , b <EOL> def cos_sin_basecase ( x , prec ) : <EOL> """<STR_LIT>""" <EOL> if prec > COS_SIN_CACHE_PREC : <EOL> return exponential_series ( x , prec , <NUM_LIT:2> ) <EOL> precs = prec - COS_SIN_CACHE_STEP <EOL> t = x >> precs <EOL> n = int ( t ) <EOL> if n not in cos_sin_cache : <EOL> w = t << ( <NUM_LIT:10> + COS_SIN_CACHE_PREC - COS_SIN_CACHE_STEP ) <EOL> cos_t , sin_t = exponential_series ( w , <NUM_LIT:10> + COS_SIN_CACHE_PREC , <NUM_LIT:2> ) <EOL> cos_sin_cache [ n ] = ( cos_t >> <NUM_LIT:10> ) , ( sin_t >> <NUM_LIT:10> ) <EOL> cos_t , sin_t = cos_sin_cache [ n ] <EOL> offset = COS_SIN_CACHE_PREC - prec <EOL> cos_t >>= offset <EOL> sin_t >>= offset <EOL> x -= t << precs <EOL> cos = MPZ_ONE << prec <EOL> sin = x <EOL> k = <NUM_LIT:2> <EOL> a = - ( ( x * x ) >> prec ) <EOL> while a : <EOL> a //= k ; cos += a ; k += <NUM_LIT:1> ; a = ( a * x ) >> prec <EOL> a //= k ; sin += a ; k += <NUM_LIT:1> ; a = - ( ( a * x ) >> prec ) <EOL> return ( ( cos * cos_t - sin * sin_t ) >> prec ) , ( ( sin * cos_t + cos * sin_t ) >> prec ) <EOL> def mpf_exp ( x , prec , rnd = round_fast ) : <EOL> sign , man , exp , bc = x <EOL> if man : <EOL> mag = bc + exp <EOL> wp = prec + <NUM_LIT> <EOL> if sign : <EOL> man = - man <EOL> if prec > <NUM_LIT> and exp >= <NUM_LIT:0> : <EOL> e = mpf_e ( wp + int ( <NUM_LIT> * mag ) ) <EOL> return mpf_pow_int ( e , man << exp , prec , rnd ) <EOL> if mag < - wp : <EOL> return mpf_perturb ( fone , sign , prec , rnd ) <EOL> if mag > <NUM_LIT:1> : <EOL> wpmod = wp + mag <EOL> offset = exp + wpmod <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> lg2 = ln2_fixed ( wpmod ) <EOL> n , t = divmod ( t , lg2 ) <EOL> n = int ( n ) <EOL> t >>= mag <EOL> else : <EOL> offset = exp + wp <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> n = <NUM_LIT:0> <EOL> man = exp_basecase ( t , wp ) <EOL> return from_man_exp ( man , n - wp , prec , rnd ) <EOL> if not exp : <EOL> return fone <EOL> if x == fninf : <EOL> return fzero <EOL> return x <EOL> def mpf_cosh_sinh ( x , prec , rnd = round_fast , tanh = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> sign , man , exp , bc = x <EOL> if ( not man ) and exp : <EOL> if tanh : <EOL> if x == finf : return fone <EOL> if x == fninf : return fnone <EOL> return fnan <EOL> if x == finf : return ( finf , finf ) <EOL> if x == fninf : return ( finf , fninf ) <EOL> return fnan , fnan <EOL> mag = exp + bc <EOL> wp = prec + <NUM_LIT> <EOL> if mag < - <NUM_LIT:4> : <EOL> if mag < - wp : <EOL> if tanh : <EOL> return mpf_perturb ( x , <NUM_LIT:1> - sign , prec , rnd ) <EOL> cosh = mpf_perturb ( fone , <NUM_LIT:0> , prec , rnd ) <EOL> sinh = mpf_perturb ( x , sign , prec , rnd ) <EOL> return cosh , sinh <EOL> wp += ( - mag ) <EOL> if mag > <NUM_LIT:10> : <EOL> if <NUM_LIT:3> * ( <NUM_LIT:1> << ( mag - <NUM_LIT:1> ) ) > wp : <EOL> if tanh : <EOL> return mpf_perturb ( [ fone , fnone ] [ sign ] , <NUM_LIT:1> - sign , prec , rnd ) <EOL> c = s = mpf_shift ( mpf_exp ( mpf_abs ( x ) , prec , rnd ) , - <NUM_LIT:1> ) <EOL> if sign : <EOL> s = mpf_neg ( s ) <EOL> return c , s <EOL> if mag > <NUM_LIT:1> : <EOL> wpmod = wp + mag <EOL> offset = exp + wpmod <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> lg2 = ln2_fixed ( wpmod ) <EOL> n , t = divmod ( t , lg2 ) <EOL> n = int ( n ) <EOL> t >>= mag <EOL> else : <EOL> offset = exp + wp <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> n = <NUM_LIT:0> <EOL> a , b = exp_expneg_basecase ( t , wp ) <EOL> cosh = a + ( b >> ( <NUM_LIT:2> * n ) ) <EOL> sinh = a - ( b >> ( <NUM_LIT:2> * n ) ) <EOL> if sign : <EOL> sinh = - sinh <EOL> if tanh : <EOL> man = ( sinh << wp ) // cosh <EOL> return from_man_exp ( man , - wp , prec , rnd ) <EOL> else : <EOL> cosh = from_man_exp ( cosh , n - wp - <NUM_LIT:1> , prec , rnd ) <EOL> sinh = from_man_exp ( sinh , n - wp - <NUM_LIT:1> , prec , rnd ) <EOL> return cosh , sinh <EOL> def mod_pi2 ( man , exp , mag , wp ) : <EOL> if mag > <NUM_LIT:0> : <EOL> i = <NUM_LIT:0> <EOL> while <NUM_LIT:1> : <EOL> cancellation_prec = <NUM_LIT:20> << i <EOL> wpmod = wp + mag + cancellation_prec <EOL> pi2 = pi_fixed ( wpmod - <NUM_LIT:1> ) <EOL> pi4 = pi2 >> <NUM_LIT:1> <EOL> offset = wpmod + exp <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> n , y = divmod ( t , pi2 ) <EOL> if y > pi4 : <EOL> small = pi2 - y <EOL> else : <EOL> small = y <EOL> if small >> ( wp + mag - <NUM_LIT:10> ) : <EOL> n = int ( n ) <EOL> t = y >> mag <EOL> wp = wpmod - mag <EOL> break <EOL> i += <NUM_LIT:1> <EOL> else : <EOL> wp += ( - mag ) <EOL> offset = exp + wp <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> n = <NUM_LIT:0> <EOL> return t , n , wp <EOL> def mpf_cos_sin ( x , prec , rnd = round_fast , which = <NUM_LIT:0> , pi = False ) : <EOL> """<STR_LIT>""" <EOL> sign , man , exp , bc = x <EOL> if not man : <EOL> if exp : <EOL> c , s = fnan , fnan <EOL> else : <EOL> c , s = fone , fzero <EOL> if which == <NUM_LIT:0> : return c , s <EOL> if which == <NUM_LIT:1> : return c <EOL> if which == <NUM_LIT:2> : return s <EOL> if which == <NUM_LIT:3> : return s <EOL> mag = bc + exp <EOL> wp = prec + <NUM_LIT:10> <EOL> if mag < <NUM_LIT:0> : <EOL> if mag < - wp : <EOL> if pi : <EOL> x = mpf_mul ( x , mpf_pi ( wp ) ) <EOL> c = mpf_perturb ( fone , <NUM_LIT:1> , prec , rnd ) <EOL> s = mpf_perturb ( x , <NUM_LIT:1> - sign , prec , rnd ) <EOL> if which == <NUM_LIT:0> : return c , s <EOL> if which == <NUM_LIT:1> : return c <EOL> if which == <NUM_LIT:2> : return s <EOL> if which == <NUM_LIT:3> : return mpf_perturb ( x , sign , prec , rnd ) <EOL> if pi : <EOL> if exp >= - <NUM_LIT:1> : <EOL> if exp == - <NUM_LIT:1> : <EOL> c = fzero <EOL> s = ( fone , fnone ) [ bool ( man & <NUM_LIT:2> ) ^ sign ] <EOL> elif exp == <NUM_LIT:0> : <EOL> c , s = ( fnone , fzero ) <EOL> else : <EOL> c , s = ( fone , fzero ) <EOL> if which == <NUM_LIT:0> : return c , s <EOL> if which == <NUM_LIT:1> : return c <EOL> if which == <NUM_LIT:2> : return s <EOL> if which == <NUM_LIT:3> : return mpf_div ( s , c , prec , rnd ) <EOL> n = ( ( man >> ( - exp - <NUM_LIT:2> ) ) + <NUM_LIT:1> ) >> <NUM_LIT:1> <EOL> man = man - ( n << ( - exp - <NUM_LIT:1> ) ) <EOL> mag2 = bitcount ( man ) + exp <EOL> wp = prec + <NUM_LIT:10> - mag2 <EOL> offset = exp + wp <EOL> if offset >= <NUM_LIT:0> : <EOL> t = man << offset <EOL> else : <EOL> t = man >> ( - offset ) <EOL> t = ( t * pi_fixed ( wp ) ) >> wp <EOL> else : <EOL> t , n , wp = mod_pi2 ( man , exp , mag , wp ) <EOL> c , s = cos_sin_basecase ( t , wp ) <EOL> m = n & <NUM_LIT:3> <EOL> if m == <NUM_LIT:1> : c , s = - s , c <EOL> elif m == <NUM_LIT:2> : c , s = - c , - s <EOL> elif m == <NUM_LIT:3> : c , s = s , - c <EOL> if sign : <EOL> s = - s <EOL> if which == <NUM_LIT:0> : <EOL> c = from_man_exp ( c , - wp , prec , rnd ) <EOL> s = from_man_exp ( s , - wp , prec , rnd ) <EOL> return c , s <EOL> if which == <NUM_LIT:1> : <EOL> return from_man_exp ( c , - wp , prec , rnd ) <EOL> if which == <NUM_LIT:2> : <EOL> return from_man_exp ( s , - wp , prec , rnd ) <EOL> if which == <NUM_LIT:3> : <EOL> return from_rational ( s , c , prec , rnd ) <EOL> def mpf_cos ( x , prec , rnd = round_fast ) : return mpf_cos_sin ( x , prec , rnd , <NUM_LIT:1> ) <EOL> def mpf_sin ( x , prec , rnd = round_fast ) : return mpf_cos_sin ( x , prec , rnd , <NUM_LIT:2> ) <EOL> def mpf_tan ( x , prec , rnd = round_fast ) : return mpf_cos_sin ( x , prec , rnd , <NUM_LIT:3> ) <EOL> def mpf_cos_sin_pi ( x , prec , rnd = round_fast ) : return mpf_cos_sin ( x , prec , rnd , <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> def mpf_cos_pi ( x , prec , rnd = round_fast ) : return mpf_cos_sin ( x , prec , rnd , <NUM_LIT:1> , <NUM_LIT:1> ) <EOL> def mpf_sin_pi ( x , prec , rnd = round_fast ) : return mpf_cos_sin ( x , prec , rnd , <NUM_LIT:2> , <NUM_LIT:1> ) <EOL> def mpf_cosh ( x , prec , rnd = round_fast ) : return mpf_cosh_sinh ( x , prec , rnd ) [ <NUM_LIT:0> ] <EOL> def mpf_sinh ( x , prec , rnd = round_fast ) : return mpf_cosh_sinh ( x , prec , rnd ) [ <NUM_LIT:1> ] <EOL> def mpf_tanh ( x , prec , rnd = round_fast ) : return mpf_cosh_sinh ( x , prec , rnd , tanh = <NUM_LIT:1> ) <EOL> def cos_sin_fixed ( x , prec , pi2 = None ) : <EOL> if pi2 is None : <EOL> pi2 = pi_fixed ( prec - <NUM_LIT:1> ) <EOL> n , t = divmod ( x , pi2 ) <EOL> n = int ( n ) <EOL> c , s = cos_sin_basecase ( t , prec ) <EOL> m = n & <NUM_LIT:3> <EOL> if m == <NUM_LIT:0> : return c , s <EOL> if m == <NUM_LIT:1> : return - s , c <EOL> if m == <NUM_LIT:2> : return - c , - s <EOL> if m == <NUM_LIT:3> : return s , - c <EOL> def exp_fixed ( x , prec , ln2 = None ) : <EOL> if ln2 is None : <EOL> ln2 = ln2_fixed ( prec ) <EOL> n , t = divmod ( x , ln2 ) <EOL> n = int ( n ) <EOL> v = exp_basecase ( t , prec ) <EOL> if n >= <NUM_LIT:0> : <EOL> return v << n <EOL> else : <EOL> return v >> ( - n ) <EOL> if BACKEND == '<STR_LIT>' : <EOL> try : <EOL> import sage . libs . mpmath . ext_libmp as _lbmp <EOL> mpf_sqrt = _lbmp . mpf_sqrt <EOL> mpf_exp = _lbmp . mpf_exp <EOL> mpf_log = _lbmp . mpf_log <EOL> mpf_cos = _lbmp . mpf_cos <EOL> mpf_sin = _lbmp . mpf_sin <EOL> mpf_pow = _lbmp . mpf_pow <EOL> exp_fixed = _lbmp . exp_fixed <EOL> cos_sin_fixed = _lbmp . cos_sin_fixed <EOL> log_int_fixed = _lbmp . log_int_fixed <EOL> except ( ImportError , AttributeError ) : <EOL> print ( "<STR_LIT>" ) </s>
