473,836 Members | 2,332 Online
Bytes | Software Development & Data Engineering Community
+ Post

Home Posts Topics Members FAQ

Sine code for ANSI C

Hello
I downloaded glibc and tried looking for the code that implements the
sine function
i couldnt find the file.
i went to the math directory and found math.h.. i guess that needs to be
included for the sine function. but which .c file implements the
sine/cosine and other trig fns
thanks
Nov 14 '05
143 8119
"Dan Pop" <Da*****@cern.c h> wrote in message
news:c8******** ***@sunnews.cer n.ch...
In <ch************ *************** ******@slb-newsm1.svr.pol. co.uk>

Christian Bau <ch***********@ cbau.freeserve. co.uk> writes:
Well, you actually _can_ find the correct answer quite well. A value of
type double represents a single real number.


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.


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


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?


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
Nov 14 '05 #101
Dan Pop wrote:
.... snip ...
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?


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 (cb********@yah oo.com) (cb********@wor ldnet.att.net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net> USE worldnet address!
Nov 14 '05 #102
In article <%2************ *******@nwrddc0 2.gnilink.net>,
"P.J. Plauger" <pj*@dinkumware .com> wrote:
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.


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.
Nov 14 '05 #103
CBFalconer wrote:
"P.J. Plauger" wrote:

... snip ...


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:jo****** **@comcast.net
"Everything should be made as simple as possible, but not simpler."
--- Albert Einstein ---
Nov 14 '05 #104
Christian Bau wrote:
.... snip ...
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 ...


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 (cb********@yah oo.com) (cb********@wor ldnet.att.net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net> USE worldnet address!
Nov 14 '05 #105
In article <40************ ***@yahoo.com>,
CBFalconer <cb********@yah oo.com> wrote:
Christian Bau wrote:

... snip ...

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


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.


When you write PI you mean a floating point approximation of the
mathematical constant pi? For large x, reduction modulo an approximation
to pi will give a result that is very much different from reduction
modulo pi. That is why some libraries store pi or 1/pi with a precision
of several thousand bits.
Nov 14 '05 #106
Joe Wright <jo********@com cast.net> wrote:
CBFalconer wrote:
"P.J. Plauger" wrote:

... snip ...
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 */


<snip>
If you see the same difference, can you explain it?


If you actually define PI somewhere and parenthesize the
replacement text for R, the code not only compiles, but
all the strangeness should be gone.

Regards
--
Irrwahn Grausewitz (ir*******@free net.de)
welcome to clc: http://www.ungerhu.com/jxh/clc.welcome.txt
clc faq-list : http://www.faqs.org/faqs/C-faq/faq/
clc OT guide : http://benpfaff.org/writings/clc/off-topic.html
Nov 14 '05 #107
Christian Bau wrote:
CBFalconer <cb********@yah oo.com> wrote:
Christian Bau wrote:

... snip ...

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


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.


When you write PI you mean a floating point approximation of the
mathematical constant pi? For large x, reduction modulo an
approximation to pi will give a result that is very much different
from reduction modulo pi. That is why some libraries store pi or
1/pi with a precision of several thousand bits.


You're right. I was thinking of the case where we start with a
number of the form "arg = 2 * PI * rotations + delta", where we
can recover delta provided we use the same value for PI and the
number has not lost precision in the first place.

--
Chuck F (cb********@yah oo.com) (cb********@wor ldnet.att.net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net> USE worldnet address!
Nov 14 '05 #108
In <%2************ *******@nwrddc0 2.gnilink.net> "P.J. Plauger" <pj*@dinkumware .com> writes:
"Dan Pop" <Da*****@cern.c h> wrote in message
news:c8******* ****@sunnews.ce rn.ch...
In <ch************ *************** ******@slb-newsm1.svr.pol. co.uk>

Christian Bau <ch***********@ cbau.freeserve. co.uk> writes:
>Well, you actually _can_ find the correct answer quite well. A value of
>type double represents a single real number.


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.


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


When I asked you for concrete examples, you provided exactly zilch.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Da*****@ifh.de
Nov 14 '05 #109
In <40************ ***@yahoo.com> CBFalconer <cb********@yah oo.com> writes:
Dan Pop wrote:

... snip ...

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?


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.


Or it may not. If I deserve 5 digits of precision, why should the
implementation waste precious CPU resources to give me 16 digits of
*bogus* precision?

To deserve 16 digits of precision, it is *my* job to call sin() with
an argument of the right magnitude. It's as simple as that.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Da*****@ifh.de
Nov 14 '05 #110

This thread has been closed and replies have been disabled. Please start a new discussion.

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.