«

Mark McIntyre wrote:

> On Sat, 03 Dec 2005 15:09:32 -0500, in comp.lang.c , Eric Sosman
> <> wrote:
>
>
>> Malcolm raises an important problem, but does not show
>>how to solve it.
>
>
> Eric raises an important problem with any solution, but fails to show
> how to solve it
>
> The solution is indeed to compare to 1+epsilon, but to define epsilon
> suitably for your application. Its likely for instance that if you're
> doing calcs that involve rounding to 5dp, then a comparison to 4dp
> will succeed. If you're doing the sum I did earlier, then epsilon
> could be 1e-10 and that'd be successful.
>
> The precise method of determining epsilon is left as an exercise for
> the reader.

What I was trying to point out is that the supposed
"solution" is no solution at all. For example,

if ( f <= 1.0 )
angle = asin(f);

is not improved by replacing 1.0 with 1.0+epsilon ...

The solution -- to the degree that one exists -- is
to use care in calculating f in the first place, and to
understand the nature of the errors that occur along the
way and the manner in which they grow and decay. That
field of study is called Numerical Analysis, and is too
rich a topic to be treated at any depth in Usenet postings.

However, the fact that floating-point numbers are usually
(but not always!) approximations and that floating-point
calculations usually (but not always!) involve error should
not be taken as meaning that all floating-point comparisions
are "dangerous!"

I once read (but cannot now find) a description of a
progression of attitudes that may apply here:

1. Initially, people trust every digit of every number
a computer calculates.

2. Having learned that some (even many) digits can be
incorrect, the once-burned-twice-shy crowd overreacts
by never trusting any floating-point calculation.

3. With further experience and study, a more rational
(pun intentional) attitude of tolerance (ditto)
eventually develops. Floating-point arithmetic is
imperfect (as is int arithmetic) but like many other
tools is useful if handled with care.

