Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > Cosine algo, was Re: sqrt algo

Reply
Thread Tools

Cosine algo, was Re: sqrt algo

 
 
pete
Guest
Posts: n/a
 
      10-15-2005
pete wrote:

> .. and now for a float sqrt algorithm:


I just wrote a cosine function called c_os, for no special reason.
It calls no functions that I haven't also written.

On my machine, it seems to work better than
my MS compiler's standard library cos function.
They don't get the zeros right. I do.

c_os(-123 * pi / 6) is 0.000000e+000
c_os(-129 * pi / 6) is 0.000000e+000

cos(-123 * pi / 6) is 7.839514e-015
cos(-129 * pi / 6) is -4.409261e-015


/* BEGIN c_os.c output */

const double pi = 4 * atan(1);

c_os(-120 * pi / 6) is 1.000000e+000
c_os(-121 * pi / 6) is 8.660254e-001
c_os(-122 * pi / 6) is 5.000000e-001
c_os(-123 * pi / 6) is 0.000000e+000
c_os(-124 * pi / 6) is -5.000000e-001
c_os(-125 * pi / 6) is -8.660254e-001
c_os(-126 * pi / 6) is -1.000000e+000
c_os(-127 * pi / 6) is -8.660254e-001
c_os(-128 * pi / 6) is -5.000000e-001
c_os(-129 * pi / 6) is 0.000000e+000
c_os(-130 * pi / 6) is 5.000000e-001
c_os(-131 * pi / 6) is 8.660254e-001

cos(-120 * pi / 6) is 1.000000e+000
cos(-121 * pi / 6) is 8.660254e-001
cos(-122 * pi / 6) is 5.000000e-001
cos(-123 * pi / 6) is 7.839514e-015
cos(-124 * pi / 6) is -5.000000e-001
cos(-125 * pi / 6) is -8.660254e-001
cos(-126 * pi / 6) is -1.000000e+000
cos(-127 * pi / 6) is -8.660254e-001
cos(-128 * pi / 6) is -5.000000e-001
cos(-129 * pi / 6) is -4.409261e-015
cos(-130 * pi / 6) is 5.000000e-001
cos(-131 * pi / 6) is 8.660254e-001

/* END c_os.c output */

/* BEGIN c_os.c */

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

double c_os(double x);
double p_i(void);
double sq_rt(double x);

int main(void)
{
const double pi = 4 * atan(1);

puts("/* BEGIN c_os.c output */\n");
puts("const double pi = 4 * atan(1);\n");
printf("c_os(-120 * pi / 6) is %e\n", c_os(-120 * pi / 6));
printf("c_os(-121 * pi / 6) is %e\n", c_os(-121 * pi / 6));
printf("c_os(-122 * pi / 6) is %e\n", c_os(-122 * pi / 6));
printf("c_os(-123 * pi / 6) is %e\n", c_os(-123 * pi / 6));
printf("c_os(-124 * pi / 6) is %e\n", c_os(-124 * pi / 6));
printf("c_os(-125 * pi / 6) is %e\n", c_os(-125 * pi / 6));
printf("c_os(-126 * pi / 6) is %e\n", c_os(-126 * pi / 6));
printf("c_os(-127 * pi / 6) is %e\n", c_os(-127 * pi / 6));
printf("c_os(-128 * pi / 6) is %e\n", c_os(-128 * pi / 6));
printf("c_os(-129 * pi / 6) is %e\n", c_os(-129 * pi / 6));
printf("c_os(-130 * pi / 6) is %e\n", c_os(-130 * pi / 6));
printf("c_os(-131 * pi / 6) is %e\n", c_os(-131 * pi / 6));
putchar('\n');
printf("cos(-120 * pi / 6) is %e\n", cos(-120 * pi / 6));
printf("cos(-121 * pi / 6) is %e\n", cos(-121 * pi / 6));
printf("cos(-122 * pi / 6) is %e\n", cos(-122 * pi / 6));
printf("cos(-123 * pi / 6) is %e\n", cos(-123 * pi / 6));
printf("cos(-124 * pi / 6) is %e\n", cos(-124 * pi / 6));
printf("cos(-125 * pi / 6) is %e\n", cos(-125 * pi / 6));
printf("cos(-126 * pi / 6) is %e\n", cos(-126 * pi / 6));
printf("cos(-127 * pi / 6) is %e\n", cos(-127 * pi / 6));
printf("cos(-128 * pi / 6) is %e\n", cos(-128 * pi / 6));
printf("cos(-129 * pi / 6) is %e\n", cos(-129 * pi / 6));
printf("cos(-130 * pi / 6) is %e\n", cos(-130 * pi / 6));
printf("cos(-131 * pi / 6) is %e\n", cos(-131 * pi / 6));
puts("\n/* END c_os.c output */");
return 0;
}

