Recursive Functions

[Eric Sosman]

Maybe it's just me but doesn't the contrived nature of the function scream
out "Homework Assignment?" Maybe I'm missing something but I just can't see
why you'd ever want to compute this value..

<off-topic>

If computing positive integer powers is of so little
interest, it's hard to see why Knuth gives the topic an
entire section of its own in "The Art of Computer Programming,
Volume II: Seminumerical Algorithms."

The method outlined is the first serious power-computing
algorithm Knuth introduces in that section. He writes that
some authors have declared the method optimal, but he is
kind enough to spare them embarrassment by omitting their
names; clearly they forgot to consider cases like n==15.

</off-topic>

 

[Glen Herrmannsfeldt]

Yes, but those typically take a floating point type for the x argument,
and an int type for the n argument.

The languages I know of supply routines for integer n, and integer, float,
double, complex, and complex double x.

The OP didn't supply the type of x, but the algorithm works for any type
where mulitply is defined.

A table dimensioned INT_MAX-INT_MIN+1 will take up a lot of memory.

That is an imaginative approach, but not what I had in mind.

#include <stdint.h>

static int8_t hbit[256] =
{-1
,0
,1,1
,2,2,2,2
,3,3,3,3,3,3,3,3
,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4
,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5
,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6
,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6
,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
};

int int_pow(int x, uint8_t n)
{
int t = 1;
if (n == 0) return 1;
if (x == 0) return 0;
switch (hbit[n]) {
case 7: if (n & 1) t *= x; n >>= 1; x *= x;
case 6: if (n & 1) t *= x; n >>= 1; x *= x;
case 5: if (n & 1) t *= x; n >>= 1; x *= x;
case 4: if (n & 1) t *= x; n >>= 1; x *= x;
case 3: if (n & 1) t *= x; n >>= 1; x *= x;
case 2: if (n & 1) t *= x; n >>= 1; x *= x;
case 1: if (n & 1) t *= x; n >>= 1; x *= x;
default: if (n & 1) t *= x;
}
return t;
}

The algorithm also works for n > 255.

-- glen

 

[Glen Herrmannsfeldt]

Maybe it's just me but doesn't the contrived nature of the function scream
out "Homework Assignment?" Maybe I'm missing something but I just can't see
why you'd ever want to compute this value..

Computing integer powers of numbers is fairly common in scientific
programming, and is one of the relatively few things missing from C used in
scientific programming.

However, requiring it as a recursive function does scream of homework. A
simple for loop should be easier and faster.

-- glen

 

