Generalized Lambert W function - premature commit

On Wed, May 23, 2012 at 2:06 AM, David Billinghurst <dbmaxima at gmail.com>wrote:

> 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.
>

Thanks for the updates.  Everything works now.

It's still a bit troubling that the bigfloat results aren't more
consistent.  Perhaps I'll try to track this down.  It could be a bug
somewhere.

>
> 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.


Looks good.

As an aside, the bigfloat package was intended so that you could use
exactly the same code for floats and bigfloats, ignoring issues with
constants and epsilon values.  If you run into a case where this is not
true, please let me know. (Of course, using the bigfloat package for floats
will be slower than just coding the same algorithm for floats.)

Ray