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 】