Dave wrote:
The only problem with fdlibm is that some funcions (especially the LOG
functions) only work properley with binary float, whereas IBM float is
Hex floating point. (that is the exponent is a power of 16, not of
two). This means that code that multiplys or divides by two by
manipulating the exponent directly does not work. I guess its fixable
though, and the code does have pretty good explainations of what is
happening, and good examples of suitable constants for use in the
polynomial approximations, which are too a greater accuracy than those
in AMS55.
Sometimes I write math functions, when I have nothing else to do.
/* BEGIN constants.c output */
sq_rt(2) is 1.414214
l_og(2) is 0.693147
e_xp(1) is 2.718282
pi is 3.141593
sq_rt(2) * sq_rt(2) - 2.0 is -4.440892e-016
sqrt(2) * sqrt(2) - 2.0 is 4.440892e-016
pow(e_xp(1), l_og(2)) - 2.0 is -2.220446e-016
pow(exp(1), log(2)) - 2.0 is -2.220446e-016
3.1415926535897932384626433832795028841971693993751 - fs_pi()
is 0.000000e+000
/* END constants.c output */
/* BEGIN constants.c */
#include <stdio.h>
#include <float.h>
#include <errno.h>
#include <math.h>
double sq_rt(double x);
double l_og(double x);
double e_xp(double x);
double fs_pi(void);
int main(void)
{
puts("\n/* BEGIN constants.c output */\n");
printf("sq_rt(2) is %f\n", sq_rt(2));
printf(" l_og(2) is %f\n", l_og(2));
printf(" e_xp(1) is %f\n", e_xp(1));
printf(" pi is %f\n\n", fs_pi( ));
printf("sq_rt(2) * sq_rt(2) - 2.0 is %e\n",
sq_rt(2) * sq_rt(2) - 2.0);
printf(" sqrt(2) * sqrt(2) - 2.0 is %e\n\n",
sqrt(2) * sqrt(2) - 2.0);
printf("pow(e_xp(1), l_og(2)) - 2.0 is %e\n",
pow(e_xp(1), l_og(2)) - 2.0);
printf(" pow(exp(1), log(2)) - 2.0 is %e\n\n",
pow(exp(1), log(2)) - 2.0);
printf("3.1415926535897932384626433832795028841971693993751"
" - fs_pi()\nis %e\n",
3.1415926535897932384626433832795028841971693993751
- fs_pi());
puts("\n/* END constants.c output */");
return 0;
}
double sq_rt(double x)
{
int n;
double a, b;
if (DBL_MAX >= x && 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 (x != 0) {
errno = EDOM;
x = HUGE_VAL;
}
}
return x;
}
double l_og(double x)
{
int n;
double a, b, c, epsilon;
static double A, B, C;
if (DBL_MAX >= x && x > 0) {
if (1 > A) {
A = sq_rt(2);
B = A / 2;
C = l_og(A);
}
for (n = 0; x > A; x /= 2) {
++n;
}
while (B > x) {
--n;
x *= 2;
}
a = (x - 1) / (x + 1);
x = C * n + a;
c = a * a;
n = 1;
epsilon = DBL_EPSILON * x;
if (0 > a) {
if (epsilon > 0) {
epsilon = -epsilon;
}
do {
n += 2;
a *= c;
b = a / n;
x += b;
} while (epsilon > b);
} else {
if (0 > epsilon) {
epsilon = -epsilon;
}
do {
n += 2;
a *= c;
b = a / n;
x += b;
} while (b > epsilon);
}
x *= 2;
} else {
errno = x != 0 ? EDOM : ERANGE;
x = -HUGE_VAL;
}
return x;
}
double e_xp(double x)
{
unsigned n, square;
double b, e;
static double x_max;
if (1 > x_max) {
x_max = l_og(DBL_MAX);
}
if (x_max >= x && x >= -x_max) {
for (square = 0; x > 1; x /= 2) {
++square;
}
while (-1 > x) {
++square;
x /= 2;
}
e = b = n = 1;
do {
b /= n++;
b *= x;
e += b;
b /= n++;
b *= x;
e += b;
} while (b > DBL_EPSILON / 4);
while (square-- != 0) {
e *= e;
}
} else {
errno = x > DBL_MAX || -DBL_MAX > x ? EDOM : ERANGE;
e = 0 > x ? 0 : HUGE_VAL;
}
return e;
}
double fs_pi(void)
{
unsigned n;
double p, a, b;
p = 0;
n = 1;
a = 3;
do {
a /= 9;
b = a / n;
n += 2;
a /= 9;
b -= a / n;
n += 2;
p += b;
} while (b > DBL_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 > DBL_EPSILON / 2);
return 4 * p;
}
/* END constants.c */
