Raymond Toy wrote:
> Now we can concentrate on why cmucl gets errors like this:
>
> ********************** Problem 143 ***************
> Input:
> test_table(lambda([z], expintegral_si(z)), 'si_2, 230, 1.95e-10)
>
>
> Result:
> error-catch
Ok. This is caused by the declaration in expintegral-e that says the
parameter is a (complex flonum). Cmucl (and sbcl) treats that as an
assertion, and expintegral_si calls expintegral-e with a fixnum. Here
is one simple fix. Not sure if this should be done in expintegral-e or
at a higher level.
diff -u -r1.5 expintegral.lisp
--- expintegral.lisp 23 Oct 2008 11:08:05 -0000 1.5
+++ expintegral.lisp 24 Oct 2008 19:47:29 -0000
@@ -531,10 +531,11 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun expintegral-e (n z)
- (declare (type integer n)
- (type (complex flonum) z))
+ (declare (type integer n))
(let ((*expint-eps* *expint-eps*)
- (*expint-maxit* *expint-maxit*))
+ (*expint-maxit* *expint-maxit*)
+ (z (coerce z '(complex flonum))))
+ (declare (type (complex flonum) z))
(when *debug-expintegral*
(format t "~&EXPINTEGRAL-E called with:~%")
Also, there are round off issues in expintegral-si and expintegral-ci.
This is caused by the expression (log (* (complex 0 -1) z)). z can be a
fixnum, so the log function must return a single-float result since the
argument is a complex rational.
One fix is to change (complex 0 1) to (complex 0d0 1d0). There are
several places where this change would need to be made. Alternatively,
before calling expintegral-ci and expintegral-si, we can make sure the
parameter is flonum or complex flonum. I don't know where the best
place would be. Maybe in the expintegral-ci/si routines themselves? At
a higher layer?
With all of these fixes, both cmucl and ecl pass rtest_expintegral. Hurray!
Ray