Commit 875c76cb authored by ulrich_y's avatar ulrich_y

Dynamic caching of bernoulli numbers

parent df3ee22c
......@@ -56,17 +56,36 @@ CONTAINS
! for m > 0 and BernoulliB[0] = 1.
! Care is taken to avoid multiple computation.
integer, intent(in) :: n
! the cache is dynamic, the first entry is the initial size, the
! second the increment
integer, parameter :: cachecontr(2) = (/ 20, 10 /)
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.
real(kind=prec), save, allocatable :: bernoulli(:)
real(kind=prec), allocatable :: buffer(:)
integer, save :: m = 2
integer k
!TODO error handling or dynamic allocation for bernoulli(:)
if (.not. allocated(bernoulli)) then
allocate(bernoulli(0:cachecontr(1)))
bernoulli(0) = 1.
bernoulli(1) = -0.5
endif
if (n > size(bernoulli)) then
! We need more buffer
k = size(bernoulli) - 1
allocate(buffer(0:k))
buffer = bernoulli
deallocate(bernoulli)
allocate(bernoulli(0:n + cachecontr(2)))
bernoulli(0:k) = buffer(0:k)
deallocate(buffer)
endif
bernoulli(0) = 1.
bernoulli(1) = -0.5
do m=m,n,2
bernoulli(m) = 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