Commit ec9e0ea2 authored by ulrich_y's avatar ulrich_y

Adding verbosity to bernoulli and harmonic subsystems

parent 875c76cb
......@@ -67,6 +67,9 @@ CONTAINS
integer k
if (.not. allocated(bernoulli)) then
#ifdef DEBUG
if(verb >= 150) print*, 'initialising Bernoulli number system for m=0,', cachecontr(1)
#endif
allocate(bernoulli(0:cachecontr(1)))
bernoulli(0) = 1.
bernoulli(1) = -0.5
......@@ -82,6 +85,12 @@ CONTAINS
allocate(bernoulli(0:n + cachecontr(2)))
bernoulli(0:k) = buffer(0:k)
#ifdef DEBUG
if(verb >= 150) then
print*, 'increasing Bernoulli number system from m=0,', k, &
' to m=0,',size(bernoulli)-1
endif
#endif
deallocate(buffer)
endif
......@@ -92,6 +101,11 @@ CONTAINS
do k=0,m-1
bernoulli(m) = bernoulli(m) - binom(m, k) * bernoulli(k) / (m - k + 1)
enddo
#ifdef DEBUG
if(verb >= 120) then
write(*,'(A,I3,A,E25.10)') "calculating B(",m,") = ",bernoulli(m)
endif
#endif
enddo
bernoullinumber = bernoulli(n)
......@@ -104,8 +118,17 @@ CONTAINS
real(kind=prec), save :: Harmonic(0:nmax) = 0
integer, save :: m = 0
if (n > nmax) then
print*,"Harmonic numbers bigger ",nmax," are not supported."
stop 1
endif
do m=m, n+1
harmonic(m+1) = harmonic(m) + 1._prec / real(m+1)
#ifdef DEBUG
if(verb >= 120) then
write(*,'(A,I3,A,E25.10)') "calculating H(",m,") = ",harmonic(m)
endif
#endif
enddo
harmonicnumber = harmonic(n)
......
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