my pow function C/C++ (last update: 2012-09-24, created: 2012-09-23) back to the list ↑
Someone on an IRC channel mentioned something about implementing a pow (power) function in C/C++, so I decided to give it a quick try.
Note: While it's trivial to do a double pow(double a, int b) function, it's a little more work to do a double pow(double a, double b) function, where b can be something like 1.234.

There are a couple of ways to go with this, so this is just one of the methods I figured out (first one, not the best one). Let's start with some ASCII-art math:

x - the big unknown
a - base
b - exponent

     b
x = a

by using one of the logarithmic identities [1] I get:

     b * ln a
x = e
                                     
this is simpler to implement, since ln a can be approximated
with e.g. a taylor series or a series based on an area hyperbolic
tangent function [2]:

           inf               2n+1
           ---    1   ( z-1 ) 
ln a = 2 * >    ----- (-----)
           ---  2*n+1 ( z+1 )
           n=0

since the power used there is just a series of multiplications
we can calculate it in a trivial loop,

                       x
the last piece is the e  function, which by definition can be
expressed as [3]:

     inf     n
 x   ---    x
e  = >    -----
     ---    n!
     n=0

again, the power is trivial, since the exponent is an integer,
and the factorial is also just a series of multiplications

[1] http://en.wikipedia.org/wiki/List_of_logarithmic_identities#Using_simpler_operations
[2] http://en.wikipedia.org/wiki/Logarithm#Power_series
[3] http://en.wikipedia.org/wiki/Exponential_function#Formal_definition

