FORTRAN program for harmonic series sum

[issacnewton]

Hi

Here is my code to get the sum of harmonic series. Harmonic series is

$$\sum_{i=1}^{\infty}\; \frac{1}{i}$$

here is the code

program harmonic

implicit none

integer :: i,n
real :: sum=0.0

write(*,*)'How many terms you want to sum ?'
read(*,*) n


do i= 1 ,n 
   

      sum=sum +(1.0/real(i))
      

end do

write(*,*) 'The sum of series is = ', sum


end program harmonic

Now if I give n= 10000000 , it gives me the answer 15.403683. But Mathematica gives me answer 16.6953. So there is something wrong happening in this simple code.

Can anybody point out the mistake ?

 

[Mark44]

I suspect that the problem is that your sum variable is of type real (or single precision), which is a four-byte type. See if changing the type of this variable to double precision, and write 1.0 as 1.0D0. You'll also need to convert your loop counter, i, to double precision, which you can do with the DBLE intrinsic function.

 

[AlephZero]

I think Mark44 is probably right.

But you might try the experiment of summing the terms smallest first, not largest, i.e. change
do i = 1,n
to
do i = n,1,-1
(and stay with real variables, not double precision).

If this "works", figuring out why it works is left as an exercise for you (and it's something important to know about floating-point calculations in general, not just for this situation).

FWIW A good approximation should be $\ln(10000000) + \gamma$ where $\gamma$ is the http://en.wikipedia.org/wiki/Euler–Mascheroni_constant which is the same as the Mathematica answer.

 

[gsal]

I would definitely start the loop with the largest value; otherwise, the little numbers are simply going to fall off.

The typical fortran integer is 4 bytes, so, no problem there for i or n to accommodate 10 million.

Actually, the problem is with the reals...they are also typically 4 bytes and maybe not large enough to accommodate little numbers; so, for the real, I would go double precision, like real*8.

Then, I would cast the loop variable right away and use the declared real variable in the division, like this

program harmonic
implicit none
integer :: i,n
real*8 :: s,sum=0.0

write(*,*)'How many terms you want to sum ?'
read(*,*) n

do i= n,1,-1 
    s = i
    sum = sum + 1.0/s
end do
write(*,*) 'The sum of series is = ', sum

end program harmonic

This way there is no need to be writing 1.0D0; instead, 1.0 would automatically be promoted to the same type as 's' in the denominator...which we already declared to be REAL*8

 

[issacnewton]

thanks everybody...

when I use the following code as suggested by alephzero ,

program harmonic

implicit none

integer :: i,n
real :: sum=0.0

write(*,*)'How many terms you want to sum ?'
read(*,*) n


do i= n ,1, -1 
   

      sum=sum +(1.0/real(i))
      

end do

write(*,*) 'The sum of series is = ', sum


end program harmonic

For n=10000000, it gives now the answer 16.686031.

Then I tried the solution suggested by gsal

program harmonic
implicit none
integer :: i,n
real*8 :: s,sum=0.0

write(*,*)'How many terms you want to sum ?'
read(*,*) n

do i= n,1,-1 
    s = i
    sum = sum + 1.0/s
end do
write(*,*) 'The sum of series is = ', sum

end program harmonic

now for n=10000000, I get the sum as 16.695311365859890 , which is very close
to the Mathematica answer, 16.6953. I don't know how to get more precise answers
in Mathematica. I used the command

Sum[1.0/i, {i,10000000}]

Next, I tried to alter gsal's code as

program harmonic
implicit none
integer :: i,n
real*8 :: s,sum=0.0

write(*,*)'How many terms you want to sum ?'
read(*,*) n

do i= 1, n 
    s = i
    sum = sum + 1.0/s
end do
write(*,*) 'The sum of series is = ', sum

end program harmonic

Now here, for n=10000000, I get the sum as 16.695311365856710 , which is different in
last 4 digits from the answer from gsal's code.

So using double precision variables certainly improves the precision , but the order in which
we sum the series seems to do some funny things...As commented by alephzero,
summing smallest terms first is probably more accurate. He has asked me to think over it
as an exercise. I found something on yahoo answers

http://in.answers.yahoo.com/question/index?qid=20110511051723AAPYwUB

Here, its suggested that

You could partly fix this by going to double precision, which is essentially pushing the problem a bit further without eliminating it.
Another way is to add only number of like magnitude, that is to divide the coefficient you are calculating in several variables, one for the thousands and above, one for the millionth, one for the billionth, and so on. Each new value in the series is split and added to the pertinent range variable, with a carry when anyone reaches a value that makes it 'spill' to the grater magnitude variable.
That is the technique that is used to compute numbers like Pi to an arbitrary large number of decimals

So how would I do this kind of splitting. Also I want to time my code to get the time the compiler takes to compute the result. I have Win XP as OS and I am using Ubuntu inside
VirtualBox. Compiler is gfortran...

thanks

 

[gsal]

If you are going to get into timing the program, that's another matter which I personally am not that good at...I am talking about setting things up according to what task takes longer, moving things in and out of the register, vectorizing loops, etc.

By the way, not every decimal number can be represented in binary with a given number of bits and so, just because you keep getting a larger and larger total sum, it does not mean it is correct...the correct number could very well be smaller than those large numbers, it is just that some of the resulting 1.0/i quantities may have had to be represented in binary in a number that when converted back to decimal is larger than 1.0/i ...do you get this?

For example, if you have a 3-bit binary system, you can exactly represent the numbers 0 thru 7...but how can you represent 5.5? You can't, you would have to chose either 101 or 110 ! Neither one of which converts back to 5.5

Also, you are going to have to choose between speed and accuracy. In Fortran, you can work with various precisions upon request...read up on "selected_real_kind".

Then, as mentioned above, there is the more elemental thing about moving things in and out of the registry for addition, leaving it there during the entire loop or what.


Example 1:

program harmonic
implicit none
integer :: i,n,d(5)
real*8 :: s,rsum(5),total

!write(*,*)'How many terms you want to sum ?'
!read(*,*) n
n=10000000

d(1) = 1*n/5
d(2) = 2*n/5
d(3) = 3*n/5
d(4) = 4*n/5
d(5) = 5*n/5

rsum=0.0
do i= 0,n/5-1
    rsum(1) = rsum(1) + 1.0/(d(1)-i)
    rsum(2) = rsum(2) + 1.0/(d(2)-i)
    rsum(3) = rsum(3) + 1.0/(d(3)-i)
    rsum(4) = rsum(4) + 1.0/(d(4)-i)
    rsum(5) = rsum(5) + 1.0/(d(5)-i)
end do

total = sum(rsum)

write(*,*) 'The sum of series is = ', total

end program harmonic

needless to say, the code above assumes that the number given can be nicely divided by 5...besides that, though:

• it has a single loop
• the loop needs to be executed only n/5 times
• BUT there are 5 different operations in the loops which require the various variables to be moved in and out of the registry.

Then, there is:

program harmonic
implicit none
integer :: i,n,d(5)
real*8 :: s,rsum(5),total

!write(*,*)'How many terms you want to sum ?'
!read(*,*) n
n=10000000

d(1) = 1*n/5
d(2) = 2*n/5
d(3) = 3*n/5
d(4) = 4*n/5
d(5) = 5*n/5

rsum=0.0
do i= 0,n/5-1
    rsum(1) = rsum(1) + 1.0/(d(1)-i)
end do
do i= 0,n/5-1
    rsum(2) = rsum(2) + 1.0/(d(2)-i)
end do
do i= 0,n/5-1
    rsum(3) = rsum(3) + 1.0/(d(3)-i)
end do
do i= 0,n/5-1
    rsum(4) = rsum(4) + 1.0/(d(4)-i)
end do
do i= 0,n/5-1
    rsum(5) = rsum(5) + 1.0/(d(5)-i)
end do

total = sum(rsum)

write(*,*) 'The sum of series is = ', total

end program harmonic

where:

• you have several loops
• but each loop is only executed n/5 times
• and the quantities inside the loop do not have to keep getting moved in and out of the re
  gistry during the entire execution of the loop

Then, of course, you could start taking advantage of Fortran 90 features about vectorization and increased precision; the following code:

program harmonic
implicit none
integer :: i,n,counter
integer, parameter :: rkind = selected_real_kind(20)
real(kind=rkind) :: total, rsum(10000000)

n=10000000
forall (i=1:n) rsum(i) = 1.0_rkind/i
total = sum( rsum((0*n/5+1):(1*n/5)) ) + &
        sum( rsum((1*n/5+1):(2*n/5)) ) + &
        sum( rsum((2*n/5+1):(3*n/5)) ) + &
        sum( rsum((3*n/5+1):(4*n/5)) ) + &
        sum( rsum((4*n/5+1):(5*n/5)) )
write(*,*) 'The sum of series is = ', total

end program harmonic

for example, takes quite a while to run, but produces the following number:

16.69531136585985181539911893954098

...then, again, I don't quite know what I am doing...using your post as an excuse, I went ahead and tried to use both "forall" and "selected_real_kind" for my first time! Be warned.

 

[uart]

Hi gsal. Did you compare the accuracy of the results obtained by splitting the code into 5 separate sums, compared with just using just the one (summing terms from smallest to largest as suggested by AlphZero of course). In my own tests I couldn't notice any improvement by splitting the sum. Increasing the real precision of course makes a big difference.

For reference, to about 80 dp the results of f(10000000) is:
16.69531136585985181539911893954045188424986975237308046278513595435628869217425467

Hi Issacn.
If you're interested in both speed and precision then it is worth considering that the highest precision that is supported natively by your CPU is probably 80 bits (10 bytes).

The syntax seems to be compiler dependent, in ftn95 use "real (kind=3) :: sum=0" for example. And in gnu g95 use "real (kind = 10) :: sum=0" for example.

As far as I can tell it seems that ftn95 treats "kind" in a simple ordinal way, with kinds 1,2,3 representing the supported precisions in order from lowest to highest (single, double, extended). Whereas the gnu "g95" compiler appears to use kind to specify the number of bytes used.

 

[issacnewton]

Thank you gsal and uart for the explanation. I am following Chapman's "Fortran 95/2003
for scientists and engineers" third edition. And I am currently on the 4th chapter about loops.
So I wanted to apply that immediately to do the sum of this famous series. As a result of this
thread , I have learned that there are various issues involved here. So may be I need to wait till I finish this book. Author doesn't discuss "kind" till , I think , chapter 11. But all the posts in this thread were very informative. The world of high computing seems very exciting...

Thanks