sums9240.txt 21 Nov. 2014
Summing a Slowly Convergent Series
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for Math. 128B
1 May 2005
Prof. W. Kahan
Roundoff obscures badly the computed sum of a slowly convergent series
unless Compensated Summation is used. The obscuration is exposed by
Recomputation with Redirected Rounding, as examples will reveal.
Presumably the limit, as K --> +infinity, of the computable sum
Sum{k = 1 to K}( 3465/( (k+1/2)^2 - 1/16 ) + 3465/( k^2 - 1/16 ) )
+ 3465/(K+1/2) + 3465/(K+1)
is 9240 . That sum exceeds it by (1155/8)/( (K+1)(K+1/2)(K+3/4) ) .
A little mathematical analysis would confirm this, but programming a
computer is easier and running it is cheaper, so ...
Term(k) := 3465/( (k+0.5)^2 - 0.0625 ) + 3465/( k^2 - 0.0625 ) ;
Tail(k) := 3465/(k + 0.5) + 3465/(k + 1.0) .
Crude Sum:
Sum := 0 ; k := 0 ; Oldsum := -1 ;
While (Sum > Oldsum) do {
Oldsum := Sum ; k = 1+k ;
Sum := Sum + Term(k) ; } ;
Sum := Sum + Tail(k) .
Compensated Sum:
Sum := 0 ; k := 0 ; Oldsum := -1 ; COMP := 0 ;
While (Sum > Oldsum) do {
Oldsum := Sum ; k = 1+k ;
COMP := COMP + Term(k) ;
Sum := Sum + COMP ;
COMP := (Oldsum - Sum) + COMP ; } ;
Sum := Sum + ( COMP + Tail(k) ) .
=======================================================================
How accurate are the sums computed below from the algorithms above?
This question has an obvious answer because the true infinite sum is
known to be 9240 . But what if this were not known? Then an error-
analysis would have to take two sources of error into account.
One source is that only finitely many ( K ) terms have been added. A
final addend Tail(K) is intended to mostly compensate for the finite
sum's omitted Term(...)s . Tail(K) is roughly Term(k)'s integral
for k > K + 1/2 , and approximates the sum of omitted Term(...)s
ever better as K increases because Term(k) decreases so slowly to
zero. Note that Tail(K) is orders of magnitude bigger than the first
Term(K+1) omitted from the While-loop's Sum . This loop determines
K by running until Term(K) contributes negligibly to its Sum . If
the error in Tail(K) were taken into account even approximately, the
loop could be stopped sooner with a far smaller K , but all that is a
long story for another day.
The second source of error is roundoff. To render roundoff in Term(k)
inconsequential, the numerator 3465 has been chosen to turn the
first three Term(k)s into small integers computed exactly. Later
Term(k)s are so much smaller that roundoff in their last digits cannot
matter. The rounding errors that matter occur in " Sum + Term(k) " ;
these accumulate to corrupt almost the last log10(K) sig.dec. of the
computed Sum . Compensated Summation fends off that corruption at an
added cost negligible compared with the cost of computing Term(k) .
Surprisingly, Compensated Summation's While-loop stops substantially
sooner (with a substantially smaller K) than the Crude Sum's loop,
thus producing a better result faster. Can you see why?
Why should the foregoing error-analysis be believed? It could easily
be mistaken. It is so tedious that it might well be omitted, leaving
nothing to be (dis)believed. Perhaps some programs to test the Sum's
accuracy will be easier to construct than an error-analysis. Programs
that choose K differently can test whether Sum is as independent of
K as it should be when K is big enough. Rerunning any such programs
with redirected roundings reveals how badly roundoff corrupts Sum .
DANGER! Rerunning Crude Sum with rounding directed upward runs it
~~~~~~~ for a predictable intolerably long time. Can you see why?
Different formulas for Term(k) and Tail(K) produce different Sums:
-----------------------------------------------------------------------
Term(k) Tail(K) Ideal Sum
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~ ~~~~~~~~~
3465/( (k+0.5)^2 - 0.0625 ) 3465/(K + 0.5) 9240
+ 3465/( k^2 - 0.0625 ) + 3465/(K + 1.0)
15/(k + k^2) (15 - 7.5/(K+0.5)) 15
/(K+0.5)
1/((k+1)sqrt(k) + k sqrt(k+1)) 2/(sqrt(K+1.5) 1
+ sqrt(K+0.5))
18/k^2 18/(K + 0.5) 3 pi^2
= 29.608813203268075856503473
=======================================================================
These sums were computed below by an IBM Thinkpad T21 under Windows
2000 running Fortran 77 programs compiled by an Intel F86 compiler
using the Floating-Point Stack's 10-byte-wide registers controlling
both the precision and direction of rounding. One of the programs'
source-texts is appended at the end.
Computing 24 sig.bit sum 9240:
Rounding to nearest 24 sig. bits:
Crude sum = 9240.26855 Time = 0.05
Termcount k = 3768
Rounding down to 24 sig. bits:
Crude sum = 9238.80371 Time = 0.05
Termcount k = 2664
Rounding up to 24 sig. bits loops near-endlessly.
Rounding 24 sig. bits to zero:
Crude sum = 9238.80371 Time = 0.05
Termcount k = 2664
Rounding to nearest 24 sig. bits:
Compensated sum = 9240.00000 Time = 0.05
Termcount k = 2698
Rounding down to 24 sig. bits:
Compensated sum = 9239.99902 Time = 0.05
Termcount k = 2711
Rounding up to 24 sig. bits:
Compensated sum = 9240.00098 Time = 0.05
Termcount k = 2682
Rounding 24 sig. bits to zero:
Compensated sum = 9239.99902 Time = 0.05
Termcount k = 2711
**** PROGRAM PAUSE. 53 sig.bits will take about 23000 times as long as 24.
Computing 53 sig.bit sum 9240:
Rounding to nearest 53 sig. bits:
Crude sum = 9240.000011475229 Time = 25.32
Termcount k = 87290410
Rounding down to 53 sig. bits:
Crude sum = 9239.999948314162 Time = 17.85
Termcount k = 61723641
Rounding up to 53 sig. bits loops near-endlessly.
Rounding 53 sig. bits to zero:
Crude sum = 9239.999948314162 Time = 17.85
Termcount k = 61723641
Rounding to nearest 53 sig. bits:
Compensated sum = 9240.000000000000 Time = 19.44
Termcount k = 61728404
Rounding down to 53 sig. bits:
Compensated sum = 9239.999999999998 Time = 19.39
Termcount k = 61730077
Rounding up to 53 sig. bits:
Compensated sum = 9240.000000000002 Time = 19.44
Termcount k = 61725293
Rounding 53 sig. bits to zero:
Compensated sum = 9239.999999999998 Time = 19.39
Termcount k = 61730077
=======================================================================
Reruns with Redirected Rounding confirm that Crude Sum loses about
log10(K) sig.dec. while Compensated Sum loses at most 1 sig.dec.
Computing 24 sig.bit sum 15:
Rounding to nearest 24 sig. bits:
Crude sum = 15.0003862 Time = 0.05
Termcount k = 5609
Rounding down to 24 sig. bits:
Crude sum = 14.9982662 Time = 0.05
Termcount k = 3966
Rounding up to 24 sig. bits loops near-endlessly.
Rounding 24 sig. bits to zero:
Crude sum = 14.9982662 Time = 0.05
Termcount k = 3966
Rounding to nearest 24 sig. bits:
Compensated sum = 15.0000000 Time = 0.05
Termcount k = 4017
Rounding down to 24 sig. bits:
Compensated sum = 14.9999990 Time = 0.05
Termcount k = 3969
Rounding up to 24 sig. bits:
Compensated sum = 15.0000010 Time = 0.05
Termcount k = 3996
Rounding 24 sig. bits to zero:
Compensated sum = 14.9999990 Time = 0.05
Termcount k = 3969
**** PROGRAM PAUSE. 53 sig.bits will take about 23000 times as long as 24.
Computing 53 sig.bit sum 15:
Rounding to nearest 53 sig. bits:
Crude sum = 15.00000001668368 Time = 13.84
Termcount k = 129955756
Rounding down to 53 sig. bits:
Crude sum = 14.99999992485424 Time = 9.83
Termcount k = 91892597
Rounding up to 53 sig. bits loops near-endlessly.
Rounding 53 sig. bits to zero:
Crude sum = 14.99999992485424 Time = 9.78
Termcount k = 91892597
Rounding to nearest 53 sig. bits:
Compensated sum = 15.00000000000000 Time = 11.70
Termcount k = 91898489
Rounding down to 53 sig. bits:
Compensated sum = 15.00000000000000 Time = 11.59
Termcount k = 91901155
Rounding up to 53 sig. bits:
Compensated sum = 15.00000000000000 Time = 11.70
Termcount k = 91894008
Rounding 53 sig. bits to zero:
Compensated sum = 15.00000000000000 Time = 11.64
Termcount k = 91901155
=======================================================================
Reruns with Redirected Rounding show Crude Sum lost about log10(K)
sig.dec. while Compensated Sum lost at most 1 sig.dec. Apparent
coincidences of four Compensated Sums with " 15 " occur because the
compiler's unformatted " print *, x " displays too few digits to
distinguish 15 from its 8-byte wide floating-point neighbors.
Computing 24 sig.bit sum 1:
Rounding to nearest 24 sig. bits:
Crude sum = 1.00036776... Time = 0.05
Termcount k = 65536
Rounding down to 24 sig. bits:
Crude sum = 0.998898983... Time = 0.05
Termcount k = 41285
Rounding up to 24 sig. bits loops near-endlessly.
Rounding 24 sig. bits to zero:
Crude sum = 0.998898983... Time = 0.05
Termcount k = 41285
Rounding to nearest 24 sig. bits:
Compensated sum = 1.000000000000000 Time = 0.05
Termcount k = 41501
Rounding down to 24 sig. bits:
Compensated sum = 1.000000000000000 Time = 0.05
Termcount k = 41421
Rounding up to 24 sig. bits:
Compensated sum = 1.000000000000000 Time = 0.05
Termcount k = 41449
Rounding 24 sig. bits to zero:
Compensated sum = 1.000000000000000 Time = 0.05
Termcount k = 41421
**** PROGRAM PAUSE. 53 sig.bits will take about 660000 times as long as 24.
Computing 53 sig.bit sum 1:
Rounding to nearest 53 sig. bits:
Crude sum = 1.000000453219023 Time = 7547
Termcount k = 43290557639
Rounding down to 53 sig. bits:
Crude sum = 0.9999986448476669 Time = 5251
Termcount k = 27271342415
Rounding up to 53 sig. bits loops near-endlessly.
Rounding 53 sig. bits to zero:
Crude sum = 0. 9999986448476669 Time = 4898
Termcount k = 27271342415
Rounding to nearest 53 sig. bits:
Compensated sum = 1.000000000000000 Time = 5650
Termcount k = 27271483790
Rounding down to 53 sig. bits:
Compensated sum = 1.000000000000000 Time = 5406
Termcount k = 27271469629
Rounding up to 53 sig. bits:
Compensated sum = 1.000000000000000 Time = 5661
Termcount k = 27271486356
Rounding 53 sig. bits to zero:
Compensated sum = 1.000000000000000 Time = 5404
Termcount k = 27271469629
=======================================================================
Reruns with Redirected Rounding show Crude Sum loses about log10(K)
sig.dec. while Compensated Sum loses not even 1 sig.dec. Unusual
coincidences of all Compensated Sums with " 1 " occur because the
quotients' directed roundings partially offset the divisors'.
Computing 24 sig.bit sum 3pi2:
3*pi**2 = 29.608813203268075856503473
Rounding to nearest 24 sig. bits:
Crude sum = 29.6094017 Time = 0.05
Termcount k = 4345
Rounding down to 24 sig. bits:
Crude sum = 29.6061382 Time = 0.05
Termcount k = 3073
Rounding up to 24 sig. bits loops near-endlessly.
Rounding 24 sig. bits to zero:
Crude sum = 29.6061382 Time = 0.05
Termcount k = 3073
Rounding to nearest 24 sig. bits:
Compensated sum = 29.6088123 Time = 0.05
Termcount k = 3111
Rounding down to 24 sig. bits:
Compensated sum = 29.6088123 Time = 0.05
Termcount k = 3124
Rounding up to 24 sig. bits:
Compensated sum = 29.6088142 Time = 0.05
Termcount k = 3076
Rounding 24 sig. bits to zero:
Compensated sum = 29.6088123 Time = 0.05
Termcount k = 3124
**** PROGRAM PAUSE. 53 sig.bits will take about 23000 times as long as 24.
Computing 53 sig.bit sum 3pi2:
3*pi**2 = 29.608813203268075856503473
Rounding to nearest 53 sig. bits:
Crude sum = 29.60881322911488 Time = 10.16
Termcount k = 100663297
Rounding down to 53 sig. bits:
Crude sum = 29.60881308685216 Time = 7.20
Termcount k = 71179700
Rounding up to 53 sig. bits loops near-endlessly.
Rounding 53 sig. bits to zero:
Crude sum = 29.60881308685216 Time = 7.25
Termcount k = 71179700
Rounding to nearest 53 sig. bits:
Compensated sum = 29.60881320326808 Time = 9.12
Termcount k = 71182173
Rounding down to 53 sig. bits:
Compensated sum = 29.60881320326807 Time = 9.17
Termcount k = 71185856
Rounding up to 53 sig. bits:
Compensated sum = 29.60881320326808 Time = 9.12
Termcount k = 71186548
Rounding 53 sig. bits to zero:
Compensated sum = 29.60881320326807 Time = 9.12
Termcount k = 71185856
=======================================================================
Reruns with Redirected Rounding confirm that Crude Sum loses about
log10(K) sig.dec. while Compensated Sum loses at most 1 sig.dec.
c==== Fortran Program sums9240 : ======================================
$ notype
program sums9240
implicit none
integer*2 icw
c t1(x) = 3465.0/x
c tail(k) = t1(dble(k) + 0.5) + t1(dble(k+1))
c t2(x) = 3465.0/(x*x - 0.0625)
c term(k) = t2(dble(k)) + t2(dble(k) + 0.5)
c
print *, ' Computing 24 sig.bit sum 9240:'
c
c Set precision to 24 sig. bits rounded to nearest:
icw = #003Ch
call ldcw87(icw)
print *, ' Rounding to nearest 24 sig. bits:'
Call crudesum
c
c Set precision to 24 sig. bits rounded down:
icw = #043Ch
call ldcw87(icw)
print *, ' Rounding down to 24 sig. bits:'
Call crudesum
c
print *, ' Rounding up to 24 sig. bits loops near-endlessly.'
c
c Set precision to 24 sig. bits rounded to zero:
icw = #0C3Ch
call ldcw87(icw)
print *, ' Rounding 24 sig. bits to zero:'
Call crudesum
c
print *, ' '
c========================================================================
c Set precision to 24 sig. bits rounded to nearest:
icw = #003Ch
call ldcw87(icw)
print *, ' Rounding to nearest 24 sig. bits:'
Call compdsum
c
c Set precision to 24 sig. bits rounded down:
icw = #043Ch
call ldcw87(icw)
print *, ' Rounding down to 24 sig. bits:'
Call compdsum
c
c Set precision to 24 sig. bits rounded up:
icw = #083Ch
call ldcw87(icw)
print *, ' Rounding up to 24 sig. bits:'
Call compdsum
c
c Set precision to 24 sig. bits rounded to zero:
icw = #0C3Ch
call ldcw87(icw)
print *, ' Rounding 24 sig. bits to zero:'
Call compdsum
c
print *, ' '
c========================================================================
PAUSE ' 53 sig.bits will take about 23000 times as long as 24.'
c
print *, ' Computing 53 sig.bit sum 9240:'
c
c Set precision to 53 sig. bits rounded to nearest:
icw = #023Ch
call ldcw87(icw)
print *, ' Rounding to nearest 53 sig. bits:'
Call crudesum
c
c Set precision to 53 sig. bits rounded down:
icw = #063Ch
call ldcw87(icw)
print *, ' Rounding down to 53 sig. bits:'
Call crudesum
c
print *, ' Rounding up to 53 sig. bits loops near-endlessly.'
c
c Set precision to 53 sig. bits rounded to zero:
icw = #0E3Ch
call ldcw87(icw)
print *, ' Rounding 53 sig. bits to zero:'
Call crudesum
c
print *, ' '
c========================================================================
c
c Set precision to 53 sig. bits rounded to nearest:
icw = #023Ch
call ldcw87(icw)
print *, ' Rounding to nearest 53 sig. bits:'
Call compdsum
c
c Set precision to 53 sig. bits rounded down:
icw = #063Ch
call ldcw87(icw)
print *, ' Rounding down to 53 sig. bits:'
Call compdsum
c
c Set precision to 53 sig. bits rounded up:
icw = #0A3Ch
call ldcw87(icw)
print *, ' Rounding up to 53 sig. bits:'
Call compdsum
c
c Set precision to 53 sig. bits rounded to zero:
icw = #0E3Ch
call ldcw87(icw)
print *, ' Rounding 53 sig. bits to zero:'
Call compdsum
c
end
c========================================================================
c
Subroutine Crudesum
c Compute crude sum 9240 :
implicit none
real*8 oldsum, sum, term, tail, x, t1, t2
real*4 sec, seconds
integer*4 k
integer*2 icw
external seconds
t1(x) = 3465.0/x
tail(k) = t1(dble(k) + 0.5) + t1(dble(k+1))
t2(x) = 3465.0/(x*x - 0.0625)
term(k) = t2(dble(k)) + t2(dble(k) + 0.5)
sec = 0.0
sec = seconds(sec)
sum = 0.0
oldsum = -1.0
k = 0
do while (sum .GT. oldsum)
oldsum = sum
k = k+1
sum = sum + term(k)
end do
sum = sum + tail(k)
sec = seconds(sec)
c
c Reset print precision to 64 sig. bits:
icw = #033Ch
call ldcw87(icw)
c
print *, ' Crude sum = ', sum, ' Time =', sec
print *, ' Termcount k =', k
c
return
end
c========================================================================
Subroutine Compdsum
c Compute compensated sum 9240 :
implicit none
real*8 comp, oldsum, sum, term, tail, x, t1, t2
real*4 sec, seconds
integer*4 k
integer*2 icw
external seconds
t1(x) = 3465.0/x
tail(k) = t1(dble(k) + 0.5) + t1(dble(k+1))
t2(x) = 3465.0/(x*x - 0.0625)
term(k) = t2(dble(k)) + t2(dble(k) + 0.5)
sec = 0.0
sec = seconds(sec)
sum = 0.0
comp = 0.0
oldsum = -1.0
k = 0
c
do while (sum .GT. oldsum)
oldsum = sum
k = k+1
comp = term(k) + comp
sum = sum + comp
comp = (oldsum - sum) + comp
end do
comp = tail(k) + comp
sum = sum + comp
sec = seconds(sec)
c
c Reset precision to 64 sig. bits:
icw = #033Ch
call ldcw87(icw)
c
print *, ' Compensated sum =', sum, ' Time =', sec
print *, ' Termcount k =', k
c
return
end
c========================================================================
c real*4 function seconds( oldsec ) is in my Math. Library
c real*4 oldsec
c +------------------------------------------------------------------+
c | seconds(oldsec) = current time - oldsec in seconds, |
c | accurate to within 0.055 sec. |
c | DANGER: Peeks at the DOS clock-count in 046Ch . |
c +------------------------------------------------------------------+
c========================================================================