Raise to the Power

Question

Challenge

Challenge is writing a program that takes a positive numbers a and a nonzero number b and outputs a^b (a raised to the power b). You can only use + - * / abs() as mathematical functions/operators. These can only be applied to scalar values, but not to whole lists or arrays.

Examples:

1.234 ^ 5.678 = 3.29980
4.5   ^ 4.5   = 869.874
4.5   ^-4.5   = 0.00114959

Relevant: http://xkcd.com/217/

Details

You can write a function or a similar construct for using in console. If you cannot use console input you can assume that both numbers are saved in variables and ouptut via standard output or writing to a file. The output has to be correct to at least 4 significat digits. You can assume that both a and b are nonzero. A runtime of significantly more than 1 minute is not acceptable. The least number of bytes will win. Please explain your program and your algorithm.

EDIT: Only positive bases have to be considered. You can assume a>0. Be aware that both numbers do not have to be integers!!!

Answer 1

JavaScript (E6) 155 174 191

Edit 2 As suggested by @bebe, using recursive function (performe worse but shorter)
Slightly changed R function to avoid 'too much recursion'
Test suite added. The function performs well for bases < 3000 and exponent in range -50..50.
Edit Golfed more and better precision

Any real number can be approximate with a rational number (and IEEE standard 'real' numbers store rationals in fact). Any rational number can be expressed as a fraction a/b with a and b integers. x^(a/b) is root b of (x^a) or (root b of x)^a. Integer exponentiation is quite easy by squaring. Integer root can be approximated using numerical methods.

Code

P=(x,e)=>(
  f=1e7,e<0&&(x=1/x,e=-e),
  F=(b,e,r=1)=>e?F(b*b,e>>1,e&1?r*b:r):r,
  R=(b,e,g=1,y=1e-30,d=(b/F(g,e-1)-g)/e)=>d>y|d<-y?R(b,e,g+d,y/.99):g,
  F(R(x,f),e*f)
)

Test In FireFox or FireBug console

for (i=0;i<100;i++)
{
  b=Math.random()*3000
  e=Math.random()*100-50
  p1=Math.pow(b,e) // standard power function, to check
  p2=P(b,e)
  d=(p1-p2)/p1 // relative difference
  if (!isFinite(p2) || d > 0.001) 
    console.log(i, b, e, p1, p2, d.toFixed(3))
}

Answer 2

Haskell, 85 90

Standard exp-log algorithm. Now with different name, shaving off a few more characters:

a%b|a>1=1/(1/a)%b|0<1=sum$scanl((/).((-b*foldr1(\n b->(1-a)*(b+1/n))c)*))1c
c=[1..99]

raise is now called (%), or % in infix notation, even making its use consume fewer bytes: 4.5%(-4.5)

The ungolfed version also uses just 172 bytes:

raise a b | a > 1     = 1 / raise (1/a) b
          | otherwise = expo (-b* ln (1-a))

ln x = foldr1 (\n a -> x*a+x/n) [1..99]

expo x = sum $ scanl ((/) . (x*)) 1 [1..99]

Answer 3

JS (ES6), 103 bytes

t=(x,m,f)=>{for(r=i=s=u=1;i<1<<7+f;r+=s/(u=i++*(f?1:u)))s*=m;return r};e=(a,b)=>t(b,t(a,1-1/a,9)*b-b,0)

Examples :

e(1.234,5.678) = 3.299798925315965
e(4.5,4.5)     = 869.8739233782269
e(4.5,-4.5)    = 0.0011495918812070608

Use Taylor series.
b^x = 1 + ln(b)*x/1! + (ln(b)*x)^2/2! + (ln(b)*x)^3/3! + (ln(b)*x)^4/4! + ...
with natural logarithm approximation :
ln(b) = (1-1/x) + (1-1/x)^2/2 + (1-1/x)^3/3 + (1-1/x)^4/4 + ...

I used 128 iterations to compute b^x (more iterations is difficult due to factorial) and 262144 iterations for ln(b)

Answer 4

Python, 77

As with some other answers this is based on log and exp. But the functions are computed by numerically solving ordinary differential equations.

def f(a,b,y=1):
 if a<1:a=1/a;b=-b
 while a>1:a/=1e-7+1;y*=b*1e-7+1
 return y

Does it satisfy the requirements? For the examples in the question, yes. For large a, it will take a very long time. For large a or b, it will become inaccurate.

Examples:

a            b            f(a, b)      pow(a, b)      <1e-5 rel error?
       1.234        5.678       3.2998       3.2998   OK
         4.5          4.5      869.873      869.874   OK
         4.5         -4.5   0.00114959   0.00114959   OK
         0.5          0.5     0.707107     0.707107   OK
         0.5         -0.5      1.41421      1.41421   OK
          80            5  3.27679e+09   3.2768e+09   OK
     2.71828      3.14159      23.1407      23.1407   OK

Update: flawr asked for more detail on the maths so here you go. I considered the following initial value problems:

• x'(t) = x(t), with x(0) = 1. The solution is exp(t).
• y'(t) = by(t), with y(0) = 1. The solution is exp(bt).

If I can find the value of t such that x(t) = a, then I'll have y(t) = exp(bt) = a^b. The simplest way to numerically solve an initial value problem is Euler's method. You compute the derivative the function is supposed to have, and then take a, step, in the direction of the derivative, and proportional to it, but scaled by a tiny constant. So that's what I do, take tiny steps until x is as big as a, and then see what y is at that time. Well, that's the way I thought of it. In my code, t is never explicitly computed (it is 1e-7 * the number of steps of the while loop), and I saved some characters by doing the calculations for x with a instead.

Answer 5

golflua 120

I use the fact that

a^b = exp(log(a^b)) = exp(b*log(a))

and wrote my own log & exp functions. The values a and b need to be entered on newlines when run in the terminal:

\L(x)g=0~@ i=1,50 c=(x-1)/x~@j=2,i c = c*(x-1)/x$g=g+c/i$~g$\E(x)g=1;R=1e7~@i=1,R g=g*(1+x/R)$~g$a=I.r()b=I.r()w(E(b*L(a)))

Sample runs:

4.5, 4.5  ==> 869.87104890175
4.5, -4.5 ==> 0.0011495904124065
3.0, 2.33 ==> 12.932794624815
9.0, 0.0  ==> 1
2.0, 2.0  ==> 3.9999996172672

An ungolfed Lua version is,

-- returns log
function L(x)
   g = 0
   for i=1,50 do
      c=(x-1)/x
      for j=2,i do
         c = c*(x-1)/x
      end
      g = g + c/i
   end
   return g
end

-- returns exp
function E(x)
   g=1;L=9999999
   for i=1,L do
      g=g*(1+x/L)
   end
   return g
end

a=io.read()
b=io.read()

print(E(b*L(a)))
print(a^b)