What is the fastest/most efficient algorithm for estimating Euler's Constant $\gamma$? What is the fastest algorithm for estimating Euler's Constant $\gamma \approx0.57721$?
Using the definition:
$$\lim_{n\to\infty} \sum_{x=1}^{n}\frac{1}{x}-\log n=\gamma$$
I finally get $2$ decimal places of accuracy when $n\geq180$.  The third correct decimal place only comes when $n \geq638$.  Clearly, this method is not very efficient (it can be expensive to compute $\log$).
What is the best method to use to numerically estimate $\gamma$ efficiently?
 A: As it turns out, the convergence of the Karatsuba series presented in Matthew's answer can be improved. This time, the geometric behavior of the error (as can be ascertained from the bounds presented) can be exploited through the use of the Shanks transformation. (Richardson can be made to work here as well, but the results are not as spectacular.)
Letting
$$\varepsilon_0^{(k)}=1-\log(k+1) \sum_{r=1}^{12k+13} \frac{ (-k)^{r+1}}{(r-1)!(r+1)} + \sum_{r=1}^{12k+13} \frac{ (-k)^{r+1} }{(r-1)!(r+1)^2}$$
Wynn's version of the Shanks transformation uses the recursion
$$\varepsilon_{k+1}^{(n)}=\varepsilon_{k-1}^{(n+1)}+\frac1{\varepsilon_{k}^{(n+1)}-\varepsilon_k^{(n)}}$$
It would seem that a two-dimensional array would be required for implementation, but one can arrange things such that only a one-dimensional array is required, through clever overwriting. Here is a Mathematica routine to demonstrate:
wynnEpsilon[seq_?VectorQ] := Module[{n = Length[seq], ep, res, v, w},
  res = {};
  Do[
   ep[k] = seq[[k]];
   w = 0;
   Do[
    v = w; w = ep[j];
    ep[j] = 
     v + (If[Abs[ep[j + 1] - w] > 10^-(Precision[w]), ep[j + 1] - w, 
         10^-(Precision[w])])^-1;
    , {j, k - 1, 1, -1}];
   res = {res, ep[If[OddQ[k], 1, 2]]};
   , {k, n}];
  Flatten[res]
  ]

(actually the same as the routine presented in this answer).
Here's a comparison of Karatsuba's series, with and without Shanks transformation:
gamprox = Table[N[1 - Log[k]*Sum[(-k)^(r + 1)/((r + 1)*(r - 1)!),
            {r, 1, 12*k + 1}] + Sum[(-k)^(r + 1)/((r + 1)^2*(r - 1)!),
            {r, 1, 12*k + 1}], 50], {k, 30}];

trans = wynnEpsilon[gamprox];

gamprox[[20]] - EulerGamma // N
1.31827*10^-7

trans[[20]] - EulerGamma // N
6.49869*10^-18

Last[gamprox] - EulerGamma // N
9.96301*10^-12

Last[trans] - EulerGamma // N
2.07059*10^-27

Not too shabby, in my humble opinion...
A: I quite like the Brent-McMillan algorithm myself (which is based on the relationships between the Euler-Mascheroni constant and modified Bessel functions):
$$\gamma=\lim_{n\to\infty}\mathscr{G}_n=\lim_{n\to\infty}\frac{\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2 (H_k-\log\,n)}{\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2}$$
where $H_k=\sum\limits_{j=1}^k \frac1{j}$ is a harmonic number.
It requires the use of a logarithm, but the algorithm is quite simple and reasonably efficient (in particular, we have the inequality $0 < \mathscr{G}_n-\gamma < \pi\exp(-4n)$).
Here's some Mathematica code for the Brent-McMillan algorithm (which should be easily translatable to your language of choice):
n = 50;
a = u = N[-Log[n], n]; b = v = 1;
i = 1;
While[True,
  k = (n/i)^2;
  a *= k; b *= k;
  a += b/i;
  If[u + a == u || v + b == v, Break[]];
  u += a; v += b;
  i++
  ];
u/v