[Arthur J. O'Dwyer]

<off-topic>
If computing positive integer powers is of so little
interest, it's hard to see why Knuth gives the topic an
entire section of its own in "The Art of Computer Programming,
Volume II: Seminumerical Algorithms."

The method outlined is the first serious power-computing
algorithm Knuth introduces in that section. He writes that
some authors have declared the method optimal, but he is
kind enough to spare them embarrassment by omitting their
names; clearly they forgot to consider cases like n==15.

What's wrong with n==15? The algorithm given is still O(lg N)
in such cases. No, probably Knuth was thinking of *better*
algorithms that could compute powers more quickly -- I don't
know of any, but perhaps a couple of exp() and log() lookup
tables, plus iteration, could do things faster than straight
recursion.

-Arthur

 

[Eric Sosman]

What's wrong with n==15? The algorithm given is still O(lg N)
in such cases. No, probably Knuth was thinking of *better*
algorithms that could compute powers more quickly -- I don't
know of any, but perhaps a couple of exp() and log() lookup
tables, plus iteration, could do things faster than straight
recursion.

The binary method starts with x and then computes x^2,
x^3, x^6, x^7, x^14, x^15 -- six multiplications.

The factor method starts with x and then computes x^2
and x^3 with two multiplications. Letting y = x^3 it then
computes y^2, y^4, y^5 (= x^15) in three more multiplications,
for five in all.

Perhaps not a big deal when `x' is something simple like
a floating-point number -- but worth paying heed to if `x'
is, say, a large matrix.

 

[James Hu]

The languages I know of supply routines for integer n, and integer, float,
double, complex, and complex double x.

Well, lets just call most of those floating point, except for integer
and complex x (which I assume you mean A+Bi with A and B integers, but
C has no such type).

The OP didn't supply the type of x, but the algorithm works for any type
where mulitply is defined.

Since the proposed algorithm does not deal with negative powers, it
is safe to assume the the type of n is unsigned.

That is an imaginative approach, but not what I had in mind.

[.. snip .. ]

The algorithm also works for n > 255.

But not without overflowing the int type. My function does not
prevent overflow from happening, but the caller can check to
make sure overflow will not before invoking the function. And
the caller can then discover the degenerate unity cases and
dispatch it without incurring the function call overhead.

However, here is an alternate implementation (now, with two
table lookups):

#include <stdint.h>
#include <limits.h>

static int32_t xpow[32] =
{0,INT_MAX,46341,1291,216,74,36,22,15,11,9,8,6,6,5,5,4,4,4,4,3,3
,3,3,3,3,3,3,3,3,3,2
};

static int8_t hbit[32] =
{-1
,0
,1,1
,2,2,2,2
,3,3,3,3,3,3,3,3
,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4
};

int32_t int_pow(int32_t x, unsigned n)
{
int32_t t = 1;
int32_t ax;
if (n == 0) return 1;
if (x == 0) return 0;
if (n == 1) return x;
if (x == INT_MIN) return 0;
if (((x<0)?-x:x) >= xpow[(n<32)?n:31])
return (x<0)?((n&1)?INT_MIN:INT_MAX):INT_MAX;
if (x == 1 || x == -1) n &= 1;
switch (hbit[n]) {
case 4: if (n & 1) t *= x; n >>= 1; x *= x;
case 3: if (n & 1) t *= x; n >>= 1; x *= x;
case 2: if (n & 1) t *= x; n >>= 1; x *= x;
case 1: if (n & 1) t *= x; n >>= 1; x *= x;
default: t *= x;
}
return t;
}


-- James

 

[James Hu]

{0,INT_MAX,46341,1291,216,74,36,22,15,11,9,8,6,6,5,5,4,4,4,4,3,3 [ ... ]
if (x == INT_MIN) return 0; [ ... ]
return (x<0)?((n&1)?INT_MIN:INT_MAX):INT_MAX;

[ ... ]

References to INT_MAX and INT_MIN should be changed to INT32_MAX
and INT32_MIN respectively.

-- James

 

[Richard Heathfield]

Computing integer powers of numbers is fairly common in scientific
programming, and is one of the relatively few things missing from C used
in scientific programming.

However, requiring it as a recursive function does scream of homework. A
simple for loop should be easier and faster.

Hmmm. I raised 7 to the power 7777 using recursion and iteration. (Since the
result occupies over 6500 decimal digits, I won't display it here!)

The recursive calculation took 0.22 seconds, and the iterative version 1.06
seconds - almost five times slower. Perhaps you could demonstrate an
iterative version that can hold a candle to the recursive technique?

 

[gc]

My take

double Power(double x, unsigned int n) {
if(n==0) return 1.;
else return (n&1) ? x*Power(x*x,(n-1)/2) : Power(x*x,n/2);
}

 

[Christian Bau]

What's wrong with n==15? The algorithm given is still O(lg N)
in such cases. No, probably Knuth was thinking of *better*
algorithms that could compute powers more quickly -- I don't
know of any, but perhaps a couple of exp() and log() lookup
tables, plus iteration, could do things faster than straight
recursion.

Knuth was thinking of the fact that this algorithm uses six
multiplications to calculate x^15, but it can be done using five
multiplications.

 

[Christian Bau]

Hmmm. I raised 7 to the power 7777 using recursion and iteration. (Since the
result occupies over 6500 decimal digits, I won't display it here!)

The recursive calculation took 0.22 seconds, and the iterative version 1.06
seconds - almost five times slower. Perhaps you could demonstrate an
iterative version that can hold a candle to the recursive technique?

double power (double x, unsigned long n) {

if (n == 0) {
return 1.0;
} else {
double result = x;
unsigned long k = 1;

while (2*k <= n) k *= 2;
while (k > 1) { k /= 2; result *= (n & k) ? result * x : result; }
return result;
}
}

But your example is much more difficult because the execution time of a
single multiplication depends on the size of the numbers involved, so it
is not the number of operations that counts, but their cost.

 

[Arthur J. O'Dwyer]

Hmmm. I raised 7 to the power 7777 using recursion and iteration. (Since
the result occupies over 6500 decimal digits, I won't display it here!)

Really. I guess you must have used 'long double', huh?

The recursive calculation took 0.22 seconds, and the iterative version 1.06
seconds - almost five times slower. Perhaps you could demonstrate an
iterative version that can hold a candle to the recursive technique?

#include <math.h>

double pow3(double x, unsigned int n)
{
do {
return exp(n*log(x));
} while (1);
}


-Arthur

 

[Richard Heathfield]

Really. I guess you must have used 'long double', huh?

No, I used bignum routines written in ISO C.

#include <math.h>

double pow3(double x, unsigned int n)
{
do {
return exp(n*log(x));
} while (1);
}

Ah, I don't think you were taking me entirely seriously. That's a shame,
because I was in fact asking a serious question.

 

[Nudge]

What's wrong with n==15? The algorithm given is still O(lg N)
in such cases. No, probably Knuth was thinking of *better*
algorithms that could compute powers more quickly -- I don't
know of any, but perhaps a couple of exp() and log() lookup
tables, plus iteration, could do things faster than straight
recursion.

(x << 4) - x

There. O(1)

/me runs for cover ^_^

 

[James Hu]

I am trying to write a recursive version of Power(x,n) that works by
[...]

... the C answer is use pow() (unless you are purposefully avoiding
floating point, in which case a table lookup is in order).

Ok, my previous versions contained minor snafus. This one is tested,
and profiled:

#include <stdint.h>
#include <limits.h>

static int32_t xpow[32] =
{0,INT32_MAX,46340,1290,215,73,35,21,14,10,8,7
,5,5,4,4,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,1
};

static int8_t hbit[32] =
{-1,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3
,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4
};

int32_t int_pow(int32_t x, unsigned n)
{
int32_t t = 1;
if (n == 0) return 1;
if (x == 0) return 0;
if (n == 1) return x;
if (x == INT32_MIN)
return (n&1)?x:INT32_MAX;
if (((x<0)?-x:x) > xpow[(n>31)?31:n])
return (x<0)?((n&1)?INT32_MIN:INT32_MAX):INT32_MAX;
if (n>31)
return (n&1)?x:1;
switch (hbit[n]) {
case 4: if (n & 1) t *= x; n >>= 1; x *= x;
case 3: if (n & 1) t *= x; n >>= 1; x *= x;
case 2: if (n & 1) t *= x; n >>= 1; x *= x;
case 1: if (n & 1) t *= x; n >>= 1; x *= x;
default: t *= x;
}
return t;
}

-- James

 

[James Hu]

Ah, I don't think you were taking me entirely seriously. That's a shame,
because I was in fact asking a serious question.

Richard, could you post your iterative version? Maybe we can help
you tweak it.

-- James

 

[Marcus Lessard]

Then maybe it is just me, but I didn't explain myself clearly. The
objective of calculating powers is all well and good and no doubt of use but
what confuses me is the highly specified nature of the algorithm: (from OP):

"...breaking n down into halves(where half of n=n/2), squaring
Power(x,n/2), and multiplying by x again if n was odd..."

Will this be more efficient? Or is it just the solution demanded by the
professor?

ML

"Eric Sosman"

 

[Tim Woodall]

Glen Herrmannsfeldt wrote:


Hmmm. I raised 7 to the power 7777 using recursion and iteration. (Since the
result occupies over 6500 decimal digits, I won't display it here!)

The recursive calculation took 0.22 seconds, and the iterative version 1.06
seconds - almost five times slower. Perhaps you could demonstrate an
iterative version that can hold a candle to the recursive technique?

I'm surprised there is much if any difference:

long mpow(long x, long e)
{
long r=1;

while(e)
{
if(e&1)
r*=x;
x*=x;
e>>=1;
}
return r;
}

Its fairly trivial to remove one multiplication at the cost of n if
tests (where n is the number of bits in e) and one multiplication can
become an assignment but other than that I don't think fewer
multiplications are possible using the square and multiply idiom. IIRC
Knuth goes into this in some detail and mentions some exponents where
this idiom does not use the minimum possible number of
multiplications. (e=15 being the smallest which can be evaluated as
mpow(mpow(x,3),5))

(Getting very OT now. squaring is particularly fast when very big
number are involved where FFT techniques become "efficient")

Tim.

 

[Wolfgang Riedel]

Hmmm. I raised 7 to the power 7777 using recursion and iteration. (Since the
result occupies over 6500 decimal digits, I won't display it here!)

The recursive calculation took 0.22 seconds, and the iterative version 1.06
seconds - almost five times slower. Perhaps you could demonstrate an
iterative version that can hold a candle to the recursive technique?

--
Richard Heathfield : (e-mail address removed)
"Usenet is a strange place." - Dennis M Ritchie, 29 July 1999.
C FAQ: http://www.eskimo.com/~scs/C-faq/top.html
K&R answers, C books, etc: http://users.powernet.co.uk/eton

6573 (!?)

 

[Arthur J. O'Dwyer]

No, I used bignum routines written in ISO C.

[I did figure as much; still, I don't think the original discussion
was meant to generalize to bignums. I thought we'd been talking
about 'foo(double, unsigned int)'.]

Ah, I don't think you were taking me entirely seriously. That's a shame,
because I was in fact asking a serious question.

Well, suppose someone goes and looks up ISO C implementations
of 'exp' and 'log', and plugs them into the above function in the
right places, then. I'd be willing to bet that there's at least
one efficient iterative method of computing exp(n) and log(n),
although maybe not to the appropriate precision for bignum work.
Certainly a bignum solution would need lots of support code.

For a serious solution, you could simulate the recursive solution
by looking at the individual bits of the exponent; but since your
own iterative solution didn't take too long, I'm guessing that
that's what you did already.

-Arthur