Sine code for ANSI C

[Dan Pop]

Re: Sine code for ANSI C

In <slrnca75nr.5gr .sam@ID-227112.user.uni-berlin.de> Sam Dennis <sam@malfunctio n.screaming.net > writes:
[color=blue]
>P.J. Plauger wrote:[color=green]
>> [Trigonometry on angles that make ordinary libraries dizzy and such]
>> Until the C Standard lets us off the hook beyond some given argument
>> magnitude, however, we have no excuse not to try.[/color]
>
>I think that the Standard lets implementors off the hook at any and all
>magnitudes, signs and, generally, values, actually, and not even _just_
>for the library... from the outset! C99, 5.2.4.2.2p4:
>
>`The accuracy of the floating-point operations (+, -, *, /) and of the
>library functions in <math.h> and <complex.h> that return floating-point
>results is implementation-defined. The implementation may state that
>the accuracy is unknown.'
>
>So if I write my implementation to always produce 42. for an expression
>involving floating-point mathematics (except in the very few situations
>where there are absolute requirements imposed) and document this, it is
>not disqualified from full conforming hosted implementation status (not
>through another (?) loophole).[/color]

True and irrelevant, as we're discussing about what a high quality
implementation should do.

If you're looking for loopholes in the C standard, it contains one so
large as to make even a *completely* useless implementation (a two line
shell script) fully conforming.

1 The implementation shall be able to translate and execute at
^^
least one program that contains at least one instance of every
^^^^^^^^^^^^^^^ ^^
one of the following limits:

Other, strictly conforming programs, must be merely accepted (whatever
that means, the standard provides no further clues).

So consider the following "implementation " for Unix systems:

echo "Program accepted."
cp /bin/true a.out

combined with an implementation of /bin/true that exercises all the
translation limits enumerated in 5.2.4.1 and the documentation required
by the standard.

Now, every translation unit, either correct or not, will receive the one
diagnostic required by the standard, all strictly conforming programs are
"accepted" (as well as all other programs ;-) and the *one* C program
mentioned above is translated and executed.

The real question is: how many users is my implementation going to have?
Answer: as many as yours ;-) This illustrates the importance of a factor
that is not even mentioned in the C standard: the quality of
implementation.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Sam Dennis]

Re: Sine code for ANSI C

Dan Pop wrote:[color=blue]
> Sam Dennis <sam@malfunctio n.screaming.net > writes:[color=green]
>>`The accuracy of the floating-point operations (+, -, *, /) and of the
>>library functions in <math.h> and <complex.h> that return floating-point
>>results is implementation-defined. The implementation may state that
>>the accuracy is unknown.'
>>
>>So if I write my implementation to always produce 42. [...][/color]
>
> True and irrelevant, as we're discussing about what a high quality
> implementation should do.[/color]

