>>>>> "Raymond" == Raymond Toy <raymond.toy at ericsson.com> writes:
Raymond> 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)))
Raymond> I was getting ready to check in all the changes we had done last week,
Raymond> but I was looking at this expression. I must be stupid. I don't see
Raymond> how f(z) = 1/2*(log(z)-log(1/z)) - log(-z) simplifies that way. For
Raymond> real z, yes, I see f(z) is -%i*%pi for positive z and 0 for negative z.
Raymond> But for complex z, not on the real axis, I always get %i*%pi (or
Raymond> -%i*%pi, depending on how you want to write -1).
As best as I can tell f(z) = -%i*%pi for all z except the negative
real axis where it is 0.
However, if we use this version, we get weird results, assuming
expintegral_ei is continuous in the right half plane:
expintegral_ei(1.0) -> 1.895 - 6.283185307179586*%i
expintegral_ei(1.0+1d-8*%i) -> 1.895 - 3.14159*%i
expintegral_ei(1.0-1d-8*%i) -> 1.895 + 3.14159*%i
If we use the actual expression f(z), we get something that seems
reasonable:
expintegral_ei(1.0) -> 1.895
expintegral_ei(1.0+1d-8*%i) -> 1.895 - 2.71828d-8*%i
expintegral_ei(1.0-1d-8*%i) -> 1.895 + 2.71828d-8*%i
But this depends on expintegral-e returning the correct complex
values, and I'm not quite sure about that yet.
Ray