eps function and numerical accuracy

[Daniel J Sebald]

[Note, I'm going to write a follow-up thread about "eps", numerical 
accuracy.]

Following up on the issues raised by the binocdf() discussion.  I'm 
wondering if the following documentation for the routine eps() is 
unclear or misleading:

 -- Built-in Function:  eps
 -- Built-in Function:  eps (X)
 -- Built-in Function:  eps (N, M)
 -- Built-in Function:  eps (N, M, K, ...)
 -- Built-in Function:  eps (..., CLASS)
     Return a scalar, matrix or N-dimensional array whose elements are
     all eps, the machine precision.  More precisely, `eps' is the
     relative spacing between any two adjacent numbers in the machine's
     floating point system.  This number is obviously system dependent.
     On machines that support IEEE floating point arithmetic, `eps' is
     approximately 2.2204e-16 for double precision and 1.1921e-07 for
     single precision.

There is a more to the documentation, but it is this first part that 
sort of imparts the meaning.  I suppose if one thinks about it long 
enough, "relative spacing between any two adjacent numbers" is kind of 
correct.  But what it doesn't impart clearly is the fact this reflects 
the nature of the mantissa.  The relative part sort of pertains to the 
exponent.  In the latter half of the above paragraph, it should probably 
state

   `eps (1.0)' is approximately 2.2204e-16...

This binocdf() discussion has got me to wondering why there is so little 
use of exp (X) in the test functions of script files for Octave.  One 
would think there should be.  For example, if I were to test

p = 1-1e-3;
pval = binopdf ([0:50]', 50, p);
assert (cumsum (pval), binocdf ([0:50]', 50, p), 2*eps)

it really isn't saying much at all:

> > assert (cumsum (pval), binocdf ([0:50]', 50, p), 2*eps)

What we really want is

assert (cumsum (pval), binocdf ([0:50]', 50, p), 2*eps(pval)/eps(1.0))

which at least sort of tracks the magnitude of the final result. 
However, I'm not sure if the assert() routine acts properly on 
vector/matrix comparisons element by element.  This test isn't 
indicating an error for the binocdf() routine prior to Rik's change, 
when I would think it should.

> > assert (cumsum (pval), binocdf ([0:50]', 50, p), 2*eps(pval)/eps(1.0))

(Again, I'm not 100% what assert() with a vector tolerance is attempting 
to do.)  I would think it should be more like the following [is there a 
CR-LF missing somewhere in this output?]:

> > assert (cumsum (pval), binocdf ([0:50]', 50, p), 2*eps(pval))error: assert 
> > (cumsum (pval),binocdf ([0:50]', 50, p),2 * eps (pval)) expected
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000001
   0.000000000000096
   0.000000000015303
   0.000000002040786
   0.000000221981888
   0.000018921655140
   0.001186482503818
   0.048794371802969
   1.000000000000000
but got
   1.00000000000007e-150
   4.99510000000035e-146
   1.22260117600006e-141
   1.95424813815765e-137
   2.29399723360430e-133
   2.10841676614637e-129
   1.57977022596915e-125
   9.92031009822021e-122
   5.32697826475548e-118
   2.48350747999872e-114
   1.01724999177793e-110
   3.69550684776946e-107
   1.19987796084329e-103
   3.50394897537906e-100
    9.25151215525190e-97
    2.21822759875234e-93
    4.84771698248160e-90
    9.68615422923810e-87
    1.77409990314262e-83
    2.98511415284120e-80
    4.62253801912459e-77
    6.59738446954312e-74
    8.68836583529609e-71
    1.05672286748147e-67
    1.18770467184289e-64
    1.23406745773931e-61
    1.18550990718313e-58
    1.05282245558334e-55
    8.64033611769894e-53
    6.54884440079011e-50
    4.58011354682620e-47
    2.95231633143266e-44
    1.75142047732177e-41
    9.54507616566026e-39
    4.76856227562884e-36
    2.17814373131045e-33
    9.06845557842062e-31
    3.42871128091703e-28
    1.17213664214074e-25
    3.60414974466429e-23
    9.90534035957436e-21
    2.41464425822922e-18
    5.17200539747901e-16
    9.61955465107593e-14
    1.53025290323127e-11
    2.04078614934109e-09
    2.21981887963310e-07
    1.89216551400609e-05
    1.18648250381793e-03
    4.87943718029691e-02
    1.00000000000000e+00
maximum absolute error 4.44089e-16 exceeds tolerance 2.71333e-166maximum 
absolute error 1.77821e-161 exceeds tolerance 2.91341e-157maximum 
absolute error 4.77334e-153 exceeds tolerance 7.82064e-149maximum 
absolute error 6.40667e-145 exceeds tolerance 5.24834e-141maximum 
absolute error 4.29944e-137 exceeds tolerance 1.76105e-133maximum 
absolute error 7.21326e-130 exceeds tolerance 2.95455e-126maximum 
absolute error 1.21018e-122 exceeds tolerance 4.95692e-119maximum 
absolute error 1.01518e-115exceeds tolerance 2.07908e-112maximum 
absolute error 8.51592e-109 exceeds tolerance 1.74406e-105maximum 
absolute error 3.57184e-102 exceeds tolerance 7.31512e-99maximum 
absolute error 7.49068e-96 exceeds tolerance 1.53409e-92maximum absolute 
error 1.57091e-89 exceeds tolerance 3.21722e-86maximum absolute error 
3.29444e-83 exceeds tolerance 3.3735e-80maximum absolute error 
3.45447e-77 exceeds tolerance 3.53737e-74maximum absolute error 
3.62227e-71 exceeds tolerance 3.70921e-68maximum absoluteerror 
1.89911e-65 exceeds tolerance 1.94469e-62maximum absolute error 
9.95682e-60 exceeds tolerance 5.09789e-57maximum absolute error 
2.61012e-54 exceeds tolerance 1.33638e-51maximum absolute error 
6.84228e-49 exceeds tolerance 3.50325e-46maximum absolute error 
8.96831e-44 exceeds tolerance 4.59177e-41maximum absolute error 
1.17549e-38 exceeds tolerance 3.00927e-36maximum absolute error 
7.70372e-34 exceeds tolerance 1.97215e-31maximum absolute error 
2.52435e-29 exceeds tolerance 6.46235e-27maximum absolute error 
8.27181e-25 exceeds tolerance 5.29396e-23maximum absolute error 
6.77626e-21 exceeds tolerance 4.33681e-19maximum absolute error 
1.38778e-17 exceeds tolerance 2.22045e-16
error: called from:
error: 
/usr/local/src/octave/octave/build-gui-11/../octave/scripts/testfun/assert.m 
at line 235, column 5
> >

Clearly, the numerical precision is being lost, and we would like to 
reflect that.

So, a couple questions:

1) Can the assert() routine be made to operate on an element-by-element 
basis?

2) How many tests in the script files are inconsequential because of 
"eps" being used rather than something like "eps(smallnumber)/eps(1.0)", 
i.e., adjusting to the scale of the expected result?

3) We probably don't want to have "eps(X)/eps(1.0)" everywhere.  Just an 
"eps(X)" would be more convenient.  Maybe I'm not thinking of eps() 
correctly in terms of its scale.

Dan