Sine code for ANSI C

[suri]

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

 

[suri]

No i knew that faq.
my question is where is the sine fn imlementation?
since u say its a native instruction set means that there is no code in
ansi C for implementing sine.?
well libm.a must surely have it. but where are the C files?

 

[suri]

I do know the series expansion of sine i was just interested to know how
its implemented in the ansi C library. like how many terms of the
infinite series are included.

I have linux and use glibc. so i could find the file in the path u mentioned

No i knew that faq.
my question is where is the sine fn imlementation?
since u say its a native instruction set means that there is no code in
ansi C for implementing sine.?

Some platforms have hardware instructions that compute sin() and the
compiler will emit them, bypassing the libm library's implementation. This
is pretty much true for ia32 and ia64.

well libm.a must surely have it. but where are the C files?

For FreeBSD, check out /usr/src/lib/libm/common/{trig.h,sincos.c} [nasty
looking code, but it works fast.] You're almost better off grabbing a
calculus book and having a look at a Taylor series expansion if you truly
want to understand the math. Or look at "Numerical Recipes in C".

 

[suri]

im sorry i meant i *could not* find the file

I do know the series expansion of sine i was just interested to know how
its implemented in the ansi C library. like how many terms of the
infinite series are included.

I have linux and use glibc. so i could find the file in the path u
mentioned

No i knew that faq.
my question is where is the sine fn imlementation?
since u say its a native instruction set means that there is no code in
ansi C for implementing sine.?

Some platforms have hardware instructions that compute sin() and the
compiler will emit them, bypassing the libm library's implementation.
This
is pretty much true for ia32 and ia64.

well libm.a must surely have it. but where are the C files?

For FreeBSD, check out /usr/src/lib/libm/common/{trig.h,sincos.c} [nasty
looking code, but it works fast.] You're almost better off grabbing a
calculus book and having a look at a Taylor series expansion if you truly
want to understand the math. Or look at "Numerical Recipes in C".

 

[Nudge]

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

http://www.eskimo.com/~scs/C-faq/q14.3.html

Do you really want an implementation for sin?

Apparently, you are running an Intel-compatible processor which
comes with a sin instruction implemented in hardware.

 

[-wombat-]

No i knew that faq.
my question is where is the sine fn imlementation?
since u say its a native instruction set means that there is no code in
ansi C for implementing sine.?

Some platforms have hardware instructions that compute sin() and the
compiler will emit them, bypassing the libm library's implementation. This
is pretty much true for ia32 and ia64.

well libm.a must surely have it. but where are the C files?

For FreeBSD, check out /usr/src/lib/libm/common/{trig.h,sincos.c} [nasty
looking code, but it works fast.] You're almost better off grabbing a
calculus book and having a look at a Taylor series expansion if you truly
want to understand the math. Or look at "Numerical Recipes in C".

 

[osmium]

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

Here is one way, but not the most advanced way. It is fairly suitable for a
float (precision wise) but not for a double.

#include <iostream.h>
#include <iomanip.h>
#include <math.h>

