Commit 813c6de9 authored by ulrich_y's avatar ulrich_y

Added algorithm for the n-th Bernoulli number

parent 819d77bd
......@@ -2,7 +2,7 @@
MODULE maths_functions
use globals
use ieps
use utils, only: factorial
use utils, only: factorial, binom
implicit none
interface polylog
module procedure polylog1, polylog2
......@@ -46,6 +46,40 @@ CONTAINS
end do
END FUNCTION naive_polylog
FUNCTION bernoullinumber(n)
! This returns the n-th Bernoulli number by computing all Bernoulli numbers
! up to the n-th recursively using the relation
! Sum[ Binomial[m+1, k] BernoulliB[k], {k,0,m} ] = 0
! for m > 0 (https://mathworld.wolfram.com/BernoulliNumber.html).
! Solving this for BernoulliB[m] results in
! BernoulliB[m] = - Sum[Binomial[m, k] BernoulliB[k] / (m-k+1), {k,0,m-1}]
! for m > 0 and BernoulliB[0] = 1.
! Care is taken to avoid multiple computation.
integer, intent(in) :: n
real(kind=prec) :: bernoullinumber
integer, parameter :: nmax = 40
! keep track of the bernoulli numbers we have already calculated
real(kind=prec), save :: bernoulli(0:nmax) = 0.
integer, save :: cmax = 2
integer m, k
!TODO error handling or dynamic allocation for bernoulli(:)
bernoulli(0) = 1.
bernoulli(1) = -0.5
do m=cmax,n,2
bernoulli(m) = 0.
do k=0,m-1
bernoulli(m) = bernoulli(m) - binom(m, k) * bernoulli(k) / (m - k + 1)
enddo
enddo
cmax=m
bernoullinumber = bernoulli(n)
END FUNCTION bernoullinumber
FUNCTION Li2(x)
!! Dilogarithm for arguments x < = 1.0
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment