Page 1 of 1
Exponential Function
Posted: Sat Jan 10, 2004 6:51 pm
by chris
Is there any way I can make this function better or is there any visible problems with it? It seems to act weird with things like 10^10 and above...
Code: Select all
int pow(int base, int exp)
{
int res = base;
if(exp == 0)
return 1;
if(exp == 1)
return base;
for (; exp-1 > 0; exp--) {
res *= base;
}
return res;
}
Re:Exponential Function
Posted: Sat Jan 10, 2004 6:55 pm
by chris
Err I guess 10^10 isn't working because I used ints?
Re:Exponential Function
Posted: Sat Jan 10, 2004 8:03 pm
by Tim
Right. Try float, double, or long long.
Re:Exponential Function
Posted: Sat Jan 10, 2004 8:41 pm
by sonneveld
Code: Select all
double pow(double base, int exp)
{
double res = 1;
while(exp--)
res *= base;
return res;
}
- Nick
Re:Exponential Function
Posted: Sat Jan 10, 2004 8:56 pm
by chris
Thanks.
Re:Exponential Function
Posted: Sun Jan 11, 2004 12:45 pm
by Schol-R-LEA
While that will work, it doesn't catch the screw cases (i.e., zero bases, zero exponents, negative exponents), it may have serious round-off errors, and won't be very efficient, either. To really do it right, you need to dig into the double-precision format and manipulate the exponent and mantissa directly. Even without that, though, a few simple checks and modifications can speed things up substantially, as shown in this example (tested for basic function but not for round-off errors):
Code: Select all
#define reciprocal(n) (1 / (n))
#define odd(x) ((x) % 2)
#define even(x) (!odd(x))
double pow(double base, int exp)
{
double accumulate = 1.0;
if (0.0 == base) // sanity check - is base zero? 0^exp == 0
{
return 0.0;
}
else if (1 == base) // sanity check - is base one? 1^exp == 1
{
return 1.0;
}
else if (0 == exp) // sanity check - is exp zero? base^0 == 1
{
return 1.0;
}
else if (0 > exp) // base^(-exp) == 1 / base^exp
{
return reciprocal(pow(base, abs(exp)));
}
else
while (0 < exp)
{
if (even(exp))
// by halving even exponents and squaring the base,
// the number of multiplications needed is reduced
// turning the algorithm from O(n) time and O(1) space
// to 0(log n) for time while remaining O(1) for space.
{
exp /= 2;
base *= base;
}
else
{
// For odd numbers, multiply the base by base^(exp-1)
// Since exp - 1 is guaranteed to be even, it will then
// use the halving formula above to reduce the total
// number of multiplications.
accumulate *= base;
exp--;
}
}
return accumulate;
}
Re:Exponential Function
Posted: Tue Jan 13, 2004 5:24 pm
by chris
Thanks alot Schol-R-LEA. Comments and everything ;D
Re:Exponential Function
Posted: Tue Jan 13, 2004 5:49 pm
by df
since exp is an integer
replace
with
might save you _1_ cycle
hehe
Re:Exponential Function
Posted: Tue Jan 13, 2004 7:02 pm
by chris
Schol-R-LEA wrote:
Code: Select all
return reciprocal(pow(base, abs(exp)));
One question: Where does the "abs" come from? I didn't notice it before
Re:Exponential Function
Posted: Tue Jan 13, 2004 8:56 pm
by Schol-R-LEA
df wrote:
since exp is an integer
replace
with
might save you _1_ cycle
hehe
I
knew I'd missed something. Actually, it might be more than one cycle, and even if it is only one, well, they do add up...
Actually, the real annoyance to me is that you're doing the division
twice, once for the modulo operation in the odd() macro, and then again here. In assembly it wouldn't be a big deal; you'd store the exponent in eax, shift, and jump on the carry flag, but trying to do the equivalent in C is too much work.
In any case, a better approach has occurred to me, at least once you've gotten the exp() and log() functions working:
Code: Select all
double pow(double base, int expt)
{
return exp((expt * log(base));
}
Writing exp() and log() are left as an exercise (hint: the log() function maps very nicely to the FPU instruction FLDL2E). ;D
Re:Exponential Function
Posted: Tue Jan 13, 2004 8:58 pm
by Schol-R-LEA
chris wrote:
Schol-R-LEA wrote:
Code: Select all
return reciprocal(pow(base, abs(exp)));
One question: Where does the "abs" come from? I didn't notice it before
Oops. It's exactly the kind of math.h function were supposed to be implementing here, isn't it? OK, then here's a quick and dirty version that should do:
Code: Select all
#define abs(x) (0 < (x)) ? (x) : (-x))
Re:Exponential Function
Posted: Wed Jan 14, 2004 7:37 am
by chris
Thanks ;D
Re:Exponential Function
Posted: Sat Jan 31, 2004 2:59 am
by Therx
I've seen it implemented before something like this (off the top of my head)
Code: Select all
double pow(double base, int expt)
{
if(expt == 0)
return 1;
else
return base * pow(base, expt - 1);
}
With the claim that its multithreaded friendly?
Pete
Re:Exponential Function
Posted: Sat Jan 31, 2004 4:04 am
by Schol-R-LEA
Pete wrote:
I've seen it implemented before something like this (off the top of my head)
Code: Select all
double pow(double base, int expt)
{
if(expt == 0)
return 1;
else
return base * pow(base, expt - 1);
}
With the claim that its multithreaded friendly?
Certainly; this is one of the simplest ways to implement this. There's no reason for it not to be multithread friendly; it's a basic recursive function, with no shared memory state. If you think about it, you'd see that
has to be at least serially reentrant, or the recursion wouldn't work properly.
I'd still do it like this, though, to speed it up:
Code: Select all
double pow(double base, int expt)
{
if(expt == 0)
return 1;
else if (expt == 1)
return base;
else if (even(base))
{
base = pow(base, (expt / 2));
return base * base;
}
else
return base * pow(base, expt - 1);
}
This is what I based my iterative version on (though I used a slightly different equivalency in the final version); like that one, this is O(n log n) in time, but unlike that version, it is also O(n log n) in space, not O(1). This one is easier to understand, however, though not as simple as yours.