arcsin 函数对小的任意浮点值的最快可能方法

Question

我需要计算 mpmath 的“mpf”浮点大数形式的小值的反正弦函数。

我所说的“小”值例如是 e/4/(10**7) = 0.000000067957045711476130884...

这是在我的机器上使用 mpmath 的内置 asin 函数进行测试的结果：

import gmpy2
from mpmath import *
from time import time

mp.dps = 10**6

val=e/4/(10**7)
print "ready"

start=time()
temp=asin(val)
print "mpmath asin: "+str(time()-start)+" seconds"

>>> 155.108999968 seconds

这是一个特殊情况：我处理的数字有点小，所以我问自己是否有一种方法可以在 python 中计算它，在这种特殊情况下实际上优于 mpmath（= 对于小值）。

Taylor 级数在这里实际上是一个不错的选择，因为它们对于小参数的收敛速度非常快。但我仍然需要以某种方式进一步加快计算速度。

其实有一些问题：
1）二进制拆分在这里无效，因为只有当您可以将参数写成一小部分时它才会发光。这里给出了一个全精度浮点数。
2) arcsin 是一个非交替序列，因此 Van Wijngaarden 或 sumalt 变换也无效（除非我不知道有一种方法可以将它们推广到非交替序列）。 https://en.wikipedia.org/wiki/Van_Wijngaarden_transformation

我能想到的唯一加速度是切比雪夫多项式。切比雪夫多项式可以应用于 arcsin 函数吗？怎么做？

Answer 1

你可以使用gmpy2中包含的mpfr类型吗？

>>> import gmpy2
>>> gmpy2.get_context().precision = 3100000
>>> val = gmpy2.exp(1)/4/10**7
>>> from time import time
>>> start=time();r=gmpy2.asin(val);print time()-start
3.36188197136

gmpy2除了支持GMP库外，还支持MPFR和MPC多精度库。

免责声明：我维护 gmpy2。

Answer 2

实际上，如果与迭代参数缩减相结合以平衡项的数量与分子和分母的大小（这称为比特突发算法），二进制拆分确实非常有效。

这是一个基于重复应用公式 atan(t) = atan(p/2^q) + atan((t*2^q-p) / (2^q+p*t) 的 mpmath 二进制拆分实现））。 Richard Brent 最近提出了这个公式（实际上 mpmath 的 atan 已经以低精度使用了这个公式的一次调用，以便从缓存中查找 atan(p/2^q)）。如果我没记错的话，MPFR 也使用比特突发算法来评估 atan，但它使用的公式略有不同，这可能更有效（而不是评估几个不同的反正切值，它使用反正切微分方程进行解析延拓）。

from mpmath.libmp import MPZ, bitcount
from mpmath import mp

def bsplit(p, q, a, b):
    if b - a == 1:
        if a == 0:
            P = p
            Q = q
        else:
            P = p * p
            Q = q * 2
        B = MPZ(1 + 2 * a)
        if a % 2 == 1:
            B = -B
        T = P
        return P, Q, B, T
    else:
        m = a + (b - a) // 2
        P1, Q1, B1, T1 = bsplit(p, q, a, m)
        P2, Q2, B2, T2 = bsplit(p, q, m, b)
        T = ((T1 * B2) << Q2) + T2 * B1 * P1
        P = P1 * P2
        B = B1 * B2
        Q = Q1 + Q2
        return P, Q, B, T

def atan_bsplit(p, q, prec):
    """computes atan(p/2^q) as a fixed-point number"""
    if p == 0:
        return MPZ(0)
    # FIXME
    nterms = (-prec / (bitcount(p) - q) - 1) * 0.5
    nterms = int(nterms) + 1
    if nterms < 1:
        return MPZ(0)
    P, Q, B, T = bsplit(p, q, 0, nterms)
    if prec >= Q:
        return (T << (prec - Q)) // B
    else:
        return T // (B << (Q - prec))

def atan_fixed(x, prec):
    t = MPZ(x)
    s = MPZ(0)
    q = 1
    while t:
        q = min(q, prec)
        p = t >> (prec - q)
        if p:
            s += atan_bsplit(p, q, prec)
            u = (t << q) - (p << prec)
            v = (MPZ(1) << (q + prec)) + p * t
            t = (u << prec) // v
        q *= 2
    return s

def atan1(x):
    prec = mp.prec
    man = x.to_fixed(prec)
    return mp.mpf((atan_fixed(man, prec), -prec))

def asin1(x):
    x = mpf(x)
    return atan1(x/sqrt(1-x**2))

使用这段代码，我得到：

>>> from mpmath import *
>>> mp.dps = 1000000
>>> val=e/4/(10**7)
>>> from time import time
>>> start = time(); y1 = asin(x); print time() - start
58.8485069275
>>> start = time(); y2 = asin1(x); print time() - start
8.26498985291
>>> nprint(y2 - y1)
-2.31674e-1000000

警告：atan1 假设 0

Answer 3

一种快速的方法是使用预先计算的查找表。

但是如果你看一下例如asin 的泰勒级数；

def asin(x):
    rv = (x + 1/3.0*x**3 + 7/30.0*x**5 + 64/315.0*x**7 + 4477/22680.0*x**9 + 
         28447/138600.0*x**11 + 23029/102960.0*x**13 + 
         17905882/70945875.0*x**15 + 1158176431/3958416000.0*x**17 + 
         9149187845813/26398676304000.0*x**19)
    return rv

您会看到，对于较小的 x 值，asin(x) ≈ x。

In [19]: asin(1e-7)
Out[19]: 1.0000000000000033e-07

In [20]: asin(1e-9)
Out[20]: 1e-09

In [21]: asin(1e-11)
Out[21]: 1e-11

In [22]: asin(1e-12)
Out[22]: 1e-12

例如对于我们使用的值：

In [23]: asin(0.000000067957045711476130884)
Out[23]: 6.795704571147624e-08

In [24]: asin(0.000000067957045711476130884)/0.000000067957045711476130884
Out[24]: 1.0000000000000016

当然，这取决于这种差异是否与您相关。