Sine code for ANSI C

[John Cochran]

Re: Sine code for ANSI C

In article <xLdpc.72825$sK 3.24627@nwrddc0 3.gnilink.net>,
P.J. Plauger <pjp@dinkumware .com> wrote:

SNIP....
[color=blue]
>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]

It seems to me that a reasonable cutoff point would be where
the difference between consecutive floating point numbers is greater
than two pi. At that point you can't even determine the *sign* of the
correct answer, yet alone determine any value that is justifiable.
The only thing that you can justify is a claim that the answer lies
somewhere between -1.0 and 1.0

[P.J. Plauger]

Re: Sine code for ANSI C

"John Cochran" <jdc@smof.fiawo l.org> wrote in message
news:c845k0$bm8 $1@smof.fiawol. org...
[color=blue]
> In article <xLdpc.72825$sK 3.24627@nwrddc0 3.gnilink.net>,
> P.J. Plauger <pjp@dinkumware .com> wrote:
>
> SNIP....
>[color=green]
> >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]
>
> It seems to me that a reasonable cutoff point would be where
> the difference between consecutive floating point numbers is greater
> than two pi. At that point you can't even determine the *sign* of the
> correct answer, yet alone determine any value that is justifiable.
> The only thing that you can justify is a claim that the answer lies
> somewhere between -1.0 and 1.0[/color]

Yes, that's a reasonable cutoff point, on the face of it. Just don't
look too close. You're falling prey to the same error in logic that
traps most people who first study this problem -- assuming that there
must be some intrinsic error, say 1/2 ulp, in the argument. If we're
going to apply that criterion to the library uniformly, then rounding
1e38 is equally suspect. (How happy would you be if round occasionally
produced a fatal error? Particularly when the answer is obvious and
easy to compute.)

But if you assume the argument is exact, as *all* library functions
really must do, the your statement is incorrect. There *is* a well
defined angle corresponding to *every* finite floating-point argument.
You (or I, to be specific) may not like the amount of work required
to compute it accurately, but the value is known and well defined.

If, OTOH, you want to let the library vendors off the hook whenever
a function is hard to compute, I've got a little list. And it doesn't
stop at sin/cos/tan.

And if, OTOOH, you want to let the library vendors off the hook
whenever a typical programmer probably doesn't know what he's doing,
*boy* do I have a list.

The trick with standards, as with library design, is to have a
rational framework that's uniformly applied. The sine function
may illustrate some of the nastiest consequences of carrying the
current framework to its logical conclusion, but I assure you
that there are plenty of other monsters lurking out there. Any
change you propose to the framework just gives you a different
set to battle.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[Christian Bau]

Re: Sine code for ANSI C

In article <c845k0$bm8$1@s mof.fiawol.org> ,
jdc@smof.fiawol .org (John Cochran) wrote:
[color=blue]
> In article <xLdpc.72825$sK 3.24627@nwrddc0 3.gnilink.net>,
> P.J. Plauger <pjp@dinkumware .com> wrote:
>
> SNIP....
>[color=green]
> >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]
>
> It seems to me that a reasonable cutoff point would be where
> the difference between consecutive floating point numbers is greater
> than two pi. At that point you can't even determine the *sign* of the
> correct answer, yet alone determine any value that is justifiable.
> The only thing that you can justify is a claim that the answer lies
> somewhere between -1.0 and 1.0[/color]

Well, you actually _can_ find the correct answer quite well. A value of
type double represents a single real number. Of course we all know that
if I assign x = a + b; then usually x is _not_ equal to the mathematical
sum of a and b, and given only x I might not draw any useful conclusions
about sin (a + b). However, sin (x) can still be calculated quite well.

[James Kanze]

Re: Sine code for ANSI C

"P.J. Plauger" <pjp@dinkumware .com> writes:

|> 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.

Maybe I'm just being naïve, but I always thought that it was a
necessary quality of a correct implementation that it give correct
results for all legal input. And that it wasn't the job of library
implementers to decide what was or was not reasonable input for my
application.

--
James Kanze
Conseils en informatique orientée objet/
Beratung in objektorientier ter Datenverarbeitu ng
9 place Sémard, 78210 St.-Cyr-l'École, France +33 (0)1 30 23 00 34

[Joe Wright]

Re: Sine code for ANSI C

