static public BigDecimal gamma ( MathContext mc ) { /* look it up if possible */ if ( mc . getPrecision ( ) < GAMMA . precision ( ) ) { return GAMMA . round ( mc ) ; } else { double eps = prec2err ( 0.577 , mc . getPrecision ( ) ) ; /* Euler-Stieltjes as shown in Dilcher, Aequat Math 48 (1) (1994) 55-85 14 */ MathContext mcloc = new MathContext ( 2 + mc . getPrecision ( ) ) ; BigDecimal resul = BigDecimal . ONE ; resul = resul . add ( log ( 2 , mcloc ) ) ; resul = resul . subtract ( log ( 3 , mcloc ) ) ; /* how many terms: zeta-1 falls as 1/2^(2n+1), so the * terms drop faster than 1/2^(4n+2). Set 1/2^(4kmax+2) < eps. * Leading term zeta(3)/(4^1*3) is 0.017. Leading zeta(3) is 1.2. Log(2) is 0.7 */ int kmax = ( int ) ( ( Math . log ( eps / 0.7 ) - 2. ) / 4. ) ; mcloc = new MathContext ( 1 + err2prec ( 1.2 , eps / kmax ) ) ; for ( int n = 1 ; ; n ++ ) { /* zeta is close to 1. Division of zeta-1 through * 4^n*(2n+1) means divion through roughly 2^(2n+1) */ BigDecimal c = zeta ( 2 * n + 1 , mcloc ) . subtract ( BigDecimal . ONE ) ; BigInteger fourn = BigInteger . valueOf ( 2 * n + 1 ) ; fourn = fourn . shiftLeft ( 2 * n ) ; c = divideRound ( c , fourn ) ; resul = resul . subtract ( c ) ; if ( c . doubleValue ( ) < 0.1 * eps ) { break ; } } return resul . round ( mc ) ; } }
Euler - Mascheroni constant .