작고 임의의 부동 소수점 값에 대해 arcsin 함수에 대해 가능한 가장 빠른 방법

Question

mpmath의 "mpf"부동 소수점 bignums 형식의 작은 값의 아크 사인 함수를 계산해야합니다. 나는 "작은"값을 부르는

예를 전자입니다/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 

이것은 특별한 경우입니다 : 작은을 위해 내가 다소 작은 숫자와 함께 작동, 그래서 실제로 특별한 경우 (대한 mpmath 뛰는 파이썬에서 그것을 계산하는 방법이 있는지 나 자신을 부탁 해요 = 값).

테일러 시리즈는 작은 인수로 매우 빠르게 수렴되기 때문에 실제로 여기에서 좋은 선택입니다. 그러나 나는 아직도 어떻게 든 계산을 더 가속화 할 필요가있다.

사실 몇 가지 문제가 있습니다.
1) 이진 분할은 인수를 작은 분수로 쓸 수있을 때만 빛나기 때문에 여기서는 효과가 없습니다. 완전한 정밀도 플로트가 여기에 주어집니다.
2) arcsin은 번갈아 사용하지 않는 시리즈이므로 Van Wijngaarden 또는 Sumalt 변환은 효과적이지 않습니다 (단, 번갈아 반복하지 않는 시리즈로 일반화하지는 않는 한). https://en.wikipedia.org/wiki/Van_Wijngaarden_transformation

내가 생각할 수있는 유일한 가속화는 Chebyshev 다항식입니다. Chebyshev 다항식이 arcsin 함수에 적용될 수 있습니까? 어떻게?

Answer 1

실제로 분자와 분모의 크기에 대한 항의 수를 균형을 맞추기 위해 반복 인수 감소와 결합하면 이진 분할이 잘 작동합니다 (비트 버스트 알고리즘이라고 함).

다음은 공식 atan (t) = atan (p/2^q) + atan (t * 2^qp)/(2^q + p *)을 반복적으로 적용하여 mpmath에 대한 이진 분할 구현입니다. 티)).이 공식은 리차드 브렌트 (Richard Brent)에 의해 최근에 제안되었습니다. 사실, mpmath의 atan은 캐시에서 atan (p/2^q)를 찾기 위해 저수준으로이 공식을 한 번 호출합니다. MPFR은 atan을 평가하기 위해 비트 버스트 알고리즘을 사용하지만, 약간 다른 공식을 사용합니다.이 방법은 몇 가지 서로 다른 아크 탄젠트 값을 평가하는 대신 arctangent 미분 방정식을 사용하여 분석 연속을 수행합니다. 이 코드

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)) 

, I 얻을 :

>>> 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 및 용어의 수의 판정을 가정 이러한 문제는 독자에게 연습 문제로 남아 있습니다.

Answer 2

빠른 계산 방법은 미리 계산 된 룩업 테이블을 사용하는 것입니다.

하지만 예를 들어 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 

물론이 차이가 귀하와 관련이 있는지에 따라 달라집니다.

Answer 3

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 

gmp2는 MPLS 및 MPC 다중 정밀도 라이브러리도 지원합니다.

면책 조항 : gmpy2를 유지합니다.