実際、バイナリ分割は、分子と分母のサイズに対して項の数をバランスさせるために反復引数削減と組み合わせると、非常にうまく機能します (これはビットバースト アルゴリズムとして知られています)。
以下は、式 atan(t) = atan(p/2^q) + atan((t*2^qp) / (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 <= x < 1/2 を想定しており、項の数の決定は最適または正しくない可能性があります (これらの問題の修正は読者の課題として残されています)。