[Jprogramming] Power for the powerless

[Dan Bron]

There's a StackExchange puzzle which challeges us to implement power (i.e.
dyad ^) using only the simple arithmetic dyads + - * % and monad | [1]. In
other words, we may not use ^ or ^. or variants.  There are still several
open questions on the puzzle, not least of which involves the domain of
the inputs (can the base be negative?) and range of the outputs (how much
precision is required?), but neverthless we can make some assumptions and
start to sketch an approach.

So, ignoring the "golf" aspect for now, and also (for the moment) any
clever but spirit-evading methods for using ^ without using ^ (including
other primitives which use ^ under the covers), how would you pursue this?
 

I took the naive approach based on the relation  x^y  <=> exp( y * ln(x) )
. Recursing into that, since we don't have ln(y) either, I took another
naive approach based on the relation ln(y) = lb(y) % lb(1x1) where lb is
2&^. . Recursing into that, since we don't have 2&^. either, I had to
resort to an algorithm I found on Wikipedia. Similar comments apply to
estimating exp(y) .  Anyway, I'm getting reasonable results, given that I
don't know what I'm doing:

	   10 pow 3     NB. Should be 1000
	999.761

	   1p1 ^ 1x1    NB. pi to the e
	22.4592
	   1p1 pow 1x1  NB. pretty close
	22.4581

	   1x1 ^ 1p1    NB. e to the pi
	23.1407
	   1x1 pow 1p1  NB. still pretty close
	23.1396

	   1x1 ^ _1p1   NB. negative exponents are ok too
	0.0432139
	   1x1 pow _1p1
	0.0432118

	   1p1 ^ _1x1
	0.0445253
	   1p1 pow _1x1
	0.0445231

So it seems like, precision questions aside, this approach can handle both
integral and rational bases and exponents, and the exponents can have any
sign.  On the other hand, while it handles positive bases ok, negative
bases are iffier (and base 0 just doesn't work):

	   _5 ^ 10      NB.  Sometimes they work ok...    
	9.76563e6
	   _5 pow 10
	9.75299e6

	   _5 ^ _10
	1.024e_7
	   _5 pow _10
	1.02267e_7

	   _10 ^ 3      NB.  Other times, 
	_1000
	   _10 pow 3    NB.  .. not so much
	999.761

This boils down to a difference between the real log, ^. and my substitute:

	   ^. _5        NB. ^. can return complex numbers
	1.60944j3.14159
	   ln _5        NB. but ln never will
	1.60944

Of course, I could start doing all sorts of fancy tricks (like switching to
a simpler definition when I detect an integral power, or something else
for a negative base) but I was hoping not to; I kind of wanted ^ to fall
out naturally from the mathematical equivalences (again, modulo
precision).  My code is below, for reference, but I'm wondering if I'm
overlooking something, or someone can suggest a potentially more fruitful
tack.

N         =:  64                            NB. Controls precision

pow       =:  exp@:(* ln)~                  NB.  x^y    = exp( y * ln(x) )
 exp      =:  1e5 spow~ 1 + %&1e5           NB.  exp(y) = lim[n -> _] (1 +
y%n)^n
 ln       =:  (lb 1x1) %~ lb                NB.  ln (y) = (lb y) % (lb 1x1)
  lb      =:  (] + lbF@:(% spow~&2)) lbI    NB.  lb (y) = n + lb(y * 2^-n)
where n =: <.2^.y
   spow   =:  */@:#~                        NB.  simple pow (positive
integer exponent)
   lbI    =:  _1 + (0,*/\N#2) I. |          NB.  log-base-2 of an integer
using lookup table
   lbF    =:  BFC +/ .* bits                NB.  (below) Binary Fraction
Coefficients
    BFC   =:  */\. ((N-1)#2) , %65536b10000
    bits  =:  (< (<0) ; 1) { bit^:(<1+N)    NB. See
http://en.wikipedia.org/wiki/Binary_logarithm#Real_number
    bit   =:  itBt@:*~@{.                   NB. Next iteration &
corresponding bit
     itBt =:  ((% {&1 2) , ]) (2&<)

-Dan

[1]  codegolf.SE "Raise to the Power":
     http://codegolf.stackexchange.com/questions/35732/raise-to-the-power