Permalink
Browse files

expintegral_e(1,z) diverged for z with imaginary part larger

than 40. This is fixed by adding region of realpart(z)<0, abs(z)>2
and phase(z)<0.9. abs(z)>2 is intentionally different from
the semicircle for realpart(z)>0, whose radius is 1.
  • Loading branch information...
1 parent 8dd1602 commit 9cce7103e12df708f539940e1d956b7268b4e32a @YasuakiHonda YasuakiHonda committed Apr 20, 2012
Showing with 126 additions and 6 deletions.
  1. +11 −6 src/expintegral.lisp
  2. +115 −0 tests/rtest_expintegral.mac
View
@@ -531,7 +531,11 @@
(format t "~& : z = ~A~%" z))
(cond
- ((and (>= (realpart z) 0) (> (abs z) 1.0))
+ ((or (and (> (abs z) 2.0) (< (abs (phase z)) (* pi 0.9)))
+ ;; (abs z)>2.0 is necessary since there is a point
+ ;; -1.700598-0.612828*%i which 1<(abs z)<2, phase z < 0.9pi and
+ ;; still c-f expansion does not converge.
+ (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 +693,18 @@
(*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)))
+ ;; The same condition as you see in expintegral-e()
+ (and (>= (realpart flz) 0) (> (abs flz) 1.0)))
;; We expand in continued fractions.
(when *debug-expintegral*
(format t "~&We expand in continued fractions.~%"))
View
@@ -251,6 +251,59 @@ expintexpand:false;
false;
/*******************************************************************************
+ Do tests for additional float evaluation
+*******************************************************************************/
+
+test_value(
+ expintegral_e(1,-1.700598-0.612828*%i),
+ 1.229787425764198*%i-3.675726471068782,15);
+true;
+
+test_value(
+ expintegral_e(1,-1.5-%i*18.0),
+ 0.181696882955049 + 0.16898654452488*%i,15);
+true;
+
+test_value(
+ expintegral_e(1,-1.5-%i*180.0),
+ .01998793885396577-.01484463667769751*%i,15);
+true;
+
+test_value(
+ expintegral_e(1,-15.0-%i*50.0),
+ 62936.65453487506*%i-462.2396897671588,15);
+true;
+
+test_value(
+ expintegral_e(1,-15.0-%i*180.0),
+ 15302.16784585461-9678.322938932864*%i,11);
+true;
+
+test_value(
+ expintegral_e(1,-150.0-%i*540.0),
+ 2.479041623323267e+62*%i+2.105063890337474e+61,15-62);
+true;
+
+test_value(
+ expintegral_e(1,1.5+%i*18),
+ .01021327940204757-.006754558125326895*%i,15);
+true;
+
+test_value(
+ expintegral_e(1,15.0+%i*180.0),
+ 1.129118870333903e-9*%i+1.261123517801268e-9,15);
+true;
+
+test_value(
+ expintegral_e(1,150.0+%i*540.0),
+ 7.503677061484707e-69-1.036533879419174e-68*%i,15-69);
+true;
+
+
+
+
+
+/*******************************************************************************
Do tests for Bigfloat evaluation (values from functions.wolfram.com)
*******************************************************************************/
@@ -442,6 +495,68 @@ test_value(
63);
true;
+
+/* restore fpprec */
+(fpprec:oldfpprec, done);
+done;
+
+/*******************************************************************************
+ Do tests for additional big float evaluation
+*******************************************************************************/
+
+/* Remember actual fpprec */
+(oldfpprec:fpprec, fpprec:64, done);
+done;
+
+test_value(
+ expintegral_e(1,-1.5B0-%i*18.0B0),
+ 1.689865445248795461022355866185975300284302677290211138169639712b-1*%i
+ +1.816968829550486868009237590887639256845378383882275967782801257b-1,63);
+true;
+
+test_value(
+ expintegral_e(1,-15.0B0-%i*50.0B0),
+ 6.293665453487508854049428479450922789578087068032898765875144175b4*%i
+ -4.62239689767167158016168999268727094583699040537958778438547305b2,59);
+true;
+
+test_value(
+ expintegral_e(1,-1.5B0-%i*180.0B0),
+ 1.998793885396577954161433015440646790143523658773070463290940391b-2
+ -1.48446366776975235341069111717359420933661819609933627048078479b-2*%i,64);
+true;
+
+test_value(
+ expintegral_e(1,-15.0B0-%i*180.0B0),
+ 1.530216784585460659290957316900806622008438558612666623998533076b4
+ -9.678322938932862392009065659301165442754946683378713459319981218b3*%i,59);
+true;
+
+test_value(
+ expintegral_e(1,-150.0B0-%i*540.0B0),
+ 2.47904162332326778788226522449358779416653162097669293288256403b62*%i
+ +2.105063890337477483794346786409993155624883726578689067616685218b61,1);
+true;
+
+test_value(
+ expintegral_e(1,1.5B0+%i*18.0B0),
+ 1.021327940204755864240745590773899499558121729025330486880547885b-2
+ -6.754558125326884112691473105665013998699325818242672495399974279b-3*%i,65);
+true;
+
+test_value(
+ expintegral_e(1,15.0B0+%i*180.0B0),
+ 1.129118870333903101943908390287222351537018836722867882019975027b-9*%i
+ +1.261123517801267593713906360450736314138456246082714062225342743b-9,72);
+true;
+
+test_value(
+ expintegral_e(1,150.0B0+%i*540.0B0),
+ 7.50367706148470796758191535417179842322266910293477644515984948b-69
+ -1.03653387941917390868865781356557208741441722016128718318801276b-68*%i,-64-69);
+true;
+
+
/* restore fpprec */
(fpprec:oldfpprec, done);
done;

0 comments on commit 9cce710

Please sign in to comment.