Velocity Reviews > cosine function returning long double (high resolution)

# cosine function returning long double (high resolution)

pete
Guest
Posts: n/a

 03-06-2008
http://www.velocityreviews.com/forums/(E-Mail Removed) wrote:
>
> Currently using cosine function in math.h
>
> Currently I get:
>
> 1 = cos(1e^-7)
>
> Is there another way for cos to return value of high accuracy say:
>
> 0.999999 = cos(1e^-7)

/* BEGIN new.c output */

LDBL_EPSILON is 2.220446e-016

fs_cosl(1.0 / 10000000) - 1.0 is -4.996004e-015
fs_cosl(1.0 / 1000000) - 1.0 is -5.000445e-013
fs_cosl(1.0 / 100000) - 1.0 is -5.000000e-011
fs_cosl(1.0 / 10000) - 1.0 is -5.000000e-009
fs_cosl(1.0 / 100) - 1.0 is -4.999958e-005

/* END new.c output */

/* BEGIN new.c */

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

#include "fs_cosl.h"

#define DEGREES_TO_RADIANS(A) ((A) * atan(1) / 45)

int main(void)
{
puts("/* BEGIN new.c output */\n");
printf("LDBL_EPSILON is %Le\n\n",
LDBL_EPSILON);
printf("fs_cosl(1.0 / 10000000) - 1.0 is %Le\n",
fs_cosl(1.0 / 10000000) - 1.0);
printf("fs_cosl(1.0 / 1000000) - 1.0 is %Le\n",
fs_cosl(1.0 / 1000000) - 1.0);
printf("fs_cosl(1.0 / 100000) - 1.0 is %Le\n",
fs_cosl(1.0 / 100000) - 1.0);
printf("fs_cosl(1.0 / 10000) - 1.0 is %Le\n",
fs_cosl(1.0 / 10000) - 1.0);
printf("fs_cosl(1.0 / 100) - 1.0 is %Le\n\n",
fs_cosl(1.0 / 100) - 1.0);
puts("/* END new.c output */");
return 0;
}

/* END new.c */

/* BEGIN fs_cosl.h */
/*
** Portable freestanding code.
*/
#ifndef H_FS_COSL_H
#define H_FS_COSL_H

long double fs_cosl(long double x);
long double fs_sqrtl(long double x);
long double fs_fmodl(long double x, long double y);

#endif

/* END fs_cosl.h */

/* BEGIN fs_cosl.c */

#include <float.h>

#include "fs_cosl.h"
/*
** pi == (atan(1.0 / 3) + atan(1.0 / 2)) * 4
*/
static long double fs_pil(void);

long double fs_cosl(long double x)
{
long unsigned n;
int negative, sine;
long double a, b, c;
static long double pi, two_pi, half_pi, third_pi;
static int initialized;

if (0 > x) {
x = -x;
}
if (LDBL_MAX >= x) {
if (!initialized) {
initialized = 1;
pi = fs_pil();
two_pi = 2 * pi;
half_pi = pi / 2;
third_pi = pi / 3;
}
if (x > two_pi) {
x = fs_fmodl(x, two_pi);
}
if (x > pi) {
x = two_pi - x;
}
if (x > half_pi) {
x = pi - x;
negative = 1;
} else {
negative = 0;
}
if (x > third_pi) {
x = half_pi - x;
sine = 1;
} else {
sine = 0;
}
c = x * x;
x = n = 0;
a = 1;
do {
b = a;
a *= c;
a /= ++n;
a /= ++n;
b -= a;
a *= c;
a /= ++n;
a /= ++n;
x += b;
} while (b > LDBL_EPSILON / 2);
if (sine) {
x = fs_sqrtl((1 - x) * (1 + x));
}
if (negative) {
x = -x;
}
} else {
x = -LDBL_MAX;
}
return x;
}

long double fs_sqrtl(long double x)
{
long int n;
long double a, b;

if (x > 0 && LDBL_MAX >= x) {
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 (x != 0) {
x = LDBL_MAX;
}
}
return x;
}

long double fs_fmodl(long double x, long double y)
{
long double a, b;
const long double c = x;

if (0 > c) {
x = -x;
}
if (0 > y) {
y = -y;
}
if (y != 0 && LDBL_MAX >= y && LDBL_MAX >= x) {
while (x >= y) {
a = x / 2;
b = y;
while (a >= b) {
b *= 2;
}
x -= b;
}
} else {
x = 0;
}
return 0 > c ? -x : x;
}

static long double fs_pil(void)
{
long unsigned n;
long double a, b;
static long double p;
static int initialized;

if (!initialized) {
initialized = 1;
n = 1;
a = 3;
do {
a /= 9;
b = a / n;
n += 2;
a /= 9;
b -= a / n;
n += 2;
p += b;
} while (b > LDBL_EPSILON / 4);
n = 1;
a = 2;
do {
a /= 4;
b = a / n;
n += 2;
a /= 4;
b -= a / n;
n += 2;
p += b;
} while (b > LDBL_EPSILON / 2);
p *= 4;
}
return p;
}

/* END fs_cosl.c */

--
pete

ratcharit@gmail.com
Guest
Posts: n/a

 03-06-2008