-
Notifications
You must be signed in to change notification settings - Fork 192
Open
@Beliavsky
Description
Rosetta Code has a section on Kahan summation, which "is a method of summing a series of numbers represented in a limited precision in a way that minimises the loss of precision in the result." I have extracted and slightly modified the sumc function from the code there. Here is an example of how it can produce a more accurate result than the intrinsic SUM. I suggest that compensated summation be added to stdlib.
module m implicit none private public :: sumc,wp integer, parameter :: wp = kind(1.0d0) contains function sumc(a) result(asum) !add elements of the array, using limited precision. real(kind=wp) , intent(in) :: a(:) real(kind=wp) :: asum real(kind=wp) :: s,c,y,t ! assistants. integer :: i ! stepper integer :: n n = size(a) if (n < 1) then asum = 0.0_wp return end if s = a(1) ! start with the first element. c = 0.0_wp ! no previous omissions to carry forward. do i = 2,n ! step through the remainder of the array. y = a(i) - c ! combine the next value with the compensation. t = s + y ! augment the sum, temporarily in t. c = (t - s) - y ! catch what part of y didn't get added to t. s = t ! place the sum. end do ! on to the next element. asum = s ! c will no longer be zero. end function sumc end module m program main ! 04/26/2021 11:39 AM driver for sumc use m, only: sumc,wp implicit none integer, parameter :: n = 100000000 real(kind=wp) :: x(n) real(kind=wp), parameter :: xadd = 10.0_wp**15 call random_seed() call random_number(x) x = x + xadd print*,([sumc(x),sum(x)] - n*xadd)/n end program main
output:
0.50331647999999996 2207252.1444556802
If you use wp = kind(1.0) (single precision) in the code above, you can see the difference even with xadd = 0.0. The output of
module m implicit none private public :: sumc,wp integer, parameter :: wp = kind(1.0) contains function sumc(a) result(asum) !add elements of the array, using limited precision. ! adapted from Rosetta Code https://rosettacode.org/wiki/Kahan_summation real(kind=wp) , intent(in) :: a(:) real(kind=wp) :: asum real(kind=wp) :: s,c,y,t ! assistants. integer :: i ! stepper integer :: n n = size(a) if (n < 1) then asum = 0.0_wp return end if s = a(1) ! start with the first element. c = 0.0_wp ! no previous omissions to carry forward. do i = 2,n ! step through the remainder of the array. y = a(i) - c ! combine the next value with the compensation. t = s + y ! augment the sum, temporarily in t. c = (t - s) - y ! catch what part of y didn't get added to t. s = t ! place the sum. end do ! on to the next element. asum = s ! c will no longer be zero. end function sumc end module m ! program main ! 04/26/2021 11:39 AM driver for sumc use m, only: sumc,wp implicit none integer, parameter :: n = 100000000 real(kind=wp) :: x(n) call random_seed() call random_number(x) print*,[sumc(x),sum(x)]/n end program main
is
0.499940693 0.167772159