183 lines
4.6 KiB
Python
183 lines
4.6 KiB
Python
#!/usr/bin/env python3
|
|
# -*- utf-8 -*-
|
|
# SPDX-FileCopyrightText: 2020 Jeff Epler <https://unpythonic.net>
|
|
#
|
|
# SPDX-License-Identifier: BSD-2-Clause
|
|
#
|
|
# Adapted from https://git.yzena.com/gavin/bc/src/branch/master/gen/lib.bc
|
|
# which states: SPDX-License-Identifier: BSD-2-Clause
|
|
#
|
|
# Copyright (c) 2018-2020 Gavin D. Howard and contributors.
|
|
#
|
|
# The algorithms use range reductions and taylor polynomaials
|
|
|
|
"""
|
|
Trig functions using jepler_udecimal
|
|
|
|
Importing this module adds the relevant methods to the Decimal object.
|
|
|
|
Generally speaking, these routines increase the precision by some amount,
|
|
perform argument range reduction followed by evaluation of a taylor polynomial,
|
|
then reduce the precision of the result to equal the origial context's
|
|
precision.
|
|
|
|
There is no guarantee that the results are correctly rounded in all cases,
|
|
however, in all but the rarest cases the digits except the last one can be
|
|
trusted.
|
|
|
|
Here are some examples of using utrig:
|
|
|
|
>>> import jepler_udecimal.utrig
|
|
>>> from jepler_udecimal import Decimal
|
|
>>> Decimal('.7').atan()
|
|
Decimal('0.6107259643892086165437588765')
|
|
>>> Decimal('.1').acos()
|
|
Decimal('1.47062890563333682288579851219')
|
|
>>> Decimal('-.1').asin()
|
|
Decimal('-0.1001674211615597963455231795')
|
|
>>> Decimal('.4').tan()
|
|
Decimal('0.4227932187381617619816354272')
|
|
>>> Decimal('.5').cos()
|
|
Decimal('0.8775825618903727161162815826')
|
|
>>> Decimal('.6').sin()
|
|
Decimal('0.5646424733950353572009454457')
|
|
|
|
"""
|
|
|
|
from . import Decimal, localcontext
|
|
|
|
__all__ = ['atan', 'sin', 'cos']
|
|
|
|
_point2 = Decimal('.2')
|
|
|
|
def atan(x, context=None):
|
|
"""Compute the arctangent of the specified value, in radians"""
|
|
if not isinstance(x, Decimal):
|
|
x = Decimal(x)
|
|
|
|
with localcontext(context) as ctx:
|
|
scale = ctx.prec
|
|
|
|
n = 1
|
|
if x < 0:
|
|
n = -1
|
|
x = -x
|
|
|
|
# Hard code values for inputs +-1 and +-.2
|
|
if scale < 65:
|
|
if x == 1:
|
|
return Decimal('.7853981633974483096156608458198757210492923498437764552437361480') / n
|
|
if x == _point2:
|
|
return Decimal('.1973955598498807583700497651947902934475851037878521015176889402') / n
|
|
|
|
if x > _point2:
|
|
ctx.prec += 5
|
|
a=atan(_point2)
|
|
else:
|
|
a=0
|
|
|
|
ctx.prec = scale + 3
|
|
|
|
# This very efficient range reduction reduces 1e300 to under .2 in
|
|
# just 6 iterations!
|
|
m = 0
|
|
while x > _point2:
|
|
m += 1
|
|
x = (x-_point2)/(1+_point2*x)
|
|
|
|
r = u = x
|
|
f = -x * x
|
|
t = Decimal(1)
|
|
i = 3
|
|
|
|
while t.logb() > -ctx.prec:
|
|
u *= f
|
|
t = u/i
|
|
i += 2
|
|
r += t
|
|
|
|
r += m*a
|
|
return r/n
|
|
|
|
def sin(x, context=None):
|
|
"""Compute the sine of the specified value, in radians"""
|
|
if not isinstance(x, Decimal):
|
|
x = Decimal(x)
|
|
|
|
with localcontext(context) as ctx:
|
|
if x < 0:
|
|
return -sin(-x)
|
|
|
|
scale = ctx.prec
|
|
|
|
ctx.prec = int(1.1*scale+2)
|
|
a = atan(1)
|
|
q = (x//a + 2)//4
|
|
x -= 4*q*a
|
|
if q % 2:
|
|
x = -x
|
|
ctx.prec = scale+2
|
|
r = a = x
|
|
q = -x*x
|
|
i=3
|
|
lim = Decimal(f"10e-{ctx.prec}")
|
|
while a.logb() > -ctx.prec:
|
|
a *= q/(i*(i-1))
|
|
r += a
|
|
i += 2
|
|
|
|
return r/1
|
|
|
|
def cos(x, context=None):
|
|
"""Compute the cosine of the specified value, in radians"""
|
|
if not isinstance(x, Decimal):
|
|
x = Decimal(x)
|
|
|
|
with localcontext(context) as ctx:
|
|
scale = ctx.prec
|
|
ctx.prec = int(scale*1.2)
|
|
r = sin(2*atan(1)+x)
|
|
return r/1
|
|
|
|
def tan(x, context=None):
|
|
"""Compute the tangent of the specified value, in radians"""
|
|
if not isinstance(x, Decimal):
|
|
x = Decimal(x)
|
|
|
|
with localcontext(context) as ctx:
|
|
scale = ctx.prec
|
|
ctx.prec += 2
|
|
s = sin(x)
|
|
r = s / (1-s*s).sqrt()
|
|
return r/1
|
|
|
|
def asin(x, context=None):
|
|
"""Compute the arcsine of the specified value, in radians"""
|
|
if not isinstance(x, Decimal):
|
|
x = Decimal(x)
|
|
|
|
with localcontext(context) as ctx:
|
|
scale = ctx.prec
|
|
ctx.prec += 2
|
|
r = atan(x / (1 - x*x).sqrt())
|
|
return r/1
|
|
|
|
def acos(x, context=None):
|
|
"""Compute the arccosine of the specified value, in radians"""
|
|
if not isinstance(x, Decimal):
|
|
x = Decimal(x)
|
|
|
|
with localcontext(context) as ctx:
|
|
scale = ctx.prec
|
|
ctx.prec += 2
|
|
r = atan((1-x*x).sqrt()/x)
|
|
if r < 0:
|
|
r += 4*atan(1)
|
|
return r
|
|
|
|
Decimal.tan = tan
|
|
Decimal.sin = sin
|
|
Decimal.cos = cos
|
|
Decimal.atan = atan
|
|
Decimal.asin = asin
|
|
Decimal.acos = acos
|