Summation Formulas

Consider the summation

On a computer, the direct evaluation

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

leads to the succesive computation of the partial sums , , , . Each step is corrupted by roundoff error

and these errors accumulate from one step to another. In our case, the effect of is present in the step

(we have ignored the terms containing ).

In general, the computed sum is

In general, the relative order can be as large as ; in single precision , so computing a long sum of terms may give a meaningless result.

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

the error made is ; we estimate this roundoff error by subtracting and then from the result

To correct the result, we need to subtract the error from the sum. Since is a small number, and the partial sum can be quite large, in general applying the correction directly to is useless: ! However, we can apply the correction to the next element (which can be much smaller). The next step reads

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

where each paranthesis contains exactly terms, and therefore evaluates to ; the exact result is .

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