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
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