The integer parameter n controls the accuracy; very roughly, the algorithm will yield n-2 or so correct digits.
The Brent-McMillan paper also presents more elaborate schemes for computing $\gamma$, such as
$$\gamma=\lim_{n\to\infty}\frac{\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2 (H_k-\log\,n)}{\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2}-\frac{\frac1{4n}\sum\limits_{k=0}^\infty \frac{(2k)!^3}{k!^4 (16n)^{2k}}}{\left(\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2\right)^2}$$
but I have no experience in using them.
A: The paper "On the computation of the Euler constant $\gamma$" by Ekatharine A. Karatsuba, in Numerical Algorithms 24(2000) 83-97, has a lot to say about this.  This link might work for you.
In particular, the author shows that for $k\ge 1$,
$$ 
\gamma= 1-\log k \sum_{r=1}^{12k+1} \frac{ (-1)^{r-1} k^{r+1}}{(r-1)!(r+1)} +
\sum_{r=1}^{12k+1} \frac{ (-1)^{r-1} k^{r+1} }{(r-1)! (r+1)^2}+\mbox{O}(2^{-k})$$
and more explicitly
$$\begin{align*}
-\frac{2}{(12k)!} - 2k^2 e^{-k} \le \gamma 
-1+&\log k \sum_{r=1}^{12k+1} \frac{ (-1)^{r-1} k^{r+1}}{(r-1)!(r+1)}
- \sum_{r=1}^{12k+1} \frac{ (-1)^{r-1} k^{r+1} }{(r-1)! (r+1)^2}\\
&\le \frac{2}{(12k)!} + 2k^2 e^{-k}\end{align*}$$
for $k\ge 1$. 
Since the series has fast convergence, you can use these to get good approximations to $\gamma$ fairly quickly.
A: I do not know about the best method, however numerically evaluating the integral
$$\gamma = - \int_0^1\!dx\,\ln \ln x^{-1}$$
seems to be pretty efficient.
A: There is another interesting formula
 $$\small 1- \gamma = \sum_{k=2}^\infty {\zeta(k)-1\over k}$$ 
found in mathworld (see eq 123) .
If we simply use approximations to the zetas by truncating their series, and write this in an array
$\small \begin{array} {lll}
 1 & 1 & 1 & 1 & \cdots & 1 \\
 {1 \over 2^2} &  {1 \over 2^3}  &  {1 \over 2^4}  &  {1 \over 2^5}  & \cdots&{1 \over 2^c}\\
 {1 \over 3^2} &  {1 \over 3^3}  &  {1 \over 3^4}  &  {1 \over 3^5}  & \cdots&{1 \over 2^c}\\
 \cdots &  \cdots &  \cdots &  \cdots &  \cdots & &\\
 {1 \over r^2} &  {1 \over r^3}  &  {1 \over r^4}  &  {1 \over r^5}  & \cdots&{1 \over r^c}\\
\hline
\zeta_r(2)&\zeta_r(3)&\zeta_r(4)&\zeta_r(5)&\cdots&\zeta_r(c)&
 \end{array} $
then we can write an approximation-formula for the Euler $\small \gamma$ 
  $$\small 1-\gamma_{r,c} = \sum_{k=2}^c {\zeta_r(k)-1\over k}$$ 
which depends on the number of rows r and the number of columns c . Now to reduce the number of coefficients needed to arrive at a good approximation 


*

*we can use the alternating (column-)sums and convert by the eta/zeta-conversion term    

*additionally we can use Eulersummation for convergence acceleration for the (now alternating) $\small \zeta_r(c) $      

*we can even introduce Euler-summation of (small) negative order to accelerate convergence of the sum of zetas (which itself is non-alternating).       
If we use all three accelerations, we get a double sum 
  $$\small 1-\gamma_{r,c} = \sum_{k=2}^c \sum_{j=1}^r a_{j,k}{ 1 \over j^k}$$ 
