Sine code for ANSI C

[P.J. Plauger]

Re: Sine code for ANSI C

"osmium" <r124c4u102@com cast.net> wrote in message
news:c75iuo$i40 p9$1@ID-179017.news.uni-berlin.de...
[color=blue]
> CBFalconer writes:
>[color=green]
> > "P.J. Plauger" wrote:[color=darkred]
> > >[/color]
> > ... snip ...[color=darkred]
> > > coefficients. But doing proper argument reduction is an open
> > > ended exercise in frustration. Just reducing the argument modulo
> > > 2*pi quickly accumulates errors unless you do arithmetic to
> > > many extra bits of precision.[/color]
> >
> > And that problem is inherent. Adding precision bits for the
> > reduction will not help, because the input value doesn't have
> > them. It is the old problem of differences of similar sized
> > quantities.[/color]
>
> Huh? If I want the phase of an oscillator after 50,000 radians are you
> saying that is not computable? Please elaborate.
>
> There was a thread hereabouts many months ago on this very subject and[/color]
AFAIK[color=blue]
> no one suggested that it was not computable, it just couldn't be done with
> doubles. And I see no inherent problems.[/color]

Right. This difference of opinion highlights two conflicting
interpretations of floating-point numbers:

1) They're fuzzy. Assume the first discarded bit is
somewhere between zero and one. With this viewpoint,
CBFalconer is correct that there's no point in trying
to compute a sine accurately for large arguments --
all the good bits get lost.

2) They are what they are. Assume that every floating-point
representation exactly represents some value, however that
representation arose. With this viewpoint, osmium is correct
that there's a corresponding sine that is worth computing
to full machine precision.

I've gone to both extremes over the past several decades.
Our latest math library, still in internal development,
can get exact function values for *all* argument values.
It uses multi-precision argument reduction that can gust
up to over 4,000 bits [sic]. "The Standard C Library"
represents an intermediate viewpoint -- it stays exact
until about half the fraction bits go away.

I still haven't decided how hard we'll try to preserve
precision for large arguments in the next library we ship.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[CBFalconer]

Re: Sine code for ANSI C

osmium wrote:[color=blue]
> CBFalconer writes:[color=green]
> > "P.J. Plauger" wrote:[color=darkred]
> > >[/color]
> > ... snip ...[color=darkred]
> > > coefficients. But doing proper argument reduction is an open
> > > ended exercise in frustration. Just reducing the argument modulo
> > > 2*pi quickly accumulates errors unless you do arithmetic to
> > > many extra bits of precision.[/color]
> >
> > And that problem is inherent. Adding precision bits for the
> > reduction will not help, because the input value doesn't have
> > them. It is the old problem of differences of similar sized
> > quantities.[/color]
>
> Huh? If I want the phase of an oscillator after 50,000 radians are
> you saying that is not computable? Please elaborate.[/color]

PI is not a rational number, so any multiple cannot be represented
exactly. In addition, for any fixed accuracy, you will lose
significant digits in normalizing. To take your value of 50,000
radians, you are spending something like 4 decimal digits just for
the multiple of 2 PI, so the accuracy of the normalized angle will
be reduced by at least those 4 places.

To illustrate, take your handy calculator and enter:

1.23456789 + 1000000000 = <something>
and follow up with:
-10000000000 = <something else>

and you will find that <something else> is not 1.23456789. You
are running into limited precision effects.

--
"I'm a war president. I make decisions here in the Oval Office
in foreign policy matters with war on my mind." - Bush.
"Churchill and Bush can both be considered wartime leaders, just
as Secretariat and Mr Ed were both horses." - James Rhodes.

[P.J. Plauger]

Re: Sine code for ANSI C

