468,512 Members | 1,390 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 468,512 developers. It's quick & easy.

replacement for 'pow()' function

Hello,

I'd like to make a simple replacement of 'pow()' function for the embedded
platform I'm working on. What is the better way, probably bit shifting?
Thanks.

Best regards, Roman Mashak.
Dec 13 '07 #1
15 6012
Roman Mashak wrote, On 14/12/07 00:43:
Hello,

I'd like to make a simple replacement of 'pow()' function for the embedded
platform I'm working on. What is the better way, probably bit shifting?
Thanks.
It all depends on the details of your requirements, including whether
speed, accuracy or space is more important and on the limits for each. I
know of code that uses a small look-up-table for trig functions, one
that leads to a lot of inaccuracy but it "good enough", fast and small.
I know other code that uses a look-up-table and interpolation of some
form (I never checked the details) because speed was a bit less
important (the processor was faster).
--
Flash Gordon
Dec 13 '07 #2
Roman Mashak wrote:
>
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 */
--
pete
Dec 13 '07 #3
pete wrote:
>
/* 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);
Interesting. What are the license requirements on the code?

--
Thad
Dec 13 '07 #4
"Roman Mashak" <mr*@tusur.ruwrote:
I'd like to make a simple replacement of 'pow()' function for the embedded
platform I'm working on. What is the better way, probably bit shifting?
Bit shifting may suffice for a function intended to produce positive
integral powers of integers, but for a complete implementation of pow()
it won't do: you try shifting a float by a double and see what the
result is...

Richard
Dec 13 '07 #5
Thad Smith wrote:
>
pete wrote:

/* 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);

Interesting. What are the license requirements on the code?
It's just trivial code.

--
pete
Dec 13 '07 #6
Thad Smith wrote:
>
pete wrote:

/* 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);

Interesting. What are the license requirements on the code?
I don't know. I think it's public domain.
Here's some more:
http://www.mindspring.com/~pfilandr/C/fs_math/
I got bored with it after writing fs_cos.

--
pete
Dec 13 '07 #7
pete wrote:
Thad Smith wrote:
>Interesting. What are the license requirements on the code?

It's just trivial code.
Licensing isn't about whether the code is trivial or not.

--
Chris "bits typed with meta-bits [typed with meta-meta-bits ...]" Dollin

Hewlett-Packard Limited registered no:
registered office: Cain Road, Bracknell, Berks RG12 1HN 690597 England

Dec 13 '07 #8
On Thu, 13 Dec 2007 16:43:45 -0800, Roman Mashak wrote:
Hello,

I'd like to make a simple replacement of 'pow()' function for the embedded
platform I'm working on. What is the better way, probably bit shifting?
Thanks.

Best regards, Roman Mashak.
You can implement exp and log via the cordic algorithm. See eg
http://en.wikipedia.org/wiki/CORDIC
and links therein

Dec 13 '07 #9
On Dec 13, 4:43 pm, "Roman Mashak" <m...@tusur.ruwrote:
Hello,

I'd like to make a simple replacement of 'pow()' function for the embedded
platform I'm working on. What is the better way, probably bit shifting?
Thanks.
Cephes by Moshier
http://www.moshier.net/#Cephes

Or libfdm by SunSoft
http://www.koders.com/noncode/fid7DF...4BEBB9931.aspx

If you already have an exp() then pow() is trivial to implement.
Dec 13 '07 #10
pete wrote:
Thad Smith wrote:
>pete wrote:
>>/* 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);
Interesting. What are the license requirements on the code?

I don't know. I think it's public domain.
If it is public domain, it is rather old, or explicitly declared such by
the author. The default in the US and most other countries is copyright
at the time of creation.

If you wrote the code and would like others to use it, I recommend that
you place a comment in the header stating that it is donated to the
public domain by the author, and include your name.
--
Thad
Dec 14 '07 #11
Thad Smith wrote:
>
pete wrote:
Thad Smith wrote:
pete wrote:
/* 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);
Interesting. What are the license requirements on the code?
I don't know. I think it's public domain.

If it is public domain, it is rather old,
or explicitly declared such by the author.
The default in the US and most other countries is copyright
at the time of creation.

If you wrote the code and would like others to use it,
I recommend that
you place a comment in the header stating that it is donated to the
public domain by the author, and include your name.
OK

--
pete
Dec 14 '07 #12
Op Thu, 13 Dec 2007 14:37:43 +0100 schreef Chris Dollin
<ch**********@hp.com>:
pete wrote:
>Thad Smith wrote:
>>Interesting. What are the license requirements on the code?

It's just trivial code.

Licensing isn't about whether the code is trivial or not.
Yes it is. Any license on something trivial cannot be enforced. Except
in an unreasonable justice system, of course...

--
Gemaakt met Opera's revolutionaire e-mailprogramma:
http://www.opera.com/mail/
Dec 14 '07 #13
pete wrote:
Thad Smith wrote:
>pete wrote:
>>Thad Smith wrote:
pete wrote:
/* 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);
Interesting. What are the license requirements on the code?
I don't know. I think it's public domain.
If it is public domain, it is rather old,
or explicitly declared such by the author.
The default in the US and most other countries is copyright
at the time of creation.

If you wrote the code and would like others to use it,
I recommend that
you place a comment in the header stating that it is donated to the
public domain by the author, and include your name.

OK
Thanks. It matters to people who program for a living and are careful
about getting legal rights to the code they use.

As a suggestion, consider contributing it to Snippets.org, which has a
public domain collection of code snippets.

--
Thad
Dec 15 '07 #14
Thad Smith wrote:
>
pete wrote:
Thad Smith wrote:
pete wrote:
Thad Smith wrote:
pete wrote:
/* 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);
Interesting. What are the license requirements on the code?
I don't know. I think it's public domain.
If it is public domain, it is rather old,
or explicitly declared such by the author.
The default in the US and most other countries is copyright
at the time of creation.

If you wrote the code and would like others to use it,
I recommend that
you place a comment in the header stating that it is donated to the
public domain by the author, and include your name.
OK

Thanks. It matters to people who program for a living and are careful
about getting legal rights to the code they use.

As a suggestion, consider contributing it to Snippets.org, which has a
public domain collection of code snippets.
OK

--
pete
Dec 16 '07 #15
"Boudewijn Dijkstra" <bo*******@indes.comwrote:
Op Thu, 13 Dec 2007 14:37:43 +0100 schreef Chris Dollin
pete wrote:
Thad Smith wrote:
>Interesting. What are the license requirements on the code?

It's just trivial code.
Licensing isn't about whether the code is trivial or not.

Yes it is. Any license on something trivial cannot be enforced. Except
in an unreasonable justice system, of course...
pete is in the USA.

Richard
Dec 17 '07 #16

This discussion thread is closed

Replies have been disabled for this discussion.

Similar topics

2 posts views Thread by Arik Peltz | last post: by
8 posts views Thread by Black Eagle | last post: by
5 posts views Thread by Magix | last post: by
1 post views Thread by fmendoza | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.