/* 

  Euler's constant (Rosetta code) in Picat.

  http://rosettacode.org/wiki/Euler%27s_constant_0.5772...
  """
  Compute the Euler constant 0.5772...

  Discovered by Leonhard Euler around 1730, it is the most ubiquitous mathematical 
  constant after pi and e, but appears more arcane than these.

  Denoted gamma (γ), it measures the amount by which the partial sums of the harmonic 
  series (the simplest diverging series) differ from the logarithmic function (its 
  approximating integral): lim n → ∞ (1 + 1/2 + 1/3 + … + 1/n − log(n)).

  The definition of γ converges too slowly to be numerically useful, but in 1735 Euler 
  himself applied his recently discovered summation formula to compute ‘the notable number’ 
  accurate to 15 places. For a single-precision implementation this is still the most economic algorithm.

  In 1961, the young Donald Knuth used Euler's method to evaluate γ to 1271 places. Knuth 
  found that the computation of the Bernoulli numbers required in the Euler-Maclaurin formula 
  was the most time-consuming part of the procedure.

  The next year Dura Sweeney computed 3566 places, using a formula based on the expansion of 
  the exponential integral which didn't need Bernoulli numbers. It's a bit-hungry method though: 
  2d digits of working precision obtain d correct places only.

  This was remedied in 1988 by David Bailey; meanwhile Richard Brent and Ed McMillan had published 
  an even more efficient algorithm based on Bessel function identities and found 30100 places in 20 hours time.

  Nowadays the old records have far been exceeded: over 6·1011 decimal places are already known. 
  These massive computations suggest that γ is neither rational nor algebraic, but this is yet to be proven.
  """  

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/


main => go.
% From mathematica
%   12345678901234567890
% 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674
% 0.577215669900188
go =>
  Gamma = 0.57721566490153286060651209008240,
  println(Gamma),
  foreach(N in 1..8)
     G = e(10**N),
     println([n=N,g=G,diff=G-Gamma])
  end,
  nl.

e(N) = [1.0/I : I in 1..N].sum-log(N).

e2(N) = E-log(N) =>
  E = 1,
  foreach(I in 2..N)
    E := E + 1/I
  end.

e3(N) = E-log(N) =>
  e3(1,N,0,E).
e3(N,N,E,E).
e3(I,N,E0,E) :-
  E1 = E0 + 1/I,
  e3(I+1,N,E1,E).

% {{trans|Rust}}
go2 =>
  Gamma = 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674,
  println(gamma=Gamma),
  member(N, 1..23),
  G = gamma(N),
  println([n=N,g=G,diff=G-Gamma]),
  fail,
  nl.

gamma(N) = Gamma =>
  Gamma = 1/2 - 1/3,
  foreach(I in 2..N)
    Power = 2**I,
    Sign = -1,
    Term = 0,
    foreach(Denominator in Power..(2*Power-1))
      Sign := Sign * -1,
      Term := Term + Sign / Denominator
    end,
    Gamma := Gamma + I*Term
  end.