> 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 .
Many power series exists for ^.z (see Abramowitz & Stegun page 68 http://people.math.sfu.ca/~cbm/aands/page_68.htm). The problem is that for some of them convergence is very slow, even slow enough that the number of required terms and the available precision works against each other. The discussion in http://www.jsoftware.com/jwiki/Essays/Extended%20Precision%20Functions#Logarithm shows a good way to compute ^.y: ^.y satisfies the identities (^.a*b) = (^.a) + ^.b (^.a^b) = b * ^.a y can be written as r * 2^m where m is an integer and r is between %%:2 and %:2 , whence ^.y can be computed as (^.r)+m*^.2 . Equation 4.1.27 in Abramowitz & Stegun is a power series for ^.y : (^.y) = 2 * +/ i %~ ((y-1)%y+1) ^ i=. 1+2*i.n The series converge for all y , rapidly if y is between %%:2 and %:2 (between approximately 0.7 and 1.4). Equation 4.1.27 can be made readily to fit the requirements of the problem. The use of ^ and i. can be changed to */\ and +/\, respectively. The value m can be computed using a loop. (m is in fact part of the internal representation of y.) ^.2 is a constant. On Tue, Aug 5, 2014 at 7:45 AM, Roger Hui <[email protected]> wrote: > 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