"CBFalconer " <cbfalconer@yah oo.com> wrote in message
news:4096913D.A 4F0D63F@yahoo.c om...
[color=blue]
> osmium wrote:[color=green]
> > CBFalconer writes:[color=darkred]
> > > "P.J. Plauger" wrote:
> > > >
> > > ... snip ...
> > > > coefficients. But doing proper argument reduction is an open
> > > > ended exercise in frustration. Just reducing the argument modulo
> > > > 2*pi quickly accumulates errors unless you do arithmetic to
> > > > many extra bits of precision.
> > >
> > > And that problem is inherent. Adding precision bits for the
> > > reduction will not help, because the input value doesn't have
> > > them. It is the old problem of differences of similar sized
> > > quantities.[/color]
> >
> > Huh? If I want the phase of an oscillator after 50,000 radians are
> > you saying that is not computable? Please elaborate.[/color]
>
> PI is not a rational number, so any multiple cannot be represented
> exactly. In addition, for any fixed accuracy, you will lose
> significant digits in normalizing. To take your value of 50,000
> radians, you are spending something like 4 decimal digits just for
> the multiple of 2 PI, so the accuracy of the normalized angle will
> be reduced by at least those 4 places.[/color]

No. The *precision* will be reduced. Whether or not the *accuracy*
is reduced depends on your conjecture about the bits that don't
exist off the right end of the fraction. Those of us who write
libraries for a living shouldn't presume that those missing bits
are garbage. We serve our customers best by assuming that they're
all zeros, and delivering up the best approximation to the correct
answer for that assumed value.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[Dik T. Winter]

Re: Sine code for ANSI C

In article <40966710$1@new s101.his.com> Rich Gibbs <rgibbs@REMOVEh is.com> writes:[color=blue]
> The other is something like, "What is the best numerical method for
> calculating the sine function?" The answer to THAT question has
> something to do with C, if that is the implementation language, but has
> more to do with numerical analysis.[/color]

It has nothing to do with the C language. The newsgroup:
sci.math.num-analysis
is more appropriate. They provide algorithms, it is only the last question
how to implement the algorithm in C. When designing these algorithms you
have quite a few contradicting objectives. (Overall relative precision,
monotony, and whatever.) Do you wish a good algorithm where sin(x) > 1
could be true? (Cody-Waite, of sixties fame provided an algorithm where
that was possible.) The algorithm could still fit in your requirements:
speed and relative precision.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/

[Dik T. Winter]

Re: Sine code for ANSI C

In article <W1wlc.78210$G_ .55362@nwrddc02 .gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:[color=blue]
> I've gone to both extremes over the past several decades.
> Our latest math library, still in internal development,
> can get exact function values for *all* argument values.
> It uses multi-precision argument reduction that can gust
> up to over 4,000 bits [sic].[/color]

Have you looked at how it was done at DEC? There is a write-up about
it in the old SigNum Notices of the sixties or seventies. Mary Decker,
I think.

It has been my opinion for a very long time that he real sine function
is not so very important. The most important in numerical mathematics
is sin2pi(x) which calculates sin(2.pi.x), and which allows easy and
exact range reduction. Alas, that function is not in C.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/

[Dik T. Winter]

Re: Sine code for ANSI C

In article <WuClc.21089$sK 3.3453@nwrddc03 .gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:[color=blue]
> No. The *precision* will be reduced. Whether or not the *accuracy*
> is reduced depends on your conjecture about the bits that don't
> exist off the right end of the fraction. Those of us who write
> libraries for a living shouldn't presume that those missing bits
> are garbage. We serve our customers best by assuming that they're
> all zeros, and delivering up the best approximation to the correct
> answer for that assumed value.[/color]

Right. The basic functions should be black boxes that assume the input
value is exact. See how the precision of these functions are defined
in Ada (and I had not nothing to do with it).
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/

[CBFalconer]

Re: Sine code for ANSI C

"Dik T. Winter" wrote:[color=blue]
> "P.J. Plauger" <pjp@dinkumware .com> writes:
>[color=green]
>> I've gone to both extremes over the past several decades.
>> Our latest math library, still in internal development,
>> can get exact function values for *all* argument values.
>> It uses multi-precision argument reduction that can gust
>> up to over 4,000 bits [sic].[/color]
>
> Have you looked at how it was done at DEC? There is a write-up
> about it in the old SigNum Notices of the sixties or seventies.
> Mary Decker, I think.
>
> It has been my opinion for a very long time that he real sine
> function is not so very important. The most important in
> numerical mathematics is sin2pi(x) which calculates sin(2.pi.x),
> and which allows easy and exact range reduction. Alas, that
> function is not in C.[/color]

