[boost] Re: test_fp_comparisons and rounding errors

From: Gennadiy Rozental (gennadiy.rozental_at_[hidden])
Date: 2003-07-05 14:47:02

"Beman Dawes" <bdawes_at_[hidden]> wrote in message
news:4.3.2.7.2.20030705093032.01ceb800_at_mailhost.esva.net...
> test_fp_comparisons has been failing for a long time. The issue has to do
> with how much tolerance to allow for rounding errors.
>
> The close_at_tolerance algorithm calculates tolerance as follows:
>
>     n*std::numeric_limits<T>::epsilon()/2   (1)
>
> where n is the number of possible floating-point rounding errors.
>
> The particular test case that is failing calls the comparison algorithm
> with a rounding error argument of 4. That allows for up to 2 epsilons
> tolerance.

Thanks, Beman for bringing this up. I actually was playing with this test
lately. It does indeed fail on following comparison:

I compare 0.11 with 0.11.

Left size calculated as: 11./100 (1 rounding error).
Right side calculated as: tmp = 11/10 (1 rounding error ); tmp*tmp - tmp (2
rounding errors).

In sum have 1+1+2 = 4 rounding errors. According to my understanding
relative error of calculation should not exceed tolerance calculated with
above formula 1. This statement does not hold.

What I've tried to do is to analyze the "relative error"/(epsilon/2) for
different compilers and FPTs. Here is the code:

----------------------------------------------------------------------------
------------

template<typename FPT>
void foo( FPT* = 0 )
{
    FPT tmp  = 11;
    tmp /= 10;
    FPT tmp1 = tmp * tmp - tmp;

    FPT tmp2 = 11./100;

    FPT diff = fpt_abs( tmp1 - tmp2 );

    FPT d1 = diff / tmp1;
    FPT d2 = diff / tmp2;

    FPT r1 = d1 / std::numeric_limits<FPT>::epsilon() * 2;
    FPT r2 = d2 / std::numeric_limits<FPT>::epsilon() * 2;

    std::cout << "diff= "    << diff << std::endl;
    std::cout << "epsilon= " << std::numeric_limits<FPT>::epsilon() <<
std::endl;
    std::cout << "r1= " << r1 << std::endl;
    std::cout << "r2= " << r2 << std::endl;
}

int
main()
{
    foo<float>();
    foo<double>();
    foo<long double>();

    return 0;
}

----------------------------------------------------------------------------
------------

Here is the result produced for different compilers:

Borland command line:
-----------------------
diff= 2.98023e-08
epsilon= 1.19209e-07
r1= 4.54545
r2= 4.54545
diff= 1.11022e-16
epsilon= 2.22045e-16
r1= 9.09091
r2= 9.09091
diff= 5.42101e-19
epsilon= 1.0842e-19
r1= 90.9091
r2= 90.9091

Metrowerks
-------------
diff= 2.98023e-08
epsilon= 1.19209e-07
r1= 4.54545
r2= 4.54545
diff= 9.71445e-17
epsilon= 2.22045e-16
r1= 7.95455
r2= 7.95455
diff= 9.71445e-17
epsilon= 2.22045e-16
r1= 7.95455
r2= 7.95455

GCC 3.2.3
-----------
diff= 2.98023e-08
epsilon= 1.19209e-07
r1= 4.54545
r2= 4.54545
diff= 1.11022e-16
epsilon= 2.22045e-16
r1= 9.09091
r2= 9.09091
diff= 5.42101e-19
epsilon= 1.0842e-19
r1= 90.9091
r2= 90.9091

MSVC 6.5
------------
diff= 0
epsilon= 1.19209e-007
r1= 0
r2= 0
diff= 9.71445e-017
epsilon= 2.22045e-016
r1= 7.95455
r2= 7.95455
diff= 9.71445e-017
epsilon= 2.22045e-016
r1= 7.95455
r2= 7.95455

MSVC6.5 + STLport
----------------------
diff= 0
epsilon= 1.19209e-007
r1= 0
r2= 0
diff= 9.71445e-017
epsilon= 2.22045e-016
r1= 7.95455
r2= 7.95455
diff= 9.71445e-017
epsilon= 2.22045e-016
r1= 7.95455
r2= 7.95455