// Source AMS 55, eqn 4.3.97. Handbook of Mathematical Functions, Pub by
U.S. Dept of Commerce
float sinx(float x)
{
static const float a[] =
{-.1666666664,.0083333315,-.0001984090,.0000027526,-.0000000239};
float xsq = x*x;
float temp = x*(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq
+a[3]*xsq*xsq*xsq*xsq
+ a[4]*xsq*xsq*xsq*xsq*xsq);
return temp;
}
// ------------------
void test()
{
float x;
while(1)
{
cin >> x;
if(x<0. || x > (3.1416/2) )
{
cout << "Argument to sinx must be in range 0>= x <= pi/2 \n";
continue;
}
cout << sinx(x) << setw(12) << (float)sin(x) << endl;
}

The test code, but not the sinx() code, is in C++ which I suppose will lead
to some hissy fits since this is a C group. It is also a pre-standard
version of C++ so I suppose there will be some more hissy fits from that
quarter. So be it.

 

[Jack Klein]

comp.lang.c:

First, please don't top-post. Material you add in a reply goes AFTER
material you are quoting. Top-posting makes technical discussions
hard to follow and is considered rude in comp.lang.c.

I do know the series expansion of sine i was just interested to know how
its implemented in the ansi C library. like how many terms of the
infinite series are included.

There is no such thing as "the" ANSI C library. The C standard
defines the functions that must be provided, their interface, and
their results when used as defined. It does not specify an
implementation at all.

The standard library functions supplied with an implementation do not
even have to be written in C.

 

[Grumble]

Some platforms have hardware instructions that compute sin() and
the compiler will emit them, bypassing the libm library's
implementation. This is pretty much true for ia32 and ia64.

<OT>

As far as I can tell, IA-64 has no such instruction

 

[Grumble]

my question is where is the sine fn imlementation?
since u say its a native instruction set means that there
is no code in ansi C for implementing sine?
well libm.a must surely have it. but where are the C files?

Why don't you ask in a forum dedicated to the GNU libc?

You might try:

comp.os.linux.questions
comp.os.linux.development.apps
http://sources.redhat.com/glibc/

(Please do not top-post.)

 

[Mark McIntyre]

since u say its a native instruction set means that there is no code in
ansi C for implementing sine.?

Quite possibly.

well libm.a must surely have it. but where are the C files?

If its native to the processor, there IS no code to do it. Whether some
(offtopic) library has an alternate implementation is offtopic here.

But if you want code, then write it yourself. Sin(x) is calculable to a
good approximation using a power series.

 

[P.J. Plauger]

Quite possibly.


If its native to the processor, there IS no code to do it. Whether some
(offtopic) library has an alternate implementation is offtopic here.

But if you want code, then write it yourself. Sin(x) is calculable to a
good approximation using a power series.

A *good* sine function is one of the hardest math functions to
write, believe it or not. For x < |pi| / 4, it it indeed easy
to approximate with a polynomial in x^2 -- a truncated power
series will do, but you can save a term or two by fiddling the
coefficients. But doing proper argument reduction is an open
ended exercise in frustration. Just reducing the argument modulo
2*pi quickly accumulates errors unless you do arithmetic to
many extra bits of precision.

Some processors have instructions that do *parts* of the job
of computing the sine, like evaluating it in the easy region
described above. Some even have functions that appear to do
the argument reduction well for you, but they don't always
deliver.

See my book, The Standard C Library, for one implementation of
sin(x) in Standard C. It's only moderately naive.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[CBFalconer]

.... snip ...
coefficients. But doing proper argument reduction is an open
ended exercise in frustration. Just reducing the argument modulo
2*pi quickly accumulates errors unless you do arithmetic to
many extra bits of precision.

And that problem is inherent. Adding precision bits for the
reduction will not help, because the input value doesn't have
them. It is the old problem of differences of similar sized
quantities.

 

[osmium]

And that problem is inherent. Adding precision bits for the
reduction will not help, because the input value doesn't have
them. It is the old problem of differences of similar sized
quantities.

Huh? If I want the phase of an oscillator after 50,000 radians are you
saying that is not computable? Please elaborate.

There was a thread hereabouts many months ago on this very subject and AFAIK
no one suggested that it was not computable, it just couldn't be done with
doubles. And I see no inherent problems.

 

[Rich Gibbs]

-wombat- said the following, on 05/02/04 18:43:

No i knew that faq.
my question is where is the sine fn imlementation?
since u say its a native instruction set means that there is no code in
ansi C for implementing sine.?

[snip]

For FreeBSD, check out /usr/src/lib/libm/common/{trig.h,sincos.c} [nasty
looking code, but it works fast.] You're almost better off grabbing a
calculus book and having a look at a Taylor series expansion if you truly
want to understand the math. Or look at "Numerical Recipes in C".

There are (at least) two possible questions that might lie behind the
OP's request. One, which is more or less on-topic, is "How might one
implement the sin() function of the C library in standard C?"

The other is something like, "What is the best numerical method for
calculating the sine function?" The answer to THAT question has
something to do with C, if that is the implementation language, but has
more to do with numerical analysis. Wombat has suggested _Numerical
Recipes in C_, which is good. The elementary N/A text I used way back
when was F.S. Acton's _Numerical Methods That [Usually] Work_. (BTW,
calculating function like this by straightforward application of the
"textbook" power series is often not the best approach.)

I'd also recommend P.J. Plauger's excellent book, _The Standard C
Library_, which includes consideration of both questions.

 

[P.J. Plauger]

Huh? If I want the phase of an oscillator after 50,000 radians are you
saying that is not computable? Please elaborate.

There was a thread hereabouts many months ago on this very subject and AFAIK
no one suggested that it was not computable, it just couldn't be done with
doubles. And I see no inherent problems.

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

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

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

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

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

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[CBFalconer]

Huh? If I want the phase of an oscillator after 50,000 radians are
you saying that is not computable? Please elaborate.

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

To illustrate, take your handy calculator and enter:

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

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

 

[P.J. Plauger]

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.

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

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[Dik T. Winter]

> The other is something like, "What is the best numerical method for
> calculating the sine function?" The answer to THAT question has
> something to do with C, if that is the implementation language, but has
> more to do with numerical analysis.

It has nothing to do with the C language. The newsgroup:
sci.math.num-analysis
is more appropriate. They provide algorithms, it is only the last question
how to implement the algorithm in C. When designing these algorithms you
have quite a few contradicting objectives. (Overall relative precision,
monotony, and whatever.) Do you wish a good algorithm where sin(x) > 1
could be true? (Cody-Waite, of sixties fame provided an algorithm where
that was possible.) The algorithm could still fit in your requirements:
speed and relative precision.

 

[Dik T. Winter]

> I've gone to both extremes over the past several decades.
> Our latest math library, still in internal development,
> can get exact function values for *all* argument values.
> It uses multi-precision argument reduction that can gust
> up to over 4,000 bits [sic].

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

It has been my opinion for a very long time that he real sine function
is not so very important. The most important in numerical mathematics
is sin2pi(x) which calculates sin(2.pi.x), and which allows easy and
exact range reduction. Alas, that function is not in C.