The power series for exp is exp(z)=sum(k=0;k<n;(z^k)%!k). This uses ^ but only for integer powers. The computation can be coded in J succinctly as:
exp=: 1 + +/ @ ($ %&(*/\) +/\@($&1)@[) 20 exp 1.2 3.32012 ^ 1.2 3.32012 (^ - 20&exp) 1.2 0 The problem statement requires that + * - % be applied only to scalar arguments. This seems gratuitous and unreasonable. If necessary, you can convert the +/, +/\, */\, $ to scalar operations using a loop. The left argument to exp is the number of terms to use, and for larger magnitude right arguments you need to use a larger number of terms. There is a way to determine the number of terms given the required precision and the magnitude of the argument, or you can just pick a great big number. (^ - 20&exp) _9 _1.51409 (^ - 40&exp) _9 _3.27306e_11 (^ - 80&exp) _9 1.12434e_14 (^ - 160&exp) _9 1.12434e_14 I'll think about the power series for ln by adapting the one I have in http://www.jsoftware.com/jwiki/Essays/Extended%20Precision%20Functions . On Mon, Aug 4, 2014 at 9:39 PM, 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
