用公式计算反三角函数
Calculating inverse trigonometric functions with formulas
我一直在尝试创建用于计算三角函数的自定义计算器。除了 Chebyshev pylonomials and/or Cordic 算法外，我还使用了泰勒级数，它精确到小数点后几位。
这是我为计算没有任何模块的简单三角函数而创建的：
from __future__ import division
def sqrt(n):
ans = n ** 0.5
return ans
def factorial(n):
k = 1
for i in range(1, n+1):
k = i * k
return k
def sin(d):
pi = 3.14159265359
n = 180 / int(d) # 180 degrees = pi radians
x = pi / n # Converting degrees to radians
ans = x - ( x ** 3 / factorial(3) ) + ( x ** 5 / factorial(5) ) - ( x ** 7 / factorial(7) ) + ( x ** 9 / factorial(9) )
return ans
def cos(d):
pi = 3.14159265359
n = 180 / int(d)
x = pi / n
ans = 1 - ( x ** 2 / factorial(2) ) + ( x ** 4 / factorial(4) ) - ( x ** 6 / factorial(6) ) + ( x ** 8 / factorial(8) )
return ans
def tan(d):
ans = sin(d) / sqrt(1 - sin(d) ** 2)
return ans
不幸的是，我找不到任何可以帮助我解释 Python 的反三角函数公式的资源。我也试过将 sin(x) 计算为 -1 (sin(x) ** -1) 的幂，但没有按预期工作。
在 Python 中执行此操作的最佳解决方案是什么（最好的，我的意思是最简单的，精度与泰勒级数相似）？这可能是幂级数还是我需要使用 cordic 算法？
这个问题的范围很广，但这里有一些简单的想法（和代码！）可以作为计算的起点 arctan。首先，好的老泰勒级数。为简单起见，我们使用固定数量的术语；在实践中，您可能希望根据 x 的大小来决定要动态使用的项数，或者引入某种收敛准则。使用固定数量的项，我们可以使用类似于霍纳方案的方法进行有效评估。
def arctan_taylor(x, terms=9):
"""
Compute arctan for small x via Taylor polynomials.
Uses a fixed number of terms. The default of 9 should give good results for
abs(x) < 0.1. Results will become poorer as abs(x) increases, becoming
unusable as abs(x) approaches 1.0 (the radius of convergence of the
series).
"""
# Uses Horner's method for evaluation.
t = 0.0
for n in range(2*terms-1, 0, -2):
t = 1.0/n - x*x*t
return x * t
上面的代码对于较小的 x（比如绝对值小于 0.1）给出了很好的结果，但是随着 x 变大，精度下降，而对于 abs(x) > 1.0，无论我们投入多少项（或多少额外精度），级数都不会收敛。所以我们需要一种更好的方法来计算更大的 x。一种解决方案是通过恒等式 arctan(x) = 2 * arctan(x / (1 + sqrt(1 + x^2))) 使用参数缩减。这给出了以下代码，它建立在 arctan_taylor 的基础上，为范围广泛的 x 提供合理的结果（但在计算 x*x 时要注意可能的溢出和下溢）。
import math
def arctan_taylor_with_reduction(x, terms=9, threshold=0.1):
"""
Compute arctan via argument reduction and Taylor series.
Applies reduction steps until x is below `threshold`,
then uses Taylor series.
"""
reductions = 0
while abs(x) > threshold:
x = x / (1 + math.sqrt(1 + x*x))
reductions += 1
return arctan_taylor(x, terms=terms) * 2**reductions
或者，给定 tan 的现有实现，您可以使用传统的求根方法简单地找到方程 tan(y) = x 的解 y。由于 arctan 已经自然地位于区间 (-pi/2, pi/2) 内，二分搜索效果很好：
def arctan_from_tan(x, tolerance=1e-15):
"""
Compute arctan as the inverse of tan, via bisection search. This assumes
that you already have a high quality tan function.
"""
low, high = -0.5 * math.pi, 0.5 * math.pi
while high - low > tolerance:
mid = 0.5 * (low + high)
if math.tan(mid) < x:
low = mid
else:
high = mid
return 0.5 * (low + high)
最后，为了好玩，这里有一个类似于 CORDIC 的实现，它确实比 Python 更适合于低级实现。这里的想法是，您一劳永逸地预先计算 1、1/2, 1/4, 等的 arctan 值的 table，然后使用它们来计算一般 arctan值，主要是通过计算真实角度的逐次近似值。值得注意的是，在预计算步骤之后，arctan 计算仅涉及加法、减法和乘以 2 的幂。（当然，这些乘法在 Python，但更接近硬件，这可能会产生很大的不同。）
cordic_table_size = 60
cordic_table = [(2**-i, math.atan(2**-i))
for i in range(cordic_table_size)]
def arctan_cordic(y, x=1.0):
"""
Compute arctan(y/x), assuming x positive, via CORDIC-like method.
"""
r = 0.0
for t, a in cordic_table:
if y < 0:
r, x, y = r - a, x - t*y, y + t*x
else:
r, x, y = r + a, x + t*y, y - t*x
return r
上述每种方法都有其优点和缺点，并且可以通过多种方式改进上述所有代码。我鼓励您尝试和探索。
总而言之，这里是对少量未经过仔细选择的测试值调用上述函数的结果，与标准库 math.atan 函数的输出进行比较：
test_values = [2.314, 0.0123, -0.56, 168.9]
for value in test_values:
print("{:20.15g} {:20.15g} {:20.15g} {:20.15g}".format(
math.atan(value),
arctan_taylor_with_reduction(value),
arctan_from_tan(value),
arctan_cordic(value),
))
我机器上的输出：
1.16288340166519 1.16288340166519 1.16288340166519 1.16288340166519
0.0122993797673 0.0122993797673 0.0122993797673002 0.0122993797672999
-0.510488321916776 -0.510488321916776 -0.510488321916776 -0.510488321916776
1.56487573286064 1.56487573286064 1.56487573286064 1.56487573286064
执行任何反函数的最简单方法是使用二进制搜索。
定义
假设函数
x = g(y)
我们要对其进行逆编码：
y = f(x) = f(g(y))
x = <x0,x1>
y = <y0,y1>
在浮点数上进行 bin 搜索
您可以在访问尾数位的整数数学运算上执行此操作，如下所示：
但是如果您在计算之前不知道结果的指数，那么您也需要使用浮点数进行 bin 搜索。
所以二分查找的思路是把y的尾数从y1一点一点的改成y0，从MSB一点一点改成LSB。然后调用直接函数 g(y)，如果结果与 x 相交，则恢复最后的位更改。
在使用浮点数的情况下，您可以使用变量来保存目标尾数位的近似值，而不是整数位访问。这将消除未知指数问题。所以在开始时将 y = y0 和实际位设置为 MSB 值，因此 b=(y1-y0)/2。每次迭代后将其减半并进行与尾数位相同的迭代次数 n... 这样您就可以在 (y1-y0)/2^n 精度内获得 n 次迭代的结果。
如果您的反函数不是单调的，请将其分解为单调区间并将每个区间作为单独的二进制搜索处理。
函数increasing/decreasing只是判断穿越条件方向（使用<或>）。
C++ acos 示例
所以 y = acos(x) 是在 x = <-1,+1> , y = <0,M_PI> 上定义的，因此递减：
double f64_acos(double x)
{
const int n=52; // mantisa bits
double y,y0,b;
int i;
// handle domain error
if (x<-1.0) return 0;
if (x>+1.0) return 0;
// x = <-1,+1> , y = <0,M_PI> , decreasing
for (y= 0.0,b=0.5*M_PI,i=0;i<n;i++,b*=0.5) // y is min, b is half of max and halving each iteration
{
y0=y; // remember original y
y+=b; // try set "bit"
if (cos(y)<x) y=y0; // if result cross x return to original y decreasing is < and increasing is >
}
return y;
}
我是这样测试的：
double x0,x1,y;
for (x0=0.0;x0<M_PI;x0+=M_PI*0.01) // cycle all angle range <0,M_PI>
{
y=cos(x0); // direct function (from math.h)
x1=f64_acos(y); // my inverse function
if (fabs(x1-x0)>1e-9) // check result and output to log if error
Form1->mm_log->Lines->Add(AnsiString().sprintf("acos(%8.3lf) = %8.3lf != %8.3lf",y,x0,x1));
}
未发现任何差异...因此实施工作正常。在 52 位尾数上进行粗二分搜索通常比多项式逼近慢...另一方面，实现非常简单...
[注释]
如果你不想处理单调区间，你可以试试
在处理测角函数时，您需要处理奇点以避免 NaN 或被零除等......
如果您对此处感兴趣，请查看更多 bin 搜索示例（主要针对整数）
我一直在尝试创建用于计算三角函数的自定义计算器。除了 Chebyshev pylonomials and/or Cordic 算法外，我还使用了泰勒级数，它精确到小数点后几位。
这是我为计算没有任何模块的简单三角函数而创建的：
from __future__ import division
def sqrt(n):
ans = n ** 0.5
return ans
def factorial(n):
k = 1
for i in range(1, n+1):
k = i * k
return k
def sin(d):
pi = 3.14159265359
n = 180 / int(d) # 180 degrees = pi radians
x = pi / n # Converting degrees to radians
ans = x - ( x ** 3 / factorial(3) ) + ( x ** 5 / factorial(5) ) - ( x ** 7 / factorial(7) ) + ( x ** 9 / factorial(9) )
return ans
def cos(d):
pi = 3.14159265359
n = 180 / int(d)
x = pi / n
ans = 1 - ( x ** 2 / factorial(2) ) + ( x ** 4 / factorial(4) ) - ( x ** 6 / factorial(6) ) + ( x ** 8 / factorial(8) )
return ans
def tan(d):
ans = sin(d) / sqrt(1 - sin(d) ** 2)
return ans
不幸的是，我找不到任何可以帮助我解释 Python 的反三角函数公式的资源。我也试过将 sin(x) 计算为 -1 (sin(x) ** -1) 的幂，但没有按预期工作。
在 Python 中执行此操作的最佳解决方案是什么（最好的，我的意思是最简单的，精度与泰勒级数相似）？这可能是幂级数还是我需要使用 cordic 算法？
这个问题的范围很广，但这里有一些简单的想法（和代码！）可以作为计算的起点 arctan。首先，好的老泰勒级数。为简单起见，我们使用固定数量的术语；在实践中，您可能希望根据 x 的大小来决定要动态使用的项数，或者引入某种收敛准则。使用固定数量的项，我们可以使用类似于霍纳方案的方法进行有效评估。
def arctan_taylor(x, terms=9):
"""
Compute arctan for small x via Taylor polynomials.
Uses a fixed number of terms. The default of 9 should give good results for
abs(x) < 0.1. Results will become poorer as abs(x) increases, becoming
unusable as abs(x) approaches 1.0 (the radius of convergence of the
series).
"""
# Uses Horner's method for evaluation.
t = 0.0
for n in range(2*terms-1, 0, -2):
t = 1.0/n - x*x*t
return x * t
上面的代码对于较小的 x（比如绝对值小于 0.1）给出了很好的结果，但是随着 x 变大，精度下降，而对于 abs(x) > 1.0，无论我们投入多少项（或多少额外精度），级数都不会收敛。所以我们需要一种更好的方法来计算更大的 x。一种解决方案是通过恒等式 arctan(x) = 2 * arctan(x / (1 + sqrt(1 + x^2))) 使用参数缩减。这给出了以下代码，它建立在 arctan_taylor 的基础上，为范围广泛的 x 提供合理的结果（但在计算 x*x 时要注意可能的溢出和下溢）。
import math
def arctan_taylor_with_reduction(x, terms=9, threshold=0.1):
"""
Compute arctan via argument reduction and Taylor series.
Applies reduction steps until x is below `threshold`,
then uses Taylor series.
"""
reductions = 0
while abs(x) > threshold:
x = x / (1 + math.sqrt(1 + x*x))
reductions += 1
return arctan_taylor(x, terms=terms) * 2**reductions
或者，给定 tan 的现有实现，您可以使用传统的求根方法简单地找到方程 tan(y) = x 的解 y。由于 arctan 已经自然地位于区间 (-pi/2, pi/2) 内，二分搜索效果很好：
def arctan_from_tan(x, tolerance=1e-15):
"""
Compute arctan as the inverse of tan, via bisection search. This assumes
that you already have a high quality tan function.
"""
low, high = -0.5 * math.pi, 0.5 * math.pi
while high - low > tolerance:
mid = 0.5 * (low + high)
if math.tan(mid) < x:
low = mid
else:
high = mid
return 0.5 * (low + high)
最后，为了好玩，这里有一个类似于 CORDIC 的实现，它确实比 Python 更适合于低级实现。这里的想法是，您一劳永逸地预先计算 1、1/2, 1/4, 等的 arctan 值的 table，然后使用它们来计算一般 arctan值，主要是通过计算真实角度的逐次近似值。值得注意的是，在预计算步骤之后，arctan 计算仅涉及加法、减法和乘以 2 的幂。（当然，这些乘法在 Python，但更接近硬件，这可能会产生很大的不同。）
cordic_table_size = 60
cordic_table = [(2**-i, math.atan(2**-i))
for i in range(cordic_table_size)]
def arctan_cordic(y, x=1.0):
"""
Compute arctan(y/x), assuming x positive, via CORDIC-like method.
"""
r = 0.0
for t, a in cordic_table:
if y < 0:
r, x, y = r - a, x - t*y, y + t*x
else:
r, x, y = r + a, x + t*y, y - t*x
return r
上述每种方法都有其优点和缺点，并且可以通过多种方式改进上述所有代码。我鼓励您尝试和探索。
总而言之，这里是对少量未经过仔细选择的测试值调用上述函数的结果，与标准库 math.atan 函数的输出进行比较：
test_values = [2.314, 0.0123, -0.56, 168.9]
for value in test_values:
print("{:20.15g} {:20.15g} {:20.15g} {:20.15g}".format(
math.atan(value),
arctan_taylor_with_reduction(value),
arctan_from_tan(value),
arctan_cordic(value),
))
我机器上的输出：
1.16288340166519 1.16288340166519 1.16288340166519 1.16288340166519
0.0122993797673 0.0122993797673 0.0122993797673002 0.0122993797672999
-0.510488321916776 -0.510488321916776 -0.510488321916776 -0.510488321916776
1.56487573286064 1.56487573286064 1.56487573286064 1.56487573286064
执行任何反函数的最简单方法是使用二进制搜索。
定义
假设函数
x = g(y)我们要对其进行逆编码：
y = f(x) = f(g(y)) x = <x0,x1> y = <y0,y1>在浮点数上进行 bin 搜索
您可以在访问尾数位的整数数学运算上执行此操作，如下所示：
但是如果您在计算之前不知道结果的指数，那么您也需要使用浮点数进行 bin 搜索。
所以二分查找的思路是把
y的尾数从y1一点一点的改成y0，从MSB一点一点改成LSB。然后调用直接函数g(y)，如果结果与x相交，则恢复最后的位更改。在使用浮点数的情况下，您可以使用变量来保存目标尾数位的近似值，而不是整数位访问。这将消除未知指数问题。所以在开始时将
y = y0和实际位设置为 MSB 值，因此b=(y1-y0)/2。每次迭代后将其减半并进行与尾数位相同的迭代次数n... 这样您就可以在(y1-y0)/2^n精度内获得n次迭代的结果。如果您的反函数不是单调的，请将其分解为单调区间并将每个区间作为单独的二进制搜索处理。
函数increasing/decreasing只是判断穿越条件方向（使用
<或>）。
C++ acos 示例
所以 y = acos(x) 是在 x = <-1,+1> , y = <0,M_PI> 上定义的，因此递减：
double f64_acos(double x)
{
const int n=52; // mantisa bits
double y,y0,b;
int i;
// handle domain error
if (x<-1.0) return 0;
if (x>+1.0) return 0;
// x = <-1,+1> , y = <0,M_PI> , decreasing
for (y= 0.0,b=0.5*M_PI,i=0;i<n;i++,b*=0.5) // y is min, b is half of max and halving each iteration
{
y0=y; // remember original y
y+=b; // try set "bit"
if (cos(y)<x) y=y0; // if result cross x return to original y decreasing is < and increasing is >
}
return y;
}
我是这样测试的：
double x0,x1,y;
for (x0=0.0;x0<M_PI;x0+=M_PI*0.01) // cycle all angle range <0,M_PI>
{
y=cos(x0); // direct function (from math.h)
x1=f64_acos(y); // my inverse function
if (fabs(x1-x0)>1e-9) // check result and output to log if error
Form1->mm_log->Lines->Add(AnsiString().sprintf("acos(%8.3lf) = %8.3lf != %8.3lf",y,x0,x1));
}
未发现任何差异...因此实施工作正常。在 52 位尾数上进行粗二分搜索通常比多项式逼近慢...另一方面，实现非常简单...
[注释]
如果你不想处理单调区间，你可以试试
在处理测角函数时，您需要处理奇点以避免 NaN 或被零除等......
如果您对此处感兴趣，请查看更多 bin 搜索示例（主要针对整数）