That just makes the problem show up elsewhere. Compare that
function with an argument of (1e18 + 0.1) against an argument of
(0.1), for example. At least on any common arithmetic system
known to me.

--
"I'm a war president. I make decisions here in the Oval Office
in foreign policy matters with war on my mind." - Bush.
"Churchill and Bush can both be considered wartime leaders, just
as Secretariat and Mr Ed were both horses." - James Rhodes.

[CBFalconer]

Re: Sine code for ANSI C

"P.J. Plauger" wrote:[color=blue]
> "CBFalconer " <cbfalconer@yah oo.com> wrote in message[color=green]
>> osmium wrote:[color=darkred]
>>> CBFalconer writes:
>>>> "P.J. Plauger" wrote:
>>>>
>>>> ... snip ...
>>>>> coefficients. But doing proper argument reduction is an open
>>>>> ended exercise in frustration. Just reducing the argument modulo
>>>>> 2*pi quickly accumulates errors unless you do arithmetic to
>>>>> many extra bits of precision.
>>>>
>>>> And that problem is inherent. Adding precision bits for the
>>>> reduction will not help, because the input value doesn't have
>>>> them. It is the old problem of differences of similar sized
>>>> quantities.
>>>
>>> Huh? If I want the phase of an oscillator after 50,000 radians are
>>> you saying that is not computable? Please elaborate.[/color]
>>
>> PI is not a rational number, so any multiple cannot be represented
>> exactly. In addition, for any fixed accuracy, you will lose
>> significant digits in normalizing. To take your value of 50,000
>> radians, you are spending something like 4 decimal digits just for
>> the multiple of 2 PI, so the accuracy of the normalized angle will
>> be reduced by at least those 4 places.[/color]
>
> No. The *precision* will be reduced. Whether or not the *accuracy*
> is reduced depends on your conjecture about the bits that don't
> exist off the right end of the fraction. Those of us who write
> libraries for a living shouldn't presume that those missing bits
> are garbage. We serve our customers best by assuming that they're
> all zeros, and delivering up the best approximation to the correct
> answer for that assumed value.[/color]

Yes, precision is a much better word. However there is an
argument for assuming the first missing bit is a 1.

--
"I'm a war president. I make decisions here in the Oval Office
in foreign policy matters with war on my mind." - Bush.
"Churchill and Bush can both be considered wartime leaders, just
as Secretariat and Mr Ed were both horses." - James Rhodes.

[P.J. Plauger]

Re: Sine code for ANSI C

"Dik T. Winter" <Dik.Winter@cwi .nl> wrote in message
news:Hx64Fv.F01 @cwi.nl...
[color=blue]
> In article <W1wlc.78210$G_ .55362@nwrddc02 .gnilink.net> "P.J. Plauger"[/color]
<pjp@dinkumware .com> writes:[color=blue][color=green]
> > I've gone to both extremes over the past several decades.
> > Our latest math library, still in internal development,
> > can get exact function values for *all* argument values.
> > It uses multi-precision argument reduction that can gust
> > up to over 4,000 bits [sic].[/color]
>
> Have you looked at how it was done at DEC? There is a write-up about
> it in the old SigNum Notices of the sixties or seventies. Mary Decker,
> I think.[/color]

I think I'm familiar with that approach, but I've studied so many
over the past ten years that I'm no longer certain. Our latest
approach is to use an array of half-precision values to represent
2/pi to as many bits as we see fit. We have a highly portable and
reasonably fast internal package that does the magic arithmetic
needed for range reduction.
[color=blue]
> It has been my opinion for a very long time that he real sine function
> is not so very important. The most important in numerical mathematics
> is sin2pi(x) which calculates sin(2.pi.x), and which allows easy and
> exact range reduction. Alas, that function is not in C.[/color]

