Sine code for ANSI C

[glen herrmannsfeldt]

Re: Sine code for ANSI C

CBFalconer wrote:

(snip regarding sin() function)
[color=blue]
> 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]


When I was in high school I knew someone with a brand
new HP-55 calculator. (You can figure out when that was
if you want.) He was so proud of it, and sure of the
answers it gave. For the sine (in degrees) of 9.999999999e99
it gave something like 0.53, which is obviously wrong
because 9.999999999e99 is a multiple of 180.

Yes, argument reduction is always a problem, because
people will expect it to be right, no matter how useless
the result is. It is a little easier in degrees than
in radians, yet few languages (PL/I being one) support
SIND() and such.

-- glen

[glen herrmannsfeldt]

Re: Sine code for ANSI C

CBFalconer wrote:[color=blue]
> "P.J. Plauger" wrote:[/color]
[color=blue]
> ... snip ...[/color]
[color=blue][color=green]
>>I agree that it's a Quality of Implementation issue just how fast
>>a library function supplies nonsense when called with nonsense
>>arguments. But I've yet to hear an objective criterion for
>>determining how much of an argument is garbage. Absent that, the
>>best way I know for library writers to satisfy customers is to
>>assume that every input value is exact, and to produce the
>>closest possible representation to the nearest internal
>>representatio n of the corresponding function value.[/color][/color]

(snip)
[color=blue]
> I am gradually coming around to your point of view here, so I am
> rewording it. To me, the argument is that the input argument
> represents a range of values (absent contrary information) with a
> known uncertainty. The most probable actual value is the exact
> value. Since d(sinx)/dx is strictly bounded, the resultant
> function is never going to fly off in all directions, and will not
> become totally meaningless unless that argument uncertainty is in
> the order of PI.[/color]

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=blue]
> We can, at the cost of some computational complexity, assume the
> input argument is exact and reduce it to a principal value with an
> accuracy governed by the precision of PI. There is no essential
> difference between this operation and forming a real value for
> 1/3. The critical thing there is the precision of our knowledge
> of 3.[/color]

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=blue]
> For a concrete example, consider a mechanism that rotates and
> detents at 45 degree intervals. The sine of the resultant angle
> can only have 5 discrete values. However the net rotation is
> described by the sum of a counter (of full turns) and one of 8
> discrete angles, in units of 45 deg. After letting that engine
> whir for a day and a half at a 1000 rpm and a half and recording
> the final angle, we want to know the sine of that angle. The
> computation machinery knows nothing about detents, all it has to
> work with is PI, the net rotation angle, and the computation
> algorithm for sin(x).[/color]

Now you multiply the angle, in multiples of 45 degrees,
by (pi/4) accurate to 53 or so bits for an IEEE double.
If the argument reduction is done with a 4000 bit accurate
pi, you find many more than 5 values for sine.
[color=blue]
> At some point the accuracy of the results will become worse than
> the accuracy of the detents, and all blows up. This is not the
> same point as that reached by simple modulo PI arithmetic.[/color]

It may have been fixed by now, but it was well known in
the past that Postscript could not do a proper 90 degree
or 180 degree rotation. It did it by calculating the sine
and cosine of angles in radians, converted from degrees.
It seems that the value of pi used was different than that
used in argument reduction, so that sin(180 degrees) was
not zero. If it is not zero, it is possible that a horizontal
line in the input will not be horizontal, in the output,
when rounded (or truncated, I forget) to pixel positions.
I have been told by professional typesetters that it was
visible in the output of high resolution phototypesetter s.

-- glen

[glen herrmannsfeldt]

Re: Sine code for ANSI C

Dik T. Winter wrote:
(snip)
(someone wrote)
[color=blue][color=green][color=darkred]
> > >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.[/color][/color][/color]
(someone else wrote)
[color=blue][color=green]
> > You *can*, but does it make physical sense to call sine with an integer
> > argument (even if represented as a float)?[/color][/color]
[color=blue]
> Must everything make physical sense? Perhaps it makes mathematical sense?[/color]

(snip)
[color=blue]
> Not in mathematical applications, where the argument to the sine function
> can very well be exact.[/color]

