c++: powering two ints correctly

[DaveNo1]

Howdy,

How do I power two integers to get the result as an int correctly and without gcc warnings?

It might be trivial or simple, but it's bothering me: I want to power two integers and the result should be an integer too. This works: i = pow(a, b)

But unfortunately gcc always gives a warning, that it is converting to int from double (i know pow returns a double). In some cases I got false results which could be corrected with i = round(pow(a,b)) both a and b were long long integers.

Thanks for advice.

David

[daniel.is.fischer]

You can get rid of the compiler warning by using an explicit cast, to avoid unnecessary errors also use round:
i = (long long int)round(pow(a,b));

[BjornEdstrom]

You can also write your own:

int pow(int b, int e);

Sometimes, if the exponent is large, a common optimization is to precompute the possible values. Say you need to compute x^11. There are only 8 possibilites that fits in an unsigned int, so you can save these in an array:

unsigned int pow11[] = {0, 1, 2048, 177147, 4194304, 48828125, 362797056, 1977326743};

[stijn263]

I always use the following power function when dealing with integer exponents. Time complexity is logarithmic in n, so it's quite fast

Code: Select all

template <class T, class U> T power(T x,U n) {
  if(n==0) return (T) (1);
  if(n%2==1) return x*power(x*x,(n-1)/2);
  return power(x*x,n/2);
}

Note that if x^n is long long, your first function variable must be long long too. If you're only using it for integers, it's much safer to use the following version:

Code: Select all

long long power(long long x, int n) {
  if(n==0) return 1;
  if(n%2==1) return x*power(x*x,(n-1)/2);
  return power(x*x,n/2);
}

It will always return the correct value given that abs(x^n) < 2^63

[DaveNo1]

Wow, thanks for the detailed replies.

Has anyone tested stijn263's implementation of the power function against the (long long int)round(pow(a,b)) alternative i.e. is one of the two faster than the other (I would say power(a,b) is faster, but I am not so shure if my testing conditions were ok)?

[daniel.is.fischer]

I'd expect (long long)round(pow(a,b)) to be rather slow, after all it is prepared to deal with nonintegral exponents. My naive testing supports that:
dafis@linux:~/EulerProblems/Testing> ./powTest 1000000
Round: 3.89s
ExpLog: 6.22s
Stijn: 2.65s
Loop: 1.12s

ExpLog is (long long)round(exp(b*log(a))), I'm a bit surprised that's so much worse than round(pow(a,b)). Stijn is Stijn's non-template function. Loop is the same binary exponentiation algorithm, but written in one loop instead of being recursive, you can see that function calls have significant overhead in C (probably in C++ too, but I don't know that language at all).

[rayfil]

Unless you have the newer 64 bit machine and appropriate compiler for it, working with simple 32-bit ints should always be faster than forcing the compiler to use 64-bit ints.

Therefore, you may find it faster to split the power up by half so that each half fits within an int and complete the power up by multiplying the two halves casting the result to a 64-bit int. (At the CPU level, this is a simple multiplication of two ints with the result in two 32-bit registers which are then stored in memory as a 64-bit int).

[stijn263]

Stijn: 2.65s
Loop: 1.12s

That's quite a big difference. Were you using a compiler that optimizes tail recursion?

Therefore, you may find it faster to split the power up by half so that each half fits within an int and complete the power up by multiplying the two halves casting the result to a 64-bit int

Did you have something like this in mind?

Code: Select all

int power2(int x, int n) {
  if(n==0) return 1;
  if(n%2==1) return x*power2(x*x,(n-1)/2);
  return power2(x*x,n/2);
}

long long power(long long x, int n) {
  if(n==0) return 1;
  if(n==1) return x;
  if(n%2==1) {long long y=power2(x,(n-1)/2); return x*y*y;}
  if(n%2==0) {long long y=power2(x,n/2); return y*y;}
}

[daniel.is.fischer]

Plain old gcc with -O, I don't know if that optimises tail recursion. But look what happens with -O3:
dafis@linux:~/EulerProblems/Testing> ./powTest3 1000000
Round: 3.73s
ExpLog: 6.00s
Stijn: 1.73s
Loop: 0.15s

w00t

[hk]

Wouldn't it be useful to mention exactly the test (what numbers are exponentiated and to what power and how often?
Then benchmarks can be run on different compilers and for several languages.

[rayfil]

stijn

Did you have something like this in mind?

I'm not entirely familiar with HLL notation but it does look like what I had in mind.

(You could also add the if(n==1) return x; instruction in your power2 proc.)

[Ikcelaks]

Stijn: 2.65s
Loop: 1.12s

That's quite a big difference. Were you using a compiler that optimizes tail recursion?

