Is there a reference for this approximation to Euler's constant, and is there a faster expression?

analysis/euler.constant New, unsolved

New entry, edited from material collected for the archive and not yet in a released version.

Problem

Is there a reference for this approximation to Euler's constant, and is
there a faster expression?

Euler's constant (0.5772...) can be approximated by the formula:

        N                  K
C(N) = sum(1/n) - ln(N) + sum(B[k]/(kN^k))
       n=1                k=1

where B[k] is the kth Bernoulli number (-1/2, 1/6, 0, -1/30, ...)
and K is chosen not larger than 2*pi*N

I discovered this experimentally (on a hand calculator).  Does anyone
have a reference?  Is there a faster expression?

-- Bob Jewett (hplabs!jewett)

Show solution

No solution was found in the collected postings.

[ Editor: the poster added a program for the UNIX calculator "bc": ]

For those of you into calculating lots of digits, the following input
to the UNIX command "bc" will produce C(N) for N=5, 10, 100, 200.  The
last is accurate to 70 digits.

Save the following in a file, then do "bc < filename"
___________________  cut here  _____________________________
scale=80
length(80)
define l(x){
        auto a, b, c, d, e, f, g, u, s, t
        if(x <=0) return(1-10^scale)
        t = scale
        scale = 0
        f = 1
        s = x
        while(s > 0){
                s = s/10
                f = f + 1
        }
        scale = t + f
        f=1
        while(x > 2){
                x = sqrt(x)
                f=f*2
        }
        while(x < .5){
                x = sqrt(x)
                f=f*2
        }
        u = (x-1)/(x+1)
        s = u*u
        b = 2*f
        c = b
        d = 1
        e = 1
        for(a=3;1==1;a=a+2){
                b=b*s
                c=c*a+d*b
                d=d*a
                g=c/d
                if(g==e){
                        scale = t
                        return(u*c/d)
                }
                e=g
        }
}
define e(x){
auto n,s,y,a,b,c
s=0
for(n=1; n<=x; n++){
  s=s+1/n
  }
x=n-1
y=x*x
a=y*y
b=a*a
c=b*b
s=s-l(x)-1/2/x+1/12/y-1/120/a+1/252/a/y
s=s-1/240/b+1/132/b/y-691/32760/b/a
s=s+1/12/b/a/y-3617/510/16/c+43867/798/18/c/y
s=s-174611/330/20/c/a+854513/138/22/c/a/y
s=s-236364091/2730/24/c/b+8553103/6/26/c/b/y
s=s-23749461029/870/28/c/b/a+8615841276005/14322/30/c/b/a/y
return(s)
}

scale=15
e(5)
scale=25
e(10)
scale=57
e(100)
scale=67
e(200)

-- Bob Jewett (hplabs!jewett)