Yes, sinq and its variants avoids the problem of needing all those
bits of pi. But many formulas actually compute radians before calling
on sin/cos/tan. For them to work in quadrants (or whole go-roundies)
would just shift the 2*pi factor to someplace in the calculation
that retains even less precision.

(If the OP is not suitably frightened about writing his/her own
sine function by now, s/he should be.)

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"CBFalconer " <cbfalconer@yah oo.com> wrote in message
news:4097068C.4 021B439@yahoo.c om...
[color=blue]
> "Dik T. Winter" wrote:[color=green]
> > "P.J. Plauger" <pjp@dinkumware .com> writes:
> >[color=darkred]
> >> I've gone to both extremes over the past several decades.
> >> Our latest math library, still in internal development,
> >> can get exact function values for *all* argument values.
> >> It uses multi-precision argument reduction that can gust
> >> up to over 4,000 bits [sic].[/color]
> >
> > Have you looked at how it was done at DEC? There is a write-up
> > about it in the old SigNum Notices of the sixties or seventies.
> > Mary Decker, I think.
> >
> > It has been my opinion for a very long time that he real sine
> > function is not so very important. The most important in
> > numerical mathematics is sin2pi(x) which calculates sin(2.pi.x),
> > and which allows easy and exact range reduction. Alas, that
> > function is not in C.[/color]
>
> That just makes the problem show up elsewhere. Compare that
> function with an argument of (1e18 + 0.1) against an argument of
> (0.1), for example. At least on any common arithmetic system
> known to me.[/color]

No, again. Working in quadrants, or go-roundies, *does* solve the
argument-reduction problem -- you just chuck the leading bits
that don't contribute to the final value. But you're *still*
conjecturing that a lack of fraction bits is a lack of accuracy,
and that's not necessarily true. Hence, it behooves the library
writer to do a good job in case it *isn't* true.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"CBFalconer " <cbfalconer@yah oo.com> wrote in message
news:4097079F.B AD3BEA1@yahoo.c om...
[color=blue]
> "P.J. Plauger" wrote:[color=green]
> > "CBFalconer " <cbfalconer@yah oo.com> wrote in message[color=darkred]
> >> osmium wrote:
> >>> CBFalconer writes:
> >>>> "P.J. Plauger" wrote:
> >>>>
> >>>> ... snip ...
> >>>>> coefficients. But doing proper argument reduction is an open
> >>>>> ended exercise in frustration. Just reducing the argument modulo
> >>>>> 2*pi quickly accumulates errors unless you do arithmetic to
> >>>>> many extra bits of precision.
> >>>>
> >>>> And that problem is inherent. Adding precision bits for the
> >>>> reduction will not help, because the input value doesn't have
> >>>> them. It is the old problem of differences of similar sized
> >>>> quantities.
> >>>
> >>> Huh? If I want the phase of an oscillator after 50,000 radians are
> >>> you saying that is not computable? Please elaborate.
> >>
> >> PI is not a rational number, so any multiple cannot be represented
> >> exactly. In addition, for any fixed accuracy, you will lose
> >> significant digits in normalizing. To take your value of 50,000
> >> radians, you are spending something like 4 decimal digits just for
> >> the multiple of 2 PI, so the accuracy of the normalized angle will
> >> be reduced by at least those 4 places.[/color]
> >
> > No. The *precision* will be reduced. Whether or not the *accuracy*
> > is reduced depends on your conjecture about the bits that don't
> > exist off the right end of the fraction. Those of us who write
> > libraries for a living shouldn't presume that those missing bits
> > are garbage. We serve our customers best by assuming that they're
> > all zeros, and delivering up the best approximation to the correct
> > answer for that assumed value.[/color]
>
> Yes, precision is a much better word. However there is an
> argument for assuming the first missing bit is a 1.[/color]

Right, and I acknowledged that argument earlier in this thread.
But when you're talking about how to write a professional sine
function, the argument is irrelevant.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[CBFalconer]

Re: Sine code for ANSI C

