''=============================================================================
#Include Once "gmp.bi"
#Include Once "mpfr.bi"
'多重定義の演算子を含める前に mpfr_digits を設定する必要があります
Dim Shared As Long mpfr_digits_precision = 50 '<- mpfr digits-precision を設定
Dim Shared As mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- まるめモードを設定
' まるめモードの選択肢は
' MPFR_RNDN 最近接偶数へのまるめ
' MPFR_RNDZ ゼロ方向へのまるめ
' MPFR_RNDU 正の無限大方向へのまるめ
' MPFR_RNDD 負の無限大方向へのまるめ
' MPFR_RNDA 小数点以下切り上げ
' MPFR_RNDF 忠実な丸め（まだ実装されていない）
' 誤差を含まない正確な計算結果を隣り合うどちらかの浮動小数点数に丸めたもの)
' MPFR_RNDNA 最も近い値への丸め（0から遠いほうへ）(mpfr_round)
' 2つの値の中間の値の場合は、正の値ならより大きいほう、負の値ならより小さいほうの値を採用。
'#include once "C:\dev\FreeBASIC-1.04.0-win64\examples\mpfr-fb2.bas" 'mpfr overloaded operators
#Include Once "mpfr-fb.bi" 'mpfr 多重定義演算子
''=============================================================================
Declare Sub gauss_jordan ( a() As mpfr, y() As mpfr, coef() As mpfr, ByVal ncol As Long, ByRef error_code As Long)
Declare Sub make_lanczos(x() As mpfr, ByVal g As String, ByRef ercode As Long)
Declare Sub gamma_lanczos( ByVal n As Integer = 6, ByRef g As String = "5", ByVal digits_prec As Long=20)
Dim Shared As mpfr e, pi, zero, one, half
Print "これは ニューメリカル・レシピ に掲載されている"
Print "ランチョス ガンンマ関数近似のサンプルです"
Print
'ガンマ関数（ランチョス近似）
gamma_lanczos(6, "5", 20)
'print "=================================="
'gamma_lanczos(14, "4.7421875",20)
gamma_lanczos(13, "6.024680040776729583740234375", 20)
Print "何かキーを押すと終了 ";
Sleep
End
'===============================================================================
Sub gamma_lanczos( ByVal n As Integer = 6, ByRef g As String = "5", ByVal digits_prec As Long=20)
Dim As Integer i, plus_prec=4*n+digits_prec
Dim As Long error_code
Dim As mpfr x(n)
Dim As String frmt, str_out
zero=0
one=1
half=".5"
e=Exp(one)
pi=pi_const()
make_lanczos(x(), g, error_code)
If Not error_code Then
'open "lanczos_gamma_"+str(n)+".txt" for output as 1
For i=0 To n
frmt="%"+Str(digits_prec+5)+"."+Str(digits_prec)+"g"
str_out=x(i).toString(frmt)
frmt=Str(i)
If i<10 Then frmt=" "+frmt
str_out="C("+frmt+") = "+str_out
Print str_out
'Print #1, str_out
Next
'close 1
Print
End If
End Sub
Sub make_lanczos(x() As mpfr, ByVal g_ As String, ByRef error_code As Long)
Dim As Long i, j
Dim As Integer sp
Dim As Integer n = UBound(x)
Dim a(n, n) As mpfr
Dim b(n) As mpfr
Dim As mpfr g, sqr2pi, pi2, fac
pi2=pi*2
sqr2pi=Sqr( pi2)
g=g_
'行列を設定します。n = 6 の場合、行列は次のようになります
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1. 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=0 To n
a(i,0)=1
For j=1 To n
a(i,j)=one/(i+j)
Next
Next
'定数項を設定します。n = 6、g = 5 の場合、
'b( 0) = 41.62443691643906820752
'b( 1) = 16.01231640525168135809
'b( 2) = 9.36473553710404851170
'b( 3) = 6.57049625931606149987
'b( 4) = 5.09521672929646781312
'b( 5) = 4.20380175313001102196
'b( 6) = 3.61487381446333490266
fac=1
For i=0 To n
If i>1 Then
fac=fac*i
End If
b(i)=Exp(i+half+g)/(i+half+g)^(i+half)*fac/sqr2pi
'print b(i).toString("%25.20f") '書式指定 <全体の幅 25>.<小数点以下の幅 20> 実数を出力
Next
'ここで線形方程式の系を解きます。n = 6、g = 5 の場合
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ] = 41.62443691643906820752
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ] = 16.01231640525168135809
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ] = 9.36473553710404851170
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ] = 6.57049625931606149987
' [ 1. 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ] = 5.09521672929646781312
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ] = 4.20380175313001102196
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]] = 3.61487381446333490266
'上記の例の解は
'C(0) = 1.00000000019001482399
'C(1) = 76.18009172947146348309
'C(2) = -86.50532032941676765250
'C(3) = 24.01409824083091049018
'C(4) = -1.23173957245015538752
'C(5) = 0.00120865097386617851
'C(6) = -0.00000539523938495313
'これは ニューメリカル・レシピ に掲載されている
'ランチョス ガンンマ関数近似です
gauss_jordan(a(), b(), x(), n, error_code)
End Sub
Sub gauss_jordan ( a() As mpfr, y() As mpfr, coef() As mpfr, ByVal ncol As Long, ByRef error_code As Long)
' ガウスの消去法（Gaussian elimination）による行列解
Dim As mpfr b(ncol, ncol), w(ncol) ' work array, ncol long
Dim As Long i,j,i1,k,l,n
Dim As mpfr hold, sm, t, ab, big
' Const TRUE = -1
' Const FALSE = Not TRUE
error_code=FALSE
n=ncol
For i=0 To n
' copy to work arrays
For j=0 To n
b(i,j)=a(i,j)
Next j
w(i)=y(i)
coef(i)=zero
Next
For i=0 To n-1
big=Abs(b(i,i))
l=i
i1=i+1
For j=i1 To n
' search for largest element
ab=Abs(b(j,i))
If ab>big Then
big=ab
l=j
End If
Next
If big=zero Then
error_code= TRUE
Else
If l<>i Then
' 対角線上で最大の要素を入れる行を交換する
For j=0 To n
hold=b(l,j)
b(l,j)=b(i,j)
b(i,j)=hold
Next
hold=w(l)
w(l)=w(i)
w(i)=hold
End If
For j=i1 To n
t=b(j,i)/b(i,i)
For k=i1 To n
b(j,k)=b(j,k)-t*b(i,k)
Next
w(j)=w(j)-t*w(i)
Next ' j-loop
End If ' if big - else
Next ' i-loop
If b(n,n)=zero Then
error_code=TRUE
Else
coef(n)=w(n)/b(n,n)
i=n-1
' back substitution
Do
Sm=zero
For j=i To n
sm=sm+b(i,j)*coef(j)
Next
coef(i)=(w(i)-sm)/b(i,i)
i=i-1
Loop Until i<0
End If
If error_code Then Print "ERROR: Matrix is singular"
End Sub
