org.nd4j.common.util.Bernoulli Maven / Gradle / Ivy
/*
* ******************************************************************************
* *
* *
* * This program and the accompanying materials are made available under the
* * terms of the Apache License, Version 2.0 which is available at
* * https://www.apache.org/licenses/LICENSE-2.0.
* *
* * See the NOTICE file distributed with this work for additional
* * information regarding copyright ownership.
* * Unless required by applicable law or agreed to in writing, software
* * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
* * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
* * License for the specific language governing permissions and limitations
* * under the License.
* *
* * SPDX-License-Identifier: Apache-2.0
* *****************************************************************************
*/
package org.nd4j.common.util;
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
/**
* Bernoulli numbers.
*/
class Bernoulli {
/*
* The list of all Bernoulli numbers as a vector, n=0,2,4,....
*/
static List a = new ArrayList();
public Bernoulli() {
if (a.isEmpty()) {
a.add(Rational.ONE);
a.add(new Rational(1, 6));
}
}
/**
* Set a coefficient in the internal table.
*
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* @param value the new value of the coefficient.
*/
protected void set(final int n, final Rational value) {
final int nindx = n / 2;
if (nindx < a.size()) {
a.set(nindx, value);
} else {
while (a.size() < nindx) {
a.add(Rational.ZERO);
}
a.add(value);
}
}
/**
* The Bernoulli number at the index provided.
*
* @param n the index, non-negative.
* @return the B_0=1 for n=0, B_1=-1/2 for n=1, B_2=1/6 for n=2 etc
*/
public Rational at(int n) {
if (n == 1) {
return (new Rational(-1, 2));
} else if (n % 2 != 0) {
return Rational.ZERO;
} else {
final int nindx = n / 2;
if (a.size() <= nindx) {
for (int i = 2 * a.size(); i <= n; i += 2) {
set(i, doubleSum(i));
}
}
return a.get(nindx);
}
}
/* Generate a new B_n by a standard double sum.
* @param n The index of the Bernoulli number.
* @return The Bernoulli number at n.
*/
private Rational doubleSum(int n) {
Rational resul = Rational.ZERO;
for (int k = 0; k <= n; k++) {
Rational jsum = Rational.ZERO;
BigInteger bin = BigInteger.ONE;
for (int j = 0; j <= k; j++) {
BigInteger jpown = BigInteger.valueOf(j).pow(n);
if (j % 2 == 0) {
jsum = jsum.add(bin.multiply(jpown));
} else {
jsum = jsum.subtract(bin.multiply(jpown));
}
/* update binomial(k,j) recursively
*/
bin = bin.multiply(BigInteger.valueOf(k - j)).divide(BigInteger.valueOf(j + 1));
}
resul = resul.add(jsum.divide(BigInteger.valueOf(k + 1)));
}
return resul;
}
} /* Bernoulli */