P.J. Plauger wrote:
[color=blue]
> "Dr Chaos" <mbkennelSPAMBE GONE@NOSPAMyaho o.com> wrote in message
> news:slrncaagm5 .gu0.mbkennelSP AMBEGONE@lyapun ov.ucsd.edu...
>
>[color=green]
>>On Fri, 14 May 2004 21:23:51 GMT, P.J. Plauger <pjp@dinkumware .com> wrote:
>>[color=darkred]
>>>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=green][color=darkred]
>>>>>>Conside r 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.
>>>
>>>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=green]
>>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?
>[/color]

P.J.
I tend to come down on your side on these things (except casting
malloc, maybe). I am not a mathematician but am very interested in
your take on the following floating point issues..

1. Accuracy vs Precision. #define Pi 3.1416 is precise to five
digits and accurate within its precision. If I do something like..

double Pi2 = Pi * 2.0;

...the constant 2.0 is accurate and precise to 16 digits. The result
of the multiplication is accurate to only five digits while it is
precise to 16. Does this make sense?

2. Large Angles. The circle is 360 degrees or '2 pi radians'. Why is
something like..

double r = 52147.3, s;
s = sin(fmod(r,2*PI ));

...not the solution for large angle argument reduction?

Keep up the good work.
--
Joe Wright mailto:joewwrig ht@comcast.net
"Everything should be made as simple as possible, but not simpler."
--- Albert Einstein ---

[CBFalconer]

Re: Sine code for ANSI C

Joe Wright wrote:[color=blue]
>[/color]
.... snip ...[color=blue]
>
> 2. Large Angles. The circle is 360 degrees or '2 pi radians'.
> Why is something like..
>
> double r = 52147.3, s;
> s = sin(fmod(r,2*PI ));
>
> ..not the solution for large angle argument reduction?[/color]

That depends highly on how you compute the fmod. Say you compute
r/(2*PI), truncate it to an integer, multiply by (2*PI), and
subtract that from r. Now you have the difference of two
comparable magnitudes, with attendant loss of significant bits.

Compare with methods of computing sigma f(n) for n = 1 ....
infinity. If you start with n=1 (using the normally available
floating point system) you will end up with something quite
divergent from the true answer, regardless of whether the series
actually converges or not. A reasonably accurate computation
requires starting with the smallest terms.

--
"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.

[Dik T. Winter]

Re: Sine code for ANSI C

In article <c82q2b$gfu$6@s unnews.cern.ch> Dan.Pop@cern.ch (Dan Pop) writes:[color=blue]
> In <HxpE44.Myq@cwi .nl> "Dik T. Winter" <Dik.Winter@cwi .nl> writes:[/color]
....[color=blue][color=green][color=darkred]
> > > 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]

I would have thought that was simple. Yes, I agree with that paragraph.
[color=blue][color=green][color=darkred]
> > > > Not in mathematical applications, where the argument to the sine
> > > > function can very well be exact.
> > >
> > > 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]

I'm so sorry.
[color=blue][color=green][color=darkred]
> > > "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.[/color]

You are forgetting a few fields: numerical algebra, computational number
theory, statistics... But in some cases (especially computational number
theory) final results can be exact while intermediate results are not.
For instance in a project I was involved in, we have shown that the first
1,200,000,000 non-trivial zeros of the Rieman zeta function have real
part 1/2. However, none of the calculated zeros was exact. (Moreover,
the calculation involved the sine function. Moreover, we had to get
reasonable precision, as it involved separating places where the sign
of the function changes. %)
---
% If you are interested, there are well known formula's that indicate in
which region the n-th non-trivial zero resides. However, these regions
overlap. But using this you can set up algorithms that will locate
groups of zero's. The problem is that these zero's will get arbitrarily
close to each other, so using floating point on the machine we did the
runs on (a CDC Cyber 205) the best separation we could get was 10**(-13).
This was not enough, so sometimes we had to resort to double precision.
However, the argument was *exact*. A precise floating-point (hence
rational) number.
--
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

"Joe Wright" <joewwright@com cast.net> wrote in message
news:Q5Sdnf1q2Z JuBzrdRVn-gg@comcast.com. ..
[color=blue]
> I tend to come down on your side on these things (except casting
> malloc, maybe). I am not a mathematician but am very interested in
> your take on the following floating point issues..
>
> 1. Accuracy vs Precision. #define Pi 3.1416 is precise to five
> digits and accurate within its precision. If I do something like..
>
> double Pi2 = Pi * 2.0;
>
> ..the constant 2.0 is accurate and precise to 16 digits. The result
> of the multiplication is accurate to only five digits while it is
> precise to 16. Does this make sense?[/color]