MSVC7.1
----------
diff= 2.98023e-008
epsilon= 1.19209e-007
r1= 4.54545
r2= 4.54545
diff= 9.71445e-017
epsilon= 2.22045e-016
r1= 7.95455
r2= 7.95455
diff= 9.71445e-017
epsilon= 2.22045e-016
r1= 7.95455
r2= 7.95455

Mingw
-------
diff= 2.98023e-08
epsilon= 1.19209e-07
r1= 4.54545
r2= 4.54545
diff= 9.71445e-17
epsilon= 1.11022e-16
r1= 15.9091
r2= 15.9091
diff= 9.71445e-17
epsilon= 1.11022e-16
r1= 15.9091
r2= 15.9091

----------------------------------------------------------------------------
------------

As you can see It's never less then 4.

One idea, that I had, is maybe "Testing tool" introduce an extra error
during relative error calculation. It take an extra 6(?) operations before
comparison is performed. It still does not explain all the results above and
I am not sure it's correct.

I attached the test program for anybody interested.

Regards,

Gennadiy.

begin 666 fpt.cpp
M(VEN8VQU9&4@/&)O;W-T+VQI;6ET<RYH<' ^#0HC:6YC;'5D92 \:6]S=')E
M86T^#0H-"G1E;7!L871E/'1Y<&5N86UE($905#X-"FEN;&EN92!&4%0-"F9P
M=%]A8G,H($905"!A<F<@*2 -"GL-"B @("!R971U<FX_at_87)G(#P@," _("UA
M<F<@.B!A<F<[#0I]#0H-"G1E;7!L871E/'1Y<&5N86UE($905#X-"G9O:60@
M9F]O*"!&4%0J(#T@," I#0I[#0H@(" @1E!4('1M<" @/2 Q,3L-"B @("!T
M;7 @+ST@,3 [#0H@(" @1E!4('1M<#$@/2!T;7 @*B!T;7 @+2!T;7 [#0H-
M"B @("!&4%0@=&UP,B ](#$Q+B\Q,# [#0H-"B @("!&4%0_at_9&EF9B ](&9P
M=%]A8G,H('1M<#$@+2!T;7 R("D[#0H-"B @("!&4%0_at_9#$@/2!D:69F("\@
M=&UP,3L-"B @("!&4%0_at_9#(@/2!D:69F("\@=&UP,CL-"@T*(" @($905"!R
M,2 ](&0Q("\@<W1D.CIN=6UE<FEC7VQI;6ET<SQ&4%0^.CIE<'-I;&]N*"D@
M*B R.PT*(" @($905"!R,B ](&0R("\@<W1D.CIN=6UE<FEC7VQI;6ET<SQ&
M4%0^.CIE<'-I;&]N*"D@*B R.PT*#0H@(" @<W1D.CIC;W5T(#P\(")D:69F
M/2 B(" @(#P\(&1I9F8@/#P@<W1D.CIE;F1L.PT*(" @('-T9#HZ8V]U=" \
M/" B97!S:6QO;CT@(B \/"!S=&0Z.FYU;65R:6-?;&EM:71S/$905#XZ.F5P
M<VEL;VXH*2 \/"!S=&0Z.F5N9&P[#0H@(" @<W1D.CIC;W5T(#P\(")R,3T@
M(B \/"!R,2 \/"!S=&0Z.F5N9&P[#0H@(" @<W1D.CIC;W5T(#P\(")R,CT@
M(B \/"!R,B \/"!S=&0Z.F5N9&P[#0I]#0H-"FEN= T*;6%I;B_at_I#0I[#0H@
M(" @9F]O/&9L;V%T/B_at_[hidden]*(" @(&9O;SQD;W5B;&4^*"D[#0H@(" @9F]O
H/&QO;F<@9&]U8FQE/B_at_[hidden]*#0H@(" @<F5T=7)N(# [#0I]#0H-"@``
`
end