A similar approach but a bit golfed, noticing that the MacLaurin series for
exp and ln are similar.
M =: 1 : (':'; '(y,(u#y)*x*{:y)') NB. The series
P =: 1 : '+/y(u M)^:99]1' NB. It's sum
epx(y) is `% P y`
ln(1+y) is ` -@(%>:) P y'
So x ^ y is:
S =: 4 :'% P y * (*(-@(%>:)P))<: x'
1.2 S 1.7
1.36335
1.2 ^ 1.7
1.36335
0.2 ^ _3.14
156.591
0.2 S _3.14
156.591
1.2 ^ 170
2.88943e13
1.2 S 170
2.88943e13
but it won't work for x >~ 2. To make it work for any x, we could evaluate
the logarithm at a power of e which is closest to x.
On Tue, Aug 5, 2014 at 7:39 AM, Dan Bron <[email protected]> wrote:
> 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
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
