Thanks!
Yes, that works, though I need it for negative numbers
very close to zero, and then I get floating point overflow.
So I've written the following version:
lambert_w1(z,decimal_digits):= block([fpprec : decimal_digits, eps],
eps : 10^(-fpprec+1),
z : bfloat(z),
block(
[w1, w2 : 0, p : bfloat(-sqrt(2*%e*z+2))],
w1 : bfloat(((((((-221/8505*p)+769/17280)*p-43/540)*p+11/72)*p-1/3)*p+1)*p-1),
while w2 - w1 >= eps do (w2 : w1, w1 : bfloat(log(z/w1))),
return(w1)))$
which works, however is very slow for the inputs I might consider it
(say, z = -2^(-100000)).
Admittedly, the whole exercise was just about to find the smallest natural
number n such that a given k we have 2^n - n >= k.
And doing that just by running through n finishes in a couple of seconds
even for k = 10^10000, while using Lambert W you have to use
z = -2^(- 10^10000),
and that likely won't work ;-).
Was an interesting exercise nevertheless, and I guess that Lambert stuff
is useful anyway.
Thanks!
Oliver
On Thu, May 22, 2008 at 09:36:36PM -0400, Viktor T. Toth wrote:
> I believe that lambert_w calculates the principal branch, which is denoted
> W_0. What you are looking for is W_{-1}. Studying the paper by Corless, I
> don't think there's a readily usable relationship between W_0 and W_{-1}.
> However, using the paper of Chapeau-Blondeau and Monir, one can easily
> construct a simple numerical implementation of W_{-1}:
>
> W1(z):=block([w1,w2,p:-sqrt(2*float(%e)*z+2)],w1:((((((-221/8505*p)+769/1728
> 0)*p-43/540)*p+11/72)*p-1/3)*p+1)*p-1,while(w1#w2) do
> (w2:w1,w1:float(log(z/w1))),w1);
>
> Do give it a try, I tested it by verifying that W1(z)*exp(W1(z)) gives back
> z for -1/%e<z<0, so it seems to do the trick.
>
>
> Viktor
>
>
> -----Original Message-----
> From: maxima-bounces at math.utexas.edu [mailto:maxima-bounces at math.utexas.edu]
> On Behalf Of Oliver Kullmann
> Sent: Thursday, May 22, 2008 7:37 PM
> To: Barton Willis
> Cc: maxima at math.utexas.edu
> Subject: Re: [Maxima] Lambert W function?
>
> Perhaps I've missed it, but what
> branch of the Lambert-W-function does
> lambert_w represent?
>
> Apparently only this one selected branch
> is available? I need the function in
> the interval [-1/e, 0[, and I need
> exactly the other real-valued branch than
> the one computed by lambert_w. How
> to get that?
>
> Oliver
>
> On Tue, May 20, 2008 at 01:11:24PM -0500, Barton Willis wrote:
> > The function "lambert_w" is new to 5.15.0 (thanks due to Raymond and
> > Stavros).
> > The solve doesn't know about lambert_w. So far, lambert_w evaluates for
> > doubles (but not big floats):
> >
> > (%i5) lambert_w(5.60);
> > (%o5) 1.392014569377103
> > (%i6) diff(lambert_w(x),x);
> > (%o6) %e^(-lambert_w(x))/(lambert_w(x)+1)
> > (%i9) taylor(lambert_w(x),x,0,4);
> > (%o9) x-x^2+(3*x^3)/2-(8*x^4)/3+...
> >
> > Barton
> >
>
> --
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima