/ jepler_udecimal / utrig.py
utrig.py
1 #!/usr/bin/env python3 2 # -*- utf-8 -*- 3 # SPDX-FileCopyrightText: 2020 Jeff Epler <https://unpythonic.net> 4 # 5 # SPDX-License-Identifier: BSD-2-Clause 6 # 7 # Adapted from https://git.yzena.com/gavin/bc/src/branch/master/gen/lib.bc 8 # which states: SPDX-License-Identifier: BSD-2-Clause 9 # 10 # Copyright (c) 2018-2020 Gavin D. Howard and contributors. 11 # 12 # The algorithms use range reductions and taylor polynomaials 13 14 # pylint: disable=invalid-name,protected-access 15 16 """ 17 Trig functions using jepler_udecimal 18 19 Importing this module adds the relevant methods to the `Decimal` object. 20 21 Generally speaking, these routines increase the precision by some amount, 22 perform argument range reduction followed by evaluation of a taylor polynomial, 23 then reduce the precision of the result to equal the origial context's 24 precision. 25 26 There is no guarantee that the results are correctly rounded in all cases, 27 however, in all but the rarest cases the digits except the last one can be 28 trusted. 29 30 Here are some examples of using utrig: 31 32 >>> import jepler_udecimal.utrig 33 >>> from jepler_udecimal import Decimal 34 >>> Decimal('.7').atan() 35 Decimal('0.6107259643892086165437588765') 36 >>> Decimal('.1').acos() 37 Decimal('1.470628905633336822885798512') 38 >>> Decimal('-.1').asin() 39 Decimal('-0.1001674211615597963455231795') 40 >>> Decimal('.4').tan() 41 Decimal('0.4227932187381617619816354272') 42 >>> Decimal('.5').cos() 43 Decimal('0.8775825618903727161162815826') 44 >>> Decimal('.6').sin() 45 Decimal('0.5646424733950353572009454457') 46 >>> Decimal('1').asin() 47 Decimal('1.570796326794896619231321692') 48 >>> Decimal('-1').acos() 49 Decimal('3.141592653589793238462643383') 50 51 52 """ 53 54 from . import Decimal, localcontext, getcontext, InvalidOperation 55 56 __all__ = ["acos", "asin", "atan", "cos", "sin", "tan"] 57 58 _point2 = Decimal(".2") 59 60 61 def atan(x, context=None): 62 """Compute the arctangent of the specified value, in radians""" 63 if not isinstance(x, Decimal): 64 x = Decimal(x) 65 66 ans = x._check_nans(context=context) 67 if ans: 68 return ans 69 70 with localcontext(context) as ctx: 71 scale = ctx.prec 72 73 n = 1 74 if x < 0: 75 n = -1 76 x = -x 77 78 # Hard code values for inputs +-1 and +-.2 79 if scale < 65: 80 if x == 1: 81 return ( 82 Decimal( 83 ".7853981633974483096156608458198757210492923498437764552437361480" 84 ) 85 / n 86 ) 87 if x == _point2: 88 return ( 89 Decimal( 90 ".1973955598498807583700497651947902934475851037878521015176889402" 91 ) 92 / n 93 ) 94 95 if x > _point2: 96 ctx.prec += 5 97 a = atan(_point2) 98 else: 99 a = 0 100 101 ctx.prec = scale + 3 102 103 # This very efficient range reduction reduces 1e300 to under .2 in 104 # just 6 iterations! 105 m = 0 106 while x > _point2: 107 m += 1 108 x = (x - _point2) / (1 + _point2 * x) 109 110 r = u = x 111 f = -x * x 112 t = Decimal(1) 113 i = 3 114 115 while t and t.logb() > -ctx.prec: 116 u *= f 117 t = u / i 118 i += 2 119 r += t 120 121 r += m * a 122 return r / n 123 124 125 def sin(x, context=None): 126 """Compute the sine of the specified value, in radians""" 127 if not isinstance(x, Decimal): 128 x = Decimal(x) 129 130 ans = x._check_nans(context=context) 131 if ans: 132 return ans 133 134 with localcontext(context) as ctx: 135 if x < 0: 136 return -sin(-x) 137 138 scale = ctx.prec 139 140 ctx.prec = int(1.1 * scale + 2) 141 a = atan(1) 142 q = (x // a + 2) // 4 143 x -= 4 * q * a 144 if q % 2: 145 x = -x 146 ctx.prec = scale + 2 147 r = a = x 148 q = -x * x 149 i = 3 150 while a and a.logb() > -ctx.prec: 151 a *= q / (i * (i - 1)) 152 r += a 153 i += 2 154 155 return r / 1 156 157 158 def cos(x, context=None): 159 """Compute the cosine of the specified value, in radians""" 160 if not isinstance(x, Decimal): 161 x = Decimal(x) 162 163 ans = x._check_nans(context=context) 164 if ans: 165 return ans 166 167 with localcontext(context) as ctx: 168 scale = ctx.prec 169 ctx.prec = int(scale * 1.2) 170 r = sin(2 * atan(1) + x) 171 return r / 1 172 173 174 def tan(x, context=None): 175 """Compute the tangent of the specified value, in radians""" 176 if not isinstance(x, Decimal): 177 x = Decimal(x) 178 179 ans = x._check_nans(context=context) 180 if ans: 181 return ans 182 183 with localcontext(context) as ctx: 184 ctx.prec += 2 185 s = sin(x) 186 r = s / (1 - s * s).sqrt() 187 return r / 1 188 189 190 def asin(x, context=None): 191 """Compute the arcsine of the specified value, in radians""" 192 if not isinstance(x, Decimal): 193 x = Decimal(x) 194 195 ans = x._check_nans(context=context) 196 if ans: 197 return ans 198 199 context = context or getcontext() 200 201 if x.compare_total_mag(Decimal(1)) > 0: 202 return context._raise_error(InvalidOperation, "asin(x), |x| > 1") 203 204 with localcontext(context) as ctx: 205 ctx.prec += 2 206 if x == 1: 207 r = atan(Decimal(1)) * 2 # pi * 1/2 radians 208 elif x == -1: 209 r = atan(Decimal(1)) * -2 # pi * -1/2 radians 210 else: 211 r = atan(x / (1 - x * x).sqrt()) 212 return r / 1 213 214 215 def acos(x, context=None): 216 """Compute the arccosine of the specified value, in radians""" 217 if not isinstance(x, Decimal): 218 x = Decimal(x) 219 220 ans = x._check_nans(context=context) 221 if ans: 222 return ans 223 224 context = context or getcontext() 225 226 if x.compare_total_mag(Decimal(1)) > 0: 227 return context._raise_error(InvalidOperation, "acos(x), |x| > 1") 228 229 with localcontext(context) as ctx: 230 ctx.prec += 2 231 if x == 1: 232 r = Decimal(0) # 0 radians 233 elif x == -1: 234 r = atan(Decimal(1)) * 4 # pi radians 235 else: 236 r = atan((1 - x * x).sqrt() / x) 237 if r < 0: 238 r += 4 * atan(1) 239 return r / 1 240 241 242 for name in __all__: 243 setattr(Decimal, name, globals()[name])