Finger exercises for Monday Sept 16
-----------------------------------

[see pp. 33-35 of Monahan]

0. Think about #2.7.  This was already a finger exercise for Fri Sept 13.

1. Think about #2.9.  Give an example, in base 10 at least.

2. Series evaluation on the computer is interesting.  Think about #'s
   2.10, 2.11, 2.12; and try to implement some of them in S.

3. Explore #2.24 using S.

======================================================================

0. Think about #2.7.  This was already a finger exercise for Fri Sept 13.

Bryce asked how one could predict the size of machine epsilon, 
 
   min (eps: 1+eps <> 1)

We found machine eps by trial and error to be somewhere around
2.220446e-16.

Recall (from the book, p. 20) that the IEEE double-precision floating
pt representation is

              (S_1, E_{11}, 1.F_{52})_2  [bias = 2^{11} - 1 = 1023 ]
   
which represents the number
   
   	       1^S \times (1.F) \times 2^{E-1023}

Therefore, 

               fl(1.00) = (0,1023, 1.000...00)

with 52 0's in F.  The smallest thing we can add to 1 and get
something different will be the smallest thing that does not underflow
when we renormalize it so that E=1023.  Consider

	       0.00000...00001 = (0,1023-52, 1.00000...00000) = 2^{-52}

(with 52 0's in F again).
  
This number must be un-normalized before it can be added to fl(1.0),
by adding 52 to E and using a leading 0 rather than a fully-normalized
leading 1, i.e.:

               fl*(2^{-52}) = (0, 1023, 0.0000...00001)

with 51 0's and a 1 to the right of the binary pt.  Adding this to 1
produces

               fl(1.0) + fl*(2^{-52}) = (0, 1023, 1.0000...0001)

which is greater than 1.  Clearly if $\epsilon<2^{-52}$ we would
expect it to "underflow" when we un-normalize for adding to 1:

               fl*(\epsilon) = (0, 1023, 0.00...0000)

(since the first significant digit of \epsilon would "fall off" the
right side of F).

Thus, in IEEE double precision, machine epsilon should be approx.
2^{-52} = 2.220446e-16.

Further exercise: Monahan says (p. 28) that machine epsilon should be
close to the machine unit U, which bounds the relative error in
translating from reals to floats:

               |fl(x) - x|
	       -----------  <= U
	           |x|

why is this?

=======================================================================

1. Think about #2.9.  Give an example, in base 10 at least.

   "Can the average of two floats (x+y)/2 be smaller than either one?"

The idea is to "lose" low-order significant digits off the right end
of F in the floating point calculation of x+y.  For example consider
the floating point representation

      (S, E_1, F_2)_{10}, bias = 5

      fl(5.1x10^{-5}) = (0, 0, 5.1)

If we add this to itself we get 10.2 x 10^{-5}, and

      fl(10.2x10^{-5}) = fl(1.02x10^{-4}) = (0,1,1.0) = fl(1.0 x 10^{-4})

(we lost the "2" since only two digits are allocated to the fraction),
and now dividing by 2 we get 0.5 x 10^{-4} = 5.0 x 10^{-5}, which is
less than what we started with.

[This doesn't work in base 2 since adding a number to itself is just a
left shift, so we don't generate an extraneous low-order digit that
can be lost off the right end of F in the sum.]

=======================================================================

2. Series evaluation on the computer is interesting.  Think about #'s
   2.10, 2.11, 2.12; and try to implement some of them in S.

2.10:
I tried the following functions:

e.1 _ function(m.max=1000) {
#
# e.1
#
# compute "e" by summing the series from largest to smallest.
#
   m _ 0
   old _ 0 
   eps _ 1  # (1/0! = 1)
   while ((m <= m.max) && ( old+eps!=old)) {
      old _ old + eps
      m _ m + 1
      oldeps _ eps
      eps _ eps/m
   }
   return(old,oldeps,m-1)
}

e.2 _ function(m.max=1000) {

# 
# e.2
#
# compute e from smallest to largest; start with a sufficiently small
# term
#

   m _ 0
   eps _ 1
   while ((eps > 0) && (m <= m.max)) {
     m _ m + 1
     oldeps _ eps 
     eps _ eps/m
   }

   m _ m - 1

   eps _ oldeps  # last nonzero epsilon above 
   old _ 0
   for (i in (m + 0 - (0:m))) {
      old _ old + eps   
      eps _ eps * i 
#      cat(i,old,eps,"\n")
#      scan()  
      } 
   return(old,oldeps,m)
}

> unlist(e.1())
      old  m
 2.718282 18
> unlist(e.2())
      old        oldeps   m
 2.822645 2.964394e-323 177
> unlist(e.2(17))
      old       oldeps  m
 2.718282 2.811457e-15 17
> unlist(e.2(50))
      old       oldeps  m
 2.718282 3.287949e-65 50
> unlist(e.2(120))
      old        oldeps   m
 2.718282 1.494879e-199 120
> unlist(e.2(150))
      old        oldeps   m
 2.718282 1.750276e-263 150
> unlist(e.2(200))
      old        oldeps   m
 2.822645 2.964394e-323 177
> e.1()$old - e.2()$old
[1] -0.1043633
> e.1()$old - e.2(150)$old
[1] -2.220446e-15

Note that the "estimate" of e becomes less accurate (way too large) as
we go "too far out in the series".  Mathematically, we know this
series converges to exactly e.  What is going wrong on the computer?



2.11
I tried the following functions:

e.1 _ function(x, m.max=1000) {
#
# e.1
#
# compute exp(x) using an alternating series
#
   sgn _ 1
   m _ 0
   old _ 0 
   eps _ 1  # (1/0! = 1)
   while ((m <= m.max) && ( old + (sgn*x^m * eps) != old)) {
      old _ old + sgn*x^m * eps
      m _ m + 1
      oldeps _ eps
      eps _ eps/m
      sgn _ - sgn
   }
   return(old,oldeps,m-1)
#   return(old)


}


e.2 _ function (x, m.max=1000) {
#
# e.2 
#
# compute exp(x) using a reciprocal series...
#
   m _ 0
   old _ 0 
   eps _ 1  # (1/0! = 1)
   while ((m <= m.max) && ( old + (x^m * eps) != old)) {
      old _ old + x^m * eps
      m _ m + 1
      oldeps _ eps
      eps _ eps/m
   }
   return(1/old,oldeps,m-1)
#   return(1/old)
}

I compared these functions as follows:

> for (k in 1:4) {
+ print (rbind(r1=e.1(k),r2=e.2(k)))
+ }
                   oldeps
r1 0.3678794 1.561921e-16 18
r2 0.3678794 2.811457e-15 17
                   oldeps
r1 0.1353353 1.611738e-24 24
r2 0.1353353 8.896791e-22 22
                    oldeps
r1 0.04978707 1.130996e-31 29
r2 0.04978707 2.479596e-27 26
                    oldeps
r1 0.01831564 1.151634e-37 33
r2 0.01831564 3.769988e-33 30

It appears that the reciprocal series method (row r2) stops sooner, by
one or two iterations.

The "oldeps" columns give the absolute error of the approximation.
The number listed for r2 is a LOWER BOUND on the error, the number
listed for r1 is an UPPER BOUND.  Hence for a couple more iterations,
method R1 seems to be getting one to four more places of accuracy in
the answer.



2.12
I tried the following function:

h.1 _ function() {
#
# h.1
#
# sum the harmonic series out to underflow for the terms.
# compare the term number for underflow with an estimate
# based on machine epsilon
#
   j _ 1
   old _ 0
#
   while (old + 1/j >  old) {

     old _ old + 1/j
     j _ j + 1
     cat(j,old, old+1/j - old,"\n")

   }
   j.empirical _ j-1
   return(j.empirical)

}

h.2 _ function() {
#
# h.2
#
# sum the harmonic series out to underflow for the terms.
# compare the term number for underflow with an estimate
# based on machine epsilon
#
   eps _ .Machine$double.eps

   n.old _ 99
   n _ 100
   while (abs(n-n.old) > 100*eps*n.old) {
     Happrox  _ .5773 + log(n)
     n.old _ n
     n _ 1/(eps*Happrox)
     # solves 1/H approx = 1/(.5773 + log(n)) = n*eps, for n
     cat(n,n.old,abs(n - n.old)/abs(n.old),"\n")
   }
   j.theoretical _ n

   return(j.theoretcal)
} 

Function h.1() seemed to take forever.  So I ran function h.2()
instead.

That yielded

> h.2()

869006374517490 100 8690063745173.9
128763778488236 869006374517490 0.851826428132094
136199076900285 128763778488236 0.0577437109981099
135968237840146 136199076900285 0.00169486508566957
135975201552675 135968237840146 5.12157297821194e-05
135974991295308 135975201552675 1.5462920017243e-06
135974997643502 135974991295308 4.66864800525202e-08
135974997451834 135974997643502 1.40958220865364e-09
135974997457621 135974997451834 4.25588351420996e-11
135974997457446 135974997457621 1.28493291610066e-12
135974997457451 135974997457446 3.87249501633409e-14
135974997457451 135974997457451 1.14910831345221e-15

This explains why (!) h.1() was taking so long!




=======================================================================

3. Explore #2.24 using S.

Which is better to calculate, exp(u-v) or exp(u)/exp(v) ?

I tried

test _ function(a=1000,delta=2e-16) {

# 
# test -  
#
# Which is better to calculate, exp(u-v) or exp(u)/exp(v) ?
#

   results _ NULL
   for (v in -a:a) {
     u _ v + delta
     results _ rbind(results, c(subtr= exp(u-v),divis =
                   (exp(u))/(exp(v))))
     }
   print(results)

   invisible(results)

}

If I do

> bozo _ test(1000)
> sink("bozo.txt")
> bozo
> sink()

then I can edit the file "bozo.txt" and see that the "division" method
breaks down from 0/0 or infinity/infinity forms at around -709 and
+711.  These numbers are larger than the ones suggested in the text,
because Monahan is probably working with single precision and we are
working with double.

[applying the "as.single" function to the exponetials, the breakdown
occurs more like at +/-90 in single precision, consistent with the
Monahan's suggestion].


=======================================================================