Hi all,
I have created a patch for expintegral.lisp for improving the behavior
of expintegral-e(n,z) when imagpart(z) is negative. Basically it adds
additional region where the continued fraction is used rather than
power series. The same modification is made on bfloat-expintegral-e(n,z).
The patch passes all the existing test as well as newly added 20 tests
in the rtest_expintegral.mac. Also, 'make check' is performed to ensure
that all the tests are OK.
So, if everyone is OK with these changes, I will submit it to the git
server.
Thanks in advance for your comments!!
Yasuaki Honda, Chiba, Japan
% git diff expintegral.lisp
diff --git a/src/expintegral.lisp b/src/expintegral.lisp
index 6bf7d82..f25866f 100644
--- a/src/expintegral.lisp
+++ b/src/expintegral.lisp
@@ -531,7 +531,8 @@
(format t "~& : z = ~A~%" z))
(cond
- ((and (>= (realpart z) 0) (> (abs z) 1.0))
+ ((or (and (> (abs z) 2) (< (abs (phase z)) (* pi 0.9)))
+ (and (>= (realpart z) 0) (> (abs z) 1.0)))
;; We expand in continued fractions.
(when *debug-expintegral*
(format t "~&We expand in continued fractions.~%"))
@@ -689,17 +690,17 @@
(*expint-maxit* 5000) ; arbitrarily chosen, we need a better choice
(bigfloattwo (add bigfloatone bigfloatone))
(bigfloat%e ($bfloat '$%e))
- (bigfloat%gamma ($bfloat '$%gamma)))
+ (bigfloat%gamma ($bfloat '$%gamma))
+ (flz (complex ($float ($realpart z)) ($float ($imagpart z)))))
(when *debug-expintegral*
(format t "~&BFLOAT-EXPINTEGRAL-E called with:~%")
(format t "~& : n = ~A~%" n)
- (format t "~& : z = ~A~%" z))
+ (format t "~& : z = ~A~%" flz))
(cond
- ((and (or (eq ($sign ($realpart z)) '$pos)
- (eq ($sign ($realpart z)) '$zero))
- (eq ($sign (sub (cabs z) bigfloatone)) '$pos))
+ ((or (and (> (abs flz) 2) (< (abs (phase flz)) (* pi 0.9)))
+ (and (>= (realpart flz) 0) (> (abs flz) 1.0)))
;; We expand in continued fractions.
(when *debug-expintegral*
(format t "~&We expand in continued fractions.~%"))
2012?4?8?15:23 Raymond Toy <toy.raymond at gmail.com>:
> On 4/7/12 9:20 PM, ???? wrote:
> > Hi Raymond san,
> >
> > > Yes, this was basically what I had in mind, except the condition
> > would have been abs(z) > 1 and,
> > > maybe, abs(phase(z)) < 3.1, so that we stay away from the negative
> axis.
> >
> > The condition (and (> (abs z) 1) (< (abs (phase z)) (* pi 0.99))) did
> > not pass the rtest_expintegral. The condition that passes the
> > rtest_expintegral is:
> > (and (> (abs z) 2) (< (abs (phase z)) (* pi 0.9)))
> How does rtest_expintegral fail for the first case? Accuracy reduced?
> >
> > If you think this is OK, then I will modify bfloat and frac versions
> > of expintegral-e() in a same manner and create a patch for everyone's
> > review.
> If you think this is the right solution, I think you should just go
> ahead. If there's a problem we'll find it. :-) And for those that are
> interested, I'm sure they'll look at your patch.
> >
> > For the z on the negative real axis, I don't know the correct behavior
> > of the expintegral_e(). Just I experimented with couple of z values
> > such as z=-100, z=-200, z=-300 with expintegral_e(1, z) in Maxima and
> > e1(z) in mpmath in Python. They coincide with each other.
> Sorry about that. I was just thinking that evaluating the power series
> with terms like 200^k would not converge very well, but I didn't check.
> And the values are reasonably close to the asymptotic value of
> exp(-z)/z, so it looks good. Might be good to add a few test cases for z
> = -100 and -200.
>
> Ray
>
>