'***************************************************
'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