Hello Ray,
you are right. I am stupid. In the complex case we have to add %i*%pi. The first
algorithm works with the testsuite. Thus, some tests for the correct complex
phase factors have to be added.
Today, I have tested this routine:
(defun expintegral-ei (z)
(+
(- (expintegral-e 1 (- z)))
; (- (* 0.5 (- (log z) (log (/ 1 z)))) (log (- z)))
(cond
((not (= (imagpart z) 0))
;; Complex value. Add phase %i*%pi.
(complex 0 (float pi)))
((> (realpart z) 0)
;; Positive real value. Add phase -%i*pi.
(complex 0 (- (float pi))))
;; Negative real value. No phase factor.
(t 0))))
To be complete here are the different cases and the results:
(%i6) expintegral_ei(1.0);
(%o6) 1.895117816355937
(%i7) expintegral_ei(-1.0);
(%o7) -0.21938393439552
(%i8) expintegral_ei(1.0+%i);
(%o8) 2.387769851510522*%i+1.764625985563854
(%i9) expintegral_ei(-1.0+%i);
(%o9) 2.962268118550434*%i-2.8162445198177954E-4
These values I have checked against wolfram.functions.com. Because the
Exponential Integral Ei has mirror symmetry the results for the arguments 1.0-%i
and -1.0-%i differ only by the sign of the imaginary part.
Later I will add this cases to the testsuite.
Dieter Kaiser
-----Urspr?ngliche Nachricht-----
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com]
Gesendet: Montag, 27. Oktober 2008 14:44
An: Dieter Kaiser
Cc: maxima at math.utexas.edu
Betreff: Re: [Maxima] Test suite failures with cvs version?
Dieter Kaiser wrote:
>
> (defun expintegral-ei (z)
> (+
> (- (expintegral-e 1 (- z)))
> ; (- (* 0.5 (- (log z) (log (/ 1 z)))) (log (- z)))
> (if (and (= (imagpart z) 0) (> (realpart z) 0))
> ;; Positive real value. Add a phase factor -%pi*%i
> (complex 0 (- (float pi)))
> ;; For all other values. No phase factor.
> 0)))
I was getting ready to check in all the changes we had done last week,
but I was looking at this expression. I must be stupid. I don't see
how f(z) = 1/2*(log(z)-log(1/z)) - log(-z) simplifies that way. For
real z, yes, I see f(z) is -%i*%pi for positive z and 0 for negative z.
But for complex z, not on the real axis, I always get %i*%pi (or
-%i*%pi, depending on how you want to write -1).
Help!
Ray