# Re: RE: st: Another precision problem with Mata

 From wgould@stata.com (William Gould, Stata) To statalist@hsphsun2.harvard.edu Subject Re: RE: st: Another precision problem with Mata Date Tue, 23 Aug 2005 09:55:31 -0500

```Ben Jann <ben.jann@soz.gess.ethz.ch> asked about the precision of abs((x-y)/z)
and in particular, abs((x-y)/z)<1 when x=y+z.  Actually, yesterday he asked
the question when y=x+z, but one way or the other, it does not make much
difference.

I replied, rather negatively and certainly in purist mode.  I said that there
was no reason to think about how one might code abs((x-y/z)<1 so that it is
guaranteed to be false in all the right circumstances because, for all
possible values, there is no way to do that.

I could have been more helpful, however, by eliminating cases when
calculating y PLUS z results in y or results in z because the other is too
small.

I have now thought about that problem, and I believe the solution to
Ben's problem is to code

if ( abs(abs(x-y)-z) < epsilon(x) ) {
...
}

The above is not guaranteed to work, but it will work in lots of
reasonable cases.

Below my signature is a do-file that Ben should run.  It provides a
demonstration of my claim and it shows exactly what the problem is
and where the precision is being lost.

-- Bill
wgould@stata.com

------------------------------------------------------------------------------
clear

/*
We will explore

abs((x-y)/z)<1

when x = y+z.  In that case, the expression becomes

abs(((y+z)-y)/z) < 1)

The correct answer is 0 because in infinite precision

abs(((y+z)-y)/z) < 1)
abs(z/z) < 1)
abs(1 < 1)

Below, we will look at two expressions:

x - y    (which is (y+z)-y in our special case)

abs((x-y)/z)
*/

mata:
real scalar f(y, z)
{
real scalar        x

x = y + z
return(abs((x-y)/z) < 1)
}

f(0, .1)                // no problem
f(1, .1)
f(1000, .1)             // no problem
f(1e+20, .1)            // problem

real scalar findprob(y, z, mult)
{
while (f(y,z)==0) y = mult*y
return(y)
}

findprob(1, .1, 2)
f(4, .1)                        // problem at 4

void fdetail(y, z)
{
real scalar        x

printf("         y = %21x\n", y)
printf("         z = %21x\n", z)
printf("-----------------------------------\n")
printf("   x = y+z = %21x\n", x = y+z)
printf("     x - y = %21x     (%s)\n", x-y,
(x-y)==z ? "okay, x-y == z" : "problem, x-y != z")
printf("\n")
printf("error:\n")
printf("         (x-y)-z = %g\n", (x-y)-z)
printf("   |(x-y)/z| - 1 = %g\n", abs((x-y)/z)-1)
printf("    note, e(y)=%g, e(x)=%g e(z)=%g e(1)=%g\n",
epsilon(y), epsilon(x), epsilon(z), epsilon(1))
}
fdetail(4, .1)
fdetail(16384, .1)
fdetail(32768, .1)
fdetail(65536, .1)
fdetail(131072, .1)
fdetail(262144, .1)

fdetail(1e+10, .1)
fdetail(1e+20, .1)                // complete loss of precision

/* ----- the point of all the above is that
if y>z, then (x-y) has error bounded by epsilon(x).
assuming numbers are such that we have any precision at all.
*/

/* ----------------------------------- */
// now consider y<z:

f(.1, 1)
f(.1, 2)
f(.1, 4)

fdetail(.1, 4)
fdetail(.1, 512)
fdetail(.1, 1024)
fdetail(.1, 8192)
fdetail(.1, 16384)
fdetail(.1, 32768)
fdetail(.1, 262144)
end
------------------------------------------------------------------------------
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```