`Consider the summation
`

`On a computer, the direct evaluation
`

s = 0.0 do i=1,n s = s + x(i) end do

`In general, the computed sum is
`

`To minimize the accumulation of roundoff, we consider 2 strategies. First, we
try to minimize at each step. For this, the partial sum
and the element must be of comparable size; therefore,
sorting the terms in increasing order before summation may help
to reduce the individual roundoffs.
`

`The second strategy is more elaborate, and is called the Kahan summation
algorithm. The idea is to account for the errors at each step and to
add the proper corrections as we go along.
At step we have
`

`Below, we give an example of a program which evaluates the sum
`

`Note the following characteristics of the implementation.
They give a flavor of how to implement large programs.
`

- each major function (terms initialization, computation of the sum using sorted up, sorted down and Kahan formula, and printing of the results) is handled by a different procedure. Each procedure is written and tested independently.
- The main program just calls the individual procedures in the appropriate order; it is therefore easy to keep track of the program flow by simply reading the main.
- Global data (summation terms, and sum values) are defined in a module, which is used in all the subroutines.

! Long summations example ! different algorithms for performing summation module terms save ! Numeration basis integer, parameter :: base=10 ! sum = 1 + base*(.) + ... + base**p_max*(.) integer, parameter :: p_max=7 ! Number of terms in the sum integer, parameter :: n_max=(base**(p_max+1)-1)/(base-1) real, dimension(n_max) :: x real :: exact, kahan, sortu, sortd end module terms subroutine init use terms implicit none integer :: i, j, k=0 ! do i = 0,p_max do j = 1,base**i x(k+j) = 1.0/base**(i) end do k = k+base**i end do ! end subroutine init subroutine sum_exact ! ! the exact value of the sum ! use terms implicit none exact = p_max+1 end subroutine sum_exact subroutine sum_sortd ! ! Performs normal summation, ! with the terms sorted in decreasing order ! use terms implicit none integer :: i sortd = 0.0 do i = 1,n_max sortd = sortd + x(i) end do end subroutine sum_sortd subroutine sum_sortu ! ! Performs normal summation, ! with the terms sorted in increasing order ! use terms implicit none integer :: i sortu = 0.0 do i = n_max,1,-1 sortu = sortu + x(i) end do end subroutine sum_sortu subroutine sum_kahan ! ! Performs the summation using Kahan's algorithm ! use terms implicit none integer :: i real :: c, y, s kahan = x(1) c = 0.0 do i = 2,n_max y = x(i) - c s = kahan + y c = (s-kahan) - y kahan = s end do end subroutine sum_kahan subroutine print_results ! ! Prints the setting, the results and their accuracies ! use terms implicit none print *, "basis = ",base,", maximal power = ",p_max print *, "no. of terms in the sum = ", n_max print *, "sort down: ",sortd," error = ",(sortd-exact)/exact print *, "sort up: ",sortu," error = ",(sortu-exact)/exact print *, "kahan: ",kahan," error = ",(kahan-exact)/exact end subroutine print_results program sum use terms implicit none ! call init call sum_exact call sum_sortd call sum_sortu call sum_kahan call print_results ! end program sum

`The results of this program show clearly the
advantage of Kahan summation.
`

! BASIS = 10, Maximal Power = 7 ! No. of terms in the Sum = 11111111 ! SORT DOWN: 6.95631695 Error = -0.130460382 ! SORT UP: 8.01876831 Error = 2.346038818E-3 ! KAHAN: 8. Error = 0.E+0