|
|
> No.4103[元記事へ]
たろささんへのお返事です。
> 『無限解析のはじまり』高瀬正仁:著 ちくま学芸文庫
> P328
> 3.34オイラーの公式の証明
>
> を眺めています。数式をプログラムにしてみました。
>
> OPTION ARITHMETIC RATIONAL
> LET format$="#." & REPEAT$("#",999)
> LET e=1
> LET f=1
> FOR i=1 TO 450
> LET f=f*i
> LET e=e+1/f ! 自然対数の底(ネピア数)1000桁
> NEXT i
> ! PRINT USING format$:e
>
> FOR r=3 TO 4
> LET n=10^r
> LET x=2
> LET x1=5
> LET m=(((1+x/n)^n)/(E^x))/(((1+x1/n)^n)/(E^x1))
> PRINT USING format$:m
> NEXT r
> END
>
> LET n=10^6 辺りから重いです。LET n=10^7 何日?
>
1000桁モードで実行すれば、r=18まで高速に計算できます。
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL FUNCTION EXP ! 外部関数宣言(必須)
LET x=2
LET x1=5
LET ex=EXP(x)
LET ex1=EXP(x1)
FOR r=3 TO 18
LET n=10^r
LET m=((1+x/n)^n/ex)/((1+x1/n)^n/ex1)
PRINT m
NEXT r
END
!
EXTERNAL FUNCTION EXP(x) ! EXP関数 (精度約1000桁 有理数モード不可)
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE NUMERIC a,lim
FUNCTION ex(t,n) ! e^a テイラー展開
IF t>lim THEN LET ex=t+ex(t*a/n,n+1) ELSE LET ex=0
END FUNCTION
LET a=ABS(x) ! 定義域 a<=2321.0057737… (10進モード a<=227.9559…)
LET lim=EPS(2.71828^INT(a))*10^(-9) ! EPS(e^a)の9桁下まで計算し加算
IF x>=0 THEN LET EXP=ex(1,1) ELSE LET EXP=1/ex(1,1)
END FUNCTION
外部関数EXP(x)は、約1000桁の精度があります。テイラー展開は関数の再帰呼出しを利用し、最後の項から足し上げることによって精度を高くしました。xは負数でも可能です。
PRINT EXP(1) で自然対数の底を小数点以下1007桁まで得られます。
【おまけ】静岡理工科大学の菅沼義昇先生が公開しているガンマ関数の近似値を求めるプログラムをBASICに移植しました。ガンマ関数が階乗の拡張と言われる理由がわかります。
!/****************************/
!/* ガンマ関数の計算 */
!/* coded by Y.Suganuma */
!/****************************/
DECLARE FUNCTION gamma
DECLARE NUMERIC x,y
INPUT PROMPT "data?":x !-69.38<=x<56.905…(10進モード)
PRINT "x ="; x
LET y=gamma(x)
PRINT " Γ(x) ="; y
PRINT "x*Γ(x) ="; x*y
PRINT "Γ(x+1) ="; gamma(x+1)
!/****************************************/
!/* Γ（ｘ）の計算（ガンマ関数，近似式） */
!/* ier : =0 : normal */
!/* =-1 : x=-n (n=0,1,2,･･･) */
!/* return : 結果 */
!/****************************************/
FUNCTION gamma(x)
LOCAL ier,err,g,s,t,v,w,y,k
LET ier=0
IF x>5 THEN
LET v=1/x
LET s=((((((-0.000592166437354*v+0.0000697281375837)*v+0.00078403922172)*v &
& -0.000229472093621)*v-0.00268132716049)*v+0.00347222222222)*v+0.0833333333333)*v+1
LET g=2.506628274631001*EXP(-x)*x^(x-0.5)*s
ELSE
LET err=1E-20
LET w=x
LET t=1
IF x<1.5 THEN
IF x<err THEN
LET k=INT(x)
LET y=k-x
IF ABS(y)<err OR ABS(1-y)<err THEN LET ier=-1
END IF
IF ier=0 THEN
DO WHILE w<1.5
LET t=t/w
LET w=w+1
LOOP
END IF
ELSE
IF w>2.5 THEN
DO WHILE w>2.5
LET w=w-1
LET t=t*w
LOOP
END IF
END IF
LET w=w-2
LET g=(((((((0.0021385778*w-0.0034961289)*w+0.0122995771)*w-0.00012513767)*w &
& +0.0740648982)*w+0.0815652323)*w+0.411849671)*w+0.422784604)*w+0.999999926
LET g=g*t
END IF
IF ier=-1 THEN PRINT "[ ier=";ier;"]"
LET gamma=g
END FUNCTION
END
|
|