|
|
> No.4674[元記事へ]
島村1243さんへのお返事です。
> その関数をBASメインプログラムからどの様に呼び出せば良いのかを知りたいのですが、上記グラフを描いたBASコードをご教示頂けないでしょうか？
specialfunction.zipファイルをダウンロードしていただくとその中にGAMMADLL?.BASというファイルが5本入っており、そのいずれかを実行すればガンマ関数のグラフを
描くことができますが、Linux版、Mac版をご利用の方は実行できません。Windows版十進BASIC用です。(更に詳しくいうとSSE2 対応のアーキテクチャを持つCPUであること)
5本あるのは、使用しているライブラリーの違いで、定義されている関数に若干違いがあります。
もし、実行時にDLLフィルがロードできない、見つからない等のエラーが出る場合は、ZIPファイルを解凍して出来たDLLフォルダごとBASIC.EXEと同じフォルダに入れてください。
若しくは、ASSIGN文の相対パスを対象のDLLファイルへの絶対パスに書き換えてください。
Linux版、Mac版をご利用の方の為に、Libraryフォルダ内のPROBDIST.LIBにガンマ関数の定義がありましたので、
そちらを使用してグラフ描いてみました。
ガンマ関数の定義域は0より大きい時 X>0 と負数では整数値以外 FP(X)<>0 の場合で描画しています。
CALL GINIT(800,800)
SET WINDOW -5,5,-5,5
DRAW GRID(1,1)
SET LINE COLOR 2
FOR X=-5 TO 5 STEP 1/128
IF FP(X)<>0 OR X>0 THEN
LET Y=GAMMA(X)
PLOT LINES:X,Y;
ELSE
PLOT LINES
END IF
NEXT X
PAUSE
CLEAR
SET WINDOW 0,5,-1,5
DRAW GRID(1,1)
SET LINE COLOR 3
FOR X=1/128 TO 5 STEP 1/128
LET Y=GammaLib.LogGamma(X) ! 負数側が計算できない(X>0 のみ)
PLOT LINES:X,Y;
NEXT X
PAUSE
CLEAR
SET WINDOW -5,5,-5,5
DRAW GRID(1,1)
SET LINE COLOR 4
FOR X=-5 TO 5 STEP 1/128
IF FP(X)<>0 OR X>0 THEN
LET Y=1/GAMMA(X) ! 1/GAMMA(0)=1/GAMMA(-1)=1/GAMMA(-2)=ゼロ ???
PLOT LINES:X,Y;
ELSE
PLOT LINES
END IF
NEXT X
PAUSE
CLEAR
SET WINDOW -1,6,-1,20
DRAW GRID(1,1)
FOR N=1 TO 6 STEP 1/4
LET COL=COL+1
SET LINE COLOR COL
FOR K=0 TO N STEP 1/8
PLOT LINES:K,BINOM(N,K);
NEXT K
PLOT LINES
NEXT N
END
MODULE GammaLib
PUBLIC FUNCTION LogGamma
SHARE NUMERIC LOG_2PI
SHARE NUMERIC N,B0,B1,B2,B4,B6,B8,B10,B12,B14,B16
LET LOG_2PI=1.83787706640934548 !/* $\LOG 2\pi$ */
LET N = 8
LET B0 = 1 !/* 以下はBernoulli数 */
LET B1 = (-1.0 / 2.0)
LET B2 = ( 1.0 / 6.0)
LET B4 = (-1.0 / 30.0)
LET B6 = ( 1.0 / 42.0)
LET B8 = (-1.0 / 30.0)
LET B10 =( 5.0 / 66.0)
LET B12 =(-691.0 / 2730.0)
LET B14 =( 7.0 / 6.0)
LET B16 =(-3617.0 / 510.0)
EXTERNAL FUNCTION LogGamma(x) !/* ガンマ関数の対数 */
DECLARE NUMERIC v, w
LET v = 1
DO WHILE x < N
LET v = v * x
LET x = x+1
LOOP
LET w = 1 / (x * x)
LET LogGamma=((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w &
& + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w &
& + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w &
& + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x &
& + 0.5 * LOG_2PI - LOG(v) - x + (x - 0.5) * LOG(x)
END FUNCTION
END MODULE
EXTERNAL FUNCTION Gamma(x) !/* ガンマ関数 */
DECLARE EXTERNAL FUNCTION GammaLib.LogGamma
IF x < 0 THEN
LET Gamma= PI / (SIN(PI * x) * EXP(LogGamma(1 - x)))
ELSE
LET Gamma= EXP(LogGamma(x))
END IF
END FUNCTION
EXTERNAL FUNCTION Beta(x, y) !/* ベータ関数 */
DECLARE EXTERNAL FUNCTION GammaLib.LogGamma
LET Beta= EXP(loggamma(x) + loggamma(y) - loggamma(x + y))
END FUNCTION
EXTERNAL FUNCTION BINOM(N,K) !'二項係数 COMB(N,K)
LET BINOM=GAMMA(N+1)/GAMMA(N-K+1)/GAMMA(K+1)
END FUNCTION
!EXTERNAL FUNCTION BINOM(N,K) !'二項係数 COMB(N,K)
!LET BINOM=1/(N+1)/BETA(N-K+1,K+1)
!END FUNCTION
EXTERNAL SUB GINIT(XSIZE,YSIZE)
SET BITMAP SIZE XSIZE,YSIZE
SET WINDOW 0,XSIZE-1,YSIZE-1,0
SET POINT STYLE 1
SET COLOR MODE "REGULAR"
FOR I=0 TO 7
SET COLOR MIX(I) BITAND(I,2)/2,BITAND(I,4)/4,BITAND(I,1)
NEXT I
CLEAR
END SUB
|
|
