Stavros Macrakis wrote:
> I've checked in some code that improves Maxima's handling of
> sign(min/max(...)), sign(min/max(...)-min/max(...)). This isn't
> complete, but I thought I'd start with the simple cases. This already
> improves simplification of min/max . I'll have to write some test cases....
>
> I also checked in a very simple LambertW simplifier, using the name
> lambertw(...). The idea is to add it to solve sometime, but it would be
> nice if we could also add some numerical evaluation routines for it....
>
> -s
Here's some simple code to numerically evaluation Lambert W.
Wikipedia also lists
W(-log(2)/2) = -log(2)
W(-1/e) = -1
W(e) = 1
which we should probably also add.
Ray
;; Initial approximation for Lambert W.
;; http://www.desy.de/~t00fri/qcdins/texhtml/lambertw/
(defun init-lambert-w (x)
(if (<= x 500)
(let ((lx1 (log (1+ x))))
(+ (* .665 (+ 1 (* .0195 lx1)) lx1)
.04))
(- (log (- x 4))
(* (- 1 (/ (log x)))
(log (log x))))))
;; Algorithm based in part on
;; http://en.wikipedia.org/wiki/Lambert's_W_function. This can also
;; be found in
;; http://www.apmaths.uwo.ca/~djeffrey/Offprints/W-adv-cm.pdf, which
;; says the iteration is just Halley's iteration applied to w*exp(w).
(defun lambert-w (z &key (maxiter 100) (prec 1d-14))
(let ((w (init-lambert-w z)))
(dotimes (k maxiter)
(let* ((we (* w (exp w)))
(w1e (* (1+ w)
(exp w)))
(delta (/ (- we z)
(- w1e (/ (* (+ w 2)
(- we z))
(+ 2 (* 2 w)))))))
(when (<= (abs (/ delta w)) prec)
(return w))
(decf w delta)))
w))