Bernoulli functions with arbitrary precision. Uses MPFloat.

Bernoulli.java

import java.io.BufferedReader;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;

import jgmp.MPFloat;

import de.luschny.math.factorial.FactorialPoorMans;

/**
 * @author jmenes
 * compute bernoulli functions using MPFloat
 *
 */
public class Bernoulli {
	
	// duplicate constants from mpfunctions here so it doesn't need to load it
	private static final MPFloat PI = 
			new MPFloat(
				"3.14159265358979323846264338327950288419716939937510582097494459230781640628620899" + 
				"8628034825342117067982148086513282306647093844609550582231725359408128481117450e+00", 
				10, 512);
	
	private static final MPFloat TWO_PI = PI.mult(2);
	private static final MPFloat HALF_PI = PI.div(2);
	
	
	private static final int[] PRIMES =
			{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,
		 	61,67,71,73,79,83,89,97,101,103,107,109,113,127,
		 	131,137,139,149,151,157,163,167,173,179,181,191,
		 	193,197,199,211,223,227,229,233,239,241,251,257,
		 	263,269,271,277,281,283,293,307,311,313,317,331,
		 	337,347,349,353,359,367,373,379,383,389,397,401,
		 	409,419,421,431,433,439,443,449,457,461,463,467,
		 	479,487,491,499,503,509,521,523,541,547,557,563,
		 	569,571,577,587,593,599,601,607,613,617,619,631,
		 	641,643,647,653,659,661,673,677,683,691,701,709,
		 	719,727,733,739,743,751,757,761,769,773,787,797,
		 	809,811,821,823,827,829,839,853,857,859,863,877,
		 	881,883,887,907,911,919,929,937,941,947,953,967,
		 	971,977,983,991,997};
	
	private static int[] primes = PRIMES;	// default values in case of file error
	
	private static int pl = primes.length;	// default values in case of file error
	
	private static String[] FACS = new String[5001]; // table of the first 5k factorials
	private static int facl = 0;
	
	private static FactorialPoorMans fpm = new FactorialPoorMans();
	
	static {
		// load primes array with PARI/GP generated file (5000 primes)
		try {
			BufferedReader pf = 
				new BufferedReader(new FileReader("primes.txt"));
			String ps = pf.readLine();
			String[] pe = ps.split(", ");
			int pel = pe.length;
			primes = new int[pel];
			for (int i=0; i<pel; i++) {
				primes[i] = Integer.parseInt(pe[i]);
			}
			pl = primes.length;
			pf.close();
			
		} catch (FileNotFoundException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		
		// load array with first 5000 factorials - n is index
		try {
			BufferedReader ff = 
				new BufferedReader(new FileReader("facs.txt"));
			String fac = ff.readLine();
			FACS = fac.split(", ");
			facl = FACS.length;
			ff.close();
			
		} catch (FileNotFoundException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		
//		FactorialPoorMans fpm = new FactorialPoorMans();
		if (facl==0) { // if facl==0 something went wrong with file load
			for (int n=0; n<5001; n++) {
				FACS[n] = fpm.factorial(n);
			}
		}
//		fpm = null;		// make it eligible for gc
		
		MPFloat.setDefaultPrec(256);
	}
	
	/**
	 * return nth Bernoulli number in rational form
	 * @param n
	 * @return
	 * ref: Computing Bernoulli Numbers Quickly, Kevin J. McGown, December 8, 2005
	 * 		http://modular.math.washington.edu/projects/168/kevin_mcgown/bernproj.pdf
	 */
	public static MPFloat BernoulliB(int n)
	{
		
		MPFloat Bn = new MPFloat(0.0); // return 0 if n is not even, >= 2
		if (n == 0)
			return new MPFloat(1.0);
		else if (n == 1)
			return new MPFloat(-0.5);
		else
			if ((n % 2) != 0)
				return Bn;
		
		MPFloat K = new MPFloat(factorial(n));
		K.mult(2, K);
		K.div(TWO_PI.pow(n), K);
		
		MPFloat d = new MPFloat(0.0);
		int prod = 1;
		for (int i = 0; i < pl; i++)
		{
			if ((n % (primes[i]-1)) == 0)
				prod *= primes[i];
		}
		d.fromString(Integer.toString(prod), 10);
		
		MPFloat N = (K.mult(d));
		N.pow(1/(n-1), N);
		N.ceil(N);
		
		//double z = 1.0;
		MPFloat z = new MPFloat(1.0);
		int ni = N.intValue();
		for (int i = 0; i < pl; i++)
		{
			if (primes[i] > ni)
				break;
			//if (i == 0)
			//{
			//	z.sub(new MPFloat(StrictMath.pow(primes[i], -n)).pow(-1), z);
			//}else {
				
				z.mult(new MPFloat(StrictMath.pow(1-primes[i], -n)).pow(-1), z);
			//}
		}
						
		//double a = StrictMath.pow((-1), (n/2)+1)*StrictMath.ceil(d*K*z);
		MPFloat dkzc = d.mult(K);
		dkzc.mult(z, dkzc);
		dkzc.ceil(dkzc);
		MPFloat a = 
			new MPFloat(StrictMath.pow((-1), (n/2)+1));
		a.mult(dkzc, a);
		
		Bn = a.div(d);
		
		return Bn;
	}
	
	/**
	 * return nth Bernoulli polynomial of x using Fourier expansion
	 * ref: Abramowitz & Stegun 23.1.16 pp 805
	 * @param n
	 * @param x
	 * @return
	 */
	public static MPFloat BernoulliP(int n, MPFloat x) {
		
		// for n>1, 1>=x>=0
		// for n=1, 1>x>0
		
		// for x=0 return the nth Bernoulli number
		if (x.doubleValue()==0.0) {
			return BernoulliB(n);
		}
		
		// for x=n return x
		if (n == x.doubleValue()) {
			return x;
		}
		
		MPFloat K = new MPFloat(-2.0);
		K.mult(new MPFloat(factorial(n)).div(TWO_PI.pow(n)), K);
		
		MPFloat S = new MPFloat();
		for (int k=1; k<=5000; k++) {
			MPFloat cosp = TWO_PI.mult(k);
			cosp.mult(x, cosp);
			cosp.sub(HALF_PI.mult(n), cosp);
			S.add(mpfunctions.cos(cosp)
					.div(new MPFloat(StrictMath.pow(k, n))), S);
		}

		return K.mult(S);
	}
	
	/**
	 * Periodic Bernoulli function - Bn({x})
	 * @param n
	 * @param x
	 * @return
	 */
	public static MPFloat BBn(int n, MPFloat x) {
		
		// pass the fractional part of x
		return BernoulliP(n, x.sub(x.floor()));
		
	}
	
	public static String factorial(int N) {
		if (N>=facl)
			return fpm.factorial(N);
		
		return FACS[N];
	}
	
	
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		for (int i=1; i<=100; i++) {
			System.out.println("B(" + i*2 + ") = " + Bernoulli.BernoulliB(i*2));
		}
		

	}

}

Dependencies:
https://github.com/PeterLuschny/Fast-Factorial-Functions/blob/master/JavaFactorial/src/de/luschny/math/factorial/FactorialPoorMans.java

For primes.txt, see: https://gist.github.com/cab1729/1318030

MPFloat: https://gist.github.com/1562150

Updated pdf link: http://wstein.org/projects/168/kevin_mcgown/bernproj.pdf

A version of this code using double is in functions.java