Yes.
[color=blue]
> 2. Large Angles. The circle is 360 degrees or '2 pi radians'. Why is
> something like..
>
> double r = 52147.3, s;
> s = sin(fmod(r,2*PI ));
>
> ..not the solution for large angle argument reduction?[/color]

People once thought it was, or should be. Indeed, that was one of the
arguments for adding the rather finicky fmod to IEEE floating point
and eventually the C Standard. But if you think about it hard enough
and long enough -- took me an afternoon and several sheets of paper --
you realize that it doesn't cut it. You effectively have to keep
subtracting 2*pi from your argument r until it's less than 2*pi.
fmod does this by subtracting the various multiples 2*pi*2^n. If
*any one* of them does not have nearly 16 good fraction digits, as
well as all the digits it needs to the left of the decimal point, it's
going to mess up the whole set of subtractions. So if you want to
reduce numbers as large as 10^38, you have to represent pi to about
log10(10^(38+16 )) or 54 digits. For 113-bit IEEE long double, you
need well over 4000 digits.

We've developed an arbitrary-precision package that represents
numbers as arrays of floating-point values, each of which uses only
half its fraction bits. So we can do adjustable precision argument
reduction fairly rapidly. Still takes a lot of storage to represent
the worst case, and a bit more logic than I wish we had to use, but
it does the job. Not only that, we need the same sort of thing
to handle several other difficult math functions, though with
nowhere near as much precision, of course. So it's not like we
indulge in heroics just for sin/cos/tan.
[color=blue]
> Keep up the good work.[/color]

Thanks.

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]
>
> We've developed an arbitrary-precision package that represents
> numbers as arrays of floating-point values, each of which uses only
> half its fraction bits. So we can do adjustable precision argument
> reduction fairly rapidly. Still takes a lot of storage to represent
> the worst case, and a bit more logic than I wish we had to use, but
> it does the job. Not only that, we need the same sort of thing
> to handle several other difficult math functions, though with
> nowhere near as much precision, of course. So it's not like we
> indulge in heroics just for sin/cos/tan.[/color]

It seems to me that the reduction could be fairly rapid if an
estimate is formed by normal division, break that up into single
bit binary portions (so that multiples of PI do not add non-zero
significant bits) to do the actual reduction. Not worked out at
all, just the glimmer of a method. The idea is to remove the
leftmost digits of the original argument first.

--
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!

[Dan Pop]

Re: Sine code for ANSI C

In <christian.ba u-3E2624.19493215 052004@slb-newsm1.svr.pol. co.uk> Christian Bau <christian.bau@ cbau.freeserve. co.uk> writes:
[color=blue]
>Well, you actually _can_ find the correct answer quite well. A value of
>type double represents a single real number.[/color]

But, when used in a real number context (as opposed to an integer number
context -- floating point can be used in both contexts) it stands for a
whole subset of the real numbers set. The real value exactly represented
is no more relevant than any other value from that set.
[color=blue]
>Of course we all know that
>if I assign x = a + b; then usually x is _not_ equal to the mathematical
>sum of a and b, and given only x I might not draw any useful conclusions
>about sin (a + b). However, sin (x) can still be calculated quite well.[/color]

The point is not whether it can be calculated, but rather how much
precision should the calculation produce? Does it makes sense to compute
sin(DBL_MAX) with 53-bit precision, ignoring the fact that DBL_MAX
stands for an interval so large as to make this function call completely
devoid of any meaning?

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

[P.J. Plauger]

Re: Sine code for ANSI C

"Dan Pop" <Dan.Pop@cern.c h> wrote in message
news:c8ahag$qfi $17@sunnews.cer n.ch...[color=blue]
> In <christian.ba u-3E2624.19493215 052004@slb-newsm1.svr.pol. co.uk>[/color]
Christian Bau <christian.bau@ cbau.freeserve. co.uk> writes:
[color=blue][color=green]
> >Well, you actually _can_ find the correct answer quite well. A value of
> >type double represents a single real number.[/color]
>
> But, when used in a real number context (as opposed to an integer number
> context -- floating point can be used in both contexts) it stands for a
> whole subset of the real numbers set. The real value exactly represented
> is no more relevant than any other value from that set.[/color]

