Consider the summation
On a computer, the direct evaluation
s = 0.0 do i=1,n s = s + x(i) end doleads to the succesive computation of the partial sums , , , . Each step is corrupted by roundoff error
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.
! 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