double c_os(double x)
{
unsigned n, negative, sine;
double a, b, c;
static double pi, two_pi, four_pi, half_pi, quarter_pi;

if (x >= -DBL_MAX && DBL_MAX >= x) {
if (1 > pi) {
pi = p_i();
two_pi = 2 * pi;
four_pi = 4 * pi;
half_pi = pi / 2;
quarter_pi = pi / 4;
}
if (0 > x) {
x = -x;
}
while (x > two_pi) {
b = x / 2;
a = two_pi;
while (b > a) {
a *= 2;
}
x -= a;
}
if (x > pi) {
x = two_pi - x;
}
if (x > half_pi) {
x = pi - x;
negative = 1;
} else {
negative = 0;
}
if (x > quarter_pi) {
x = half_pi - x;
sine = 1;
} else {
sine = 0;
}
c = x * x;
x = n = 0;
a = 1;
do {
b = a;
a *= c / ++n;
a /= ++n;
b -= a;
a *= c / ++n;
a /= ++n;
x += b;
} while (b > DBL_EPSILON / 2);
if (sine) {
x = sq_rt(1 - x * x);
}
if (negative) {
x = -x;
}
} else {
x = -HUGE_VAL;
errno = EDOM;
}
return x;
}

double p_i(void)
{
unsigned n;
double p, a, b;

p = 0;
a = 3;
n = 1;
do {
a /= 9;
b = a / n;
a /= 9;
n += 2;
b -= a / n;
n += 2;
p += b;
} while (b > DBL_EPSILON / 4);
a = 2;
n = 1;
do {
a /= 4;
b = a / n;
a /= 4;
n += 2;
b -= a / n;
n += 2;
p += b;
} while (b > DBL_EPSILON / 2);
return 4 * p;
}

double sq_rt(double x)
{
int n;
double a, b;

if (DBL_MAX >= x) {
if (x > 0) {
for (n = 0; x > 2; x /= 4) {
++n;
}
while (0.5 > x) {
--n;
x *= 4;
}
a = x;
b = (1 + x) / 2;
do {
x = b;
b = (a / x + x) / 2;
} while (x > b);
while (n > 0) {
x *= 2;
--n;
}
while (0 > n) {
x /= 2;
++n;
}
} else {
if (0 > x) {
errno = EDOM;
x = HUGE_VAL;
}
}
} else {
errno = ERANGE;
}
return x;
}

/* END c_os.c */

--
pete
 
Reply With Quote
 
 
 
 
Dik T. Winter
Guest
Posts: n/a
 
      10-15-2005
In article <(E-Mail Removed)> http://www.velocityreviews.com/forums/(E-Mail Removed) writes:
> pete wrote:

....
> They don't get the zeros right. I do.
>
> c_os(-123 * pi / 6) is 0.000000e+000
> c_os(-129 * pi / 6) is 0.000000e+000
>
> cos(-123 * pi / 6) is 7.839514e-015
> cos(-129 * pi / 6) is -4.409261e-015


The values you deliver are correct if and only if "-123 * pi / 6" is
exact, which it is not. There is a rounding error, so your function
recevies is an actual value the exact value plus some small value.
The cosinus of that is approximately that small value.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/
 
