Created
November 14, 2013 01:21
-
-
Save cab1729/7459622 to your computer and use it in GitHub Desktop.
Bernoulli functions with arbitrary precision. Uses MPFloat.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)); | |
} | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
A version of this code using double is in functions.java