[Maxima] Generalized Lambert W function - premature commit
toy.raymond at gmail.com
Wed May 23 14:47:19 CDT 2012
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.
>>>> (%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
> 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
> 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.
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.)
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the Maxima