How can I calculate non-integer exponents?

Question

I can calculate the result of $x^y$ provided that $y \in\mathbb{N}, x \neq 0$ using a simple recursive function:

$$ f(x,y) = \begin {cases} 1 & y = 0 \\ (x)f(x, y-1) & y > 0 \end {cases} $$

or, perhaps more simply stated, by multiplying $x$ by itself $y$ times.

Unfortunately, I am unsure how I can numerically approximate $x^y$ for non-integer rationals.

For example, what method can I use to approximate 3^3.3?

If possible, I would like to be able to do this using only elementary arithmetic operations, i.e. addition, subtraction, multiplication, and division.

Answer 1

There is a concept using "fractional continued fractions" invented by D. Gomez Morin , see http://domingogomez.web.officelive.com/gcf.aspx . I've fiddled some time ago with this and think I recall right, that you can find roots of m'th order and that the apparent complexitiy of the "fractional continued fraction" reduces to rational fractions recursively applied. I'll try next days to recover the algorithm again, but perhaps the link to the page is already helpful.
Just found the hint of D.G.Morin to an excerpt of Steven Finch's book "mathematical constants", where S.F. mentions also this method in a much compacted way (see page 4 of http://assets.cambridge.org/052181/8052/sample/0521818052ws.pdf ).

---

[update] For the first broken link there is an entry in the wayback internet archive: Domingo Fomez Morin [/update]

Answer 2

I'll consider the problem of computing $x^\frac1{q}, \; q > 0$; as I've already mentioned in the comments, one can decompose any positive rational number as $m+\dfrac{p}{q}$, where $m,p$ are nonnegative integers, $q$ is a positive integer, and $p < q$. Thus for computing $x^{m+\frac{p}{q}}$, one could use binary exponentiation on $x^m$ and $\left(x^\frac1{q}\right)^p$ and multiply the results accordingly.

A.N. Khovanskiĭ, in his book on continued fractions, displays a continued fraction representation for the binomial function:

$$(1+z)^\alpha=1+\cfrac{2\alpha z}{2+(1-\alpha)z-\cfrac{(1-\alpha^2)z^2}{3(z+2)-\cfrac{(4-\alpha^2)z^2}{5(z+2)-\cfrac{(9-\alpha^2)z^2}{7(z+2)-\cdots}}}}$$

which converges for $|\arg(z+1)| < \pi$.

Letting $z=x-1$ and $\alpha=\dfrac1{q}$, one can then evaluate this continued fraction (with, say, Lentz-Thompson-Barnett) to generate a "seed" that can be subsequently polished with Newton-Raphson, Halley, or any of a number of iterations with high-order convergence. You'll have to experiment with how accurate a seed you need to start up the iteration, by picking a not-too-small tolerance when evaluating the continued fraction.

---

Here's some Mathematica code demonstrating what I've been saying earlier, for computing $\sqrt[3]{55}$:

With[{q = 3, t = 55, prec = 30},
 y = N[2 + (1 - 1/q) (t - 1), prec];
 c = y; d = 0; k = 1;
 While[True,
  u = (k^2 - q^-2) (t - 1)^2; v = (2 k + 1) (t + 1);
  c = v - u/c; d = 1/(v - u d);
  h = c*d; y *= h;
  If[Abs[h - 1] <= 10^-4, Break[]];
  k++];
 FixedPoint[
  Function[x, x ((1 + q) t - x^q (1 - q))/(x^q (1 + q) - (1 - q) t)], 
  1 + 2 (t - 1)/q/y]]

Here, I've arbitrarily chosen to stop when the continued fraction has already converged to $\approx 4$ digits, and then polished the result with Halley's method. The result here is good to $\approx 28$ digits. Again, you'll have to experiment on the accuracy versus expense of evaluating the "seed", as well as picking the appropriate iteration method for polishing the seed.