Roman said:
Hello,
I'd like to make a simple replacement of 'pow()'
function for the embedded platform I'm working on.
/* BEGIN fs_pow_test.c */
#include <stdio.h>
/*
** fs_pow is implemented in portable freestanding code
*/
#include "fs_pow.h"
int main(void)
{
puts("/* BEGIN fs_pow_test.c output */\n");
printf("fs_pow(2, 4) is %f\n\n",
fs_pow(2, 4));
printf("fs_pow(0.0001, -0.25) - 10 is %e\n\n",
fs_pow(0.0001, -0.25) - 10);
puts("/* END fs_pow_test.c output */");
return 0;
}
/* END fs_pow_test.c */
/* BEGIN fs_pow.h */
/*
** Portable code for either freestanding or hosted
** implementations of C.
*/
#ifndef H_FS_POW_H
#define H_FS_POW_H
double fs_pow(double x, double y);
double fs_fmod(double x, double y);
double fs_exp(double x);
double fs_log(double x);
double fs_sqrt(double x);
#endif
/* END fs_pow.h */
/* BEGIN fs_pow.c */
#include <float.h>
#include "fs_pow.h"
double fs_pow(double x, double y)
{
double p;
if (0 > x && fs_fmod(y, 1) == 0) {
if (fs_fmod(y, 2) == 0) {
p = fs_exp(fs_log(-x) * y);
} else {
p = -fs_exp(fs_log(-x) * y);
}
} else {
if (x != 0 || 0 >= y) {
p = fs_exp(fs_log( x) * y);
} else {
p = 0;
}
}
return p;
}
double fs_fmod(double x, double y)
{
double a, b;
const double c = x;
if (0 > c) {
x = -x;
}
if (0 > y) {
y = -y;
}
if (y != 0 && DBL_MAX >= y && DBL_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;
}
double fs_exp(double x)
{
unsigned n, square;
double b, e;
static double x_max, x_min;
static int initialized;
if (!initialized) {
initialized = 1;
x_max = fs_log(DBL_MAX);
x_min = fs_log(DBL_MIN);
}
if (x_max >= x && x >= x_min) {
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 {
e = x > 0 ? DBL_MAX : 0;
}
return e;
}
double fs_log(double x)
{
int n;
double a, b, c, epsilon;
static double A, B, C;
static int initialized;
if (x > 0 && LDBL_MAX >= x) {
if (!initialized) {
initialized = 1;
A = fs_sqrt(2);
B = A / 2;
C = fs_log(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 {
x = -DBL_MAX;
}
return x;
}
double fs_sqrt(double x)
{
int n;
double a, b;
if (x > 0 && DBL_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 / b + b) / 2;
} while (x > b);
while (n > 0) {
x *= 2;
--n;
}
while (0 > n) {
x /= 2;
++n;
}
} else {
if (x != 0) {
x = DBL_MAX;
}
}
return x;
}
/* END fs_pow.c */
