lambdaway :: euler

4 | list | login |load

| |

_h1 [[euler]] | euler's constant 
_h3 C  =  0.57721566490153286...

_ul [[https://rosettacode.org/wiki/Euler%27s_constant_0.5772...|https://rosettacode.org/wiki/Euler%27s_constant_0.5772...]]
_ul [[http://numbers.computation.free.fr/Constants/Gamma/gamma.html|http://numbers.computation.free.fr/Constants/Gamma/gamma.html]]

;; { require lib_BN}

_h2 1) from the definition

_h3 1.1) wren version
{pre
var eps = 1e-6

System.print("From the definition, err. 3e-10")
var n = 400
var h = 1
for (k in 2..n) h = h + 1/k
//faster convergence: Negoi, 1997
var a = (n + 0.5 + 1/(24*n)).log

Fmt.print("Hn    $0.14f", h)
Fmt.print("gamma $0.14f\nk = $d\n", h - a, n)

From the definition, err. 3e-10
Hn    6.56992969117651
gamma 0.57721566457657
k = 400
}

_h2 1.2) lambdatalk version

{pre
'{def negoi
 {lambda {:n}
  {let { {:n :n}
         {:h {+ {S.map {lambda {:k} {/ 1 :k}} {S.serie 1 :n}}} }
         {:a {log {+ :n 0.5 {/ 1 {* 24 :n}}}}}   // Negoi, 1997
       } {div}-> Hn :h
         {div}gamma {- :h :a}
         {div}k :n
}}}
-> {def negoi
 {lambda {:n}
  {let { {:n :n}
         {:h {+ {S.map {lambda {:k} {/ 1 :k}} {S.serie 1 :n}}} }
         {:a {log {+ :n 0.5 {/ 1 {* 24 :n}}}}}
       } {div}-> Hn :h
         {div}gamma {- :h :a} with k = :n
}}}

'{negoi 400}
{negoi 400}
     (0.57721566457657)
}

_h2 2) Sweeney, 1963

_h3 2.1) wren version

{pre
System.print("Sweeney, 1963, err. idem")
n = 21
var s = [0, n]
var r = n
var k = 1
while (true) '{
    k = k + 1
    r = r * n / k
    s[k & 1] = s[k & 1] + r/k
    if (r <= eps) break
}
Fmt.print("gamma $0.14f\nk = $d\n", s[1] - s[0] - n.log, k)

Sweeney, 1963, err. idem
gamma 0.57721566456363
k = 68
}

_h3 2.2) lambdatalk version

{pre
'{def sweeney
 {def sweeney.set!
  {lambda {:s :r :k :i}
   {A.set! :i {+ {A.get :i :s} {/ :r :k}} :s}
 }} 
 {def sweeney.loop
  {lambda {:n :s :r :k}
   {if {<= :r 1.e-10}
    then gamma = {- {A.get 1 :s} {A.get 0 :s} {log :n}} with k=:k
    else {sweeney.loop :n          
                       {sweeney.set! :s {* :r {/ :n :k}} :k {% :k 2}}      
                       {* :r {/ :n :k}}
                       {+ :k 1} }
 }}}
 {lambda {:n}
  {sweeney.loop :n {A.new 0 :n} :n 2} }} 
-> {def sweeney
 {def sweeney.set!
  {lambda {:s :r :k :i}
   {A.set! :i {+ {A.get :i :s} {/ :r :k}} :s}
 }} 
 {def sweeney.loop
  {lambda {:n :s :r :k}
   {if {<= :r 1.e-10}
    then gamma = {- {A.get 1 :s} {A.get 0 :s} {log :n}} with k=:k
    else {sweeney.loop :n          
                       {sweeney.set! :s {* :r {/ :n :k}} :k {% :k 2}}      
                       {* :r {/ :n :k}}
                       {+ :k 1} }
 }}}
 {lambda {:n}
  {sweeney.loop :n {A.new 0 :n} :n 2} }}

'{sweeney 21}
-> {sweeney 21}  
          (0.57721566456363)
}

_h2 3) Bailey, 1968

_h2 3.1) freeBasic version

{pre
n = 5
a = 1
h = 1
n2 = 2^n
r = 1
k = 1
do
   k += 1
   r *= n2 / k
   h += 1 / k
   b = a: a += r * h
loop until abs(b - a) < eps
a *= n2 / exp(n2)

? "gamma"; a - n * log(2); !"\nk ="; k
?
}

_h3 3.2) lambdatalk version

{pre

'{def bailey
 {def bailey.exit
  {lambda {:a :b :k}
   {or {< {abs {- :b :a}} 1.e-10} {> :k 100}}}}
 {def bailey.loop
  {lambda {:n :a :h :r :k}
   {if {bailey.exit {+ :a {* :r :h}} :a :k}
    then {- {* :a {/ {pow 2 :n} {exp {pow 2 :n}}}} {* :n {log 2}}}  with k=:k
    else {bailey.loop :n
                      {+ :a {* :r :h}} 
                      {+ :h {/ 1 :k}}
                      {* :r {/ {pow 2 :n} :k}}
                      {+ :k 1}} }}}
 {lambda {:n}
  {bailey.loop :n 1 1 1 2}}}
-> {def bailey
 {def bailey.exit
  {lambda {:a :b :k}
   {or {< {abs {- :b :a}} 1.e-10} {> :k 100}}}}
 {def bailey.loop
  {lambda {:n :a :h :r :k}
   {if {bailey.exit {+ :a {* :r :h}} :a :k}
    then {- {* :a {/ {pow 2 :n} {exp {pow 2 :n}}}} {* :n {log 2}}}  with k=:k
    else {bailey.loop :n
                      {+ :a {* :r :h}} 
                      {+ :h {/ 1 :k}}
                      {* :r {/ {pow 2 :n} :k}}
                      {+ :k 1}} }}}
 {lambda {:n}
  {bailey.loop :n 1 1 1 2}}}

'{bailey 5} 
-> {bailey 5}
  (0.57721566490153286...)
}

_h2 4) Brent-McMillan, 1980

_h3 4.1) freeBasic version

{pre
n = 13

a = -log(n)
b = 1
u = a
v = b
n2 = n * n
k2 = 0
k = 0
do
   k2 += 2*k + 1
   k += 1
   a *= n2 / k
   b *= n2 / k2
   a = (a + b) / k
   u += a
   v += b
loop until abs(a) < eps

? "gamma"; u / v; !"\nk ="; k
?
}

_h3 4.2) lambdatalk version

{pre
...
}

_h2 5) how Euler did it in 1735

_h3 5.1) wren version

{pre
System.print("How Euler did it in 1735")
// Bernoulli numbers with even indices
var b2 = [1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510, 43867/798]
var m = 7
n = 10
// n'th harmonic number
h = 1
for (k in 2..n) h = h + 1/k
Fmt.print("Hn    $0.14f", h)
h = h - n.log
Fmt.print("  -ln $0.14f", h)
// expansion C = -digamma(1)
a = -1 / (2*n)
n2 = n * n
r = 1
for (k in 1..m) '{
    r = r * n2
    a = a + b2[k] / (2*k*r)
}
Fmt.print("err  $0.14f\ngamma $0.14f\nk = $d", a, h + a, n + m) 
System.print("\nC = 0.57721566490153286...")

How Euler did it in 1735
Hn    2.92896825396825
  -ln 0.62638316097421
err  -0.04916749607268
gamma 0.57721566490153
k = 17

}

_h3 lambdatalk version

{pre
...
}