Show me a useful, real life, mathematical problem that requires the
evaluation of 167873 radians. Now, I can almost imagine problems
involving large integer numbers of degrees. If those degrees
are multiplied by (pi/180), possibly with a different value
of pi than is used in the argument reduction, all useful bits
are lost.

Yes, radians are nice in many cases, but this is not one.

-- glen

[glen herrmannsfeldt]

Re: Sine code for ANSI C

P.J. Plauger wrote:

(snip, someone wrote)
[color=blue][color=green][color=darkred]
>>>Huh? If I want the phase of an oscillator after 50,000 radians are
>>>you saying that is not computable? Please elaborate.[/color][/color][/color]

(someone else wrote)
[color=blue][color=green]
>>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][/color]
[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]

Why assume they are all zero? In what cases would someone
want the sine of a number that was a large integer number of
radians? Well, 50,000 isn't so large, but say that it
was 1e50?

The first math library that I used would refuse to do
double precision sine at pi*2**50, pretty much no fraction
bits left. I can almost imagine some cases for an integer
number of degrees, but that would almost only make sense
in decimal floating point arithmetic. (Many calculators
but few computers.)

-- glen

[glen herrmannsfeldt]

Re: Sine code for ANSI C

P.J. Plauger wrote:
(snip)
[color=blue]
> Yes, it's concrete enough to show that you still don't get it.[/color]

(snip)
[color=blue]
> You may think that it's unlikely that all the bits of that value are
> meaningful in a given application, and you're probably right. But
> you're not *definitely* right. So as a library vendor, best engineering
> practice is to supply all the bits that *might* make sense, in case
> they actually do. A good quality implementation will continue to
> compute the sine of small arguments quickly, but if it has to take
> progressively longer to reduce ever larger arguments, then so be it.
> You don't want the cost, reduce the arguments quickly, by your own
> crude methods that are good enough for your application, before calling
> sine.[/color]

By returning a value you are making a claim that there is some
sense in calculating that value. There is no sense in calculating
such large values of radians such that the uncertainty in the
angle is greater than pi.
[color=blue]
> Note also that even sin(1e-38) *might not* be good to full precision,
> because the argument happens not to be good to full precision, but
> you have no criterion for judging how many bits are worth computing.
> Since it's so easy to compute them, and since they're often all
> meaningful, nobody wastes much time debating the cost of producing
> bits that are potentially garbage.[/color]

Well, small arguments it radians are easy to compute. One
could, for example, do a numeric derivative at this point,
say (sin(2e-38)-sin(1e-38))/(2e38-1e38) and likely get an
exact 1.00.

(snip)
[color=blue]
> You also didn't get that you're still not giving a sufficiently concrete
> criterion for when the sine function should stop trying. If you don't
> like that a sine function wastes your program microseconds computing
> what you shouldn't have told it to compute, then you need to be more
> precise about where it can and should give up. You state above that
> the function result is definitely meaningless once 1 ulp in the argument
> weighs more than 2*pi, but why go that far? Aside from a phase factor,
> you've lost all angle information at pi/2. But then how meaningful is it
> to have just a couple of bits of fraction information? To repeat what
> you stated above:[/color]

(snip)

Any of those are close enough for me.
[color=blue]
> So you've taken on a serious responsibility here. You have to tell us
> library vendors just how small "not so large" is, so we know where to
> produce quick garbage instead of slower correct answers. If you don't,
> you have no right to deem our products unacceptable, or even simply
> wasteful.[/color]

At that point, all you are saying is that the user should do
their own argument reduction, using the appropriate method
which only the user can know, before calling sin().

Java, at least, defines a standard value for pi, so that programs
can used that in generating arguments. Should you assume that
the argument, in radians, is a multiple of that value of pi,
or a much more accurate one?

-- glen

[glen herrmannsfeldt]

Re: Sine code for ANSI C

osmium wrote:

(snip)
[color=blue]
> You might mention to Mr. Pop the notion of a revolution counter, such as a
> Veeder-Root counter on the propeller shaft of an ocean going ship. My
> guess is that there are a lot of radians between England and New Jersey. If
> Dan Pop needs the name of a particular ship to make the point sufficiently
> concrete, I expect some one could provide the name of such a ship.[/color]

Design a rotating shaft counter that can count in exact radians
and I will figure out how to calculate the sine of it.

I can easily design one that will count revolutions, degrees,
or most any other integer multiple of revolutions, but not
radians.

There is no point in doing argument reduction with an exact,
to thousands of bits, representation of pi when the user can't
generate arguments with such pi, and has no source of such
arguments.

Go over to comp.lang.pl1, and you will find the sind() and
cosd() functions, which could do exact argument reduction
on large integers. (I have no idea whether they do or not.)

-- glen

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:GKGoc.8085 6$Ik.6105956@at tbi_s53...
[color=blue]
> P.J. Plauger wrote:
>
> (snip, someone wrote)
>[color=green][color=darkred]
> >>>Huh? If I want the phase of an oscillator after 50,000 radians are
> >>>you saying that is not computable? Please elaborate.[/color][/color]
>
> (someone else wrote)
>[color=green][color=darkred]
> >>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][/color]
>[color=green]
> > 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]
>
> Why assume they are all zero?[/color]

For the same reason that the library assumes 0.5 (which really looks
like 0.5000000 or thereabouts) is 0.5. The value of a floating-point
number is defined by its bits. Part of that definition effectively
requires you to treat all the bits you don't see as zeros.
[color=blue]
> In what cases would someone
> want the sine of a number that was a large integer number of
> radians? Well, 50,000 isn't so large, but say that it
> was 1e50?[/color]

One of the interesting challenges for us library writers is that
we have to serve all possible customers. We don't get to rule on
which ones are being silly -- or overly demanding. But if we
slight so much as a single case, however extreme or apparently
unlikely, customers rightly castigate us, and reviewers emphasize
the failure all out of proportion to the number of cases we get
right. So *you* can decide that 50,000 radians is worth reducing
properly but 1e50 is not. Mazeltov. Until the C Standard lets us
off the hook beyond some given argument magnitude, however, we
have no excuse not to try.

In fact, one of the commercial test suites we currently use tests
sin(x) up to about 1e18. Is that silly? I can't find any place in
the C Standard that says it isn't. So far, no major potential
customer has seized on this suite as an acceptance test, but *you*
try to explain to some high-level decision maker why it's okay
to pay us for a product that fails a handful of tests.
[color=blue]
> The first math library that I used would refuse to do
> double precision sine at pi*2**50, pretty much no fraction
> bits left.[/color]

That's nice. Who got to pick the cutoff point?
[color=blue]
> I can almost imagine some cases for an integer
> number of degrees, but that would almost only make sense
> in decimal floating point arithmetic. (Many calculators
> but few computers.)[/color]

Really? I can't imagine why.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:HZGoc.7992 7$0H1.7555079@a ttbi_s54...
[color=blue]
> Dik T. Winter wrote:
> (snip)
> (someone wrote)
>[color=green][color=darkred]
> > > >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.[/color][/color]
> (someone else wrote)
>[color=green][color=darkred]
> > > You *can*, but does it make physical sense to call sine with an[/color][/color][/color]
integer[color=blue][color=green][color=darkred]
> > > argument (even if represented as a float)?[/color][/color]
>[color=green]
> > Must everything make physical sense? Perhaps it makes mathematical[/color][/color]
sense?[color=blue]
>
> (snip)
>[color=green]
> > Not in mathematical applications, where the argument to the sine[/color][/color]
function[color=blue][color=green]
> > can very well be exact.[/color]
>
> Show me a useful, real life, mathematical problem that requires the
> evaluation of 167873 radians.[/color]

Now you're playing Dan Pop's game. Somebody contrives a problem
that's arguably practical and you get to rule, from the catbird
seat, that it's not really a real-life problem. (And yet both
you and Dan Pop toss out equally academic problems as though
they're pragmatic engineering issues.)

But it doesn't *matter* whether or not I can score points at your
game, because that's not the game I'm in (as I keep repeating).
My company endeavors to make high quality libraries that satisfy
programming language standards and meet the needs of demanding
customers. We don't have the luxury of demanding in turn that our
customers prove that everything they do makes sense. Our treaty
point is the C Standard, in this case. We don't match it, they
have a reason to be pissed. We match it, they don't have a case.
(Well, there are also QOI issues, but that's not at issue here.)
[color=blue]
> Now, I can almost imagine problems
> involving large integer numbers of degrees.[/color]

Why is that better? It's just easier to comprehend the truncations
involved, and easier to compute the associated nonsense, or valid
reduced argument. Still the same problem.
[color=blue]
> If those degrees
> are multiplied by (pi/180), possibly with a different value
> of pi than is used in the argument reduction, all useful bits
> are lost.[/color]

Indeed. As a library vendor, you use the best possible value of
pi. If that's not what the customer used, it's the customer's
problem. (I once had a customer who was sure we had a bad tangent
function, until I pointed out that he was approximating pi in
his program with 22/7.)
[color=blue]
> Yes, radians are nice in many cases, but this is not one.[/color]

Only because it's harder to reduce the argument accurately.
In either case, you can't give any more meaning to an argument
reduced from a large magnitude, *or any less.*

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:CcHoc.3405 $1q3.396588@att bi_s51...
[color=blue]
> osmium wrote:
>
> (snip)
>[color=green]
> > You might mention to Mr. Pop the notion of a revolution counter, such as[/color][/color]
a[color=blue][color=green]
> > Veeder-Root counter on the propeller shaft of an ocean going ship. My
> > guess is that there are a lot of radians between England and New Jersey.[/color][/color]
If[color=blue][color=green]
> > Dan Pop needs the name of a particular ship to make the point[/color][/color]
sufficiently[color=blue][color=green]
> > concrete, I expect some one could provide the name of such a ship.[/color]
>
> Design a rotating shaft counter that can count in exact radians
> and I will figure out how to calculate the sine of it.
>
> I can easily design one that will count revolutions, degrees,
> or most any other integer multiple of revolutions, but not
> radians.
>
> There is no point in doing argument reduction with an exact,
> to thousands of bits, representation of pi when the user can't
> generate arguments with such pi, and has no source of such
> arguments.[/color]

Your logic is specious. What if the customer is computing
directly in radians? It's an artifact of the library function
that requires the use of many bits of pi. It turns out that
the easiest way to compute pi is to convert it to quadrants
and compute it over the interval [-0.5, 0.5]. So we have to
work in high precision to give the customer the best answer.
The customer *does not* have to work in equally high precision
to give us an input worthy of computing a sine.

You've once again fallen into the trap of making up stories
about the validity of arguments presented to a library function.
That's *not* a luxury we library writers can indulge.
[color=blue]
> Go over to comp.lang.pl1, and you will find the sind() and
> cosd() functions, which could do exact argument reduction
> on large integers. (I have no idea whether they do or not.)[/color]

Yes, it's easier to do. Why do you think the resulting reduced
arguments are any worthier of evaluating, or any less?

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:PnHoc.8119 9$Ik.6124808@at tbi_s53...
[color=blue]
> P.J. Plauger wrote:
> (snip)
>[color=green]
> > Yes, it's concrete enough to show that you still don't get it.[/color]
>
> (snip)
>[color=green]
> > You may think that it's unlikely that all the bits of that value are
> > meaningful in a given application, and you're probably right. But
> > you're not *definitely* right. So as a library vendor, best engineering
> > practice is to supply all the bits that *might* make sense, in case
> > they actually do. A good quality implementation will continue to
> > compute the sine of small arguments quickly, but if it has to take
> > progressively longer to reduce ever larger arguments, then so be it.
> > You don't want the cost, reduce the arguments quickly, by your own
> > crude methods that are good enough for your application, before calling
> > sine.[/color]
>
> By returning a value you are making a claim that there is some
> sense in calculating that value. There is no sense in calculating
> such large values of radians such that the uncertainty in the
> angle is greater than pi.[/color]

You've made up a story again, about the "uncertaint y" in the reduced
angle. All the library writer knows is that there's little *precision*
in the representation of that reduced angle. But it still has a well
defined value and it may well be meaningful to the customer.
[color=blue][color=green]
> > Note also that even sin(1e-38) *might not* be good to full precision,
> > because the argument happens not to be good to full precision, but
> > you have no criterion for judging how many bits are worth computing.
> > Since it's so easy to compute them, and since they're often all
> > meaningful, nobody wastes much time debating the cost of producing
> > bits that are potentially garbage.[/color]
>
> Well, small arguments it radians are easy to compute. One
> could, for example, do a numeric derivative at this point,
> say (sin(2e-38)-sin(1e-38))/(2e38-1e38) and likely get an
> exact 1.00.[/color]

You missed the point. Of *course* they're easy to compute. So
too is the sine of 1e38 quadrants -- it's zero. But should we
decide how much garbage to give back to customers based on how
easy it is for us to compute it? Suppose I make up a story that
most people who call sine with tiny arguments don't know what
the hell they're doing. I could modulate your claim above:

: By returning a value you are making a claim that there is some
: sense in calculating that value. There is no sense in calculating
: such tiny values of radians such that the magnitude of the
: angle is much less than than pi.

Who knows, maybe I could even cite studies. But in any case, I
wouldn't last long in front of a customer waving a bug report
under my nose. And I would have no defense in the C Standard.

Conjectures about the validity of input arguments just don't
cut it.
[color=blue]
> (snip)
>[color=green]
> > You also didn't get that you're still not giving a sufficiently concrete
> > criterion for when the sine function should stop trying. If you don't
> > like that a sine function wastes your program microseconds computing
> > what you shouldn't have told it to compute, then you need to be more
> > precise about where it can and should give up. You state above that
> > the function result is definitely meaningless once 1 ulp in the argument
> > weighs more than 2*pi, but why go that far? Aside from a phase factor,
> > you've lost all angle information at pi/2. But then how meaningful is it
> > to have just a couple of bits of fraction information? To repeat what
> > you stated above:[/color]
>
> (snip)
>
> Any of those are close enough for me.[/color]

That's nice. Would you mind writing that up as a DR? I'd love it
if the C Standard was changed to make less work for me.
[color=blue][color=green]
> > So you've taken on a serious responsibility here. You have to tell us
> > library vendors just how small "not so large" is, so we know where to
> > produce quick garbage instead of slower correct answers. If you don't,
> > you have no right to deem our products unacceptable, or even simply
> > wasteful.[/color]
>
> At that point, all you are saying is that the user should do
> their own argument reduction, using the appropriate method
> which only the user can know, before calling sin().[/color]

Yes, that's exactly what I'm saying. And I said it in response to
Dan Pop's claim. He alleges that taking a long time to compute the
sine of a large angle is "wasteful" of CPU cycles, even "unacceptab le."

Groucho Marks: Doctor, doctor, it hurts when I do this. What should I do?

Doctor: Don't do that.
[color=blue]
> Java, at least, defines a standard value for pi, so that programs
> can used that in generating arguments. Should you assume that
> the argument, in radians, is a multiple of that value of pi,
> or a much more accurate one?[/color]

Whatever value pi may have when truncated to machine precision
is useless for reducing arguments properly. The error shows up
after just a few radians. It's hard to defend a conjecture that
the sine of 5*pi is never worth computing properly.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"glen herrmannsfeldt" <gah@ugcs.calte ch.edu> wrote in message
news:AzHoc.3903 8$xw3.2533208@a ttbi_s04...
[color=blue]
> CBFalconer wrote:[color=green]
> > "P.J. Plauger" wrote:[/color]
>[color=green]
> > ... snip ...[/color]
>[color=green][color=darkred]
> >>I agree that it's a Quality of Implementation issue just how fast
> >>a library function supplies nonsense when called with nonsense
> >>arguments. But I've yet to hear an objective criterion for
> >>determining how much of an argument is garbage. Absent that, the
> >>best way I know for library writers to satisfy customers is to
> >>assume that every input value is exact, and to produce the
> >>closest possible representation to the nearest internal
> >>representatio n of the corresponding function value.[/color][/color]
>
> (snip)
>[color=green]
> > I am gradually coming around to your point of view here, so I am
> > rewording it. To me, the argument is that the input argument
> > represents a range of values (absent contrary information) with a
> > known uncertainty. The most probable actual value is the exact
> > value. Since d(sinx)/dx is strictly bounded, the resultant
> > function is never going to fly off in all directions, and will not
> > become totally meaningless unless that argument uncertainty is in
> > the order of PI.[/color]
>
> 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]

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.

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=blue][color=green]
> > We can, at the cost of some computational complexity, assume the
> > input argument is exact and reduce it to a principal value with an
> > accuracy governed by the precision of PI. There is no essential
> > difference between this operation and forming a real value for
> > 1/3. The critical thing there is the precision of our knowledge
> > of 3.[/color]
>
> 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]

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?

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[CBFalconer]

