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
|