It seems to me that someone whose alarm bells clang at `f<=1.0'
may not have completed the journey from stage 2 to 3. Malcolm
is probably at 2+epsilon or 3-delta, but afraid to calculate the
remaining distance

Recommended (strongly recommended) reading: "What Every
Computer Scientist Should Know About Floating-Point Arithmetic"
by David Goldberg, widely archived on the Web.

--
Eric Sosman
lid

On Sat, 03 Dec 2005 18:29:30 -0500, in comp.lang.c , Eric Sosman
<> wrote:

>
> What I was trying to point out is that the supposed
>"solution" is no solution at all.

I realise that, but my point is that you're wrong.

>For example,
>
> if ( f <= 1.0 )
> angle = asin(f);
>
>is not improved by replacing 1.0 with 1.0+epsilon ...

Of course it is. You just need to choose epsilon carefully, and decide
how to handle failures. The choice of both is highly dependent on how
you arrived at the FP value you want to examine, and how you want to
handle the problem cases.

#define EPS -1e-12

if ( f <= 1.0+EPS )
angle = asin(f);
else
{
puts("f ->1");
f = asin(1.0);
}

> Recommended (strongly recommended) reading: "What Every
>Computer Scientist Should Know About Floating-Point Arithmetic"
>by David Goldberg, widely archived on the Web.

Absolutely.
--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>

----== Posted via Newsfeeds.Com - Unlimited-Unrestricted-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! 120,000+ Newsgroups
----= East and West-Coast Server Farms - Total Privacy via Encryption =----

Mark McIntyre wrote:
>
> On Sat, 03 Dec 2005 15:09:32 -0500, in comp.lang.c , Eric Sosman
> <> wrote:
>
> > Malcolm raises an important problem, but does not show
> >how to solve it.
>
> Eric raises an important problem with any solution, but fails to show
> how to solve it
>
> The solution is indeed to compare to 1+epsilon, but to define epsilon
> suitably for your application.

It depends on what you're doing.

In the very few times that I have written code
which compares a double variable against a double constant
with either a relational operator or an equality operator,
it would be inappropriate to make use of DBL_EPSILON.
The (DBL_MAX >= x && x > 0) expression, is correct as is,
and so is the (0 != x) expression.

double sssqrt(double x)
{
if (DBL_MAX >= x && x > 0) {
const double a = x;
double b = x / 2 + 0.5;

do {
x = b;
b = (a / x + x) / 2;
} while (x > b);
} else {
if (0 != x) {
errno = EDOM;
x = HUGE_VAL;
}
}
return x;
}

In reviewing my toy math library, I've noticed that
I don't have any code which compares two double variables
with an equality operator.

All of my epsilon usage is kind of like this:

static double p_i(void)
{
unsigned n;
double p, a, b;

p = 0;
a = 3;
n = 1;
do {
a /= 9;
b = a / n;
a /= 9;
n += 2;
b -= a / n;
n += 2;
p += b;
} while (b > DBL_EPSILON / 4);
a = 2;
n = 1;
do {
a /= 4;
b = a / n;
a /= 4;
n += 2;
b -= a / n;
n += 2;
p += b;
} while (b > DBL_EPSILON / 2);
return 4 * p;
}

--
pete

pete wrote:
> == is not a relational operator.

?


August

--
I am the "ILOVEGNU" signature virus. Just copy me to your
signature. This email was infected under the terms of the GNU
General Public License.

August Karlstrom <> writes:

> pete wrote:
>> == is not a relational operator.
>
> ?

6.5.8 Relational operators
Syntax
1 relational-expression:
shift-expression
relational-expression < shift-expression
relational-expression > shift-expression
relational-expression <= shift-expression
relational-expression >= shift-expression

[...]

6.5.9 Equality operators
Syntax
1 equality-expression:
relational-expression
equality-expression == relational-expression
equality-expression != relational-expression

--
int main(void){char p[]="ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuv wxyz.\
\n",*q="kl BIcNBFr.NKEzjwCIxNJC";int i=sizeof p/2;char *strchr();int putchar(\
);while(*q){i+=strchr(p,*q++)-p;if(i>=(int)sizeof p)i-=sizeof p-1;putchar(p[i]\
);}return 0;}

On Sat, 03 Dec 2005 23:59:57 GMT, in comp.lang.c , pete
<> wrote:

>Mark McIntyre wrote:
>
>> The solution is indeed to compare to 1+epsilon, but to define epsilon
>> suitably for your application.
>
>it would be inappropriate to make use of DBL_EPSILON.

Note carefully that I did not say DBL_EPSILON.
--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>

----== Posted via Newsfeeds.Com - Unlimited-Unrestricted-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! 120,000+ Newsgroups
----= East and West-Coast Server Farms - Total Privacy via Encryption =----

Mark McIntyre wrote:
>
> On Sat, 03 Dec 2005 23:59:57 GMT, in comp.lang.c , pete
> <> wrote:
>
> >Mark McIntyre wrote:
> >
> >> The solution is indeed to compare
> >> to 1+epsilon, but to define epsilon
> >> suitably for your application.
> >
> >it would be inappropriate to make use of DBL_EPSILON.
>
> Note carefully that I did not say DBL_EPSILON.

I would say that in cases
where a suitable definition of epsilon, is zero,
that the use of epsilon isn't required.

--
pete

Ben Pfaff wrote:
> August Karlstrom <> writes:
>
>
>>pete wrote:
>>
>>>== is not a relational operator.
>>
>>?
>
>
> 6.5.8 Relational operators
> Syntax
> 1 relational-expression:
> shift-expression
> relational-expression < shift-expression
> relational-expression > shift-expression
> relational-expression <= shift-expression
> relational-expression >= shift-expression
>
> [...]
>
> 6.5.9 Equality operators
> Syntax
> 1 equality-expression:
> relational-expression
> equality-expression == relational-expression
> equality-expression != relational-expression
>

OK, that's weird. Still `==' works as a relational operator in the
common (mathematical) sense.


August

--
I am the "ILOVEGNU" signature virus. Just copy me to your
signature. This email was infected under the terms of the GNU
General Public License.

August Karlstrom wrote:
> Ben Pfaff wrote:
>> August Karlstrom <> writes:
>>
>>
>>> pete wrote:
>>>
>>>> == is not a relational operator.
>>>
>>> ?
>>
>>
>> 6.5.8 Relational operators
>> Syntax
>> 1 relational-expression:
>> shift-expression
>> relational-expression < shift-expression
>> relational-expression > shift-expression
>> relational-expression <= shift-expression
>> relational-expression >= shift-expression
>>
>> [...]
>>
>> 6.5.9 Equality operators
>> Syntax
>> 1 equality-expression:
>> relational-expression
>> equality-expression == relational-expression
>> equality-expression != relational-expression
>>
>
> OK, that's weird. Still `==' works as a relational operator in the
> common (mathematical) sense.
>
Yes, it's a nice gotcha, isn't it? Like "byte", this is one of those terms C
appropriates in its own way and doesn't quite mean what you'd expect it to mean.

S.

August Karlstrom wrote:
>
> Ben Pfaff wrote:
> > August Karlstrom <> writes:
> >
> >
> >>pete wrote:
> >>
> >>>== is not a relational operator.
> >>
> >>?

> > 6.5.8 Relational operators

> > 6.5.9 Equality operators

> OK, that's weird. Still `==' works as a relational operator in the
> common (mathematical) sense.

I wrote what I wrote, plainly:

"== is an equality operator.
== is not a relational operator.

<= is a relational operator.
<= is not an equality operator."

--
pete