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========================================================================