where the $\small a_{j,k} $ contain the factors due to the denominator in the $\small \gamma$-formula and due to the threefold convergence-acceleration.        
I did actually implement this in Pari/GP and the surprising result was, that the best approximations were (using order 0.5 in the Eulersummation for the columns and -0.25 for the Eulersummation of the approximated zetas), if roughly r=c . Then the number of correct digits were about r/2; so with r=64 and c=64 we get $\small \gamma$ to 31 digits accuracy.
So the effort comes out to be 
$$\small  \text{ # of correct digits} \sim r/2 \qquad \text{ if } r \sim c $$
The cost of computation of the complete array of zeta-terms is thus in principle quadratic in d (the required number of correct digits); for the Euler-sums a vector for the column-acceleration and another vector for the row-acceleration is required whose values can recursively be computed and are thus linear with the number of rows resp the number of columns and thus also linear with d. (The convergence-acceleration (1.) by using the alternating sums costs nearly nothing)
A: No heavy machinery or special series manipulation required for this. Simple enough to do on a standard scientific calculator.
By applying the Euler-Maclaurin summation formula, one finds that
$$\lim_{b\to\infty}H_b-H_{n-1}-\ln(b)+\ln(n)=\frac1{2n}+\left(\sum_{j=1}^p\frac{B_{2j}}{2j\times n^{2j}}\right)+R_{n,p}$$
and thus, we may derive
$$\gamma=H_{n-1}-\ln(n)+\frac1{2n}+\left(\sum_{j=1}^p\frac{B_{2j}}{2j\times n^{2j}}\right)+R_{n,p}$$
where
$$H_n=\sum_{k=1}^n\frac1k$$
$$|R_{n,p}|\le\frac{2\zeta(2p)(2p-1)!}{(2\pi n)^{2p}}$$
where $\zeta$ is the Riemann zeta function. It suffices to use
$$\zeta(s)<\frac1{1-2^{1-s}}\left(1-\frac1{2^s}\right)$$
Or more simply,
$$\zeta(s)<2$$
both for $s\ge2$.
For example, with $n=5$ and $p=4$, we can approximate out $7$ places:
$$\gamma=H_4-\ln(5)+\frac1{10}+\frac1{300}-\frac1{75000}+\frac1{3937500}-\frac1{93750000}+R_{5,4}$$
$$\gamma\approx0.577215664+R_{5,4}$$
where
$$|R_{5,4}|<1.1\times10^{-8}$$
A: A good place for fast evaluation of constants is Gourdon and Sebah's 'Numbers, constants and computation'. 
They got $108\cdot 10^6$ digits for $\gamma$ in 1999 (see the end of their 2004 article 'The Euler constant') and propose a free program for high precision evaluation of various constants 'PiFast'.
On his page of constants Simon Plouffe has Euler's constants to 10^6 digits (the file looks much smaller sorry...) using Brent's splitting algorithm (see the 1980 paper of Brent 'Some new algorithms for high-precision computation of Euler’s constant' or more recently 3.1 in Haible and Papanikolaou's 'Fast multiprecision evaluation of series of rational numbers'). 
It seems that the 1999 record was broken in 2009 by A. Yee & R. Chan with 29,844,489,545 digits 'Mathematical Constants - Billions of Digits' (warning: the torrent file proposed there is more than 11Gb large! An earlier 52Mb file of 'only' 116 million digits is available here using the method proposed by Gourdon and Sebah).
A: (N.B. The previous version of this answer featured both the Brent-McMillan algorithm and the acceleration of Macys's series; I have decided to move the Brent-McMillan material into a new answer in the interest of having only one method per answer.)

The convergence properties of Macys's series in Will's answer can be improved a fair bit, if you're willing to devote some amount of computational effort; due to the $n^{-2}$ behavior of the error, one obvious choice for a convergence acceleration method is Richardson extrapolation.
Skipping some hairy details (which I might include later if I find time, but see Marchuk/Shaidurov if you must), the working formula is
$$\gamma=\lim_{n\to\infty} G_n=2\lim_{n\to\infty} \sum_{i=1}^{n+1} \frac{(-1)^{n-i} i^{2n+2}}{(n+i+1)!(n-i+1)!}\left(\sum_{k=i+1}^{i(i+1)} \frac1{k}-\sum_{k=1}^i \frac1{k}\right)$$
Here are some sample results:
$$\begin{array}{ccc}n&G_n&\gamma-G_n\\10&0.577210083083&5.581818\times10^{-6}\\50&0.577215659731&5.170456\times10^{-9}\\100&0.577215664665&2.362333\times10^{-10}\\200&0.577215664891&1.061648\times10^{-11}\\250&0.577215664898&3.902515\times10^{-12}\\300&0.577215664900&1.721878\times10^{-12}\\350&0.577215664901&8.618620\times10^{-13}\end{array}$$
For higher precision, there isn't much of an improvement; I would still recommend Brent-McMillan if one needs many digits of $\gamma$.
A: Finch's Mathematical Constants mentions these papers:


*

*D. W. DeTemple, A quicker convergence to Euler's constant, Amer. Math. Monthly (1993) 468-470.

*T. Negoi, A faster convergence to Euler's constant, Math. Gazette 83 (1999) 487-489.

