Python

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 的回答

的帮助下更改代码以获得精确结果

将无限级数计算封装到另一个函数中 asin():

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

使用 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]")

结果：

 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

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

来近似求和项

如果有人能进行精确的微积分，我将非常高兴。

Answer 2

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 or fdlibm。

Python - 计算反正弦值时不准确

Python - Inaccuracy while computing arcsin

math

trigonometry

python-3.x

结果

问题