stijn263, your code is not fully tail-recursive (look at how you treat the odd exponent case), so I doubt that any non-genius compiler would optimize it to a loop.

It's not too hard to make it tail-recursive with an accumulator. Here is the sml code:

Code: Select all

fun pow' (a,1,e) = a*e
  | pow' (a,b,e) = pow' (a*a,b div 2,if b mod 2 = 0 then e else a*e)

fun pow (_,0) = 1
  | pow (a,b) = pow' (a,b,1)

And C (modified stijn263's code, but possibly incorrect, as I don't have a C compiler handy):

Code: Select all

int power2(int x, int n, int e) {
  if(n==1) return (x*e);
  if(n%2==1) return power2(x*x,(n-1)/2,x*e);
  return power2(x*x,n/2,e);
}

int power(int x, int n) {
  if(n==0) return 1;
  return power2(x,n,1);
}

[hk]

And the winner is Raymond.
For Raymond's proposal I did really simple.

Code: Select all

function raymond(x,n:integer):int64;
var h,i:integer;
begin
if n=0 then begin result:=1;exit end;
if n=1 then begin result:=x;exit end;
h:=x;
for i:=2 to n shr 1 do h:=h*x;
result:=int64(h)*h;
if odd(n) then result:=x*result;
end;

I invented my own benchmark and programmed several of the proposals in Delphi:

Code: Select all

procedure TForm1.Button1Click(Sender: TObject);
var i,j:integer;sum:int64;
start,stop,freq:int64;


function stijn(x:int64;n:integer):int64;
begin
 if n=0 then begin result:=1;exit end;
 if n=1 then begin result:=x;exit end;
 if n and 1=1 then begin result:=x*stijn(x*x,n shr 1);exit end;
 result:=stijn(x*x,n shr 1);
end;

function loop(x,n:integer):int64;
 var x2:int64;
 begin
   if n=0 then begin result:=1;exit end;
   if n=1 then begin result:=x;exit end;
   x2:=x;
   result:=1;
   while n>1 do
   begin
    if n and 1=1 then  result:=result*x2;
     x2:=x2*x2;
     n:=n shr 1
   end;
   result:=x2*result
 end;

function raymond(x,n:integer):int64;
var h,i:integer;
begin
if n=0 then begin result:=1;exit end;
if n=1 then begin result:=x;exit end;
h:=x;
for i:=2 to n shr 1 do h:=h*x;
result:=int64(h)*h;
if odd(n) then result:=x*result;
end;

function ikcelaks(x:int64;n:integer):int64;
  function ikcelaks2(x:int64;n:integer;e:int64):int64;
  begin
  if n=1 then  begin result:=x*e;exit end;
  if n and 1=1 then begin result:=ikcelaks2(x*x,n shr 1,x*e);exit end;
  result:=ikcelaks2(x*x,n shr 1,e)
  end;
begin
 if n=0 then begin result:=1 ;exit end;
 result:=ikcelaks2(x,n,1);
end;

begin
queryperformancecounter(start);
queryperformancefrequency(freq);
sum:=0;
for i:=1 to 100000 do
  for j:=0 to 39 do
 sum:=sum+trunc(power(3,j));
queryperformancecounter(stop);
memo1.lines.add('round '+inttostr(sum)+' '+floattostr((stop-start)/freq*1000));


queryperformancecounter(start);
queryperformancefrequency(freq);
sum:=0;
for i:=1 to 100000 do
  for j:=0 to 39 do
 sum:=sum+stijn(3,j);
queryperformancecounter(stop);
memo1.lines.add('stijn '+inttostr(sum)+' '+floattostr((stop-start)/freq*1000));


queryperformancecounter(start);
queryperformancefrequency(freq);
sum:=0;
for i:=1 to 100000 do
 for j:=0 to 39 do
 sum:=sum+loop(3,j);
queryperformancecounter(stop);
memo1.lines.add('loop '+inttostr(sum)+' '+floattostr((stop-start)/freq*1000));
queryperformancecounter(start);
queryperformancefrequency(freq);
sum:=0;
for i:=1 to 100000 do
 for j:=0 to 39 do
 sum:=sum+raymond(3,j);
queryperformancecounter(stop);
memo1.lines.add('raymond '+inttostr(sum)+' '+floattostr((stop-start)/freq*1000));

queryperformancecounter(start);
queryperformancefrequency(freq);
sum:=0;
for i:=1 to 100000 do
 for j:=0 to 39 do
 sum:=sum+ikcelaks(3,j);
queryperformancecounter(stop);
memo1.lines.add('ikcelaks '+inttostr(sum)+' '+floattostr((stop-start)/freq*1000));
end;