"P.J. Plauger" wrote:[color=blue]
>
> "CBFalconer " <cbfalconer@yah oo.com> wrote in message
> news:4097068C.4 021B439@yahoo.c om...
>[color=green]
> > "Dik T. Winter" wrote:[color=darkred]
> > > "P.J. Plauger" <pjp@dinkumware .com> writes:
> > >
> > >> I've gone to both extremes over the past several decades.
> > >> Our latest math library, still in internal development,
> > >> can get exact function values for *all* argument values.
> > >> It uses multi-precision argument reduction that can gust
> > >> up to over 4,000 bits [sic].
> > >
> > > Have you looked at how it was done at DEC? There is a write-up
> > > about it in the old SigNum Notices of the sixties or seventies.
> > > Mary Decker, I think.
> > >
> > > It has been my opinion for a very long time that he real sine
> > > function is not so very important. The most important in
> > > numerical mathematics is sin2pi(x) which calculates sin(2.pi.x),
> > > and which allows easy and exact range reduction. Alas, that
> > > function is not in C.[/color]
> >
> > That just makes the problem show up elsewhere. Compare that
> > function with an argument of (1e18 + 0.1) against an argument of
> > (0.1), for example. At least on any common arithmetic system
> > known to me.[/color]
>
> No, again. Working in quadrants, or go-roundies, *does* solve the
> argument-reduction problem -- you just chuck the leading bits
> that don't contribute to the final value. But you're *still*
> conjecturing that a lack of fraction bits is a lack of accuracy,
> and that's not necessarily true. Hence, it behooves the library
> writer to do a good job in case it *isn't* true.[/color]

In the example I gave above, using any binary representation, it
is true. Maybe I should have used + (1.0/3.0) to ensure problems
on any non-trinary system :-)

At any rate, the ideal system produces instantaneous results with
no range limitations and the full accuracy of the arithmetic
representation. So the designer should pick a compromise that
will suffice for most anticipated users. Here there be dragons.

--
"I'm a war president. I make decisions here in the Oval Office
in foreign policy matters with war on my mind." - Bush.
"Churchill and Bush can both be considered wartime leaders, just
as Secretariat and Mr Ed were both horses." - James Rhodes.

[P.J. Plauger]

Re: Sine code for ANSI C

"CBFalconer " <cbfalconer@yah oo.com> wrote in message
news:409764B2.4 E0914FA@yahoo.c om...
[color=blue][color=green][color=darkred]
> > > > It has been my opinion for a very long time that he real sine
> > > > function is not so very important. The most important in
> > > > numerical mathematics is sin2pi(x) which calculates sin(2.pi.x),
> > > > and which allows easy and exact range reduction. Alas, that
> > > > function is not in C.
> > >
> > > That just makes the problem show up elsewhere. Compare that
> > > function with an argument of (1e18 + 0.1) against an argument of
> > > (0.1), for example. At least on any common arithmetic system
> > > known to me.[/color]
> >
> > No, again. Working in quadrants, or go-roundies, *does* solve the
> > argument-reduction problem -- you just chuck the leading bits
> > that don't contribute to the final value. But you're *still*
> > conjecturing that a lack of fraction bits is a lack of accuracy,
> > and that's not necessarily true. Hence, it behooves the library
> > writer to do a good job in case it *isn't* true.[/color]
>
> In the example I gave above, using any binary representation, it
> is true.[/color]

What is true? A typical IEEE double representation will render
1e18 + 0.1 as 1e18, if that's what you intend to signify.
When fed as input to a function, the function can't know that
some poor misguided programmer had aspirations for some of the
bits not retained. But that has nothing to do with how hard
the sine function should work to deliver the best representation
of sin(x) when x == 1e18. You're once again asking the library
writer to make up stories about bits dead and gone -- that's
the responsibility of the programmer in designing the program
and judging the meaningfulness of the result returned by sin(x)
*for that particular program.* The library writer can't and
shouldn't care about what might have been.
[color=blue]
> Maybe I should have used + (1.0/3.0) to ensure problems
> on any non-trinary system :-)[/color]