Reply With Quote
 
 
 
 
Mark McIntyre
Guest
Posts: n/a
 
      10-15-2005
On Sat, 15 Oct 2005 13:07:20 GMT, in comp.lang.c , pete
<(E-Mail Removed)> wrote:

>pete wrote:
>
>> .. and now for a float sqrt algorithm:

>
>I just wrote a cosine function called c_os, for no special reason.
>It calls no functions that I haven't also written.
>
>On my machine, it seems to work better than
>my MS compiler's standard library cos function.
>They don't get the zeros right. I do.


Its equally possible that you're getting the cos of a number close to
n pi/2 wrong....

Seriously tho, both answers are within the possible accuracy of a
double. One might be mathematically better than the other, but neither
is more right than the other in computing terms.

Its even possible, (I've not read your algo), that you have
special-case handling.
--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>

----== Posted via Newsfeeds.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! 120,000+ Newsgroups
----= East and West-Coast Server Farms - Total Privacy via Encryption =----
 
Reply With Quote
 
Christian Bau
Guest
Posts: n/a
 
      10-15-2005
In article <(E-Mail Removed)>,
pete <(E-Mail Removed)> wrote:

> pete wrote:
>
> > .. and now for a float sqrt algorithm:

>
> I just wrote a cosine function called c_os, for no special reason.
> It calls no functions that I haven't also written.
>
> On my machine, it seems to work better than
> my MS compiler's standard library cos function.
> They don't get the zeros right. I do.
>
> c_os(-123 * pi / 6) is 0.000000e+000
> c_os(-129 * pi / 6) is 0.000000e+000
>
> cos(-123 * pi / 6) is 7.839514e-015
> cos(-129 * pi / 6) is -4.409261e-015


Well, they might get it "righter" than you do.

Whatever value "pi" is in your C program, it is _not_ the mathematical
constant Pi, but it is equal to Pi + epsilon for some very small
epsilon. Pi is an irrational number, and floating point numbers in your
C implementation are most likely all rational numbers.

When you calculate (-129 * pi / 6), that error epsilon is multiplied by
-129/6 = -21.5, so your argument is off by -21.5 epsilon from a zero of
the cosine function (plus whatever rounding errors the multiplication
and division added).

So calculating cos (-129 * pi / 6) using a _good_ implementation is very
unlikely to produce a value of zero.
 
Reply With Quote
 
pete
Guest
Posts: n/a
 
      10-16-2005
Christian Bau wrote:
>
> In article <(E-Mail Removed)>,
> pete <(E-Mail Removed)> wrote:
>
> > pete wrote:
> >
> > > .. and now for a float sqrt algorithm:

> >
> > I just wrote a cosine function called c_os, for no special reason.
> > It calls no functions that I haven't also written.
> >
> > On my machine, it seems to work better than
> > my MS compiler's standard library cos function.
> > They don't get the zeros right. I do.
> >
> > c_os(-123 * pi / 6) is 0.000000e+000
> > c_os(-129 * pi / 6) is 0.000000e+000
> >
> > cos(-123 * pi / 6) is 7.839514e-015
> > cos(-129 * pi / 6) is -4.409261e-015

>
> Well, they might get it "righter" than you do.


They might, but ...
I subtracted the theoretical value from the calculated value,
in the cases where the values were rational numbers
which could most likely be represented exactly in floating point,
like 1.0 and 0.5.

c_os(-120 * pi / 6) - 1.0 is 0.000000e+000
c_os(-122 * pi / 6) - 0.5 is 1.221245e-015
c_os(-123 * pi / 6) + 0.0 is 0.000000e+000
c_os(-124 * pi / 6) + 0.5 is 1.998401e-015
c_os(-126 * pi / 6) + 1.0 is 0.000000e+000
c_os(-128 * pi / 6) + 0.5 is -4.329870e-015
c_os(-129 * pi / 6) + 0.0 is 0.000000e+000
c_os(-130 * pi / 6) - 0.5 is 1.221245e-015

cos(-120 * pi / 6) - 1.0 is 0.000000e+000
cos(-122 * pi / 6) - 0.5 is 3.182025e-015
cos(-123 * pi / 6) + 0.0 is 7.839514e-015
cos(-124 * pi / 6) + 0.5 is 4.242944e-015
cos(-126 * pi / 6) + 1.0 is 0.000000e+000
cos(-128 * pi / 6) + 0.5 is -6.364809e-015
cos(-129 * pi / 6) + 0.0 is -4.409261e-015
cos(-130 * pi / 6) - 0.5 is -1.272257e-015

In all cases, c_os had the smaller or same difference.

I tried it with a different value for pi,
pi = 3.141592653589793238462643383279502884197169399375 1;
and got the same results.

I think the c_os function is OK.

--
pete
 
Reply With Quote
 
Christian Bau
Guest
Posts: n/a
 
      10-16-2005
In article <(E-Mail Removed)>,
pete <(E-Mail Removed)> wrote:

> Christian Bau wrote:
> >
> > In article <(E-Mail Removed)>,
> > pete <(E-Mail Removed)> wrote:
> >
> > > pete wrote:
> > >
> > > > .. and now for a float sqrt algorithm:
> > >
> > > I just wrote a cosine function called c_os, for no special reason.
> > > It calls no functions that I haven't also written.
> > >
> > > On my machine, it seems to work better than
> > > my MS compiler's standard library cos function.
> > > They don't get the zeros right. I do.
> > >
> > > c_os(-123 * pi / 6) is 0.000000e+000
> > > c_os(-129 * pi / 6) is 0.000000e+000
> > >
> > > cos(-123 * pi / 6) is 7.839514e-015
> > > cos(-129 * pi / 6) is -4.409261e-015

> >
> > Well, they might get it "righter" than you do.

>
> They might, but ...
> I subtracted the theoretical value from the calculated value,
> in the cases where the values were rational numbers
> which could most likely be represented exactly in floating point,
> like 1.0 and 0.5.
>
> c_os(-120 * pi / 6) - 1.0 is 0.000000e+000
> c_os(-122 * pi / 6) - 0.5 is 1.221245e-015
> c_os(-123 * pi / 6) + 0.0 is 0.000000e+000
> c_os(-124 * pi / 6) + 0.5 is 1.998401e-015
> c_os(-126 * pi / 6) + 1.0 is 0.000000e+000
> c_os(-128 * pi / 6) + 0.5 is -4.329870e-015
> c_os(-129 * pi / 6) + 0.0 is 0.000000e+000
> c_os(-130 * pi / 6) - 0.5 is 1.221245e-015
>
> cos(-120 * pi / 6) - 1.0 is 0.000000e+000
> cos(-122 * pi / 6) - 0.5 is 3.182025e-015
> cos(-123 * pi / 6) + 0.0 is 7.839514e-015
> cos(-124 * pi / 6) + 0.5 is 4.242944e-015
> cos(-126 * pi / 6) + 1.0 is 0.000000e+000
> cos(-128 * pi / 6) + 0.5 is -6.364809e-015
> cos(-129 * pi / 6) + 0.0 is -4.409261e-015
> cos(-130 * pi / 6) - 0.5 is -1.272257e-015
>
> In all cases, c_os had the smaller or same difference.


Seems you didn't get the point.

Whatever value your variable "pi" has, it is not the mathematical
constant Pi. Therefore, cos ((2k + 1/2) pi) is not _supposed_ to be
zero. The larger you choose k, the further away from zero it _should_
be.

Your c_os function does its argument reduction with an approximation to
the mathematical Pi rounded to double, most likely with 53 bits of
mantissa. If the library you used uses the built-in FCOS instruction of
the P4 chip, that uses an approximation to Pi with 66 bits of mantissa.
It will _correctly_ give non-zero results when you try to calculate cos
(-129 * pi / 6) using your incorrect approximation pi.
 
Reply With Quote
 
pete
Guest
Posts: n/a
 
      10-16-2005
Christian Bau wrote:
>
> In article <(E-Mail Removed)>,
> pete <(E-Mail Removed)> wrote:
>
> > Christian Bau wrote:
> > >
> > > In article <(E-Mail Removed)>,
> > > pete <(E-Mail Removed)> wrote:
> > >
> > > > pete wrote:
> > > >
> > > > > .. and now for a float sqrt algorithm:
> > > >
> > > > I just wrote a cosine function called c_os,
> > > > for no special reason.
> > > > It calls no functions that I haven't also written.
> > > >
> > > > On my machine, it seems to work better than
> > > > my MS compiler's standard library cos function.
> > > > They don't get the zeros right. I do.


> > > Well, they might get it "righter" than you do.


> > c_os(-120 * pi / 6) - 1.0 is 0.000000e+000
> > c_os(-126 * pi / 6) + 1.0 is 0.000000e+000


> > cos(-120 * pi / 6) - 1.0 is 0.000000e+000
> > cos(-126 * pi / 6) + 1.0 is 0.000000e+000


> Seems you didn't get the point.


I didn't.
I've never seen a pocket calculator that didn't yield 0 exactly
as the cosine of the product of an odd integer
and the calculator's approximation of pi / 2.

Both cos and c_os hit 1.0 and -1.0 exactly,
for cos(-120 * pi / 6) and cos(-120 * pi / 6).
That makes me think that c_os and cos are normalizing
large x the same way as each other on my system.

c_os normalizes x down into the range of between 0 an pi / 4.
I assume that most implementations of cos
normalize x down to at least 2 * pi.
I think that if c_os's approximation of pi,
matches the implementation's approximation of pi,
then I shouldn't see any difference in the normalized value for x,
the point not being that that the value of x is exact,
(I know it usually isn't)
but rather that it is possible for cos and c_os to see large x,
the same.

--
pete
 
Reply With Quote
 
pete
Guest
Posts: n/a
 
      10-16-2005
Dik T. Winter wrote:
>
> In article <(E-Mail Removed)> (E-Mail Removed) writes:
> > pete wrote:

> ...
> > They don't get the zeros right. I do.
> >
> > c_os(-123 * pi / 6) is 0.000000e+000
> > c_os(-129 * pi / 6) is 0.000000e+000
> >
> > cos(-123 * pi / 6) is 7.839514e-015
> > cos(-129 * pi / 6) is -4.409261e-015

>
> The values you deliver are correct if and only if "-123 * pi / 6" is
> exact, which it is not.


Couldn't 0.0 also be delivered if the true cosine
of the actual argument was positive
and also closer to 0.0 than to DBL_MIN ?

--
pete
 
Reply With Quote
 
pete
Guest
Posts: n/a
 
      10-16-2005
Mark McIntyre wrote:
>
> On Sat, 15 Oct 2005 13:07:20 GMT, in comp.lang.c , pete
> <(E-Mail Removed)> wrote:
>
> >pete wrote:
> >
> >> .. and now for a float sqrt algorithm:

> >
> >I just wrote a cosine function called c_os, for no special reason.
> >It calls no functions that I haven't also written.
> >
> >On my machine, it seems to work better than
> >my MS compiler's standard library cos function.
> >They don't get the zeros right. I do.

>
> Its equally possible that you're getting the cos of a number close to
> n pi/2 wrong....


> Its even possible, (I've not read your algo), that you have
> special-case handling.


There's just a smidgeon of special case handling in sq_rt.

In these cases:

c_os(-123 * pi / 6) + 0.0 is 0.000000e+000
c_os(-129 * pi / 6) + 0.0 is 0.000000e+000

cos(-123 * pi / 6) + 0.0 is 7.839514e-015
cos(-129 * pi / 6) + 0.0 is -4.409261e-015

For values of x, which are equal to the product
of an odd integer and c_os's approximation of pi / 2,
the c_os function, normalizes x all the way down to 0.0,
and calculates the sine of 0.0 from the cosine.
c_os then reports the sine of 0.0, as the cosine of pi / 2.
The special case is that that involves using sq_rt
to find the root of zero. sq_rt(0) is a special case.

When I was writing c_os, it didn't originally
normalize x any lower than pi / 2,
but I found that I was having accuaracy problems all around
the vicinity of pi / 2.

The infinite series that c_os approximates,
converges rapidly for small x,
but not rapidly when x is greater than 1.

The two generally recognized sources of error
when summing a series, are:
1 the error in each term
2 the error from only summing a finite number of terms.
A more rapidly converging series, reduces both errors.

I was wondering if my implementation's cos function
doesn't normalized x, when x is about pi / 2.

c_os agrees with cos on my system for some input.

c_os(-120 * pi / 6) - 1.0 is 0.000000e+000
c_os(-126 * pi / 6) + 1.0 is 0.000000e+000

cos(-120 * pi / 6) - 1.0 is 0.000000e+000
cos(-126 * pi / 6) + 1.0 is 0.000000e+000

I'm am suspecting that it's because the series summing do loop,
converges after just one iteration when x == 0.0,
and also because c_os and cos are using a value for pi
which is effectively the same as each other's.

--
pete
 
Reply With Quote
 
Christian Bau
Guest
Posts: n/a
 
      10-16-2005
In article <(E-Mail Removed)>,
pete <(E-Mail Removed)> wrote:

> Christian Bau wrote:
> >
> > In article <(E-Mail Removed)>,
> > pete <(E-Mail Removed)> wrote:
> >
> > > Christian Bau wrote:
> > > >
> > > > In article <(E-Mail Removed)>,
> > > > pete <(E-Mail Removed)> wrote:
> > > >
> > > > > pete wrote:
> > > > >
> > > > > > .. and now for a float sqrt algorithm:
> > > > >
> > > > > I just wrote a cosine function called c_os,
> > > > > for no special reason.
> > > > > It calls no functions that I haven't also written.
> > > > >
> > > > > On my machine, it seems to work better than
> > > > > my MS compiler's standard library cos function.
> > > > > They don't get the zeros right. I do.

>
> > > > Well, they might get it "righter" than you do.

>
> > > c_os(-120 * pi / 6) - 1.0 is 0.000000e+000
> > > c_os(-126 * pi / 6) + 1.0 is 0.000000e+000

>
> > > cos(-120 * pi / 6) - 1.0 is 0.000000e+000
> > > cos(-126 * pi / 6) + 1.0 is 0.000000e+000

>
> > Seems you didn't get the point.

>
> I didn't.
> I've never seen a pocket calculator that didn't yield 0 exactly
> as the cosine of the product of an odd integer
> and the calculator's approximation of pi / 2.
>
> Both cos and c_os hit 1.0 and -1.0 exactly,
> for cos(-120 * pi / 6) and cos(-120 * pi / 6).
> That makes me think that c_os and cos are normalizing
> large x the same way as each other on my system.


The values 1.0 and -1.0 are much harder to miss than 0.0.

For x = 2*k*Pi + eps, cos (x) is about 1 - eps^2 / 2. With typical 53
mantissa for double numbers, you should get a result of 1.0 as long as
eps < 2^-27. That's quite a large number
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
cosine function returning long double (high resolution) ratcharit@gmail.com C Programming 11 03-06-2008 08:36 PM
function/process to generate sine and cosine wave FPGA VHDL 15 02-09-2008 02:21 PM
sine and cosine wave generation FPGA VHDL 8 01-14-2008 09:11 PM
cosine calcs Niv VHDL 9 08-16-2006 05:25 PM
having trouble with Discrete Cosine Transform program khaleel.alyasini@gmail.com C++ 1 03-24-2006 01:21 PM



Advertisments