The code in C++ (please note it was written in a couple of minutes, and it basically is low quality - I treated doubles as if they were real numbers, which they are not, and that alone introduces an error in the calculations; it's OK to look the code, but please please please do not use it in calculations, the results might be some poor Mars rover getting stuck on Mars again):

#include <stdio.h>
#include <math.h>

#define LN_STEP 20
#define EXP_STEP 13

double
my_spow(double a, int cnt) {
  double res = 1.0;
  while(cnt --> 0)
    res *= a;

  return res;
}

double
my_ln(double z) {
  double res = 0;
  int i;

  for(i = 0; i < LN_STEP; i++) {
    double n = (double)i;
    double p0 = 1.0 / (2.0 * n + 1.0);
    double p1 = (z - 1.0) / (z + 1.0);

    res += p0 * my_spow(p1, 2 * i + 1);
  }

  return 2.0 * res;
}

double
my_fact(int n) {
  double res = 1;
  int i;
  for(i = 2; i <= n; i++) {
    res *= (double)i;
  }

  return res;
}

double
my_exp(double x) {
  double res = 1.0 + x;
  int i;

  for(i = 2; i < EXP_STEP; i++) {
    // this is where the error happens for "big" numbers
    res += my_spow(x, i) / my_fact(i); 
  }

  return res;
}

double
my_pow(double a, double b) {
  return my_exp(b * my_ln(a));
}


int
main(void) {

  double nums[] = { 0.4, 2.2, 2.452, 1.23, 4.0, 2.0, 8.32, 0.1 };
  size_t nums_sz = sizeof(nums) / sizeof(nums[0]);
  size_t i, j;

  for(j = 0; j < nums_sz; j++) {
    for(i = 0; i < nums_sz; i++) {

      printf("%8.3f**%-8.3f ---> libc: %15.5f ||| my: %15.5f\n",
          nums[i], nums[j],
          pow(nums[i], nums[j]),
          my_pow(nums[i], nums[j]));
    }
  }

  return 0;
}


The approximation of power is quite decent for small numbers, and gets really bad for bigger. Here's the output:

   0.400**0.400    ---> libc:         0.69314 ||| my:         0.69314
   2.200**0.400    ---> libc:         1.37078 ||| my:         1.37078
   2.452**0.400    ---> libc:         1.43156 ||| my:         1.43156
   1.230**0.400    ---> libc:         1.08633 ||| my:         1.08633
   4.000**0.400    ---> libc:         1.74110 ||| my:         1.74110
   2.000**0.400    ---> libc:         1.31951 ||| my:         1.31951
   8.320**0.400    ---> libc:         2.33372 ||| my:         2.33372
   0.100**0.400    ---> libc:         0.39811 ||| my:         0.39811
   0.400**2.200    ---> libc:         0.13321 ||| my:         0.13321
   2.200**2.200    ---> libc:         5.66670 ||| my:         5.66670
   2.452**2.200    ---> libc:         7.19358 ||| my:         7.19358
   1.230**2.200    ---> libc:         1.57685 ||| my:         1.57685
   4.000**2.200    ---> libc:        21.11213 ||| my:        21.11172
   2.000**2.200    ---> libc:         4.59479 ||| my:         4.59479
   8.320**2.200    ---> libc:       105.74779 ||| my:       105.63012
   0.100**2.200    ---> libc:         0.00631 ||| my:         0.17599
   0.400**2.452    ---> libc:         0.10574 ||| my:         0.10575
   2.200**2.452    ---> libc:         6.91227 ||| my:         6.91226
   2.452**2.452    ---> libc:         9.01787 ||| my:         9.01787
   1.230**2.452    ---> libc:         1.66130 ||| my:         1.66130
   4.000**2.452    ---> libc:        29.93995 ||| my:        29.93825
   2.000**2.452    ---> libc:         5.47174 ||| my:         5.47174
   8.320**2.452    ---> libc:       180.36090 ||| my:       179.85576
   0.100**2.452    ---> libc:         0.00353 ||| my:         0.67730
   0.400**1.230    ---> libc:         0.32399 ||| my:         0.32399
   2.200**1.230    ---> libc:         2.63742 ||| my:         2.63742
   2.452**1.230    ---> libc:         3.01377 ||| my:         3.01377
   1.230**1.230    ---> libc:         1.28998 ||| my:         1.28998
   4.000**1.230    ---> libc:         5.50217 ||| my:         5.50217
   2.000**1.230    ---> libc:         2.34567 ||| my:         2.34567
   8.320**1.230    ---> libc:        13.54415 ||| my:        13.54400
   0.100**1.230    ---> libc:         0.05888 ||| my:         0.05899
   0.400**4.000    ---> libc:         0.02560 ||| my:         0.02833
   2.200**4.000    ---> libc:        23.42560 ||| my:        23.42497
   2.452**4.000    ---> libc:        36.14780 ||| my:        36.14430
   1.230**4.000    ---> libc:         2.28887 ||| my:         2.28887
   4.000**4.000    ---> libc:       256.00000 ||| my:       254.78292
   2.000**4.000    ---> libc:        16.00000 ||| my:        15.99989
   8.320**4.000    ---> libc:      4791.74066 ||| my:      4363.29401
   0.100**4.000    ---> libc:         0.00010 ||| my:       328.63746
   0.400**2.000    ---> libc:         0.16000 ||| my:         0.16000
   2.200**2.000    ---> libc:         4.84000 ||| my:         4.84000
   2.452**2.000    ---> libc:         6.01230 ||| my:         6.01230
   1.230**2.000    ---> libc:         1.51290 ||| my:         1.51290
   4.000**2.000    ---> libc:        16.00000 ||| my:        15.99989
   2.000**2.000    ---> libc:         4.00000 ||| my:         4.00000
   8.320**2.000    ---> libc:        69.22240 ||| my:        69.18925
   0.100**2.000    ---> libc:         0.01000 ||| my:         0.06040
   0.400**8.320    ---> libc:         0.00049 ||| my:        30.27515
   2.200**8.320    ---> libc:       706.24714 ||| my:       694.15292
   2.452**8.320    ---> libc:      1741.04612 ||| my:      1669.14017
   1.230**8.320    ---> libc:         5.59771 ||| my:         5.59771
   4.000**8.320    ---> libc:    102126.65979 ||| my:     64247.95490
   2.000**8.320    ---> libc:       319.57262 ||| my:       317.49769
   8.320**8.320    ---> libc:  45229946.57243 ||| my:   4805712.01950
   0.100**8.320    ---> libc:         0.00000 ||| my:   3096498.25791
   0.400**0.100    ---> libc:         0.91244 ||| my:         0.91244
   2.200**0.100    ---> libc:         1.08204 ||| my:         1.08204
   2.452**0.100    ---> libc:         1.09384 ||| my:         1.09384
   1.230**0.100    ---> libc:         1.02092 ||| my:         1.02092
   4.000**0.100    ---> libc:         1.14870 ||| my:         1.14870
   2.000**0.100    ---> libc:         1.07177 ||| my:         1.07177
   8.320**0.100    ---> libc:         1.23598 ||| my:         1.23598
   0.100**0.100    ---> libc:         0.79433 ||| my:         0.79433

Another way to go would be to focus on a being equal to b-th root of x, but that's for another time.

Update: Another time

A random idea by Arashi was to separate the integer part in the exponent from the rest. This allows one to trivially calculate the integer-part, and use another way to approximate the result of a to the non-integer power (aka the root).

Some ASCII-art math again:

x - the big unknown
a - base
b - exponent

let b_i and b_ni be the integer and the rest of the exponent,

the final result will be:

     b_i    b_ni
x = a    * a

                           b_i
calculating the result of a    is trivial,

     b_ni
the a     part can be expressed as:

                                                      P
 b_ni    b_ni1 / 10    b_ni2 / 100          b_niP / 10
a     = a           * a            * ... * a

where the b_ni1 ... b_niP are the P-th decimal digit after the radix,
going from the radix onwards, P being the desired precision (integer),
(please note that this can be easilly calculated by multiplying b_ni
times 10 and taking the integer part as b_niP, and then removing the
integer part by subtracting; see code below)

                    P
          b_niP / 10                      P             b_niP
now, the a            is actually a the 10 -th root of a
                                            P
(short note: calculating gcd of b_niP and 10 , and dividing both
by it, might speed things up a little, especially if you code is
as bad as my is)

     b_niP
the a      is trivial to calculate (since b_niP is integer)
the n-th root of A (x being the result) can be approximated e.g. using
the Newton's method approach [1]:

x  = whatever, take a guess
 0
        1  (                 A    )
x    = --- ( (n - 1) x  + ------- )
 k+1    n  (          k     n - 1 )
                           x
                            k

(the steps you make, the better precision you get)

[1] http://en.wikipedia.org/wiki/Nth_root_algorithm

Now the code and the results (some are OK, and other are totally off; I'm blaming the floats or/and Newton, but it might be an error in my code too):

#include <stdio.h>
#include <math.h>

#define DOWN_STEP 4
#define ROOT_STEP 8

double
my_spow(double a, int cnt) {

  double res = 1.0;
  while(cnt --> 0)
    res *= a;

  return res;
}

int
my_gcd(int a, int b) {
  if(a == b) return a;
  if(b < a) return my_gcd(a - b, b);
  return my_gcd(a, b - a);
}

double
my_sroot(double A, int n_int) {
  double n = (double)n_int;
  double x = A;

  for(int i = 0; i < ROOT_STEP; i++) { 
    x = (1.0 / n) *
        ( (n - 1.0) * x +
          (A / my_spow(x, n_int - 1)));
  }

  return x;
}


double
my_srootpow(double a, int p, int n) {
  int gcd = my_gcd(p, n);
  p /= gcd;
  n /= gcd;

  a = my_spow(a, p);
  a = my_sroot(a, n);

  return a;
}

double
my_pow(double a, double b) {

  int b_int = (int)b;
  b -= (double)b_int;

  double res = my_spow(a, b_int);

  int n = 10;
  for(int i = 0; i < DOWN_STEP; i++, n *= 10) {
    b *= 10.0;
    b_int = (int)b;
    b -= (double)b_int;

    if(b_int == 0)
      continue;

    res *= my_srootpow(a, b_int, n);
  }

  return res;
}


int
main(void) {

  double nums[] = { 0.4, 2.2, 2.452, 1.23, 4.0, 2.0, 8.32, 0.1 };
  size_t nums_sz = sizeof(nums) / sizeof(nums[0]);
  size_t i, j;

  for(j = 0; j < nums_sz; j++) {
    for(i = 0; i < nums_sz; i++) {

      printf("%8.3f**%-8.3f ---> libc: %15.5f ||| my: %15.5f\n",
          nums[i], nums[j],
          pow(nums[i], nums[j]),
          my_pow(nums[i], nums[j]));
    }
  }

  return 0;
}


   0.400**0.400    ---> libc:         0.69314 ||| my:        10.26685
   2.200**0.400    ---> libc:         1.37078 ||| my:         1.37095
   2.452**0.400    ---> libc:         1.43156 ||| my:         1.43717
   1.230**0.400    ---> libc:         1.08633 ||| my:         1.08633
   4.000**0.400    ---> libc:         1.74110 ||| my:         2.72139
   2.000**0.400    ---> libc:         1.31951 ||| my:         1.31951
   8.320**0.400    ---> libc:         2.33372 ||| my:        11.61406
   0.100**0.400    ---> libc:         0.39811 ||| my:     41943.04168
   0.400**2.200    ---> libc:         0.13321 ||| my:         0.13868
   2.200**2.200    ---> libc:         5.66670 ||| my:         5.66670
   2.452**2.200    ---> libc:         7.19358 ||| my:         7.19358
   1.230**2.200    ---> libc:         1.57685 ||| my:         1.57685
   4.000**2.200    ---> libc:        21.11213 ||| my:        21.11214
   2.000**2.200    ---> libc:         4.59479 ||| my:         4.59479
   8.320**2.200    ---> libc:       105.74779 ||| my:       112.20871
   0.100**2.200    ---> libc:         0.00631 ||| my:         0.41960
   0.400**2.452    ---> libc:         0.10574 ||| my:         1.#INF0
   2.200**2.452    ---> libc:         6.91227 ||| my:     25497.44499
   2.452**2.452    ---> libc:         9.01787 ||| my:    109456.41639
   1.230**2.452    ---> libc:         1.66130 ||| my:        13.04653
   4.000**2.452    ---> libc:        29.93995 ||| my: 120097870.12520
   2.000**2.452    ---> libc:         5.47174 ||| my:      7109.78994
   8.320**2.452    ---> libc:       180.36090 ||| my: 6991211176864.52250
   0.100**2.452    ---> libc:         0.00353 ||| my:         1.#INF0
   0.400**1.230    ---> libc:         0.32399 ||| my:         1.#INF0
   2.200**1.230    ---> libc:         2.63742 ||| my:   6965033.52228
   2.452**1.230    ---> libc:         3.01377 ||| my:  62272810.74654
   1.230**1.230    ---> libc:         1.28998 ||| my:        55.29688
   4.000**1.230    ---> libc:         5.50217 ||| my: 1223480976646.81910
   2.000**1.230    ---> libc:         2.34567 ||| my:   1015759.68243
   8.320**1.230    ---> libc:        13.54415 ||| my: 3453277764497018900.00000
   0.100**1.230    ---> libc:         0.05888 ||| my:         1.#INF0
   0.400**4.000    ---> libc:         0.02560 ||| my:         0.02560
   2.200**4.000    ---> libc:        23.42560 ||| my:        23.42560
   2.452**4.000    ---> libc:        36.14780 ||| my:        36.14780
   1.230**4.000    ---> libc:         2.28887 ||| my:         2.28887
   4.000**4.000    ---> libc:       256.00000 ||| my:       256.00000
   2.000**4.000    ---> libc:        16.00000 ||| my:        16.00000
   8.320**4.000    ---> libc:      4791.74066 ||| my:      4791.74066
   0.100**4.000    ---> libc:         0.00010 ||| my:         0.00010
   0.400**2.000    ---> libc:         0.16000 ||| my:         0.16000
   2.200**2.000    ---> libc:         4.84000 ||| my:         4.84000
   2.452**2.000    ---> libc:         6.01230 ||| my:         6.01230
   1.230**2.000    ---> libc:         1.51290 ||| my:         1.51290
   4.000**2.000    ---> libc:        16.00000 ||| my:        16.00000
   2.000**2.000    ---> libc:         4.00000 ||| my:         4.00000
   8.320**2.000    ---> libc:        69.22240 ||| my:        69.22240
   0.100**2.000    ---> libc:         0.01000 ||| my:         0.01000
   0.400**8.320    ---> libc:         0.00049 ||| my: 24404541157229269000000.00000
   2.200**8.320    ---> libc:       706.24714 ||| my:      4707.83158
   2.452**8.320    ---> libc:      1741.04612 ||| my:     17297.93556
   1.230**8.320    ---> libc:         5.59771 ||| my:         5.84202
   4.000**8.320    ---> libc:    102126.65979 ||| my:   6144245.73927
   2.000**8.320    ---> libc:       319.57262 ||| my:      1500.06304
   8.320**8.320    ---> libc:  45229946.57243 ||| my: 40292906546.66707
   0.100**8.320    ---> libc:         0.00000 ||| my: 8304435027255038600000000000000000000000000000000000000000000.00000
   0.400**0.100    ---> libc:         0.91244 ||| my:        73.15450
   2.200**0.100    ---> libc:         1.08204 ||| my:         1.08793
   2.452**0.100    ---> libc:         1.09384 ||| my:         1.13043
   1.230**0.100    ---> libc:         1.02092 ||| my:         1.02092
   4.000**0.100    ---> libc:         1.14870 ||| my:         1.72365
   2.000**0.100    ---> libc:         1.07177 ||| my:         1.07212
   8.320**0.100    ---> libc:         1.23598 ||| my:         3.58149
   0.100**0.100    ---> libc:         0.79433 ||| my:   4782969.04305

【 design & art by Xa / Gynvael Coldwind 】 【 logo font (birdman regular) by utopiafonts / Dale Harris 】