So what if you did? That still doesn't change the analysis.
There exists a best internal approximation for 1/3 in any given
floating-point representation. There also exists a best internal
approximation for the sine *of that representation. * The library
can't and shouldn't surmise that 0.333333 (speaking decimal)
is the leading part of 0.33333333333.. . It could just as well
be the leading part of 0.3333334.
[color=blue]
> At any rate, the ideal system produces instantaneous results with
> no range limitations and the full accuracy of the arithmetic
> representation. So the designer should pick a compromise that
> will suffice for most anticipated users. Here there be dragons.[/color]

I mostly agree with that sentiment, except for the "most anticipated
users" part. If I followed that maxim to its logical extreme, my
company would produce very fast math functions that yield results
of mediocre accuracy, with spurious overflows and traps for extreme
arguments. An incredibly large number of math libraries meet that
de facto spec, *and the users almost never complain.*

We aim to please more exacting users, however, to the point that
they often don't notice the pains we've taken to stay out of their
way. So, to review your checklist:

-- We don't aim for instantaneous results, even though a few of
the better designed modern CPUs can now evaluate the common
math functions in single hardware instructions. Our compromise
is to keep the code as fast as possible while still writing in
highly portable C.

-- We *do* aim for no range limitations. We think through every
possible input representation, and produce either a suitable
special value (Inf, -Inf, 0, -0, or NaN) or a finite result close
to the best internal approximation of the function value for that
input, treating the input as exact.

-- We provide full accuracy only for a handful of the commonest
math functions. But we generally aim for worst-case errors of
3 ulp (units in the least-significant place) for float, 4 for
double, and 5 for 113-bit long double. The rare occasions
where we fail to achieve this goal are around the zeros of
messy functions.

Most important, we have tools to *measure* how well we do,
and how well other libraries do. At some point in the not
too distant future, we'll publish some objective comparisons
at our web site. We find that programmers are more likely to
consider upgrading an existing library when they have a
better understanding of its limitations.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[CBFalconer]

Re: Sine code for ANSI C

"P.J. Plauger" wrote:[color=blue]
> "CBFalconer " <cbfalconer@yah oo.com> wrote in message
>[/color]
.... snip ...[color=blue]
>[color=green]
> > At any rate, the ideal system produces instantaneous results with
> > no range limitations and the full accuracy of the arithmetic
> > representation. So the designer should pick a compromise that
> > will suffice for most anticipated users. Here there be dragons.[/color]
>
> I mostly agree with that sentiment, except for the "most anticipated
> users" part. If I followed that maxim to its logical extreme, my
> company would produce very fast math functions that yield results
> of mediocre accuracy, with spurious overflows and traps for extreme
> arguments. An incredibly large number of math libraries meet that
> de facto spec, *and the users almost never complain.*[/color]

They probably know no better. However I would not castigate traps
for extremes arguments too much - it is surely better for the user
to know something has happened than to use values with zero
significant digits. In a way it is unfortunate that hardware has
replaced software floating point, because it makes it
impracticable to fix problems.

Not all applications require the same FP library, and the usual
library/linker etc construction of C applications makes such
customization fairly easy.
[color=blue]
>
> We aim to please more exacting users, however, to the point that
> they often don't notice the pains we've taken to stay out of their
> way. So, to review your checklist:
>
> -- We don't aim for instantaneous results, even though a few of
> the better designed modern CPUs can now evaluate the common
> math functions in single hardware instructions. Our compromise
> is to keep the code as fast as possible while still writing in
> highly portable C.
>
> -- We *do* aim for no range limitations. We think through every
> possible input representation, and produce either a suitable
> special value (Inf, -Inf, 0, -0, or NaN) or a finite result close
> to the best internal approximation of the function value for that
> input, treating the input as exact.[/color]

A precision measure (count of significant bits) would be handy.
[color=blue]
>
> -- We provide full accuracy only for a handful of the commonest
> math functions. But we generally aim for worst-case errors of
> 3 ulp (units in the least-significant place) for float, 4 for
> double, and 5 for 113-bit long double. The rare occasions
> where we fail to achieve this goal are around the zeros of
> messy functions.[/color]

