Subject: Generalized Lambert W function - premature commit
From: Raymond Toy
Date: Tue, 22 May 2012 20:30:17 -0700
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