That's just one interpretation of floating-point representation. It is
akin to the oft-repeated presumption in this thread that every value
carries a +/- 1/2 ulp error with it. But that's not always the case.
Many a floating-point calculation is *exact*. So the real value exactly
represented is way more relevant than any other value from the set.
In particular, it is the value that library functions should assume
when deciding what function value to return. Do that and you will
please a large number of your customers. Don't do it and you won't
even please those who really believe the 1/2 ulp fuzz. (They should
be doing their own error propagation analysis no matter what, and
they'll almost always benefit from a result in the center of the range.)
[color=blue][color=green]
> >Of course we all know that
> >if I assign x = a + b; then usually x is _not_ equal to the mathematical
> >sum of a and b, and given only x I might not draw any useful conclusions
> >about sin (a + b). However, sin (x) can still be calculated quite well.[/color]
>
> The point is not whether it can be calculated, but rather how much
> precision should the calculation produce? Does it makes sense to compute
> sin(DBL_MAX) with 53-bit precision, ignoring the fact that DBL_MAX
> stands for an interval so large as to make this function call completely
> devoid of any meaning?[/color]

You have to ignore the possibility that DBL_MAX stands for a large range
because *it might not.* Granted, it probably does have that much fuzz
practically all the time; but you, the library vendor, don't know that.
And as we've discussed at length here, it's hard to say just where along
the road to possible perdition it's okay to stop trying. Once you tool
up for doing the job right along a good part of that road, it ain't that
much more work to go all the way. Then you don't have to worry about
some customer coming along and saying, "Hey, why did you give up on
computing that value properly? I happen to know what I'm doing." And
that customer might even be right.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[CBFalconer]

Re: Sine code for ANSI C

Dan Pop wrote:[color=blue]
>[/color]
.... snip ...[color=blue]
>
> The point is not whether it can be calculated, but rather how much
> precision should the calculation produce? Does it makes sense to
> compute sin(DBL_MAX) with 53-bit precision, ignoring the fact that
> DBL_MAX stands for an interval so large as to make this function
> call completely devoid of any meaning?[/color]

PJPs point is that that interval is the range of values specified,
with the actual value being the most likely value in the range.
The shape of the probability curve in the interval is known only
to the end user. It may be an impulse function precisely at the
specified value.

If I ask for the value of sin(12345678*PI ) I can confidently
assert that the value is 0. Using that value will never cause an
unexplained fault. The following illustrates the problem:

[1] c:\c\junk>a
1 3.1415926535897 931 0.0000000000000 001
2 6.2831853071795 862 -0.0000000000000 002
4 12.566370614359 1725 -0.0000000000000 005
8 25.132741228718 3449 -0.0000000000000 010
16 50.265482457436 6899 -0.0000000000000 020
32 100.53096491487 33797 -0.0000000000000 039
64 201.06192982974 67594 -0.0000000000000 078
128 402.12385965949 35188 -0.0000000000000 157
256 804.24771931898 70377 -0.0000000000000 313
512 1608.4954386379 740754 -0.0000000000000 627
1024 3216.9908772759 481508 -0.0000000000001 254
2048 6433.9817545518 963016 -0.0000000000002 508
4096 12867.963509103 7926031 -0.0000000000005 016
8192 25735.927018207 5852063 -0.0000000000010 032
16384 51471.854036415 1704125 -0.0000000000020 064
32768 102943.70807283 03408250 -0.0000000000040 128
65536 205887.41614566 06816500 -0.0000000000080 256
131072 411774.83229132 13633001 -0.0000000000160 512
262144 823549.66458264 27266002 -0.0000000000321 023
524288 1647099.3291652 854532003 -0.0000000000642 046
1048576 3294198.6583305 709064007 -0.0000000001284 093
2097152 6588397.3166611 418128014 -0.0000000002568 186
4194304 13176794.633322 2836256027 -0.0000000005136 371
8388608 26353589.266644 5672512054 -0.0000000010272 743
16777216 52707178.533289 1345024109 -0.0000000020545 485
33554432 105414357.06657 82690048218 -0.0000000041090 971
67108864 210828714.13315 65380096436 -0.0000000082181 941
134217728 421657428.26631 30760192871 -0.0000000164363 883
268435456 843314856.53262 61520385742 -0.0000000328727 765
536870912 1686629713.0652 523040771484 -0.0000000657455 530
1073741824 3373259426.1305 046081542969 -0.0000001314911 060

[1] c:\c\junk>cat sines.c
#include <stdio.h>
#include <limits.h>
#include <math.h>
int main(void)
{
double PI = 4.0 * atan(1.0);
unsigned int i;

for (i = 1; i < INT_MAX; i += i) {
printf("%11d %28.16f %20.16f\n", i, PI * i, sin(PI * i));
}
return 0;
}

--
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!

[Christian Bau]

Re: Sine code for ANSI C

In article <%27qc.135575$G _.78424@nwrddc0 2.gnilink.net>,
"P.J. Plauger" <pjp@dinkumware .com> wrote:
[color=blue]
> You have to ignore the possibility that DBL_MAX stands for a large range
> because *it might not.* Granted, it probably does have that much fuzz
> practically all the time; but you, the library vendor, don't know that.
> And as we've discussed at length here, it's hard to say just where along
> the road to possible perdition it's okay to stop trying. Once you tool
> up for doing the job right along a good part of that road, it ain't that
> much more work to go all the way. Then you don't have to worry about
> some customer coming along and saying, "Hey, why did you give up on
> computing that value properly? I happen to know what I'm doing." And
> that customer might even be right.[/color]

Somehow your mentioning of sin (DBL_MAX) sparked a little idea that
might help you calculating sin (x) for arbitrarily large values quite
quickly, especially useful if you use long double with 15 exponent bits
and an enormous range:

Instead of storing huge numbers of bits of 1/2pi and doing huge
precision argument reduction, just store the values of sin (65536 ^ N)
and cos (65536 ^ N) with lets say 20 or 24 extra bits of precision.
Every large floating-point number can be written as x = +/- k * 2^n
where k is for example a 53 bit or 64 bit integer. So take a value of x,
write it as the sum of very few small multiples of (65536 ^ N) which
should be quite simple, then use some simple theorems for calculating
sin (a+b) when sin and cos of a and b are known. That should work quite
quickly for arbitrary large x.

[Joe Wright]

Re: Sine code for ANSI C

CBFalconer wrote:
[color=blue]
> "P.J. Plauger" wrote:
>
> ... snip ...[/color]

I got really interested in P.J.'s argument reduction to sin() et.al.
and while constructing something to test with I have found something
very strange. A very simple program..

#include <stdio.h>
#include <math.h>

#define R PI*2/360 /* Radians per Degree */

int main(void) {
double d, r = R;
int i, lin;
for (lin = 0, i = 30; i < 6000000; i += 360000) {
d = i * R;
printf("\n%2d. Deg %8d, Rad %.16e Sin %+.16e", ++lin, i, d,
sin(d));
}
puts("");
for (lin = 0, i = 30; i < 6000000; i += 360000) {
d = i * r;
printf("\n%2d. Deg %8d, Rad %.16e Sin %+.16e", ++lin, i, d,
sin(d));
}
puts("");
return 0;
}

Sorry about the line wraps.

I use the DJGPP versions of gcc.

The point of the exercise was to simulate fairly large integral
(degree) angles, convert the (int) degrees to (double) radians and
present the result to sin() and see what we get. Because I know the
sine of 30 degrees is 0.5 and I suppose the sine 390 and 750 degrees
(and so on) is the same.

The program consists of two virtually identical loops. The
difference is only the assignment of d.

The loops produce 17 lines of output, identical except (in My case)
for line 14. In that case, the calculation of d yields a different
result, off by one in the least significant bit of the mantissa of
d. Passing these two to sin() result in a difference of 18 lsb bits
of the result's mantissa.

If you see the same difference, can you explain it?
--
Joe Wright mailto:joewwrig ht@comcast.net
"Everything should be made as simple as possible, but not simpler."
--- Albert Einstein ---

[CBFalconer]

Re: Sine code for ANSI C

Christian Bau wrote:[color=blue]
>[/color]
.... snip ...[color=blue]
>
> Somehow your mentioning of sin (DBL_MAX) sparked a little idea
> that might help you calculating sin (x) for arbitrarily large
> values quite quickly, especially useful if you use long double
> with 15 exponent bits and an enormous range:
>
> Instead of storing huge numbers of bits of 1/2pi and doing huge
> precision argument reduction, just store the values of
> sin (65536 ^ N) and cos (65536 ^ N) with lets say 20 or 24
> extra bits of precision. ... snip ...[/color]

That isn't necessary. Assuming the argument is precise (i.e.
extended by zero bits on the right) all we need to compute (arg
mod PI) to the precision of PI is a couple of guard bits.

--
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!