173 lines
2.9 KiB
Plaintext
173 lines
2.9 KiB
Plaintext
'***************************************************
|
|
'Subject: Computation of Euler's constant 0.5772...
|
|
' with the Brent-McMillan algorithm B1,
|
|
' Math. Comp. 34 (1980), 305-312
|
|
'tested : FreeBasic 1.08.1 with gmp 6.2.0
|
|
'---------------------------------------------------
|
|
#include "gmp.bi"
|
|
|
|
'multi-precision float pointers
|
|
Dim as mpf_ptr a, b
|
|
Dim shared as mpf_ptr k2, u, v
|
|
'unsigned long integers
|
|
Dim as ulong k, n, n2, r, s, t
|
|
'precision parameters
|
|
Dim shared as ulong e10, e2
|
|
Dim shared e as clong
|
|
Dim shared f as double
|
|
Dim as double tim = TIMER
|
|
CLS
|
|
|
|
a = allocate(len(__mpf_struct))
|
|
b = allocate(len(__mpf_struct))
|
|
u = allocate(len(__mpf_struct))
|
|
v = allocate(len(__mpf_struct))
|
|
k2 = allocate(len(__mpf_struct))
|
|
|
|
'log(x/y) with the Taylor series for atanh(x-y/x+y)
|
|
Sub ln (byval s as mpf_ptr, byval x as ulong, byval y as ulong)
|
|
Dim as mpf_ptr d = u, q = v
|
|
Dim k as ulong
|
|
'Möbius transformation
|
|
k = x: x -= y: y += k
|
|
|
|
If x <> 1 Then
|
|
Print "ln: illegal argument x - y <> 1"
|
|
End
|
|
End If
|
|
|
|
's = 1 / (x + y)
|
|
mpf_set_ui (s, y)
|
|
mpf_ui_div (s, 1, s)
|
|
'k2 = s * s
|
|
mpf_mul (k2, s, s)
|
|
mpf_set (d, s)
|
|
|
|
k = 1
|
|
Do
|
|
k += 2
|
|
'd *= k2
|
|
mpf_mul (d, d, k2)
|
|
'q = d / k
|
|
mpf_div_ui (q, d, k)
|
|
's += q
|
|
mpf_add (s, s, q)
|
|
|
|
f = mpf_get_d_2exp (@e, q)
|
|
Loop until abs(e) > e2
|
|
|
|
's *= 2
|
|
mpf_mul_2exp (s, s, 1)
|
|
End Sub
|
|
|
|
'Main
|
|
|
|
'n = 2^i * 3^j * 5^k
|
|
|
|
'log(n) = r * log(16/15) + s * log(25/24) + t * log(81/80)
|
|
|
|
'solve linear system for r, s, t
|
|
' 4 -3 -4| i
|
|
'-1 -1 4| j
|
|
'-1 2 -1| k
|
|
|
|
'examples
|
|
t = 1
|
|
select case t
|
|
case 1
|
|
n = 60
|
|
r = 41
|
|
s = 30
|
|
t = 18
|
|
'100 digits
|
|
case 2
|
|
n = 4800
|
|
r = 85
|
|
s = 62
|
|
t = 37
|
|
'8000 digits, 0.6 s
|
|
case 3
|
|
n = 9375
|
|
r = 91
|
|
s = 68
|
|
t = 40
|
|
'15625 digits, 2.5 s
|
|
case else
|
|
n = 18750
|
|
r = 98
|
|
s = 73
|
|
t = 43
|
|
'31250 digits, 12 s. @2.00GHz
|
|
end select
|
|
|
|
'decimal precision
|
|
e10 = n / .6
|
|
'binary precision
|
|
e2 = (1 + e10) / .30103
|
|
|
|
'initialize mpf's
|
|
mpf_set_default_prec (e2)
|
|
mpf_inits (a, b, u, v, k2, Cptr(mpf_ptr, 0))
|
|
|
|
'Compute log terms
|
|
|
|
ln b, 16, 15
|
|
|
|
'a = r * b
|
|
mpf_mul_ui (a, b, r)
|
|
|
|
ln b, 25, 24
|
|
|
|
'a += s * b
|
|
mpf_mul_ui (u, b, s)
|
|
mpf_add (a, a, u)
|
|
|
|
ln b, 81, 80
|
|
|
|
'a += t * b
|
|
mpf_mul_ui (u, b, t)
|
|
mpf_add (a, a, u)
|
|
|
|
''gmp_printf (!"log(%lu) %.*Ff\n", n, e10, a)
|
|
|
|
'B&M, algorithm B1
|
|
|
|
'a = -a, b = 1
|
|
mpf_neg (a, a)
|
|
mpf_set_ui (b, 1)
|
|
mpf_set (u, a)
|
|
mpf_set (v, b)
|
|
|
|
k = 0
|
|
n2 = n * n
|
|
'k2 = k * k
|
|
mpf_set_ui (k2, 0)
|
|
do
|
|
'k2 += 2k + 1
|
|
mpf_add_ui (k2, k2, (k shl 1) + 1)
|
|
k += 1
|
|
|
|
'b = b * n2 / k2
|
|
mpf_div (b, b, k2)
|
|
mpf_mul_ui (b, b, n2)
|
|
'a = (a * n2 / k + b) / k
|
|
mpf_div_ui (a, a, k)
|
|
mpf_mul_ui (a, a, n2)
|
|
mpf_add (a, a, b)
|
|
mpf_div_ui (a, a, k)
|
|
|
|
'u += a, v += b
|
|
mpf_add (u, u, a)
|
|
mpf_add (v, v, b)
|
|
|
|
f = mpf_get_d_2exp (@e, a)
|
|
Loop until abs(e) > e2
|
|
|
|
mpf_div (u, u, v)
|
|
gmp_printf (!"gamma %.*Ff (maxerr. 1e-%lu)\n", e10, u, e10)
|
|
|
|
gmp_printf (!"k = %lu\n\n", k)
|
|
|
|
gmp_printf (!"time: %.7f s\n", TIMER - tim)
|
|
end
|