Re: Sine code for ANSI C

"P.J. Plauger" wrote:[color=blue]
>[/color]
.... snip ...[color=blue]
>
> Whatever value pi may have when truncated to machine precision
> is useless for reducing arguments properly. The error shows up
> after just a few radians. It's hard to defend a conjecture that
> the sine of 5*pi is never worth computing properly.[/color]

This I have to strongly disagree with. What you have to do is
effectively extend the significant bits of the argument while
reducing it. As long as you assume the argument is exact, you can
do this from the unlimited supply of zeroes. The accuracy of 10
mod 3 is dependant on the accuracy of 3. It cannot exceed the
accuracy of either 10 or 3, but it need lose nothing in its
computation.

--
A: Because it fouls the order in which people normally read text.
Q: Why is top-posting such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

[CBFalconer]

Re: Sine code for ANSI C

glen herrmannsfeldt wrote:[color=blue]
> CBFalconer wrote:[color=green]
>> "P.J. Plauger" wrote:[/color]
>[color=green]
>> ... snip ...[/color]
>[color=green][color=darkred]
>>> I agree that it's a Quality of Implementation issue just how fast
>>> a library function supplies nonsense when called with nonsense
>>> arguments. But I've yet to hear an objective criterion for
>>> determining how much of an argument is garbage. Absent that, the
>>> best way I know for library writers to satisfy customers is to
>>> assume that every input value is exact, and to produce the
>>> closest possible representation to the nearest internal
>>> representation of the corresponding function value.[/color][/color]
>
> (snip)
>[color=green]
>> I am gradually coming around to your point of view here, so I am
>> rewording it. To me, the argument is that the input argument
>> represents a range of values (absent contrary information) with a
>> known uncertainty. The most probable actual value is the exact
>> value. Since d(sinx)/dx is strictly bounded, the resultant
>> function is never going to fly off in all directions, and will not
>> become totally meaningless unless that argument uncertainty is in
>> the order of PI.[/color]
>
> 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]

Consider how you compute 1e30 DIV 3 to any desired accuracy. The
result depends on the accuracy of 3 alone if you assume that 1e30
is exact. Similarly, consider the results on division of 1e30 by
0.33333. (exactly 5 3s)

--
A: Because it fouls the order in which people normally read text.
Q: Why is top-posting such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

[Grumble]

Re: Sine code for ANSI C

Tim Prince wrote:
[color=blue]
> Grumble wrote:
>[color=green]
>> AMD64 supports x87. Thus FSIN is available. What did you mean by
>> "not really true of AMD64 either" ?[/color]
>
> Possibly referring to compilers complying with ABIs disallowing
> x87, or taking advantage of higher performance of SSE parallel
> libraries. Use of fsin on IA64 is extremely unlikely, even though
> it's still there.[/color]

I assume you meant AMD64, not IA64?

Where did you read that the AMD64 ABI disallowed x87?

http://www.x86-64.org/abi.pdf

[Sam Dennis]

Re: Sine code for ANSI C

P.J. Plauger wrote:[color=blue]
> [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).

I think that it'd be particularly ironic to combine this with unlimited
precision for floating-point constants.

--
++acr@,ka"