main
prec = 1024 ' bits precision, about 300 decimal digits
n = 6
DIM x
(n
) AS mpfr_struct
, g
AS mpfr_struct
CALL mpfr_init2
(x
(i
), prec
) c = mpfr_set_si(g, 5, MPFR_RNDN)
CALL make_lanczos
(x
(), g
, error_code
) PRINT "lanczos approximation for the gamma function" CALL print_mpfr_f
(x
(i
), 27, 20)
PRINT "===============================================" PRINT "lanczos approximation for the gamma function" PRINT " n = 14, g = 4.7421875" n = 14
CALL mpfr_init2
(x
(i
), prec
) CALL valmpfr
(g
, "4.7421875", 10) CALL make_lanczos
(x
(), g
, error_code
) CALL print_mpfr_f
(x
(i
), 27, 20)
SUB make_lanczos
(x
() AS mpfr_struct
, g
AS mpfr_struct
, error_code
AS LONG)
DIM a
(n
, n
) AS mpfr_struct
DIM sqr2pi
AS mpfr_struct
, pi2
AS mpfr_struct
, fac
AS mpfr_struct
DIM one
AS mpfr_struct
, tmp
AS mpfr_struct
, tmp2
AS mpfr_struct
CALL mpfr_init2
(a
(i
, j
), prec
) CALL mpfr_init2
(b
(i
), prec
) CALL mpfr_init2
(sqr2pi
, prec
) CALL mpfr_init2
(pi2
, prec
) CALL mpfr_init2
(fac
, prec
) CALL mpfr_init2
(one
, prec
) CALL mpfr_init2
(tmp
, prec
) CALL mpfr_init2
(tmp2
, prec
) c = mpfr_set_si(one, 1, MPFR_RNDN)
c = mpfr_const_pi(pi2, MPFR_RNDN)
c = mpfr_mul_si(pi2, pi2, 2, MPFR_RNDN)
c = mpfr_sqrt(sqr2pi, pi2, MPFR_RNDN)
'setup the matrix, for n=6 the matrix looks like this
' [[ 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 ]]
c = mpfr_set_si(a(i, 0), 1, MPFR_RNDN)
c = mpfr_div_si(a(i, j), one, i + j, MPFR_RNDN)
'setup the constant terms, for n=6 and g=5 they are
'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
c = mpfr_set_si(fac, 1, MPFR_RNDN)
c = mpfr_mul_si(fac, fac, i, MPFR_RNDN)
'b(i)=exp(i+half+g)/(i+half+g)^(i+half)*fac/sqr2pi
c = mpfr_add_d(tmp, g, 0.5, MPFR_RNDN)
c = mpfr_add_si(tmp, tmp, i, MPFR_RNDN)
c = mpfr_exp(b(i), tmp, MPFR_RNDN)
c = mpfr_set_d(tmp2, i + 0.5, MPFR_RNDN)
c = mpfr_pow(tmp, tmp, tmp2, MPFR_RNDN)
c = mpfr_div(b(i), b(i), tmp, MPFR_RNDN)
c = mpfr_mul(b(i), b(i), fac, MPFR_RNDN)
c = mpfr_div(b(i), b(i), sqr2pi, MPFR_RNDN)
'CALL print_mpfr_f(b(i), 27, 20)
'now solve the system of linear equations, for n=6 and g=5 they are
' [[ 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
'and the solution for the above example is
'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
'which is the Lanczos gamma function approximation
'as published in Numerical recipes
CALL gauss_jordan
(a
(), b
(), x
(), n
, error_code
)
SUB gauss_jordan
(a
() AS mpfr_struct
, y
() AS mpfr_struct
, coef
() AS mpfr_struct
, ncol
AS LONG, error_code
AS LONG)
' matrix solution by Gaussian Elimination
DIM b
(ncol
, ncol
) AS mpfr_struct
, w
(ncol
) AS mpfr_struct
' work array, ncol long
DIM hold
AS mpfr_struct
, sm
AS mpfr_struct
, tmp
AS mpfr_struct
DIM t
AS mpfr_struct
, ab
AS mpfr_struct
, big
AS mpfr_struct
error_code = FALSE
n = ncol
CALL mpfr_init2
(b
(i
, j
), prec
) CALL mpfr_init2
(w
(i
), prec
) CALL mpfr_init2
(hold
, prec
) CALL mpfr_init2
(sm
, prec
) CALL mpfr_init2
(ab
, prec
) CALL mpfr_init2
(big
, prec
) CALL mpfr_init2
(tmp
, prec
)
' copy to work arrays
'b(i,j)=a(i,j)
c = mpfr_set(b(i, j), a(i, j), MPFR_RNDN)
'w(i)=y(i)
c = mpfr_set(w(i), y(i), MPFR_RNDN)
'coef(i)=zero
c = mpfr_set_si(coef(i), 0, MPFR_RNDN)
c = mpfr_abs(big, b(i, i), MPFR_RNDN)
l = i
i1 = i + 1
' search for largest element
c = mpfr_abs(ab, b(j, i), MPFR_RNDN)
c = mpfr_set(big, ab, MPFR_RNDN)
l = j
IF mpfr_cmp_si
(big
, 0) = 0 THEN error_code = TRUE
' interchange rows to put largest element on diagonal
c = mpfr_set(hold, b(l, j), MPFR_RNDN)
c = mpfr_set(b(l, j), b(i, j), MPFR_RNDN)
c = mpfr_set(b(i, j), hold, MPFR_RNDN)
c = mpfr_set(hold, w(l), MPFR_RNDN)
c = mpfr_set(w(l), w(i), MPFR_RNDN)
c = mpfr_set(w(i), hold, MPFR_RNDN)
c = mpfr_div(t, b(j, i), b(i, i), MPFR_RNDN)
c = mpfr_mul(tmp, t, b(i, k), MPFR_RNDN)
c = mpfr_sub(b(j, k), b(j, k), tmp, MPFR_RNDN)
c = mpfr_mul(tmp, t, w(i), MPFR_RNDN)
c = mpfr_sub(w(j), w(j), tmp, MPFR_RNDN)
IF mpfr_cmp_si
(b
(n
, n
), 0) = 0 THEN error_code = TRUE
c = mpfr_div(coef(n), w(n), b(n, n), MPFR_RNDN)
i = n - 1
' back substitution
c = mpfr_set_si(sm, 0, MPFR_RNDN)
c = mpfr_mul(tmp, b(i, j), coef(j), MPFR_RNDN)
c = mpfr_add(sm, sm, tmp, MPFR_RNDN)
c = mpfr_sub(tmp, w(i), sm, MPFR_RNDN)
c = mpfr_div(coef(i), tmp, b(i, i), MPFR_RNDN)
i = i - 1