Ridiculous accuracy for silly arguments is unnecessary for high quality
programs, which should accommodate libraries coded by mere mortals that
can furthermore calculate logarithms in less time than it takes for the
operator to look up the answer himself. (Admittedly, such programs can
not operate on my hypothetical implementation if they actually make use
of the floating-point facilities provided by C itself; the portable way
to get guaranteed accuracy is to use a library and its types for all of
one's non-integral calculations, written in strictly conforming C.)
[color=blue]
> If you're looking for loopholes in the C standard, it contains one so
> large as to make even a *completely* useless implementation (a two line
> shell script) fully conforming.
>
> 1 The implementation shall be able to translate and execute at
> least one program that contains at least one instance of every
> one of the following limits:[/color]

I've decided that that was probably a political decision to allow buggy
compilers to claim conformance. Specifying accuracy for floating-point
operations such that enough current implementations can conform, on the
other hand, is a genuinely hellish task.
[color=blue]
> [A conforming implementation, once documentation is attached]
> echo "Program accepted."
> cp /bin/true a.out[/color]

(The extension below can be ignored if you don't consider `accepted' to
require even an appropriate diagnostic from a one line file that starts
with `#error'.)

Well, obviously there's the #error thing, which means that #include has
to be handled, which means...

cat "$@"
echo "Program accepted."
cp /bin/true a.out

....and the documentation must state that #include source file searching
is not supported and that all output is diagnostics (or similar words).
But, yes, that'd be a conforming implementation; and, because any input
is `acceptable' to it, any data is a conforming program; that's no news
to you, of course.

Yes, the Standard has huge loopholes. Those two are best ignored, with
useful definitions of `conforming' implementation and programs mentally
substituted (in the latter case, by `strictly conforming program').
[color=blue]
> This illustrates the importance of a factor that is not even mentioned
> in the C standard: the quality of implementation.[/color]

Oh, but it would be a perfectly good implementation in most other ways.

Actually, I've been thinking for a long time about writing this, except
that the quality of implementation, as it's normally understood, should
be variable, for conformance testing of programs.

Catching all the potential problems isn't feasible, but a useful number
is probably just about manageable.

--
++acr@,ka"

[glen herrmannsfeldt]

Re: Sine code for ANSI C

P.J. Plauger wrote:
(I wrote)
[color=blue][color=green]
>>You say that, and then your argument seems to indicate the
>>opposite. If the argument was generated using floating
>>point arithmetic that either truncated our rounded the
>>result (do you know which?), the uncertainty is at least
>>0.5ulp.[/color][/color]
[color=blue]
> Another conjecture. And what if it wasn't? Then you'd have no
> excuse other than to take the argument at face value and do
> the best you can with it. I say again -- it is up to the programmer
> to do error-propagation analysis and decide how much precision to
> trust in the final answer. The library writer doesn't have the
> luxury of giving back garbage bits because they *might* not be
> meaningful.[/color]

Consider this program fragment, as no real examples have
been presented so far:

for(x=0;x+60!=x ;x+=60) printf("%20.15g \n",sin(x*A)) ;

A should have the appropriate value so that x is
expressed in degrees. Now, consider the effect
of ordinary argument reduction, and unusually
accurate argument reduction.

If you do argument reduction using an appropriate
approximation to pi/4 to machine precision, you get
the answer expected by the programmer. If you use
a pi/4 much more accurate than machine precision, the
result will slowly diverge from the expected value.

[color=blue]
> But I kinda like the freedom you offer us. We could write *all*
> library functions on the basis that the return value can be
> anywhere between f(x - 1/2 ulp) and f(x + 1/2 ulp) -- speaking
> symbolically, of course. Boy would that make it easier to compute
> exponentials, powers, gamma functions, etc. But I bet you wouldn't
> like the results.[/color]

(snip)
[color=blue][color=green]
>>Which value of pi do you use? Rounded to the precision of
>>the argument? (As the user might have used in generating it).
>>Much more accurate than the argument? (Unlikely used by
>>the user.)[/color][/color]
[color=blue]
> Again, you're *presuming* that the user did some sloppy calculation
> to arrive at a large argument to sine. If that's what happened, then
> the result will be correspondingly inaccurate. But that's *always*
> true of every calculation, whether it involves sin(x) or any other
> math function. It's true even if no math functions are called at
> all. The authors of the math functions can only take what they're
> given and do the best possible job with them. Why is that such
> a hard principle to grasp?[/color]

Why did people buy expensive computers that Cray built,
even though they couldn't even divide to machine precision?

The IBM 360/91 was documented as not following the S/360
architecture because its floating point divide gave rounded
results (0.5ulp) instead of the truncated results (1ulp)
that the architecture specified.

Consider someone doing a single precision sine. Most
likely they use single precision instead of double
because they don't need so much accuracy and hope that
the result will be generated faster.

As there are plenty of multiple precision libraries
around, users wanting more than machine precision
know where to find it.

The OS/360 math library writers had to work extra hard
to get accurate results given the base 16 floating
point of S/360. Sometimes it took a little extra work
to get good answers, which is fine. (The last iteration
on the square root algorithm is done slightly different
to avoid precision loss in a halve operation.) A lot of
extra work to get useless answers is not. (The sine
routine gives a fatal error when the argument is greater
than pi*2**18 (single), pi*2**50 (double), or pi*2**100
(extended precision).

Giving a fatal error tells the user that something is
wrong and should be fixed. Supplying answers using bits
that don't really exist allows the user to believe that
useful answers are being generated, possibly wasting
much calculation and money.

Actually, the time when this happened to me, probably
when I was in ninth grade, (due to an undefined variable
if I remember right), someone carefully explained to me
why the library would refuse such a calculation. It made
sense at the time, and it still does.

-- glen

[Richard Bos]

Re: Sine code for ANSI C

glen herrmannsfeldt <gah@ugcs.calte ch.edu> wrote:
[color=blue]
> P.J. Plauger wrote:
>[color=green]
> > Again, you're *presuming* that the user did some sloppy calculation
> > to arrive at a large argument to sine. If that's what happened, then
> > the result will be correspondingly inaccurate. But that's *always*
> > true of every calculation, whether it involves sin(x) or any other
> > math function. It's true even if no math functions are called at
> > all. The authors of the math functions can only take what they're
> > given and do the best possible job with them. Why is that such
> > a hard principle to grasp?[/color]
>
> Why did people buy expensive computers that Cray built,
> even though they couldn't even divide to machine precision?[/color]

Because Crays are _very_ specialised machines, which can make
assumptions that a library writer for a normal computer cannot.

I must agree with Mr. Plauger, here. If a library gives me a value with
more precision than I need or expect, I can always ignore that extra
precision. If it gives me a value with less precision than it could, as
some people seem to be advocating, I can _never_ add that precision
myself without doing the library's job for it, which should not be the
aim of a good quality library.

Richard

[Dan Pop]

Re: Sine code for ANSI C

In <slrnca7h8o.5gr .sam@ID-227112.user.uni-berlin.de> Sam Dennis <sam@malfunctio n.screaming.net > writes:
[color=blue]
>Dan Pop wrote:[color=green]
>> If you're looking for loopholes in the C standard, it contains one so
>> large as to make even a *completely* useless implementation (a two line
>> shell script) fully conforming.
>>
>> 1 The implementation shall be able to translate and execute at
>> least one program that contains at least one instance of every
>> one of the following limits:[/color]
>
>I've decided that that was probably a political decision to allow buggy
>compilers to claim conformance.[/color]

It's actually a pragmatic decision: for any given implementation, it is
possible to construct a strictly conforming program exceeding the
implementation' s resources, while still staying within the translation
limits mentioned in the standard.

One year or so ago I explored the ways of avoiding this wording by
adding a few more translation limits. Then, the standard could reasonably
require the translation and execution of *any* program staying within the
translation limits.
[color=blue]
>Specifying accuracy for floating-point
>operations such that enough current implementations can conform, on the
>other hand, is a genuinely hellish task.[/color]

I don't think so. If the standard could specify minimal accuracy
requirements for the floating point types, it could also specify
minimal accuracy requirements for the floating point operations. I have
already debated this issue in comp.std.c and my opponents came with bogus
counterargument s (like the unfeasibility of documenting the accuracy of
all the <math.h> functions, although I was strictly and explicitly
talking about the floating point operations).
[color=blue][color=green]
>> [A conforming implementation, once documentation is attached]
>> echo "Program accepted."
>> cp /bin/true a.out[/color]
>
>(The extension below can be ignored if you don't consider `accepted' to
>require even an appropriate diagnostic from a one line file that starts
>with `#error'.)[/color]

Who says that the one and only correctly translated program has to
contain the #error pragma?
[color=blue][color=green]
>> This illustrates the importance of a factor that is not even mentioned
>> in the C standard: the quality of implementation.[/color]
>
>Oh, but it would be a perfectly good implementation in most other ways.[/color]

I.e. for people with no need for floating point. This is actually the
case with certain implementations for embedded control applications:
floating point support is too expensive and irrelevant to the typical
application domain of the processor, yet the standard says that it must
be provided. A much better solution would have been to make floating
point support optional for freestanding implementations .

In the case of hosted implementations , people would choose the ones
providing proper support for floating point even if they don't need it.

OTOH, considering that practically every hosted implementation in current
use today is for systems with native floating point support, the issue is
moot: the implementation delivers whatever the hardware delivers.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Dan Pop]

Re: Sine code for ANSI C

In <40a46d80.65302 293@news.indivi dual.net> rlb@hoekstra-uitgeverij.nl (Richard Bos) writes:
[color=blue]
>glen herrmannsfeldt <gah@ugcs.calte ch.edu> wrote:
>[color=green]
>> P.J. Plauger wrote:
>>[color=darkred]
>> > Again, you're *presuming* that the user did some sloppy calculation
>> > to arrive at a large argument to sine. If that's what happened, then
>> > the result will be correspondingly inaccurate. But that's *always*
>> > true of every calculation, whether it involves sin(x) or any other
>> > math function. It's true even if no math functions are called at
>> > all. The authors of the math functions can only take what they're
>> > given and do the best possible job with them. Why is that such
>> > a hard principle to grasp?[/color]
>>
>> Why did people buy expensive computers that Cray built,
>> even though they couldn't even divide to machine precision?[/color]
>
>Because Crays are _very_ specialised machines, which can make
>assumptions that a library writer for a normal computer cannot.
>
>I must agree with Mr. Plauger, here. If a library gives me a value with
>more precision than I need or expect, I can always ignore that extra
>precision. If it gives me a value with less precision than it could, as
>some people seem to be advocating, I can _never_ add that precision
>myself without doing the library's job for it, which should not be the
>aim of a good quality library.[/color]

You have failed to understand the substance of this discussion. Nobody
argued for less precission than could be *meaningfully* computed.

Back to my extreme example, one could consider DBL_MAX as an *exact*
value and spend significant CPU resources computing sin(DBL_MAX) with
double precision accuracy, but when would be that precision *meaningful*?

Give the user as much precision as it makes sense from the magnitude of
the argument and let *him* do the argument reduction if he needs more
precision, because only he can decide what is the proper way of performing
the argument reduction for his application. More often than not, trying
to outsmart the user results in dumb behaviour.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Dik T. Winter]

Re: Sine code for ANSI C

In article <c7tf3u$iub$1@s unnews.cern.ch> Dan.Pop@cern.ch (Dan Pop) writes:[color=blue]
> In <HxLvHG.5J@cwi. nl> "Dik T. Winter" <Dik.Winter@cwi .nl> writes:
>[color=green]
> >In article <c7t8i0$1sv$1@s unnews.cern.ch> Dan.Pop@cern.ch (Dan Pop) writes:[color=darkred]
> > > In <1dpoc.118917$G _.72884@nwrddc0 2.gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:[/color]
> >...[color=darkred]
> > > >> Concrete examples, please.
> > > >
> > > >What is the sine of 162,873 radians? If you're working in radians,
> > > >you can represent this input value *exactly* even in a float.
> > >
> > > You *can*, but does it make physical sense to call sine with an integer
> > > argument (even if represented as a float)?[/color]
> >
> >Must everything make physical sense? Perhaps it makes mathematical sense?[/color]
>
> And does it? Last time I checked, mathematics define the sine as having
> a real argument and the C programming language provides zilch support for
> real numbers.[/color]

Yup.
[color=blue][color=green][color=darkred]
> > > In real life applications, the argument of sine is computed using
> > > floating point arithmetic (on non-integer values), so it *is* a fuzzy
> > > value, with the degree of fuzziness implied by its magnitude.[/color]
> >
> >Not in mathematical applications, where the argument to the sine function
> >can very well be exact.[/color]
>
> Please elaborate, again, with concrete examples.[/color]

You want a concrete example, I just do think that such examples are possible.
[color=blue]
> Where do these
> "mathematic al applications" get their input data from? How do they
> handle the *incorrect* (from the mathematics POV) result of calling sine
> for pretty much *any* argument except 0?[/color]

Pretty much as in every such case. Careful error analysis.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:V7Xoc.6287 $6f5.499223@att bi_s54...
[color=blue]
> Consider this program fragment, as no real examples have
> been presented so far:
>
> for(x=0;x+60!=x ;x+=60) printf("%20.15g \n",sin(x*A)) ;
>
> A should have the appropriate value so that x is
> expressed in degrees. Now, consider the effect
> of ordinary argument reduction, and unusually
> accurate argument reduction.
>
> If you do argument reduction using an appropriate
> approximation to pi/4 to machine precision, you get
> the answer expected by the programmer. If you use
> a pi/4 much more accurate than machine precision, the
> result will slowly diverge from the expected value.[/color]

Okay, so what? The programmer has unrealistic expectations.
It's the library vendor's responsibility to meet possibly
realistic expectations. The vendor can't possibly guess
all the ways that library functions will be called with
arguments that need "correcting " so as not to hurt the
programmer's feelings.

What does this tell the designer of a sine function about
where it's okay to stop delivering accurate results?
[color=blue][color=green]
> > Again, you're *presuming* that the user did some sloppy calculation
> > to arrive at a large argument to sine. If that's what happened, then
> > the result will be correspondingly inaccurate. But that's *always*
> > true of every calculation, whether it involves sin(x) or any other
> > math function. It's true even if no math functions are called at
> > all. The authors of the math functions can only take what they're
> > given and do the best possible job with them. Why is that such
> > a hard principle to grasp?[/color]
>
> Why did people buy expensive computers that Cray built,
> even though they couldn't even divide to machine precision?[/color]

Not sure what this has to do with the current discussion, but
I'll answer anyway. There was a market for machines that were
fast -- at least by the standards of those days -- and a prevailing
view that inaccurate floating-point was safe enough if you just
carried a few extra sacrificial precision bits. That view was
correct often enough that not all of Cray's customers went
bankrupt.


What does this tell the designer of a sine function about
where it's okay to stop delivering accurate results?
[color=blue]
> The IBM 360/91 was documented as not following the S/360
> architecture because its floating point divide gave rounded
> results (0.5ulp) instead of the truncated results (1ulp)
> that the architecture specified.[/color]

As before.


What does this tell the designer of a sine function about
where it's okay to stop delivering accurate results?
[color=blue]
> Consider someone doing a single precision sine. Most
> likely they use single precision instead of double
> because they don't need so much accuracy and hope that
> the result will be generated faster.[/color]

Most likely.

What does this tell the designer of a sine function about
where it's okay to stop delivering accurate results?
[color=blue]
> As there are plenty of multiple precision libraries
> around, users wanting more than machine precision
> know where to find it.[/color]


Okay.

What does this tell the designer of a sine function about
where it's okay to stop delivering accurate results?
[color=blue]
> The OS/360 math library writers had to work extra hard
> to get accurate results given the base 16 floating
> point of S/360. Sometimes it took a little extra work
> to get good answers, which is fine. (The last iteration
> on the square root algorithm is done slightly different
> to avoid precision loss in a halve operation.) A lot of
> extra work to get useless answers is not. (The sine
> routine gives a fatal error when the argument is greater
> than pi*2**18 (single), pi*2**50 (double), or pi*2**100
> (extended precision).[/color]

Nice criteria. I don't see them anywhere in the C Standard.
Are you going to propose them as a DR? Be fun to see you
defend why you pick the exponent 18 for 24-bit (IEEE)
precision, 50 for 53-bit and 100 for 64-bit(!) or 113-bit.

And what "fatal" error are you going to introduce, given that
the math library currently has no such concept. Even IEEE 654
always provides for continuing execution with some error
code. Looks like you're bravely guiding the standards community
into new and uncharted territory. Brave of you.
[color=blue]
> Giving a fatal error tells the user that something is
> wrong and should be fixed.[/color]

Can you quantify the "wrongness" in terms of your criteria?
Should they not also be applied to the various rounding
functions? If not, why not?
[color=blue]
> Supplying answers using bits
> that don't really exist allows the user to believe that
> useful answers are being generated, possibly wasting
> much calculation and money.[/color]

Happens all the time, when users call math functions with
arguments having bits that don't really exist. Just whose
responsibility is that, anyway?
[color=blue]
> Actually, the time when this happened to me, probably
> when I was in ninth grade, (due to an undefined variable
> if I remember right), someone carefully explained to me
> why the library would refuse such a calculation. It made
> sense at the time, and it still does.[/color]

Things that make sense for a ninth-grade class don't necessarily
scale to professional libraries licensed to professional
programmers, both working to language standards developed
by international organizations.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[Dan Pop]

Re: Sine code for ANSI C

In <HxpE44.Myq@cwi .nl> "Dik T. Winter" <Dik.Winter@cwi .nl> writes:
[color=blue]
>In article <c7tf3u$iub$1@s unnews.cern.ch> Dan.Pop@cern.ch (Dan Pop) writes:[color=green]
> > In <HxLvHG.5J@cwi. nl> "Dik T. Winter" <Dik.Winter@cwi .nl> writes:
> >[color=darkred]
> > >In article <c7t8i0$1sv$1@s unnews.cern.ch> Dan.Pop@cern.ch (Dan Pop) writes:
> > > > In <1dpoc.118917$G _.72884@nwrddc0 2.gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:
> > >...
> > > > >> Concrete examples, please.
> > > > >
> > > > >What is the sine of 162,873 radians? If you're working in radians,
> > > > >you can represent this input value *exactly* even in a float.
> > > >
> > > > You *can*, but does it make physical sense to call sine with an integer
> > > > argument (even if represented as a float)?
> > >
> > >Must everything make physical sense? Perhaps it makes mathematical sense?[/color]
> >
> > And does it? Last time I checked, mathematics define the sine as having
> > a real argument and the C programming language provides zilch support for
> > real numbers.[/color]
>
>Yup.[/color]

Yup what?!? Please elaborate.
[color=blue][color=green][color=darkred]
> > > > In real life applications, the argument of sine is computed using
> > > > floating point arithmetic (on non-integer values), so it *is* a fuzzy
> > > > value, with the degree of fuzziness implied by its magnitude.
> > >
> > >Not in mathematical applications, where the argument to the sine function
> > >can very well be exact.[/color]
> >
> > Please elaborate, again, with concrete examples.[/color]
>
>You want a concrete example, I just do think that such examples are possible.[/color]

This is not good enough.
[color=blue][color=green]
> > "mathematic al applications" get their input data from? How do they
> > handle the *incorrect* (from the mathematics POV) result of calling sine
> > for pretty much *any* argument except 0?[/color]
>
>Pretty much as in every such case. Careful error analysis.[/color]

With the exception of numerical analysis, mathematical results are exact,
by definition.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Christian Bau]

Re: Sine code for ANSI C

In article <C43pc.127380$G _.63467@nwrddc0 2.gnilink.net>,
"P.J. Plauger" <pjp@dinkumware .com> wrote:

<snipped lots of stuff>

Apart from the precision of an individual call to sin(x), would people
think that fabs (sin (x)) <= 1.0 must be true for _all_ values of x, no
matter how large? And should sin (x)*sin(x) + cos(x)*cos(x) be close to
1 for all x, no matter how large?

[glen herrmannsfeldt]

Re: Sine code for ANSI C

P.J. Plauger wrote:[color=blue]
> "glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
> news:V7Xoc.6287 $6f5.499223@att bi_s54...[/color]
[color=blue][color=green]
>>Consider this program fragment, as no real examples have
>>been presented so far:[/color][/color]
[color=blue][color=green]
>>for(x=0;x+60! =x;x+=60) printf("%20.15g \n",sin(x*A)) ;[/color][/color]

(snip)
[color=blue]
> Okay, so what? The programmer has unrealistic expectations.
> It's the library vendor's responsibility to meet possibly
> realistic expectations. The vendor can't possibly guess
> all the ways that library functions will be called with
> arguments that need "correcting " so as not to hurt the
> programmer's feelings.[/color]

That example is more realistic than any you have posted
so far.
[color=blue]
> What does this tell the designer of a sine function about
> where it's okay to stop delivering accurate results?[/color]

(snip)
[color=blue][color=green]
>>Consider someone doing a single precision sine. Most
>>likely they use single precision instead of double
>>because they don't need so much accuracy and hope that
>>the result will be generated faster.[/color][/color]
[color=blue]
> Most likely.[/color]
[color=blue]
> What does this tell the designer of a sine function about
> where it's okay to stop delivering accurate results?[/color]

When they use a single precision function they expect
less accurate answers than a double precision function.
[color=blue][color=green]
>>As there are plenty of multiple precision libraries
>>around, users wanting more than machine precision
>>know where to find it.[/color][/color]
[color=blue]
> Okay.[/color]
[color=blue]
> What does this tell the designer of a sine function about
> where it's okay to stop delivering accurate results?[/color]

It says that when the operations require exceptional
ability to perform, far above and beyond the call of
duty, it is time to give up. When thousands of
digits of pi are required to compute the result on
a 64 bit operand, something is wrong.
[color=blue][color=green]
>>The OS/360 math library writers had to work extra hard
>>to get accurate results given the base 16 floating
>>point of S/360. Sometimes it took a little extra work
>>to get good answers, which is fine. (The last iteration
>>on the square root algorithm is done slightly different
>>to avoid precision loss in a halve operation.) A lot of
>>extra work to get useless answers is not. (The sine
>>routine gives a fatal error when the argument is greater
>>than pi*2**18 (single), pi*2**50 (double), or pi*2**100
>>(extended precision).[/color][/color]
[color=blue]
> Nice criteria. I don't see them anywhere in the C Standard.
> Are you going to propose them as a DR? Be fun to see you
> defend why you pick the exponent 18 for 24-bit (IEEE)
> precision, 50 for 53-bit and 100 for 64-bit(!) or 113-bit.[/color]

The requirements of the C standard have already been
discussed, and that didn't seem to bother you any.

The numbers above are not for IEEE, as it hadn't been
invented yet. Extended precision has 112 bit fraction.
[color=blue]
> And what "fatal" error are you going to introduce, given that
> the math library currently has no such concept. Even IEEE 654
> always provides for continuing execution with some error
> code. Looks like you're bravely guiding the standards community
> into new and uncharted territory. Brave of you.[/color]

It seems that on most implementations sqrt() has the
ability to generate a fatal error when given a negative
argument.
[color=blue][color=green]
>>Giving a fatal error tells the user that something is
>>wrong and should be fixed.[/color][/color]
[color=blue]
> Can you quantify the "wrongness" in terms of your criteria?
> Should they not also be applied to the various rounding
> functions? If not, why not?[/color]

There are no problems of any use to anyone where they
are useful.

(snip)
[color=blue]
> Happens all the time, when users call math functions with
> arguments having bits that don't really exist. Just whose
> responsibility is that, anyway?[/color]
[color=blue][color=green]
>>Actually, the time when this happened to me, probably
>>when I was in ninth grade, (due to an undefined variable
>>if I remember right), someone carefully explained to me
>>why the library would refuse such a calculation. It made
>>sense at the time, and it still does.[/color][/color]
[color=blue]
> Things that make sense for a ninth-grade class don't necessarily
> scale to professional libraries licensed to professional
> programmers, both working to language standards developed
> by international organizations.[/color]

It wasn't in a class. They didn't teach it in ninth grade,
and as far as I know they don't now. It was explained to
me by a professional physicist.

-- glen

[Christian Bau]

Re: Sine code for ANSI C

In article <a59pc.1495$gr. 83718@attbi_s52 >,
glen herrmannsfeldt <gah@ugcs.calte ch.edu> wrote:
[color=blue]
> It says that when the operations require exceptional
> ability to perform, far above and beyond the call of
> duty, it is time to give up. When thousands of
> digits of pi are required to compute the result on
> a 64 bit operand, something is wrong.[/color]

Just for a thought: There is this thing called "inverse error analysis".
You try to calculate y = f (x). Usually you assume that your
floating-point math will produce a value y' instead which is not quite
equal to y, and you want to keep the difference between y and y' small.
But instead you might assume that you calculate y = f (x'), that is your
floating-point arithmetic produced the exact function value for some
different x', hopefully close to x.

Now consider y = cos (x). For x = 0, cos (x) = 1, and it becomes smaller
very slowly; for small x we have cos (x) approximately 1 - x^2/2. Now
take 64 bit IEEE floating point. The largest number smaller than 1 is 1
- 2^-53. Choose x such that cos x is about 1 - 2^-54, x = pow (2.0,
-53.0/2). For this x, the best that your implementation of cos (x) can
do is to return a result of either 1 or 1 - 2^-53. 1 = cos (x') for x' =
0, 1 - 2^-53 = cos (x') for x' about pow (2.0, -26). So the best you can
possibly get is cos (x'), with x' either zero or 41 percent larger than
x. And there is nothing you can do about it.

No matter how good your cos () implementation is, the worst case with
inverse error analysis is horrendously bad.

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:a59pc.1495 $gr.83718@attbi _s52...
[color=blue]
> P.J. Plauger wrote:[color=green]
> > "glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
> > news:V7Xoc.6287 $6f5.499223@att bi_s54...[/color]
>[color=green][color=darkred]
> >>Consider this program fragment, as no real examples have
> >>been presented so far:[/color][/color]
>[color=green][color=darkred]
> >>for(x=0;x+60! =x;x+=60) printf("%20.15g \n",sin(x*A)) ;[/color][/color]
>
> (snip)
>[color=green]
> > Okay, so what? The programmer has unrealistic expectations.
> > It's the library vendor's responsibility to meet possibly
> > realistic expectations. The vendor can't possibly guess
> > all the ways that library functions will be called with
> > arguments that need "correcting " so as not to hurt the
> > programmer's feelings.[/color]
>
> That example is more realistic than any you have posted
> so far.[/color]

Again, so what? We're talking about the requirements placed on
a vendor of high quality math functions. How various innocents
misuse the library doesn't give the vendor any more latitude.
It's what the *professionals* expect, and the C Standard
indicates, that matter. Sadly, the C Standard gives *no*
latitude for copping out once an argument to sine gets large
in magnitude.
[color=blue][color=green][color=darkred]
> >>Consider someone doing a single precision sine. Most
> >>likely they use single precision instead of double
> >>because they don't need so much accuracy and hope that
> >>the result will be generated faster.[/color][/color]
>[color=green]
> > Most likely.[/color]
>[color=green]
> > What does this tell the designer of a sine function about
> > where it's okay to stop delivering accurate results?[/color]
>
> When they use a single precision function they expect
> less accurate answers than a double precision function.[/color]

No, they expect less *precise* answers. There's a difference,
and until you understand it you're not well equipped to
critique the design of professional math libraries.
[color=blue][color=green][color=darkred]
> >>As there are plenty of multiple precision libraries
> >>around, users wanting more than machine precision
> >>know where to find it.[/color][/color]
>[color=green]
> > Okay.[/color]
>[color=green]
> > What does this tell the designer of a sine function about
> > where it's okay to stop delivering accurate results?[/color]
>
> It says that when the operations require exceptional
> ability to perform, far above and beyond the call of
> duty, it is time to give up. When thousands of
> digits of pi are required to compute the result on
> a 64 bit operand, something is wrong.[/color]

"Far above and beyond" is once again qualitative arm waving.
You may hear a weaker call to duty than others, but there's
nowhere in the C Standard that justifies an arbitrary
point beyond which it's okay to give up on sine, or any other
function. You think sine is hard to get right for some
arguments, try lgamma, tgamma, and erfc (all in C99).
[color=blue][color=green][color=darkred]
> >>The OS/360 math library writers had to work extra hard
> >>to get accurate results given the base 16 floating
> >>point of S/360. Sometimes it took a little extra work
> >>to get good answers, which is fine. (The last iteration
> >>on the square root algorithm is done slightly different
> >>to avoid precision loss in a halve operation.) A lot of
> >>extra work to get useless answers is not. (The sine
> >>routine gives a fatal error when the argument is greater
> >>than pi*2**18 (single), pi*2**50 (double), or pi*2**100
> >>(extended precision).[/color][/color]
>[color=green]
> > Nice criteria. I don't see them anywhere in the C Standard.
> > Are you going to propose them as a DR? Be fun to see you
> > defend why you pick the exponent 18 for 24-bit (IEEE)
> > precision, 50 for 53-bit and 100 for 64-bit(!) or 113-bit.[/color]
>
> The requirements of the C standard have already been
> discussed, and that didn't seem to bother you any.[/color]

What, you mean that a conforming implementation has no explicit
precision requirements? I understand that; in fact I was one
of the principal proponents of that (lack of) requirement in
the C Standard. Dan Pop already responded to that issue better
than I can. It is well understood that every implementation
has certain "quality of implementation" features. (I also
pioneered use of that term.) If getc takes two days per character,
nobody will buy your library. Similarly, if math functions
trap out on certain unspecified values, word will get out and
your customers will stay away in droves.

Math libraries are tougher to quantify, because they require so
many measurements. But once you do the measurements and publish
them, you'd better be prepared to justify the results. You can't
just say, "I decided not to implement pow for negative exponents."
Or more to the point, "at some point the results of sine begin
to really suck, but I'd rather not say where that is." And if
you *do* choose to publish your threshold of suckiness, you'd
better have a rationale for choosing it. Otherwise, some third
party tester will report that your sine function occasionally
is off by 40 billion ulp, and not provide the rationale for why
that might be.
[color=blue]
> The numbers above are not for IEEE, as it hadn't been
> invented yet. Extended precision has 112 bit fraction.[/color]

Oh, I see. You were talking about criteria accompanying a
40-year-old architecture, with software technology to match.
[color=blue][color=green]
> > And what "fatal" error are you going to introduce, given that
> > the math library currently has no such concept. Even IEEE 654
> > always provides for continuing execution with some error
> > code. Looks like you're bravely guiding the standards community
> > into new and uncharted territory. Brave of you.[/color]
>
> It seems that on most implementations sqrt() has the
> ability to generate a fatal error when given a negative
> argument.[/color]

Really? The C Standard says it should set errno to EDOM and
return some failure value. If your implementation also obeys
IEC 60559 (aka IEEE 754) then you should return a NaN. I
find very few customers who want a program to terminate rather
than produce a NaN. I know even fewer who would *ever* expect
sine to return a NaN for any finite argument.
[color=blue][color=green][color=darkred]
> >>Giving a fatal error tells the user that something is
> >>wrong and should be fixed.[/color][/color]
>[color=green]
> > Can you quantify the "wrongness" in terms of your criteria?
> > Should they not also be applied to the various rounding
> > functions? If not, why not?[/color]
>
> There are no problems of any use to anyone where they
> are useful.[/color]

Is that a Zen koan or just the nonsense it appears to be?
[color=blue][color=green]
> > Things that make sense for a ninth-grade class don't necessarily
> > scale to professional libraries licensed to professional
> > programmers, both working to language standards developed
> > by international organizations.[/color]
>
> It wasn't in a class. They didn't teach it in ninth grade,
> and as far as I know they don't now. It was explained to
> me by a professional physicist.[/color]

As a one-time professional physicist myself, I now begin to
understand the nature of your education.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[Dr Chaos]

Re: Sine code for ANSI C

On Fri, 14 May 2004 21:23:51 GMT, P.J. Plauger <pjp@dinkumware .com> wrote:[color=blue]
> Again, so what? We're talking about the requirements placed on
> a vendor of high quality math functions. How various innocents
> misuse the library doesn't give the vendor any more latitude.
> It's what the *professionals* expect, and the C Standard
> indicates, that matter. Sadly, the C Standard gives *no*
> latitude for copping out once an argument to sine gets large
> in magnitude.[/color]

Then that's a likely thought-bug in the C Standard.

[color=blue]
>[color=green][color=darkred]
>> >>Consider someone doing a single precision sine. Most
>> >>likely they use single precision instead of double
>> >>because they don't need so much accuracy and hope that
>> >>the result will be generated faster.[/color]
>>[color=darkred]
>> > Most likely.[/color]
>>[color=darkred]
>> > What does this tell the designer of a sine function about
>> > where it's okay to stop delivering accurate results?[/color]
>>
>> When they use a single precision function they expect
>> less accurate answers than a double precision function.[/color]
>
> No, they expect less *precise* answers. There's a difference,
> and until you understand it you're not well equipped to
> critique the design of professional math libraries.[/color]

The design of the professional math libraries is not the issue, it's
whether the effort is worthwhile, as opposed to accomodating likely
poorly-thought out algorithms by the user.

I think accumulation of rotations is probably best done
with complex multiplication.

[P.J. Plauger]

Re: Sine code for ANSI C

"Dr Chaos" <mbkennelSPAMBE GONE@NOSPAMyaho o.com> wrote in message
news:slrncaagm5 .gu0.mbkennelSP AMBEGONE@lyapun ov.ucsd.edu...
[color=blue]
> On Fri, 14 May 2004 21:23:51 GMT, P.J. Plauger <pjp@dinkumware .com> wrote:[color=green]
> > Again, so what? We're talking about the requirements placed on
> > a vendor of high quality math functions. How various innocents
> > misuse the library doesn't give the vendor any more latitude.
> > It's what the *professionals* expect, and the C Standard
> > indicates, that matter. Sadly, the C Standard gives *no*
> > latitude for copping out once an argument to sine gets large
> > in magnitude.[/color]
>
> Then that's a likely thought-bug in the C Standard.[/color]

Likely, but it's there, and some of us have to live with it.
[color=blue][color=green][color=darkred]
> >> >>Consider someone doing a single precision sine. Most
> >> >>likely they use single precision instead of double
> >> >>because they don't need so much accuracy and hope that
> >> >>the result will be generated faster.
> >>
> >> > Most likely.
> >>
> >> > What does this tell the designer of a sine function about
> >> > where it's okay to stop delivering accurate results?
> >>
> >> When they use a single precision function they expect
> >> less accurate answers than a double precision function.[/color]
> >
> > No, they expect less *precise* answers. There's a difference,
> > and until you understand it you're not well equipped to
> > critique the design of professional math libraries.[/color]
>
> The design of the professional math libraries is not the issue, it's
> whether the effort is worthwhile, as opposed to accomodating likely
> poorly-thought out algorithms by the user.[/color]

The design of professional math libraries *is* the issue. Until
such time as standards quantify what calculations are "worthwhile "
and what merely accommodate poorly thought out algorithms, we
have an obligation to assume that whatever is specified might
be considered worthwhile to some serious users.

The other side of the coin is knowing where to stop once the
"worthwhile " police get empowered. Several people who have
contributed to this thread are convinced that computing the
sine of a sufficiently large angle is not worthwhile, but *nobody*
has ventured a cutoff point that has any defensible logic behind
it. And I assure you that as soon as any such defense is mounted,
I and others can apply it to a variety of other math functions.
You will then hear the usual howls, "but *that's* different."
[color=blue]
> I think accumulation of rotations is probably best done
> with complex multiplication.[/color]

And why do you think this can be done with any more accuracy,
or precision, than the techniques cited (and sneered at) so far
for generating large angles?

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com