Those objectives are generally broad enough to eliminate any worry
over whether the next missing bit is 0 or 1.
[color=blue]
>
> Most important, we have tools to *measure* how well we do,
> and how well other libraries do. At some point in the not
> too distant future, we'll publish some objective comparisons
> at our web site. We find that programmers are more likely to
> consider upgrading an existing library when they have a
> better understanding of its limitations.[/color]

Constructing those tools, in itself, is an application requiring a
much different math library. I suspect Gauss has no small hand.

--
Chuck F (cbfalconer@yah oo.com) (cbfalconer@wor ldnet.att.net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net> USE worldnet address!

[P.J. Plauger]

Re: Sine code for ANSI C

"CBFalconer " <cbfalconer@yah oo.com> wrote in message
news:40979D8B.C 3EEBDEC@yahoo.c om...
[color=blue][color=green]
> > I mostly agree with that sentiment, except for the "most anticipated
> > users" part. If I followed that maxim to its logical extreme, my
> > company would produce very fast math functions that yield results
> > of mediocre accuracy, with spurious overflows and traps for extreme
> > arguments. An incredibly large number of math libraries meet that
> > de facto spec, *and the users almost never complain.*[/color]
>
> They probably know no better. However I would not castigate traps
> for extremes arguments too much - it is surely better for the user
> to know something has happened than to use values with zero
> significant digits.[/color]

Only sin/cos/tan suffers from catastrophic significance loss, and
these functions rarely trap. I'm talking about things like

hypot(DBL_MAX / 2, DBL_MAX / 2)

which has a well defined and precise value, but which often
causes an intermediate overflow and/or trap in naive implementations .
[color=blue]
> In a way it is unfortunate that hardware has
> replaced software floating point, because it makes it
> impracticable to fix problems.[/color]

I don't understand that comment.
[color=blue]
> Not all applications require the same FP library, and the usual
> library/linker etc construction of C applications makes such
> customization fairly easy.[/color]

Indeed. That's how we can sell replacement and add-on libraries.
[color=blue][color=green]
> > We aim to please more exacting users, however, to the point that
> > they often don't notice the pains we've taken to stay out of their
> > way. So, to review your checklist:
> >
> > -- We don't aim for instantaneous results, even though a few of
> > the better designed modern CPUs can now evaluate the common
> > math functions in single hardware instructions. Our compromise
> > is to keep the code as fast as possible while still writing in
> > highly portable C.
> >
> > -- We *do* aim for no range limitations. We think through every
> > possible input representation, and produce either a suitable
> > special value (Inf, -Inf, 0, -0, or NaN) or a finite result close
> > to the best internal approximation of the function value for that
> > input, treating the input as exact.[/color]
>
> A precision measure (count of significant bits) would be handy.[/color]

That's the ulp measure described below.
[color=blue][color=green]
> > -- We provide full accuracy only for a handful of the commonest
> > math functions. But we generally aim for worst-case errors of
> > 3 ulp (units in the least-significant place) for float, 4 for
> > double, and 5 for 113-bit long double. The rare occasions
> > where we fail to achieve this goal are around the zeros of
> > messy functions.[/color]
>
> Those objectives are generally broad enough to eliminate any worry
> over whether the next missing bit is 0 or 1.[/color]

No, they're orthogonal to that issue, as I keep trying to explain
to you. We do the best we can with the function arguments given.
It is the inescapable responsibility of the programmer to
understand the effect of uncertainties in those argument values on
the uncertainties in the function results.
[color=blue][color=green]
> > Most important, we have tools to *measure* how well we do,
> > and how well other libraries do. At some point in the not
> > too distant future, we'll publish some objective comparisons
> > at our web site. We find that programmers are more likely to
> > consider upgrading an existing library when they have a
> > better understanding of its limitations.[/color]
>
> Constructing those tools, in itself, is an application requiring a
> much different math library. I suspect Gauss has no small hand.[/color]

His spirit indeed lives on. We mostly use Maplesoft these
days, having long since given up on Mathematica as too
buggy.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com