110 lines
3.3 KiB
Java
110 lines
3.3 KiB
Java
|
/*******************************************************************************
|
||
|
|
* Copyright (c) 2015-2018 Skymind, Inc.
|
||
|
|
*
|
||
|
|
* 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.
|
||
|
|
*
|
||
|
|
* 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.linalg.util;
|
||
|
|
|
||
|
|
/*
|
||
|
|
* To change this template, choose Tools | Templates and open the template in the editor.
|
||
|
|
*/
|
||
|
|
|
||
|
|
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<Rational> a = new ArrayList<Rational>();
|
||
|
|
|
||
|
|
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 */
|