Calculating Non-Integer Exponent
Start with:
$$2^{3.1} = 2^3 2^{0.1} = 2^3 e^{0.1 \log{2}}$$
Now use a Taylor expansion, so that the above is approximately
$$2^3 \left [1+0.1 \log{2} + \frac{1}{2!} (0.1 \log{2})^2 + \frac{1}{3!} (0.1 \log{2})^3+\cdots + \frac{1}{n!} (0.1 \log{2})^n\right ] $$
wheer $n$ depends on the tolerance you require. In this case, if this error tolerance is $\epsilon$, then we want
$$\frac{1}{(n+1)!} (0.1 \log{2})^{n+1} \lt \epsilon$$
For example, if $\epsilon=5 \cdot 10^{-7}$, then $n=4$.
To answer a question that me from the past had :P The basic idea is that you could actually calculate both of these directly.
$$ 2^{3.1} = e^{2 \log\left( 3.1 \right)} $$
So really the procedure is simple:
- Take the natural log of 3.1
- Multiply by the base (2)
- Use the result as the argument to $\exp(x) $ or $e^x$
So this problem is really reduced to, how do calculators calculate logarithms and exponentials? The answer to that isn't so hard. I'm not giving the fastest methods by any means, there are some series that converge to these values much more quickly than the ones I've suggested, but this should get anybody who sees this in the future on track pretty quickly :).
Calculating the Log
To calculate the logarithm, a numerical method like the trapezoidal rule or Simpsons rule would be sufficient; because:
$$ \log (x) = \int_1^x { 1 \over t } dt $$
Calculating the Exponential
Here we could just use the direct limit definition:
$$ \exp (x) = \lim_{n\rightarrow\infty} \left(1 + x \over n \right)^n $$
Where n is just some large integer (so you just multiply directly), or you could again, approximate $e^x $ with a taylor series, which should converge to the actual answer a little more quickly, since:
$$ \exp (x) = \lim_{n\rightarrow \infty}\left( 1 + x + {x^2 \over 2!} + {x^3 \over 3!} + ... + {x^n \over n!} \right) $$
However, most calculators store some of these values in a LUT (look up table), since they will be (at least partially) pre-calculated with more accuracy than you are likely to need. This allows modern computers to find the exponential and logarithm of a function very quickly [although X86 - which is the instruction set used by most common computers today actually calculates $2^x$ and $log_2 (x)$, because its easy to calculate them in binary (base 2)].
I just wanted to directly calculate the value of the number $2^{3.1}$ as I was wondering how a computer would do it
(emphasis mine)
What follows is quite possibly one of the most detailed answers I've ever written on this network. I believe this thoroughly answers your question about "how a computer would do it", so it warrants a read.
This answer is broken into four main sections, titled:
- Figuring out the problem domain - where I show how to work out the problem by hand
- Getting to $8\sqrt[10]2$ programmatically - where I show how to get a computer to work out the problem
- Implementing root(exp, x) - describing how to implement nth root approximations
- An actual implementation - a full C implementation culminating in the calculation of $2^{3.1}$, among other things.
This isn't StackOverflow, so I'll try to respect that and refrain from going into too much technical detail. That said, any non-pseudocode example code I include will be in C, since C has a relatively low level of abstraction and thus an algorithm implemented in C is more concrete, and thus more easily ported to other languages.
Figuring out the problem domain
That said, the way I would approach this is by first representing the exponent as a fraction, factoring and simplifying that, and solving the simplified radical.
$ 2^{3.1} =\\ 2^{\frac{31}{10}} =\\ \sqrt[10]{2^{31}} =\\ \sqrt[10] {2^{10} * 2^{10} * 2^{10} * 2^1} =\\ 2 * 2 * 2 * \sqrt[10] 2 =\\ 8\sqrt[10] 2 $
Now we have something that can be represented in a program:
8 * root(10, 2)
Getting to $8\sqrt[10]2$ programmatically
Of course, going from pow(2, 3.1) to 8 * root(10, 2) with a machine that only understands how to perform addition and multiplication isn't exactly trivial, so you would probably want to separate the integer and decimal components; i.e. what you did.
#define GET_INT(n) \
( (double)((int)(n)) )
#define GET_DEC(n) \
( (n) - GET_INT(n) )
#define IS_INT(n) \
( GET_INT(n) == (n) )
Thus, $2^{3.1} =~...$
pow(a, b):
IS_INT(b):
pow(a, b)
otherwise:
bi = GET_INT(b)
bf = GET_DEC(b)
pow(a, bi) * root( get_denom(bf), pow(a, get_numer(bf)) )
pow(2, 3.0):
IS_INT(3.0)
8
pow(2, 3.2):
otherwise
bi = 3.0
bf = 0.2
pow(2, 3.0) * root(10, pow(2, 2))
where get_numer and get_denom are functions that take e.g. 0.2 and return 2 and 10, respectively. I'm not sure how to implement these two functions, but I know it's doable. For example, you might multiply bf by 10 until IS_INT(bf) is true, and then return bf for the numerator and $10^{the~number~of~times~you~multiplied~bf~by~10}$ for the denominator, but does this satisfactorily cover the problem domain and perform without unexpected floating-point related bugs, rounding errors or over/underflow? I cannot answer that. All I can say is that the approach should work in theory for a decent range of values. The implementation of our naive approach looks like this:
double *dec_to_frac(double x){
double magnitude = 1;
static double ret[2];
while(IS_INT(x) == 0){
x *= 10;
magnitude *= 10;
}
ret[0] = x;
ret[1] = magnitude;
return ret;
}
double get_numer(double x){
double *frac = dec_to_frac(x);
return frac[0];
}
double get_denom(double x){
double *frac = dec_to_frac(x);
return frac[1];
}
That said, we have now established an algorithm for calculating $a^b$ for float-compatible values > 0.
So now, for an input of $2^{3.1}$, you have the following C expression:
pow(a, bi) * root(get_denom(bf), pow(a, get_numer(bf)))
where a = 2.0, bf = 0.1, bi = 3.0
We have everything we need to calculate a good approximation, except for one function, root, which I will now go over.
Implementing root(exp, x)
First, we need to establish based on an example, so we'll start with cube roots.
Put simply, the goal is to turn the problem into an equation, turn the equation into a function $f(n)$, find the derivative of that function $f'(n)$, and then repeatedly calculate $n - \frac{f(n)}{f'(n)}$ where $n$ is the seed value or "guess". This will give you an approximation.
$ \sqrt[3]x = n~where~n^3 = x\\ n^3 = x\\ n^3 - x = 0\\ f(n) = n^3 - x\\ f'(n) = 3n^2\\ g(n) = n - \frac{f(n)}{f'(n)}\\ \sqrt[3]x \approx g^P(Q)~where~P~is~the~desired~precision~and~Q~is~the~initial~guess $
Generalizing this, we get
$ \sqrt[I]x = n~where~n^I = x\\ n^I = x\\ n^I - x = 0\\ f(n) = n^I - x\\ f'(n) = In^{I - 1}\\ g(n) = n - \frac{f(n)}{f'(n)}\\ \sqrt[I]x \approx g^P(Q) $
Which I believe is properly expressed as...
$ \sqrt[I]x \approx g^P(Q) ~where...\\ g(n) = n - \frac{f(n)}{f'(n)}\\ f(n) = n^I - x\\ f'(n) = In^{I - 1}\\ $
Now, we can implement
#define PRECISION 10
#define MAX(a, b) ( (a) > (b) ? (a) : (b) )
double root(int exp, double x){
double guess = MAX(1.0, x / exp);
for(int i = 0; i < PRECISION; ++i)
guess -= (pow(guess, (double)exp) - x) / ( exp * pow(guess, (double)(exp - 1)) );
return guess;
}
Putting all of that together into a specialized pow(double, double) function, with an input of (2.0, 3.1), you should get an output of roughly 8.574187700290345.
An actual implementation
What follows is my implementation, which I make no guarantees about being robust. I've also added some extra details to account for some complicated idiosyncrasies regarding how the computer represents floats.
#include <stdio.h>
#include <float.h>
#define PRECISION 20
#define MAX(a, b) \
( (a) > (b) ? (a) : (b) )
#define MIN(a, b) \
( (a) < (b) ? (a) : (b) )
#define GET_INT(n) \
( (double)((int)(n)) )
#define GET_DEC(n) \
( (n) - GET_INT(n) )
#define IS_INT(n) \
( MAX((n), GET_INT(n)) - MIN((n), GET_INT(n)) <= FLT_EPSILON )
int *dec_to_frac(double);
int get_numer(double);
int get_denom(double);
double fpow(double, double);
double ipow(double, int);
double root(int exp, double x);
int *
dec_to_frac(double x)
{
int magnitude = 1;
static int ret[2];
while(IS_INT(x) == 0){
x *= 10;
magnitude *= 10;
}
ret[0] = (int)x;
ret[1] = magnitude;
return ret;
}
int
get_numer(double x)
{
return (dec_to_frac(x))[0];
}
int
get_denom(double x)
{
return (dec_to_frac(x))[1];
}
double
fpow(double a, double b)
{
double bi;
double bf;
if(IS_INT(b))
return ipow(a, (int)b);
bi = GET_INT(b);
bf = GET_DEC(b);
return ipow(a, (int)bi) * root( get_denom(bf), ipow(a, get_numer(bf)) );
}
double
ipow(double a, int b)
{
double total = 1.0;
while(b--)
total *= a;
return total;
}
double
root(int exp, double x)
{
double guess = MAX(1.0, x / exp);
for(int i = 0; i < PRECISION; i++)
guess -= (ipow(guess, exp) - x) / ( exp * ipow(guess, exp - 1) );
return guess;
}
int
main(int argc, char **argv)
{
printf("pow(2, 3) = %.1f\n", fpow(2.0, 3.0));
printf("pow(2, 3.1) = %.15f\n", fpow(2.0, 3.1));
printf("root(2, 2) = %.15lf\n", root(2, 2.0));
printf("root(3, 3) = %.15lf\n", root(3, 3.0));
printf("root(2, 9) = %.15lf\n", root(2, 9.0));
printf("root(4, 81) = %.15lf\n", root(4, 81.0));
return 0;
}
Running this program generates the following output:
pow(2, 3) = 8.0
pow(2, 3.1) = 8.574187700290345
root(2, 2) = 1.414213562373095
root(3, 3) = 1.442249570307408
root(2, 9) = 3.000000000000000
root(4, 81) = 3.000000000000000