round 7715491895585597952 579.834587915503ms
stijn 7715491895585597952 472.62375525381ms
loop 7715491895585597952 398.608202997867ms
raymond 7715491895585597952 216.679773546638ms
ikcelaks 7715491895585597952 554.623003761651ms

[Ikcelaks]

Oh! I didn't know that this was a contest for computing large powers. Otherwise, I would have added the trick of limiting the 64-bit arrithmatic to the final squaring. As was, I just wanted to show how to make the function fully tail-recursive.

Also, does the Delphi compiler eliminate tail-recursive calls? I'm thinking that it doesn't. And why should it, when code of that style looks so ugly in that syntax?

[hk]

I don't think it does.
But it's not a big problem to write this recursion out by hand, doesn't it?
In "Data Structures Using Pascal" (1981) there is paid much attention to writing out recursion by hand.
A good exercise to do it yourself, instead of hoping that the compiler will.

Well, it's not a contest but, but working in 32 bits and doing the last step by a squaring makes for fast code.

About ugliness, that's a matter of taste, I presume.
BTW Delphi is not a C compiler, but a Pascal compiler, and Pacal is another language.

[Ikcelaks]

I wasn't criticizing Delphi's syntax in general. Pascal was my first real programming language, and I still have some fond memories. My point was that it isn't natural to write recursive functions in it. Nobody in their right mind would write a tail recursive Pascal function, because it would be much more clear to write the equivalent loops explicitly (not to mention that the tail recursion wouldn't be optimized).

In "Data Structures Using Pascal" (1981) there is paid much attention to writing out recursion by hand.
A good exercise to do it yourself, instead of hoping that the compiler will.

Absolutely, in an imperative language like Pascal, it would even be more natural to write it as a loop. However, since I typically write programs in a functional programming language, it's much more natural to use the tail-recursive definition (SML doesn't even have a primative for-loop) and know that it's a fundamental feature of the compiler to produce the proper machine code. It doesn't make sense to benchmark a tail-recursive function in a language that isn't designed for tail-recursion.

For example, my pow' function is essentailly equivalent (and my compiler knows that) to the following code using a while-loop (pseudo code):

Code: Select all

Function pow' (a,b,e)
  Begin
    While b <> 1 do
      if (b mod 2 = 1) then e := a*e;
      a := a*a;
      b := b div 2;
    endWhile
    return (a*e);
  End

[rayfil]

But unfortunately gcc always gives a warning, that it is converting to int from double (i know pow returns a double). In some cases I got false results which could be corrected with i = round(pow(a,b)) both a and b were long long integers.

One thing I forgot to mention. If you really needed a and b being long long integers (i.e. I assume you meant 64-bit integers), it would mean that the exponent could possibly exceed an unsigned 32-bit integer, in other words greater than 4294967295. Even 2 raised to that power would be an enormous number and could not be returned as an integer.

I would also mean that the number requiring to be raised to that power could also exceed an unsigned 32-bit integer. Such a number, even only squared, would immediately exceed 64 bits and could not be returned as an integer. Your compiler was thus well justified to issue a warning.

In addition, many compilers will automaticly use floats to compute a^b. If the compiler uses doubles (i.e. 32-bit floats), it can only return an integer with 23 bits of precision. Even using a "long double" (i.e. 64-bit float), you can only have 53 bits of precision. Thus, you should never rely entirely on your compiler's internal functions to return a precise integer with integers as input to a power function.

You must always take into account the limits of your declared data types and their implication.

[BjornEdstrom]

Some (verbose) special cases that may be useful:

Code: Select all

inline unsigned long long powx2(unsinged long long b) { return b*b; }
inline unsigned long long powx3(unsinged long long b) { return b*b*b; }
inline unsigned long long powx4(unsinged long long b) { register unsigned long long bb = b*b; return bb*bb; }
inline unsigned long long powx5(unsinged long long b) { register unsigned long long bb = b*b; return bb*bb*b; }

inline unsigned long long pow2x(unsinged long long e) { return 1ULL<<e; }
inline unsigned long long pow3x(unsigned int e) { const unsigned long long table[] = {1, 3, 9, 27, 81, 243, 729, ..., 12157665459056928801ULL}; return table[e]; }
inline unsigned long long pow4x(unsinged long long e) { return 1ULL<<(e<<1ULL); }
inline unsigned long long pow5x(unsigned int e) { const unsigned long long table[] = {1, 5, 25, 125, 625, 3125, ..., 1490116119384765625ULL}; return table[e]; }

Should work okay on 64 bit processors at least.

Edit: Code above also demonstrates nicely why C++ is so annoying. :-)