Subject: Generalized Lambert W function - premature commit
From: David Billinghurst
Date: Wed, 23 May 2012 19:06:47 +1000
On 23/05/2012 1:30 PM, Raymond Toy wrote:
> On 5/22/12 6:12 AM, David Billinghurst wrote:
>> On 22/05/2012 12:20 AM, David Billinghurst wrote:
>>> I have checked in a fix. It works for gcl under windows. Found
>>> another problem.
>>>
>>> fpprec:80$
>>> z:bfloat(-1/%e)-1b-20$
>>> (%i130) w:generalized_lambert_w(-1,z);
>>>
>>> Maxima encountered a Lisp error:
>>>
>>> Error in SETQ [or a callee]: Zero divisor.
>>>
>> Fixed this one too. The bigfloat guess for the Halley iteration is
>> generated in floating point. This doesn't work when z is -1/%e to
>> machine precision. Need to use a bigfloat power series around the
>> branch point.
> Pretty cool stuff!
>
> I found a few test failures with cmucl. I think they all stem for the
> use of (exp 1) in 3 places. This produces a single-float approximation
> to %e, but we really want a double-float value. I just used (exp 1.0)
> in those 3 places. (I'm not surprised that this works for gcl (because
> single-float = double-float there), but it should fail for other lisps.)
>
> Now the only test that fails is 44, where cmucl only achieves an
> accuracy of 3.89b-69 instead of 1b-77. Not sure why there should be
> this discrepancy because bigfloats should produce the same answer
> everywhere.
>
> Ray
Thanks for the report. I have replaced (exp 1) with %e-val.
I've had a look at the test 44 failure. The root finder is
ill-conditioned so close to the branch point, and loses about 15 digits
of accuracy. I have relaxed the test tolerance significantly.
I have also used the bigfloat package to remove some duplicate code and
there is scope for further improvement. I needed a bigfloat %e so I
added the method to numeric.lisp, based on the code for the %pi method.