Python - 计算arcsin时不准确

Question

我正在尝试在 Python 中实现 arcsin，而不使用任何外部库。

这是我的代码：

from time import process_time as pt

class TrigoCalc(metaclass=__readonly):
    # This class evaluates various Trigonometric functions 
    # including Inverse Trigonometric functions
    def __setattr__(self, name, value):
        raise Exception("Value can't be changed")
    

    @staticmethod
    def asin(x):
        '''Implementation from Taylor series
        asin(x) => summation[(2k)! * x^(2k + 1) / (2^(2k) * (k!)^2 * (2k + 1))]
                  k = [0, inf)
        x should be real
        '''
        # a0 = 1                                                                           
        # a1 = 1/(2*3)                                                                     
        # a2 = 1/2 * 3/(4*5) 
        # a3 = 1/2 * 3/4 * 5/(6*7)
        # a4 = 1/2 * 3/4 * 5/6 * 7/(8*9)
        # a5 = 1/2 * 3/4 * 5/6 * 7/8 * 9/(10*11)
        # a6 = 1/2 * 3/4 * 5/6 * 7/8 * 9/10 * 11/(12*13)
        # a7 = 1/2 * 3/4 * 5/6 * 7/8 * 9/10 * 11/12 * 13/(14*15)
        # a8 = 1/2 * 3/4 * 5/6 * 7/8 * 9/10 * 11/12 * 13/14 * 15/(16*17)
        # a9 = 1/2 * 3/4 * 5/6 * 7/8 * 9/10 * 11/12 * 13/14 * 15/16 * 17/(18*19)
        # a10 = 1/2 * 3/4 * 5/6 * 7/8 * 9/10 * 11/12 * 13/14 * 15/16 * 17/18 * 19/(20*21)
        
        # taking 10 coefficients for arriving at a common sequence
        
        # N = n, D = n + 1; (N/D) --> Multiplication, number of times the coefficient number, n >= 1
        
        start_time = pt()

        coeff_list = []
        NUM_ITER = 10000
        for k in range(NUM_ITER):
            if k == 0:
                coeff_list.append(1)
            else:
                N = 1
                D = N + 1
                C = N/D
                if k >= 2:
                    for i in range(k-1):
                        N += 2; D += 2
                        C = C * N/D
                coeff_list.append(C)
        
        __sum = 0
        for k in range(NUM_ITER):
            n = coeff_list[k] * math_utils.power(x, 2*k + 1) / (2*k + 1)
            __sum += n
    
        # Radian conversion to degrees
        __sum = __sum/TrigoCalc.pi * 180
        
        end_time = pt()

        print(f'Execution time: {end_time - start_time} seconds')

        return __sum

结果

当NUM_ITER 为60（无限级数迭代60 次）时，x = 1 极点处的计算存在明显的不准确性，而x = 1/2 给出了 14 点的精度。

In [2]: TrigoCalc.asin(0.5)
Execution time: 0.0 seconds
Out[2]: 30.000000000000007

In [3]: TrigoCalc.asin(1)
Execution time: 0.0 seconds
Out[3]: 85.823908877692

两次运行的执行时间都不明显。

当NUM_ITER 为10000 时，在x = 1 极点处，结果比上次运行更准确，但在x = 1/2 处，精度完全相同。

In [4]: TrigoCalc.asin(0.5)
Execution time: 19.109375 seconds
Out[4]: 30.000000000000007

In [5]: TrigoCalc.asin(1)
Execution time: 19.109375 seconds
Out[5]: 89.67674183336727 

对于这种类型的计算，这 2 次运行的执行时间非常长。

问题

如何平衡代码，使其在较小的NUM_ITER 极点x = 1 处至少提供 1 点精度？

请随时对代码提出建议或更新。

Python 版本： 3.7.7

编辑：在@Joni 的回答的帮助下更改代码以获得精确的结果

1. 将无限级数计算封装到asin()内的另一个函数中：

    def asin(x):     
        def __arcsin_calc(x):                                        
            # ....                                                 
            # Computations                                               
            # ....                                                             
            # Removing the radian to degree conversion from this function             
            return __sum            
   

2. 使用asin() 中的新函数将限制添加到x 以避免收敛缓慢：

    if -1.0 <= x < -0.5:
        return -(TrigoCalc.pi/2 - __arcsin_calc(math_utils.power((1 - x*x), 0.5))) / TrigoCalc.pi * 180    # Radian to Degree conversion
   
    elif -0.5 <= x <= 0.5:
        return __arcsin_calc(x)/TrigoCalc.pi * 180
   
    elif 0.5 < x <= 1.0:          
        return (TrigoCalc.pi/2 - __arcsin_calc(math_utils.power((1 - x*x), 0.5))) / TrigoCalc.pi * 180
   
    else:      
        raise ValueError("x should be in range of [-1, 1]")
   

3. 结果：

    In [2]: TrigoCalc.asin(0.99)
    Execution time: 0.0 seconds
    Out[2]: 81.89022502527023
   
    In [3]: math.asin(0.99)/TrigoCalc.pi*180
    Out[3]: 81.89038554400582
   
    In [4]: TrigoCalc.asin(1)
    Execution time: 0.0 seconds
    Out[4]: 90.0
   
    In [5]: math.asin(1)/TrigoCalc.pi*180
    Out[5]: 90.0
   

Answer 1

arcsin(1) 的问题是arcsin(x) 在 x=1 处垂直（导数无限增长）。像泰勒级数这样的多项式逼近跟不上它。你会得到一个非常缓慢的收敛，并且需要大量的项来获得一个不错的近似值。您需要改变处理问题的方式。

例如，对于小 x，y = sin(pi/2 - x) 近似为 1 - x^2/2，您可以从中推导出近似值 asin(y) = pi/2 - sqrt(2 - 2*y)。这种近似值适用于非常接近 1 的值 - 您可以直接使用它。

如果你再努力一点，你就能证明准确身份

asin(x) = pi/2 - 2*asin( sqrt( (1-x)/2 ) )

使用此恒等式，您可以使用现有的 asin 函数计算接近 1 的 x asin(x)，该函数适用于接近 0 的 x。

例如：要计算 asin(0.99)，您需要计算：

asin(0.99) = pi/2 - 2*asin( sqrt( (1-.99)/2 ) )
           = pi/2 - 2*asin( sqrt(.005) )
           = pi/2 - 2*asin(0.07071067811865475)

...然后您将使用现有算法来获得asin(0.07071067811865475) 的高质量近似值。

这是生产质量数学库实现中使用的技术 - 参见例如 OpenLibm 或 fdlibm。

Answer 2

一个非常基本的近似值将给出sum from 0 to N 近似于arcsin 1e(-N)（以弧度表示）。 在这里，您以度数给出结果，因为在度数和弧度之间大致存在1e2 的比率，您需要将NUM_ITER = 1e(N+2) 设置为近似于1e(-N) 的arcsin。

因此，对于您的具体问题，您需要使用 N = 1（大约 1 分）进行测试，因此是 NUM_ITER = 1e(1+2) = 1,000。这一点都不准确，但可以让您了解您正在寻找的价值。

然后，如果您想查找确切的值，我看不到每次都可以使用的精确数学方法（无论 x.point 精度如何）。 但是，如果这是您算法的目标，您可以使用二分法算法来查找 NUM_ITER。第一个近似值将减少您的计算时间。

精确的近似值来自比率或x^O(n) 和4^O(n)，4^O(n) 更大。我们可以用O(1/10^n)来近似求和。

如果有人可以进行精确的微积分，我会非常高兴地看到它。