1;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2;;; 3;;; Exponential Integrals 4;;; 5;;; This file contains the following Maxima User functions: 6;;; 7;;; $expintegral_e (n,z) - Exponential Integral En(z) 8;;; $expintegral_e1 (z) - Exponential Integral E1(z) 9;;; $expintegral_ei (z) - Exponential Integral Ei(z) 10;;; 11;;; $expintegral_li (z) - Logarithmic Integral Li(z) 12;;; 13;;; $expintegral_si (z) - Exponential Integral Si(z) 14;;; $expintegral_ci (z) - Exponential Integral Ci(z) 15;;; 16;;; $expintegral_shi (z) - Exponential Integral Shi(z) 17;;; $expintegral_chi (z) - Exponential Integral Chi(z) 18;;; 19;;; $expint (x) - Exponential Integral E1(x) (depreciated) 20;;; 21;;; Global variables for the Maxima User: 22;;; 23;;; $expintrep - Change the representation of the Exponential Integral to 24;;; gamma_incomplete, expintegral_e1, expintegral_ei, 25;;; expintegral_li, expintegral_trig, expintegral_hyp 26;;; 27;;; $expintexpand - Expand the Exponential Integral E[n](z) 28;;; for half integral values in terms of Erfc or Erf and 29;;; for positive integers in terms of Ei 30;;; 31;;; The following features are implemented: 32;;; 33;;; 1. Numerical evaluation for complex Flonum and Bigfloat numbers 34;;; using an expansion in a power series or continued fractions. 35;;; The numerical support is fully implemented for the E[n](z) function. 36;;; All other functions call E[n](z) for numerical evaluation. 37;;; 38;;; 2. For a negative integer parameter E[n](z) is automatically expanded in 39;;; a finite series in terms of powers and the Exponential function. 40;;; 41;;; 3. When $expintexpand is set to TRUE or ERF E[n](z) expands 42;;; a) for n a half integral number in terms of Erfc (TRUE) or Erf (ERF) 43;;; b) for n a positive integer number in terms of Ei 44;;; 45;;; 3. Simplifications for special values: Ev(0), E[0](z), Li(0), Li(1),... 46;;; 47;;; 4. Derivatives of the Exponential Integrals 48;;; 49;;; 5. Change the representation of every Exponential Integral through other 50;;; Exponential Integrals or the Incomplete Gamma function. 51;;; 52;;; 6. Mirror symmetry for all functions and reflection symmetry for 53;;; the Exponential Inegral Si and Shi are implemented. 54;;; 55;;; 7. Handling of taylor expansions as argument. 56;;; 57;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 58;;; This library is free software; you can redistribute it and/or modify it 59;;; under the terms of the GNU General Public License as published by the 60;;; Free Software Foundation; either version 2 of the License, or (at 61;;; your option) any later version. 62;;; 63;;; This library is distributed in the hope that it will be useful, but 64;;; WITHOUT ANY WARRANTY; without even the implied warranty of 65;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 66;;; Library General Public License for more details. 67;;; 68;;; You should have received a copy of the GNU General Public License along 69;;; with this library; if not, write to the Free Software 70;;; Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 71;;; 72;;; Copyright (C) 2008 Dieter Kaiser 73;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 74 75(in-package :maxima) 76 77;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 78;;; Globals to help debugging the code 79;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 80 81(defvar *debug-expintegral* nil 82 "When enabled print debug information.") 83 84(defvar *debug-expint-maxit* 0 85 "When in debug mode count the maximum of iterations needed by the algorithm.") 86 87(defvar *debug-expint-fracmaxit* 0 88 "When in debug mode count the maximum of iterations needed by the algorithm.") 89 90(defvar *debug-expint-bfloatmaxit* 0 91 "When in debug mode count the maximum of iterations needed by the algorithm.") 92 93(defvar *debug-expint-fracbfloatmaxit* 0 94 "When in debug mode count the maximum of iterations needed by the algorithm.") 95 96;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 97;;; Globals for the Maxima Users 98;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 99 100(defvar $expintexpand nil 101 "When not nil we expand for a half integral parameter of the Exponetial 102 Integral in a series in terms of the Erfc or Erf function and for positive 103 integer in terms of the Ei function.") 104 105(defvar $expintrep nil 106 "Change the representation of the Exponential Integral. 107 Values are: gamma_incomplete, expintegral_e1, expintegral_ei, 108 expintegral_li, expintegral_trig, expintegral_hyp.") 109 110;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 111;;; Global to this file 112;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 113 114(defvar *expintflag* '(%expintegral_e1 %expintegral_ei %expintegral_li 115 $expintegral_trig $expintegral_hyp %gamma_incomplete) 116 "Allowed flags to transform the Exponential Integral.") 117 118(defun simp-domain-error (&rest args) 119 (if errorsw 120 (throw 'errorsw t) 121 (apply #'merror args))) 122 123;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 124;;; 125;;; Part 1: The implementation of the Exponential Integral En 126;;; 127;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 128 129(defmfun $expintegral_e (v z) 130 (simplify (list '(%expintegral_e) v z))) 131 132;;; Set properties to give full support to the parser and display 133 134(defprop $expintegral_e %expintegral_e alias) 135(defprop $expintegral_e %expintegral_e verb) 136 137(defprop %expintegral_e $expintegral_e reversealias) 138(defprop %expintegral_e $expintegral_e noun) 139 140;;; Exponential Integral E is a simplifying function 141 142(defprop %expintegral_e simp-expintegral-e operators) 143 144;;; Exponential Integral E distributes over bags 145 146(defprop %expintegral_e (mlist $matrix mequal) distribute_over) 147 148;;; Exponential Integral E has mirror symmetry, 149;;; but not on the real negative axis. 150 151(defprop %expintegral_e conjugate-expintegral-e conjugate-function) 152 153(defun conjugate-expintegral-e (args) 154 (let ((n (first args)) 155 (z (second args))) 156 (cond ((off-negative-real-axisp z) 157 ;; Definitely not on the negative real axis for z. Mirror symmetry. 158 (take '(%expintegral_e) 159 (take '($conjugate) n) 160 (take '($conjugate) z))) 161 (t 162 ;; On the negative real axis or no information. Unsimplified. 163 (list '($conjugate simp) 164 (take '(%expintegral_e) n z)))))) 165 166;;; Differentiation of Exponential Integral E 167 168(defprop %expintegral_e 169 ((n z) 170 ;; The derivative wrt the parameter n is expressed in terms of the 171 ;; Regularized Hypergeometric function 2F2 (see functions.wolfram.com) 172 ((mplus) 173 ((mtimes) -1 174 (($hypergeometric_regularized) 175 ((mlist) 176 ((mplus) 1 ((mtimes) -1 n)) 177 ((mplus) 1 ((mtimes) -1 n))) 178 ((mlist) 179 ((mplus) 2 ((mtimes) -1 n)) 180 ((mplus) 2 ((mtimes) -1 n))) 181 ((mtimes) -1 z)) 182 ((mexpt) 183 ((%gamma) ((mplus) 1 ((mtimes) -1 n))) 2)) 184 ((mtimes) 185 ((%gamma) ((mplus) 1 ((mtimes) -1 n))) 186 ((mexpt) z ((mplus) -1 n)) 187 ((mplus) 188 ((mtimes) -1 189 ((mqapply) 190 (($psi array) 0) 191 ((mplus) 1 ((mtimes) -1 n)))) 192 ((%log) z)))) 193 194 ;; The derivative wrt the argument of the function 195 ((mtimes) -1 ((%expintegral_e) ((mplus) -1 n) z))) 196 grad) 197 198;;; Integral of Exponential Integral E 199 200(defprop %expintegral_e 201 ((n z) 202 nil 203 ((mtimes) -1 ((%expintegral_e) ((mplus) 1 n) z))) 204 integral) 205 206;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 207 208;;; We support a simplim%function. The function is looked up in simplimit and 209;;; handles specific values of the function. 210 211(defprop %expintegral_e simplim%expintegral_e simplim%function) 212 213(defun simplim%expintegral_e (expr var val) 214 ;; Look for the limit of the arguments. 215 (let ((a (limit (cadr expr) var val 'think)) 216 (z (limit (caddr expr) var val 'think))) 217 (cond ((and (onep1 a) 218 (or (eq z '$zeroa) 219 (eq z '$zerob) 220 (zerop1 z))) 221 ;; Special case order a=1 222 '$inf) 223 224 ((member ($sign (add ($realpart a) -1)) '($neg $nz $zero)) 225 ; realpart of order < 1 226 (cond ((eq z '$zeroa) 227 ;; from above, always inf 228 '$inf) 229 ((eq z '$zerob) 230 ;; this can be further improved to give a directed infinity 231 '$infinity) 232 ((zerop1 z) 233 ;; no direction, return infinity 234 '$infinity) 235 (t 236 ($expintegral_e a z)))) 237 (t 238 ;; All other cases are handled by the simplifier of the function. 239 ($expintegral_e a z))))) 240 241;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 242 243(defun simp-expintegral-e (expr ignored z) 244 (declare (ignore ignored)) 245 (twoargcheck expr) 246 (let ((order (simpcheck (cadr expr) z)) 247 (arg (simpcheck (caddr expr) z)) 248 (ratorder)) 249 (cond 250 ;; Check for special values 251 ((or (eq arg '$inf) 252 (alike1 arg '((mtimes) -1 $minf))) 253 ;; arg is inf or -minf, return zero 254 0) 255 256 ((zerop1 arg) 257 (let ((sgn ($sign (add ($realpart order) -1)))) 258 (cond 259 ((eq sgn '$pos) 260 ;; we handle the special case E[v](0) = 1/(v-1), for realpart(v)>1 261 (inv (add order -1))) 262 ((member sgn '($neg $nz $zero)) 263 (simp-domain-error 264 (intl:gettext 265 "expintegral_e: expintegral_e(~:M,~:M) is undefined.") 266 order arg)) 267 (t (eqtest (list '(%expintegral_e) order arg) expr))))) 268 269 ((or (and (symbolp order) (member order infinities)) 270 (and (symbolp arg) (member arg infinities))) 271 ;; order or arg is one of the infinities, we return a noun form, 272 ;; but we have already handled the special value inf for arg. 273 (eqtest (list '(%expintegral_e) order arg) expr)) 274 275 ((and (numberp order) (integerp order)) 276 ;; The parameter of the Exponential integral is an integer. For this 277 ;; case we can do further simplifications or numerical evaluation. 278 (cond 279 ((= order 0) 280 ;; Special case E[0](z) = %e^(-z)/z, z<>0 281 ;; The case z=0 is already handled. 282 (div (power '$%e (mul -1 arg)) arg)) 283 284 ((= order -1) 285 ;; Special case E[-1](0) = ((z+1)*%e^(-z))/z^2, z<>0 286 ;; The case z=0 is already handled. 287 (div (mul (power '$%e (mul -1 arg)) (add arg 1)) (mul arg arg))) 288 289 ((< order -1) 290 ;; We expand in a series, z<>0 291 (mul 292 (factorial (- order)) 293 (power arg (+ order -1)) 294 (power '$%e (mul -1 arg)) 295 (let ((index (gensumindex))) 296 (dosum 297 (div (power arg index) 298 (take '(mfactorial) index)) 299 index 0 (- order) t)))) 300 301 ((and (> order 0) 302 (complex-float-numerical-eval-p arg)) 303 ;; Numerical evaluation for double float real or complex arg 304 ;; order is an integer > 0 and arg <> 0 for order < 2 305 (let ((carg (complex ($float ($realpart arg)) 306 ($float ($imagpart arg))))) 307 (complexify (expintegral-e order carg)))) 308 309 ((and (> order 0) 310 (complex-bigfloat-numerical-eval-p arg)) 311 ;; Numerical evaluation for Bigfloat real or complex arg. 312 (let* (($ratprint nil) 313 (carg (add ($bfloat ($realpart arg)) 314 (mul '$%i ($bfloat ($imagpart arg))))) 315 (result (bfloat-expintegral-e order carg))) 316 (add ($realpart result) (mul '$%i ($imagpart result))))) 317 318 ((and $expintexpand (> order 0)) 319 ;; We only expand in terms of the Exponential Integral Ei 320 ;; if the expand flag is set. 321 (sub (mul -1 322 (power (mul -1 arg) (- order 1)) 323 (inv (factorial (- order 1))) 324 (add (take '(%expintegral_ei) (mul -1 arg)) 325 (mul (inv 2) 326 (sub (take '(%log) (mul -1 (inv arg))) 327 (take '(%log) (mul -1 arg)))) 328 (take '(%log) arg))) 329 (mul (power '$%e (mul -1 arg)) 330 (let ((index (gensumindex))) 331 (dosum 332 (div (power arg (add index -1)) 333 (take '($pochhammer) (- 1 order) index)) 334 index 1 (- order 1) t))))) 335 336 ((eq $expintrep '%gamma_incomplete) 337 ;; We transform to the Incomplete Gamma function. 338 (mul (power arg (- order 1)) 339 (take '(%gamma_incomplete) (- 1 order) arg))) 340 341 (t 342 (eqtest (list '(%expintegral_e) order arg) expr)))) 343 344 ((complex-float-numerical-eval-p order arg) 345 (cond 346 ((and (setq z (integer-representation-p order)) 347 (plusp z)) 348 ;; We have a pure real positive order and the realpart is a float 349 ;; representation of an integer value. 350 ;; We call the routine for an integer order. 351 (let ((carg (complex ($float ($realpart arg)) 352 ($float ($imagpart arg))))) 353 (complexify (expintegral-e z carg)))) 354 (t 355 ;; The general case, order and arg are complex or real. 356 (let ((corder (complex ($float ($realpart order)) 357 ($float ($imagpart order)))) 358 (carg (complex ($float ($realpart arg)) 359 ($float ($imagpart arg))))) 360 (complexify (frac-expintegral-e corder carg)))))) 361 362 ((complex-bigfloat-numerical-eval-p order arg) 363 (cond 364 ((and (setq z (integer-representation-p order)) 365 (plusp z)) 366 ;; We have a real positive order and the realpart is a Float or 367 ;; Bigfloat representation of an integer value. 368 ;; We call the routine for an integer order. 369 (let* (($ratprint nil) 370 (carg (add ($bfloat ($realpart arg)) 371 (mul '$%i ($bfloat ($imagpart arg))))) 372 (result (bfloat-expintegral-e z carg))) 373 (add ($realpart result) 374 (mul '$%i ($imagpart result))))) 375 (t 376 ;; the general case, order and arg are bigfloat or complex bigfloat 377 (let* (($ratprint nil) 378 (corder (add ($bfloat ($realpart order)) 379 (mul '$%i ($bfloat ($imagpart order))))) 380 (carg (add ($bfloat ($realpart arg)) 381 (mul '$%i ($bfloat ($imagpart arg))))) 382 (result (frac-bfloat-expintegral-e corder carg))) 383 (add ($realpart result) 384 (mul '$%i ($imagpart result))))))) 385 386 ((and $expintexpand 387 (setq ratorder (max-numeric-ratio-p order 2))) 388 ;; We have a half integral order and $expintexpand is not NIL. 389 ;; We expand in a series in terms of the Erfc or Erf function. 390 (let ((func (cond ((eq $expintexpand '%erf) 391 (sub 1 ($erf (power arg '((rat simp) 1 2))))) 392 (t 393 ($erfc (power arg '((rat simp) 1 2))))))) 394 (cond 395 ((= ratorder 1/2) 396 (mul (power '$%pi '((rat simp) 1 2)) 397 (power arg '((rat simp) -1 2)) 398 func)) 399 ((= ratorder -1/2) 400 (add 401 (mul 402 (power '$%pi '((rat simp) 1 2)) 403 (inv (mul 2 (power arg '((rat simp) 3 2)))) 404 func) 405 (div (power '$%e (mul -1 arg)) arg))) 406 (t 407 (let ((n (- ratorder 1/2))) 408 (mul 409 (power arg (sub n '((rat simp) 1 2))) 410 (add 411 (mul func (take '(%gamma) (sub '((rat simp) 1 2) n))) 412 (mul 413 (power '$%e (mul -1 arg)) 414 (let ((index (gensumindex))) 415 (dosum 416 (div 417 (power arg (add index '((rat simp) 1 2))) 418 (take '($pochhammer) 419 (sub '((rat simp) 1 2) n) 420 (add index n 1))) 421 index 0 (mul -1 (add n 1)) t))) 422 (mul -1 423 (power '$%e (mul -1 arg)) 424 (let ((index (gensumindex))) 425 (dosum 426 (div 427 (power arg (add index '((rat simp) 1 2))) 428 (take '($pochhammer) 429 (sub '((rat simp) 1 2) n) 430 (add index n 1))) 431 index (- n) -1 t)))))))))) 432 433 ((eq $expintrep '%gamma_incomplete) 434 ;; We transform to the Incomplete Gamma function. 435 (mul (power arg (sub order 1)) 436 (take '(%gamma_incomplete) (sub 1 order) arg))) 437 438 (t 439 (eqtest (list '(%expintegral_e) order arg) expr))))) 440 441;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 442;;; Numerical evaluation of the Exponential Integral En(z) 443;;; 444;;; The following numerical routines are implemented: 445;;; 446;;; expintegral-e - n positive integer, z real or complex 447;;; frac-expintegral-e - n,z real or complex; n not a positive integer 448;;; bfloat-expintegral-e - n positive integer, 449;;; z Bigfloat or Complex Bigfloat 450;;; frac-bfloat-expintegral-e - n Bigfloat, z Bigfloat or Complex Bigfloat 451;;; 452;;; The algorithm are implemented for full support of Flonum and Bigfloat real 453;;; or Complex parameter and argument of the Exponential Integral. 454;;; Because we have no support for Complex Bigfloat arguments of the Gamma 455;;; function the evaluation for a Complex Bigfloat parameter don't give 456;;; the desiered accuracy. 457;;; 458;;; The flonum versions return a CL complex number. The Bigfloat versions 459;;; a Maxima Complex Bigfloat number. It is assumed that the calling routine 460;;; check the values. We don't handle any special case. This has to be done by 461;;; the calling routine. 462;;; 463;;; The evaluation uses an expansion in continued fractions for arguments with 464;;; realpart(z) > 0 and abs(z)> 1.0 (A&S 5.1.22). This expansion works for 465;;; every Real or Complex numbers including Bigfloat numbers for the parameter n 466;;; and the argument z: 467;;; 468;;; 1 n 1 n+1 2 469;;; En(z) = e^(-z) * ( --- --- --- --- --- ... ) 470;;; z+ 1+ z+ 1+ z+ 471;;; 472;;; The continued fraction is evaluated by the modified Lentz's method 473;;; for the more rapidly converging even form. 474;;; 475;;; For the parameter n an positive integer we do an expansion in a power series 476;;; (A&S 5.1.12): 477;;; inf 478;;; === 479;;; (-z)^(n-1) \ (-z)^m 480;;; En(z) = --------- * (-log(z)+psi(n)) * > ---------- ; n an integer 481;;; (n-1)! / (m-n+1)*m! 482;;; === 483;;; m=0 (m <> n-1) 484;;; 485;;; For an parameter n not an integer we expand in the following series 486;;; (functions.wolfram.com ): 487;;; inf 488;;; === 489;;; \ (-1)^m * z^m 490;;; Ev(z) = gamma(1-v) * z^(v-1) * > ------------- ; n not an integer 491;;; / (m-v+1)*m! 492;;; === 493;;; m=0 494;;; 495;;; The evaluation stops if an accuracy better than *expint-eps* is achived. 496;;; If the expansion don't converge within *expint-maxit* steps a Maxima 497;;; Error is thrown. 498;;; 499;;; The algorithm is based on technics desribed in Numerical Recipes, 2nd Ed. 500;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 501 502;;; Constants to terminate the numerical evaluation 503;;; 504;;; The accuracy *expint-eps* is fixed to 1.0e-15 for double float calculations. 505;;; The variable is declared global, so we can later give the Maxima User access 506;;; to the variable. The routine for Bigfloat numerical evaluation change this 507;;; value to the desired precision of the global $fpprec. 508;;; The maximum number of iterations is arbitrary set to 1000. For Bigfloat 509;;; evaluation this number is for very Big numbers too small. 510;;; 511;;; The maximum iterations counted for the test file rtest-expintegral.mac are 512;;; 101 for Complex Flonum and 1672 for Complex Bigfloat evaluation. 513 514(defvar *expint-eps* 1.0e-15) 515(defvar *expint-maxit* 1000) 516 517;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 518 519(defun expintegral-e (n z) 520 (declare (type integer n)) 521 (let ((*expint-eps* *expint-eps*) 522 (*expint-maxit* *expint-maxit*) 523 ;; Add (complex) 0 to get rid of any signed zeroes, and make z 524 ;; be a complex number. 525 (z (+ (coerce 0 '(complex flonum)) z))) 526 (declare (type (complex flonum) z)) 527 528 (when *debug-expintegral* 529 (format t "~&EXPINTEGRAL-E called with:~%") 530 (format t "~& : n = ~A~%" n) 531 (format t "~& : z = ~A~%" z)) 532 533 (cond 534 ((or (and (> (abs z) 2.0) (< (abs (phase z)) (* pi 0.9))) 535 ;; (abs z)>2.0 is necessary since there is a point 536 ;; -1.700598-0.612828*%i which 1<(abs z)<2, phase z < 0.9pi and 537 ;; still c-f expansion does not converge. 538 (and (>= (realpart z) 0) (> (abs z) 1.0))) 539 ;; We expand in continued fractions. 540 (when *debug-expintegral* 541 (format t "~&We expand in continued fractions.~%")) 542 (let* ((b (+ z n)) 543 (c (/ 1.0 (* *expint-eps* *expint-eps*))) 544 (d (/ 1.0 b)) 545 (n1 (- n 1)) 546 (h d) 547 (e 0.0)) 548 (do* ((i 1 (+ i 1)) 549 (a (* -1 n) (* (- i) (+ n1 i)))) 550 ((> i *expint-maxit*) 551 (merror 552 (intl:gettext "expintegral_e: continued fractions failed."))) 553 554 (setq b (+ b 2.0)) 555 (setq d (/ 1.0 (+ (* a d) b))) 556 (setq c (+ b (/ a c))) 557 (setq e (* c d)) 558 (setq h (* h e)) 559 560 (when (< (abs (- e 1.0)) *expint-eps*) 561 (when *debug-expintegral* 562 (setq *debug-expint-maxit* (max *debug-expint-maxit* i))) 563 (return (* h (exp (- z)))))))) 564 (t 565 ;; We expand in a power series. 566 (when *debug-expintegral* 567 (format t "~&We expand in a power series.~%")) 568 (let* ((n1 (- n 1)) 569 (euler (mget '$%gamma '$numer)) 570 (r (if (= n1 0) (- (- euler) (log z)) (/ 1.0 n1))) 571 (f 1.0) 572 (e 0.0)) 573 (do ((i 1 (+ i 1))) 574 ((> i *expint-maxit*) 575 (merror (intl:gettext "expintegral_e: series failed."))) 576 (setq f (* -1 f (/ z i))) 577 (cond 578 ((= i n1) 579 (let ((psi (- euler))) 580 (dotimes (ii n1) 581 (setq psi (+ psi (/ 1.0 (+ ii 1))))) 582 (setq e (* f (- psi (log z)))))) 583 (t 584 (setq e (/ (- f) (- i n1))))) 585 (setq r (+ r e)) 586 (when (< (abs e) (* (abs r) *expint-eps*)) 587 (when *debug-expintegral* 588 (setq *debug-expint-maxit* (max *debug-expint-maxit* i))) 589 (return r)))))))) 590 591;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 592;;; Numerical evaluation for a real or complex parameter. 593;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 594 595(defun frac-expintegral-e (n z) 596 (declare (type (complex flonum) n) 597 (type (complex flonum) z)) 598 599 (let ((*expint-eps* *expint-eps*) 600 (*expint-maxit* *expint-maxit*)) 601 602 (when *debug-expintegral* 603 (format t "~&FRAC-EXPINTEGRAL-E called with:~%") 604 (format t "~& : n = ~A~%" n) 605 (format t "~& : z = ~A~%" z)) 606 607 (cond 608 ((and (> (realpart z) 0) (> (abs z) 1.0)) 609 ;; We expand in continued fractions. 610 (when *debug-expintegral* 611 (format t "~&We expand in continued fractions.~%")) 612 (let* ((b (+ z n)) 613 (c (/ 1.0 (* *expint-eps* *expint-eps*))) 614 (d (/ 1.0 b)) 615 (n1 (- n 1)) 616 (h d) 617 (e 0.0)) 618 (do* ((i 1 (+ i 1)) 619 (a (* -1 n) (* (- i) (+ n1 i)))) 620 ((> i *expint-maxit*) 621 (merror 622 (intl:gettext "expintegral_e: continued fractions failed."))) 623 624 (setq b (+ b 2.0)) 625 (setq d (/ 1.0 (+ (* a d) b))) 626 (setq c (+ b (/ a c))) 627 (setq e (* c d)) 628 (setq h (* h e)) 629 630 (when (< (abs (- e 1.0)) *expint-eps*) 631 (when *debug-expintegral* 632 (setq *debug-expint-fracmaxit* (max *debug-expint-fracmaxit* i))) 633 (return (* h (exp (- z)))))))) 634 635 ((and (= (imagpart n) 0) 636 (> (realpart n) 0) 637 (= (nth-value 1 (truncate (realpart n))) 0)) 638 ;; We have a positive integer n or an float representation of an 639 ;; integer. We call expintegral-e which does this calculation. 640 (when *debug-expintegral* 641 (format t "~&We call expintegral-e.~%")) 642 (expintegral-e (truncate (realpart n)) z)) 643 644 (t 645 ;; At this point the parameter n is a real (not an float representation 646 ;; of an integer) or complex. We expand in a power series. 647 (when *debug-expintegral* 648 (format t "~&We expand in a power series.~%")) 649 (let* ((n1 (- n 1)) 650 ;; It would be possible to call the numerical implementation 651 ;; gamm-lanczos directly. But then the code would depend on the 652 ;; details of the implementation. 653 (gm (let ((tmp (take '(%gamma) (complexify (- 1 n))))) 654 (complex ($realpart tmp) ($imagpart tmp)))) 655 (r (- (* (expt z n1) gm) (/ 1.0 (- 1 n)))) 656 (f 1.0) 657 (e 0.0)) 658 (do ((i 1 (+ i 1))) 659 ((> i *expint-maxit*) 660 (merror (intl:gettext "expintegral_e: series failed."))) 661 (setq f (* -1 f (/ z (float i)))) 662 (setq e (/ (- f) (- (float i) n1))) 663 (setq r (+ r e)) 664 (when (< (abs e) (* (abs r) *expint-eps*)) 665 (when *debug-expintegral* 666 (setq *debug-expint-fracmaxit* (max *debug-expint-fracmaxit* i))) 667 (return r)))))))) 668 669;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 670;;; Helper functions for Bigfloat numerical evaluation. 671;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 672 673(defun cmul (x y) ($rectform (mul x y))) 674 675(defun cdiv (x y) ($rectform (div x y))) 676 677(defun cpower (x y) ($rectform (power x y))) 678 679;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 680;;; We have not changed the above algorithm, but generalized it to handle 681;;; complex and real Bigfloat numbers. By carefully examination of the 682;;; algorithm some of the additional calls to $rectform can be eliminated. 683;;; But the algorithm works and so we leave the extra calls for later work 684;;; in the code. 685;;; The accuracy of the result is determined by *expint-eps*. The value is 686;;; chosen to correspond to the value of $fpprec. We don't give any extra 687;;; digits to fpprec, so we loose 1 to 2 digits of precision. 688;;; One problem is to chose a sufficient big *expint-maxit*. 689;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 690 691(defun bfloat-expintegral-e (n z) 692 (let ((*expint-eps* (power ($bfloat 10.0) (- $fpprec))) 693 (*expint-maxit* 5000) ; arbitrarily chosen, we need a better choice 694 (bigfloattwo (add bigfloatone bigfloatone)) 695 (bigfloat%e ($bfloat '$%e)) 696 (bigfloat%gamma ($bfloat '$%gamma)) 697 (flz (complex ($float ($realpart z)) ($float ($imagpart z))))) 698 699 (when *debug-expintegral* 700 (format t "~&BFLOAT-EXPINTEGRAL-E called with:~%") 701 (format t "~& : n = ~A~%" n) 702 (format t "~& : z = ~A~%" flz)) 703 704 (cond 705 ((or (and (> (abs flz) 2) (< (abs (phase flz)) (* pi 0.9))) 706 ;; The same condition as you see in expintegral-e() 707 (and (>= (realpart flz) 0) (> (abs flz) 1.0))) 708 ;; We expand in continued fractions. 709 (when *debug-expintegral* 710 (format t "~&We expand in continued fractions.~%")) 711 (let* ((b (add z n)) 712 (c (div bigfloatone (mul *expint-eps* *expint-eps*))) 713 (d (cdiv bigfloatone b)) 714 (n1 (- n 1)) 715 (h d) 716 (e 0.0)) 717 (do* ((i 1 (+ i 1)) 718 (a (* -1 n) (* (- i) (+ n1 i)))) 719 ((> i *expint-maxit*) 720 (merror 721 (intl:gettext "expintegral_e: continued fractions failed."))) 722 723 (setq b (add b bigfloattwo)) 724 (setq d (cdiv bigfloatone (add (mul a d) b))) 725 (setq c (add b (cdiv a c))) 726 (setq e (cmul c d)) 727 (setq h (cmul h e)) 728 729 (when (eq ($sign (sub (cabs (sub e bigfloatone)) *expint-eps*)) 730 '$neg) 731 (when *debug-expintegral* 732 (setq *debug-expint-bfloatmaxit* 733 (max *debug-expint-bfloatmaxit* i))) 734 (return (cmul h (cpower bigfloat%e (mul -1 z)))))))) 735 (t 736 ;; We expand in a power series. 737 (when *debug-expintegral* 738 (format t "~&We expand in a power series.~%")) 739 (let* ((n1 (- n 1)) 740 (meuler (mul -1 bigfloat%gamma)) 741 (r (if (= n1 0) (sub meuler ($log z)) (div bigfloatone n1))) 742 (f bigfloatone) 743 (e bigfloatzero)) 744 (do* ((i 1 (+ i 1))) 745 ((> i *expint-maxit*) 746 (merror (intl:gettext "expintegral_e: series failed."))) 747 (setq f (mul -1 (cmul f (cdiv z i)))) 748 (cond 749 ((= i n1) 750 (let ((psi meuler)) 751 (dotimes (ii n1) 752 (setq psi (add psi (cdiv bigfloatone (+ ii 1))))) 753 (setq e (cmul f (sub psi ($log z)))))) 754 (t 755 (setq e (cdiv (mul -1 f) (- i n1))))) 756 (setq r (add r e)) 757 (when (eq ($sign (sub (cabs e) (cmul (cabs r) *expint-eps*))) 758 '$neg) 759 (when *debug-expintegral* 760 (setq *debug-expint-bfloatmaxit* 761 (max *debug-expint-bfloatmaxit* i))) 762 (return r)))))))) 763 764;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 765;;; Numerical Bigfloat evaluation for a real (Bigfloat) parameter. 766;;; The algorithm would work for a Complex Bigfloat parameter too. But we 767;;; need the values of Gamma for Complex Bigfloats. This is at this time (2008) 768;;; not implemented in Maxima. 769;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 770 771(defun frac-bfloat-expintegral-e (n z) 772 (let ((*expint-eps* (power ($bfloat 10.0) (- $fpprec))) 773 (*expint-maxit* 5000) ; arbitrarily chosen, we need a better choice 774 (bigfloattwo (add bigfloatone bigfloatone)) 775 (bigfloat%e ($bfloat '$%e)) 776 (bigfloat%gamma ($bfloat '$%gamma))) 777 778 (when *debug-expintegral* 779 (format t "~&FRAC-BFLOAT-EXPINTEGRAL-E called with:~%") 780 (format t "~& : n = ~A~%" n) 781 (format t "~& : z = ~A~%" z)) 782 783 (cond 784 ((and (or (eq ($sign ($realpart z)) '$pos) 785 (eq ($sign ($realpart z)) '$zero)) 786 (eq ($sign (sub (cabs z) bigfloatone)) '$pos)) 787 ;; We expand in continued fractions. 788 (when *debug-expintegral* 789 (format t "We expand in continued fractions.~%")) 790 (let* ((b (add z n)) 791 (c (div bigfloatone (mul *expint-eps* *expint-eps*))) 792 (d (cdiv bigfloatone b)) 793 (n1 (sub n 1)) 794 (h d) 795 (e 0.0)) 796 (do* ((i 1 (+ i 1)) 797 (a (mul -1 n) (cmul (- i) (add n1 i)))) 798 ((> i *expint-maxit*) 799 (merror 800 (intl:gettext "expintegral_e: continued fractions failed."))) 801 802 (setq b (add b bigfloattwo)) 803 (setq d (cdiv bigfloatone (add (mul a d) b))) 804 (setq c (add b (cdiv a c))) 805 (setq e (cmul c d)) 806 (setq h (cmul h e)) 807 808 (when (eq ($sign (sub (cabs (sub e bigfloatone)) *expint-eps*)) 809 '$neg) 810 (when *debug-expintegral* 811 (setq *debug-expint-fracbfloatmaxit* 812 (max *debug-expint-fracbfloatmaxit* i))) 813 (return (cmul h (cpower bigfloat%e (mul -1 z)))))))) 814 815 ((or (and (numberp n) 816 (= ($imagpart n) 0) 817 (> ($realpart n) 0) 818 (= (nth-value 1 (truncate ($realpart n))) 0)) 819 (and ($bfloatp n) 820 (eq ($sign n) '$pos) 821 (equal (sub (mul 2 ($fix n)) (mul 2 n)) 822 bigfloatzero))) 823 ;; We have a Float or Bigfloat representation of positive integer. 824 ;; We call bfloat-expintegral-e. 825 (when *debug-expintegral* 826 (format t "frac-Bigfloat with integer ~A~%" n)) 827 (bfloat-expintegral-e ($fix ($realpart n)) z)) 828 829 (t 830 ;; At this point the parameter n is a real (not a float representation 831 ;; of an integer) or complex. We expand in a power series. 832 (when *debug-expintegral* 833 (format t "We expand in a power series.~%")) 834 (let* ((n1 (sub n bigfloatone)) 835 (n2 (sub bigfloatone n)) 836 (gm (take '(%gamma) n2)) 837 (r (sub (cmul (cpower z n1) gm) (cdiv bigfloatone n2))) 838 (f bigfloatone) 839 (e bigfloatzero)) 840 (do ((i 1 (+ i 1))) 841 ((> i *expint-maxit*) 842 (merror (intl:gettext "expintegral_e: series failed."))) 843 (setq f (cmul (mul -1 bigfloatone) (cmul f (cdiv z i)))) 844 (setq e (cdiv (mul -1 f) (sub i n1))) 845 (setq r (add r e)) 846 (when (eq ($sign (sub (cabs e) (cmul (cabs r) *expint-eps*))) 847 '$neg) 848 (when *debug-expintegral* 849 (setq *debug-expint-fracbfloatmaxit* 850 (max *debug-expint-fracbfloatmaxit* i))) 851 (return r)))))))) 852 853;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 854;;; 855;;; Part 2: The implementation of the Exponential Integral E1 856;;; 857;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 858 859(defmfun $expintegral_e1 (z) 860 (simplify (list '(%expintegral_e1) z))) 861 862;;; Set properties to give full support to the parser and display 863 864(defprop $expintegral_e1 %expintegral_e1 alias) 865(defprop $expintegral_e1 %expintegral_e1 verb) 866 867(defprop %expintegral_e1 $expintegral_e1 reversealias) 868(defprop %expintegral_e1 $expintegral_e1 noun) 869 870;;; Exponential Integral E1 is a simplifying function 871 872(defprop %expintegral_e1 simp-expintegral_e1 operators) 873 874;;; Exponential Integral E1 distributes over bags 875 876(defprop %expintegral_e1 (mlist $matrix mequal) distribute_over) 877 878;;; Exponential Integral E1 has mirror symmetry, 879;;; but not on the real negative axis. 880 881(defprop %expintegral_e1 conjugate-expintegral-e1 conjugate-function) 882 883(defun conjugate-expintegral-e1 (args) 884 (let ((z (first args))) 885 (cond ((off-negative-real-axisp z) 886 ;; Definitely not on the negative real axis for z. Mirror symmetry. 887 (take '(%expintegral_e1) (take '($conjugate) z))) 888 (t 889 ;; On the negative real axis or no information. Unsimplified. 890 (list '($conjugate simp) (take '(%expintegral_e1) z)))))) 891 892;;; Differentiation of Exponential Integral E1 893 894(defprop %expintegral_e1 895 ((x) 896 ((mtimes) -1 897 ((mexpt) x -1) 898 ((mexpt) $%e ((mtimes) -1 x)))) 899 grad) 900 901;;; Integral of Exponential Integral E1 902 903(defprop %expintegral_e1 904 ((z) 905 ((mtimes) -1 ((%expintegral_e) 2 z))) 906 integral) 907 908;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 909 910;;; We support a simplim%function. The function is looked up in simplimit and 911;;; handles specific values of the function. 912 913(defprop %expintegral_e1 simplim%expintegral_e1 simplim%function) 914 915(defun simplim%expintegral_e1 (expr var val) 916 ;; Look for the limit of the argument. 917 (let ((z (limit (cadr expr) var val 'think))) 918 (cond 919 ;; Handle an argument 0 at this place 920 ((or (zerop1 z) 921 (eq z '$zeroa) 922 (eq z '$zerob)) 923 ;; limit is inf from both sides 924 '$inf) 925 (t 926 ;; All other cases are handled by the simplifier of the function. 927 (take '(%expintegral_e1) z))))) 928 929;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 930 931(defun simp-expintegral_e1 (expr ignored z) 932 (declare (ignore ignored)) 933 (oneargcheck expr) 934 (let ((arg (simpcheck (cadr expr) z))) 935 (cond 936 ;; Check for special values 937 ((or (eq arg '$inf) 938 (alike1 arg '((mtimes) -1 $minf))) 939 0) 940 ((zerop1 arg) 941 (simp-domain-error 942 (intl:gettext "expintegral_e1: expintegral_e1(~:M) is undefined.") arg)) 943 944 ;; Check for numerical evaluation 945 ((complex-float-numerical-eval-p arg) 946 ;; For E1 we call En(z) with n=1 directly. 947 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 948 (complexify (expintegral-e 1 carg)))) 949 950 ((complex-bigfloat-numerical-eval-p arg) 951 ;; For E1 we call En(z) with n=1 directly. 952 (let* (($ratprint nil) 953 (carg (add ($bfloat ($realpart arg)) 954 (mul '$%i ($bfloat ($imagpart arg))))) 955 (result (bfloat-expintegral-e 1 carg))) 956 (add ($realpart result) 957 (mul '$%i ($imagpart result))))) 958 959 ;; Check argument simplifications and transformations 960 ((taylorize (mop expr) (second expr))) 961 962 ((and $expintrep 963 (member $expintrep *expintflag* :test #'eq) 964 (not (eq $expintrep '%expintegral_e1))) 965 (case $expintrep 966 (%gamma_incomplete 967 (take '(%gamma_incomplete) 0 arg)) 968 (%expintegral_ei 969 (add (mul -1 (take '(%expintegral_ei) (mul -1 arg))) 970 (mul (inv 2) 971 (sub (take '(%log) (mul -1 arg)) 972 (take '(%log) (mul -1 (inv arg))))) 973 (mul -1 (take '(%log) arg)))) 974 (%expintegral_li 975 (add (mul -1 (take '(%expintegral_li) (power '$%e (mul -1 arg)))) 976 (mul -1 (take '(%log) arg)) 977 (mul (inv 2) 978 (sub (take '(%log) (mul -1 arg)) 979 (take '(%log) (mul -1 (inv arg))))))) 980 ($expintegral_trig 981 (add (mul -1 '$%i (take '(%expintegral_si) (mul '$%i arg))) 982 (mul -1 (take '(%expintegral_ci) (mul '$%i arg))) 983 (take '(%log) (mul '$%i arg)) 984 (mul -1 (take '(%log) arg)))) 985 ($expintegral_hyp 986 (sub (take '(%expintegral_shi) arg) 987 (take '(%expintegral_chi) arg))) 988 (t 989 (eqtest (list '(%expintegral_e1) arg) expr)))) 990 991 (t 992 (eqtest (list '(%expintegral_e1) arg) expr))))) 993 994;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 995;;; 996;;; Part 3: The implementation of the Exponential Integral Ei 997;;; 998;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 999 1000(defmfun $expintegral_ei (z) 1001 (simplify (list '(%expintegral_ei) z))) 1002 1003;;; Set properties to give full support to the parser and display 1004 1005(defprop $expintegral_ei %expintegral_ei alias) 1006(defprop $expintegral_ei %expintegral_ei verb) 1007 1008(defprop %expintegral_ei $expintegral_ei reversealias) 1009(defprop %expintegral_ei $expintegral_ei noun) 1010 1011;;; Exponential Integral Ei is a simplifying function 1012 1013(defprop %expintegral_ei simp-expintegral-ei operators) 1014 1015;;; Exponential Integral Ei distributes over bags 1016 1017(defprop %expintegral_ei (mlist $matrix mequal) distribute_over) 1018 1019;;; Exponential Integral Ei has mirror symmetry 1020 1021(defprop %expintegral_ei t commutes-with-conjugate) 1022 1023;;; Differentiation of Exponential Integral Ei 1024 1025(defprop %expintegral_ei 1026 ((x) 1027 ((mtimes) ((mexpt) x -1) ((mexpt) $%e x))) 1028 grad) 1029 1030;;; Integral of Exponential Ei 1031 1032(defprop %expintegral_ei 1033 ((x) 1034 ((mplus) 1035 ((mtimes) -1 ((mexpt) $%e x)) 1036 ((mtimes) x ((%expintegral_ei) x)))) 1037 integral) 1038 1039;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1040 1041;;; We support a simplim%function. The function is looked up in simplimit and 1042;;; handles specific values of the function. 1043 1044(defprop %expintegral_ei simplim%expintegral_ei simplim%function) 1045 1046(defun simplim%expintegral_ei (expr var val) 1047 ;; Look for the limit of the arguments. 1048 (let ((z (limit (cadr expr) var val 'think))) 1049 (cond 1050 ;; Handle an argument 0 at this place 1051 ((or (zerop1 z) 1052 (eq z '$zeroa) 1053 (eq z '$zerob)) 1054 '$minf) 1055 (t 1056 ;; All other cases are handled by the simplifier of the function. 1057 (take '(%expintegral_ei) z))))) 1058 1059;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1060 1061(defun simp-expintegral-ei (expr ignored z) 1062 (declare (ignore ignored)) 1063 (oneargcheck expr) 1064 (let ((arg (simpcheck (cadr expr) z))) 1065 (cond 1066 ;; Check special values 1067 ((zerop1 arg) 1068 (simp-domain-error 1069 (intl:gettext "expintegral_ei: expintegral_ei(~:M) is undefined.") 1070 arg)) 1071 ((or (eq arg '$inf) 1072 (alike1 arg '((mtimes) -1 $minf))) 1073 '$inf) 1074 ((or (eq arg '$minf) 1075 (alike1 arg '((mtimes) -1 $inf))) 1076 0) 1077 ((or (alike1 arg '((mtimes) $%i $inf)) 1078 (alike1 arg '((mtimes) -1 $%i $minf))) 1079 (mul '$%i '$%pi)) 1080 ((or (alike1 arg '((mtimes) $%i $minf)) 1081 (alike1 arg '((mtimes) -1 $%i $inf))) 1082 (mul -1 '$%i '$%pi)) 1083 1084 ;; Check numerical evaluation 1085 ((complex-float-numerical-eval-p arg) 1086 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 1087 (complexify (expintegral-ei carg)))) 1088 1089 ((complex-bigfloat-numerical-eval-p arg) 1090 (let* (($ratprint nil) 1091 (carg (add ($bfloat ($realpart arg)) 1092 (mul '$%i ($bfloat ($imagpart arg))))) 1093 (result (bfloat-expintegral-ei carg))) 1094 (add ($realpart result) 1095 (mul '$%i ($imagpart result))))) 1096 1097 ;; Check argument simplifications and transformations 1098 ((taylorize (mop expr) (second expr))) 1099 1100 ((and $expintrep 1101 (member $expintrep *expintflag*) 1102 (not (eq $expintrep '%expintegral_ei))) 1103 (case $expintrep 1104 (%gamma_incomplete 1105 (add (mul -1 (take '(%gamma_incomplete) 0 (mul -1 arg))) 1106 (mul (inv 2) 1107 (sub (take '(%log) arg) 1108 (take '(%log) (inv arg)))) 1109 (mul -1 (take '(%log) (mul -1 arg))))) 1110 (%expintegral_e1 1111 (add (mul -1 (take '(%expintegral_e1) (mul -1 arg))) 1112 (mul (inv 2) 1113 (sub (take '(%log) arg) 1114 (take '(%log) (inv arg)))) 1115 (mul -1 (take '(%log) (mul -1 arg))))) 1116 (%expintegral_li 1117 (take '(%expintegral_li) (power '$%e arg))) 1118 ($expintegral_trig 1119 (add (take '(%expintegral_ci) (mul '$%i arg)) 1120 (mul -1 '$%i (take '(%expintegral_si) (mul '$%i arg))) 1121 (mul (inv -2) 1122 (sub (take '(%log) (inv arg)) 1123 (take '(%log) arg))) 1124 (mul -1 (take '(%log) (mul '$%i arg))))) 1125 ($expintegral_hyp 1126 (add (take '(%expintegral_chi) arg) 1127 (take '(%expintegral_shi) arg) 1128 (mul (inv -2) 1129 (add (take '(%log) (inv arg)) 1130 (take '(%log) arg))))))) 1131 1132 (t 1133 (eqtest (list '(%expintegral_ei) arg) expr))))) 1134 1135;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1136;;; Numerical evaluation of the Exponential Integral Ei(z): 1137;;; 1138;;; We use the following representation (see functions.wolfram.com): 1139;;; 1140;;; Ei(z) = -E1(-z) + 0.5*(log(z)-log(1/z))-log(-z) 1141;;; 1142;;; z is a CL Complex number. Because we evaluate for Complex values we have to 1143;;; take into account the complete Complex phase factors. 1144;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1145 1146(defun expintegral-ei (z) 1147 (+ (- (expintegral-e 1 (- z))) 1148 ;; Carefully compute 1/2*(log(z)-log(1/z))-log(-z), using the 1149 ;; branch cuts that we want, not the one that Lisp wants. 1150 ;; (Mostly an issue with Lisps that support signed zeroes.) 1151 (cond 1152 ((> (imagpart z) 0) 1153 ;; Positive imaginary part. Add phase %i*%pi. 1154 (complex 0 (float pi))) 1155 ((< (imagpart z) 0) 1156 ;; Negative imaginary part. Add phase -%i*%pi. 1157 (complex 0 (- (float pi)))) 1158 ((> (realpart z) 0) 1159 ;; Positive real value. Add phase -%i*pi. 1160 (complex 0 (- (float pi)))) 1161 ;; Negative real value. No phase factor. 1162 (t 0)))) 1163 1164;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1165;;; We have not modified the algorithm for Bigfloat numbers. It is only 1166;;; generalized for Bigfloats. The calcualtion of the complex phase factor 1167;;; can be simplified to conditions about the sign of the realpart and 1168;;; imagpart. We leave this for further work to optimize the speed of the 1169;;; calculation. 1170;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1171 1172(defun bfloat-expintegral-ei (z) 1173 (let ((mz (mul -1 z))) 1174 (add (cmul (mul -1 bigfloatone) 1175 (bfloat-expintegral-e 1 mz)) 1176 (sub (cmul (div bigfloatone 2) 1177 (sub (take '(%log) z) 1178 (take '(%log) (cdiv bigfloatone z)))) 1179 (take '(%log) mz))))) 1180 1181;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1182;;; 1183;;; Part 4: The implementation of the Logarithmic integral li(z) 1184;;; 1185;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1186 1187(defmfun $expintegral_li (z) 1188 (simplify (list '(%expintegral_li) z))) 1189 1190;;; Set properties to give full support to the parser and display 1191 1192(defprop $expintegral_li %expintegral_li alias) 1193(defprop $expintegral_li %expintegral_li verb) 1194 1195(defprop %expintegral_li $expintegral_li reversealias) 1196(defprop %expintegral_li $expintegral_li noun) 1197 1198;;; Exponential Integral Li is a simplifying function 1199 1200(defprop %expintegral_li simp-expintegral-li operators) 1201 1202;;; Exponential Integral Li distributes over bags 1203 1204(defprop %expintegral_li (mlist $matrix mequal) distribute_over) 1205 1206;;; Exponential Integral Li has mirror symmetry, 1207;;; but not on the real negative axis. 1208 1209(defprop %expintegral_li conjugate-expintegral-li conjugate-function) 1210 1211(defun conjugate-expintegral-li (args) 1212 (let ((z (first args))) 1213 (cond ((off-negative-real-axisp z) 1214 ;; Definitely not on the negative real axis for z. Mirror symmetry. 1215 (take '(%expintegral_li) (take '($conjugate) z))) 1216 (t 1217 ;; On the negative real axis or no information. Unsimplified. 1218 (list '($conjugate simp) (take '(%expintegral_li) z)))))) 1219 1220;;; Differentiation of Exponential Integral Li 1221 1222(defprop %expintegral_li 1223 ((x) 1224 ((mtimes) ((mexpt) ((%log) x) -1))) 1225 grad) 1226 1227;;; Integral of Exponential Li 1228 1229(defprop %expintegral_li 1230 ((x) 1231 ((mplus) 1232 ((mtimes) x ((%expintegral_li) x)) 1233 ((mtimes) -1 ((%expintegral_ei) ((mtimes) 2 ((%log) x)))))) 1234 integral) 1235 1236;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1237 1238;;; We support a simplim%function. The function is looked up in simplimit and 1239;;; handles specific values of the function. 1240 1241(defprop %expintegral_li simplim%expintegral_li simplim%function) 1242 1243(defun simplim%expintegral_li (expr var val) 1244 ;; Look for the limit of the argument. 1245 (let ((z (limit (cadr expr) var val 'think))) 1246 (cond 1247 ;; Handle an argument 1 at this place 1248 ((onep1 z) '$minf) 1249 (t 1250 ;; All other cases are handled by the simplifier of the function. 1251 (take '(%expintegral_li) z))))) 1252 1253;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1254 1255(defun simp-expintegral-li (expr ignored z) 1256 (declare (ignore ignored)) 1257 (oneargcheck expr) 1258 (let ((arg (simpcheck (cadr expr) z))) 1259 (cond 1260 ((zerop1 arg) arg) 1261 ((onep1 arg) 1262 (simp-domain-error 1263 (intl:gettext "expintegral_li: expintegral_li(~:M) is undefined.") arg)) 1264 ((or (eq arg '$inf) 1265 (alike1 arg '((mtimes) -1 $minf))) 1266 '$inf) 1267 ((eq arg '$infinity) '$infinity) 1268 1269 ((complex-float-numerical-eval-p arg) 1270 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 1271 (complexify (expintegral-li carg)))) 1272 1273 ((complex-bigfloat-numerical-eval-p arg) 1274 (let* (($ratprint nil) 1275 (carg (add ($bfloat ($realpart arg)) 1276 (mul '$%i ($bfloat ($imagpart arg))))) 1277 (result (bfloat-expintegral-li carg))) 1278 (add (mul '$%i ($imagpart result)) 1279 ($realpart result)))) 1280 1281 ;; Check for argument simplifications and transformations 1282 ((taylorize (mop expr) (second expr))) 1283 1284 ((and $expintrep 1285 (member $expintrep *expintflag*) 1286 (not (eq $expintrep '%expintegral_li))) 1287 (let ((logarg (take '(%log) arg))) 1288 (case $expintrep 1289 (%gamma_incomplete 1290 (add (mul -1 (take '(%gamma_incomplete) 0 (mul -1 logarg))) 1291 (mul (inv 2) 1292 (sub (take '(%log) logarg) 1293 (take '(%log) (inv logarg)))) 1294 (mul -1 (take '(%log) (mul -1 logarg))))) 1295 (%expintegral_e1 1296 (add (mul -1 (take '(%expintegral_e1) (mul -1 logarg))) 1297 (mul (inv 2) 1298 (sub (take '(%log) logarg) 1299 (take '(%log) (inv logarg)))) 1300 (mul -1 (take '(%log) (mul -1 logarg))))) 1301 (%expintegral_ei 1302 ($expintegral_ei logarg)) 1303 ($expintegral_trig 1304 (add (take '(%expintegral_ci) (mul '$%i logarg)) 1305 (mul -1 '$%i (take '(%expintegral_si) (mul '$%i logarg))) 1306 (mul (inv -2) 1307 (sub (take '(%log) (inv logarg)) 1308 (take '(%log) logarg))) 1309 (mul -1 (take '(%log) (mul '$%i logarg))))) 1310 ($expintegral_hyp 1311 (add (take '(%expintegral_chi) logarg) 1312 (take '(%expintegral_shi) logarg) 1313 (mul (inv -2) 1314 (add (take '(%log) (inv logarg)) 1315 (take '(%log) logarg)))))))) 1316 1317 (t 1318 (eqtest (list '(%expintegral_li) arg) expr))))) 1319 1320;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1321;;; Numerical evaluation of the Expintegral Li 1322;;; 1323;;; We use the representation: 1324;;; 1325;;; Li(z) = Ei(log(z)) 1326;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1327 1328(defun expintegral-li (z) 1329 (expintegral-ei (log z))) 1330 1331(defun bfloat-expintegral-li (z) 1332 (bfloat-expintegral-ei ($log z))) 1333 1334;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1335;;; 1336;;; Part 5: The implementation of the Exponential Integral Si 1337;;; 1338;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1339 1340(defmfun $expintegral_si (z) 1341 (simplify (list '(%expintegral_si) z))) 1342 1343;;; Set properties to give full support to the parser and display 1344 1345(defprop $expintegral_si %expintegral_si alias) 1346(defprop $expintegral_si %expintegral_si verb) 1347 1348(defprop %expintegral_si $expintegral_si reversealias) 1349(defprop %expintegral_si $expintegral_si noun) 1350 1351;;; Exponential Integral Si is a simplifying function 1352 1353(defprop %expintegral_si simp-expintegral-si operators) 1354 1355;;; Exponential Integral Si distributes over bags 1356 1357(defprop %expintegral_si (mlist $matrix mequal) distribute_over) 1358 1359;;; Exponential Integral Si has mirror symmetry 1360 1361(defprop %expintegral_si t commutes-with-conjugate) 1362 1363;;; Exponential Integral Si is a odd function 1364 1365(defprop %expintegral_si odd-function-reflect reflection-rule) 1366 1367;;; Differentiation of Exponential Integral Si 1368 1369(defprop %expintegral_si 1370 ((x) 1371 ((mtimes) ((%sin) x) ((mexpt) x -1))) 1372 grad) 1373 1374;;; Integral of Exponential Si 1375 1376(defprop %expintegral_si 1377 ((x) 1378 ((mplus) 1379 ((%cos) x) 1380 ((mtimes) x ((%expintegral_si) x)))) 1381 integral) 1382 1383;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1384 1385;;; We support a simplim%function. 1386 1387(defprop %expintegral_si simplim%expintegral_si simplim%function) 1388 1389(defun simplim%expintegral_si (expr var val) 1390 ;; Look for the limit of the argument. 1391 (let ((z (limit (cadr expr) var val 'think))) 1392 ;; All cases are handled by the simplifier of the function. 1393 (take '(%expintegral_si) z))) 1394 1395;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1396 1397(defun simp-expintegral-si (expr ignored z) 1398 (declare (ignore ignored)) 1399 (oneargcheck expr) 1400 (let ((arg (simpcheck (cadr expr) z))) 1401 (cond 1402 ;; Check for special values 1403 ((zerop1 arg) arg) 1404 ((or (eq arg '$inf) 1405 (alike1 arg '((mtimes) -1 $minf))) 1406 (div '$%pi 2)) 1407 ((or (eq arg '$minf) 1408 (alike1 arg '((mtimes) -1 $inf))) 1409 (mul -1 (div '$%pi 2))) 1410 1411 ;; Check for numerical evaluation 1412 ((complex-float-numerical-eval-p arg) 1413 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 1414 (complexify (expintegral-si carg)))) 1415 1416 ((complex-bigfloat-numerical-eval-p arg) 1417 (let* (($ratprint nil) 1418 (carg (add ($bfloat ($realpart arg)) 1419 (mul '$%i ($bfloat ($imagpart arg))))) 1420 (result (bfloat-expintegral-si carg))) 1421 (add (mul '$%i ($imagpart result)) 1422 ($realpart result)))) 1423 1424 ;; Check for argument simplifications and transformations 1425 ((taylorize (mop expr) (second expr))) 1426 ((apply-reflection-simp (mop expr) arg $trigsign)) 1427 1428 ((and $expintrep 1429 (member $expintrep *expintflag*) 1430 (not (eq $expintrep '$expintegral_trig))) 1431 (case $expintrep 1432 (%gamma_incomplete 1433 (mul (div '$%i 2) 1434 (add (take '(%gamma_incomplete) 0 (mul -1 '$%i arg)) 1435 (mul -1 (take '(%gamma_incomplete) 0 (mul '$%i arg))) 1436 (take '(%log) (mul -1 '$%i arg)) 1437 (mul -1 (take '(%log) (mul '$%i arg)))))) 1438 (%expintegral_e1 1439 (mul (div '$%i 2) 1440 (add (take '(%expintegral_e1) (mul -1 '$%i arg)) 1441 (mul -1 (take '(%expintegral_e1) (mul '$%i arg))) 1442 (take '(%log) (mul -1 '$%i arg)) 1443 (mul -1 (take '(%log) (mul '$%i arg)))))) 1444 (%expintegral_ei 1445 (mul (div '$%i 4) 1446 (add (mul 2 1447 (sub (take '(%expintegral_ei) (mul -1 '$%i arg)) 1448 (take '(%expintegral_ei) (mul '$%i arg)))) 1449 (take '(%log) (div '$%i arg)) 1450 (mul -1 (take '(%log) (mul -1 (div '$%i arg)))) 1451 (mul -1 (take '(%log) (mul -1 '$%i arg))) 1452 (take '(%log) (mul '$%i arg))))) 1453 (%expintegral_li 1454 (mul (inv (mul 2 '$%i)) 1455 (add (take '(%expintegral_li) (power '$%e (mul '$%i arg))) 1456 (mul -1 1457 (take '(%expintegral_li) 1458 (power '$%e (mul -1 '$%e arg)))) 1459 (mul (div '$%pi -2) 1460 (take '(%signum) ($realpart arg)))))) 1461 ($expintegral_hyp 1462 (mul -1 '$%i (take '(%expintegral_shi) (mul '$%i arg)))))) 1463 1464 (t 1465 (eqtest (list '(%expintegral_si) arg) expr))))) 1466 1467;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1468;;; Numerical evaluation of the Exponential Integral Si 1469;;; 1470;;; We use the representation: 1471;;; 1472;;; Si(z) = %i/2 * (E1(-%i*z) - E1(*%i*z) + log(%i*z) - log(-%i*z)) 1473;;; 1474;;; For the Sin, Cos, Sinh and Cosh Exponential Integrals we have to call the 1475;;; numerical evaluation twice. In principle we could use a direct expansion 1476;;; in a power series or continued fractions to optimize the speed of the code. 1477;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1478 1479(defun expintegral-si (z) 1480 (let ((z (coerce z '(complex flonum)))) 1481 (* (complex 0 0.5) 1482 (+ (expintegral-e 1 (* (complex 0 -1) z)) 1483 (- (expintegral-e 1 (* (complex 0 1) z))) 1484 (log (* (complex 0 -1) z)) 1485 (- (log (* (complex 0 1) z))))))) 1486 1487(defun bfloat-expintegral-si (z) 1488 (let ((z*%i (cmul '$%i z)) 1489 (mz*%i (cmul (mul -1 '$%i) z))) 1490 (cmul 1491 (mul 0.5 '$%i) 1492 (add 1493 (bfloat-expintegral-e 1 mz*%i) 1494 (mul -1 (bfloat-expintegral-e 1 z*%i)) 1495 ($log mz*%i) 1496 (mul -1 ($log z*%i)))))) 1497 1498;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1499;;; 1500;;; Part 6: The implementation of the Exponential Integral Shi 1501;;; 1502;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1503 1504(defmfun $expintegral_shi (z) 1505 (simplify (list '(%expintegral_shi) z))) 1506 1507;;; Set properties to give full support to the parser and display 1508 1509(defprop $expintegral_shi %expintegral_shi alias) 1510(defprop $expintegral_shi %expintegral_shi verb) 1511 1512(defprop %expintegral_shi $expintegral_shi reversealias) 1513(defprop %expintegral_shi $expintegral_shi noun) 1514 1515;;; Exponential Integral Shi is a simplifying function 1516 1517(defprop %expintegral_shi simp-expintegral-shi operators) 1518 1519;;; Exponential Integral Shi distributes over bags 1520 1521(defprop %expintegral_shi (mlist $matrix mequal) distribute_over) 1522 1523;;; Exponential Integral Shi has mirror symmetry 1524 1525(defprop %expintegral_si t commutes-with-conjugate) 1526 1527;;; Exponential Integral Shi is a odd function 1528 1529(defprop %expintegral_si odd-function-reflect reflection-rule) 1530 1531;;; Differentiation of Exponential Integral Shi 1532 1533(defprop %expintegral_shi 1534 ((x) 1535 ((mtimes) ((%sinh) x) ((mexpt) x -1))) 1536 grad) 1537 1538;;; Integral of Exponential Shi 1539 1540(defprop %expintegral_shi 1541 ((x) 1542 ((mplus) 1543 ((mtimes) -1 ((%cosh) x)) 1544 ((mtimes) x ((%expintegral_shi) x)))) 1545 integral) 1546 1547;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1548 1549;;; We support a simplim%function. The function is looked up in simplimit and 1550;;; handles specific values of the function. 1551 1552(defprop %expintegral_shi simplim%expintegral_shi simplim%function) 1553 1554(defun simplim%expintegral_shi (expr var val) 1555 ;; Look for the limit of the argument. 1556 (let ((z (limit (cadr expr) var val 'think))) 1557 (cond 1558 ;; Handle infinities at this place 1559 ((eq z '$inf) 1560 '$inf) 1561 ((eq z '$minf) 1562 '$minf) 1563 (t 1564 ;; All other cases are handled by the simplifier of the function. 1565 (take '(%expintegral_shi) z))))) 1566 1567;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1568 1569(defun simp-expintegral-shi (expr ignored z) 1570 (declare (ignore ignored)) 1571 (oneargcheck expr) 1572 (let ((arg (simpcheck (cadr expr) z))) 1573 (cond 1574 ;; Check for special values 1575 ((zerop1 arg) arg) 1576 ((or (alike1 arg '((mtimes) $%i $inf)) 1577 (alike1 arg '((mtimes) -1 $%i $minf))) 1578 (div (mul '$%i '$%pi) 2)) 1579 ((or (alike1 arg '((mtimes) $%i $minf)) 1580 (alike1 arg '((mtimes) -1 $%i $inf))) 1581 (div (mul -1 '$%i '$%pi) 2)) 1582 1583 ;; Check for numrical evaluation 1584 ((float-numerical-eval-p arg) 1585 (realpart (expintegral-shi arg))) 1586 1587 ((complex-float-numerical-eval-p arg) 1588 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 1589 (complexify (expintegral-shi carg)))) 1590 1591 ((complex-bigfloat-numerical-eval-p arg) 1592 (let* (($ratprint nil) 1593 (carg (add ($bfloat ($realpart arg)) 1594 (mul '$%i ($bfloat ($imagpart arg))))) 1595 (result (bfloat-expintegral-shi carg))) 1596 (add (mul '$%i ($imagpart result)) 1597 ($realpart result)))) 1598 1599 ;; Check for argument simplifications and transformations 1600 ((taylorize (mop expr) (second expr))) 1601 ((apply-reflection-simp (mop expr) arg $trigsign)) 1602 1603 ((and $expintrep 1604 (member $expintrep *expintflag*) 1605 (not (eq $expintrep '$expintegral_hyp))) 1606 (case $expintrep 1607 (%gamma_incomplete 1608 (mul (inv 2) 1609 (add (take '(%gamma_incomplete) 0 arg) 1610 (mul -1 (take '(%gamma_incomplete) 0 (mul -1 arg))) 1611 (mul -1 (take '(%log) (mul -1 arg))) 1612 (take '(%log) arg)))) 1613 (%expintegral_e1 1614 (mul (inv 2) 1615 (add (take '(%expintegral_e1) arg) 1616 (mul -1 (take '(%expintegral_e1) (mul -1 arg))) 1617 (mul -1 (take '(%log) (mul -1 arg))) 1618 (take '(%log) arg)))) 1619 (%expintegral_ei 1620 (mul (inv 4) 1621 (add (mul 2 1622 (sub (take '(%expintegral_ei) arg) 1623 (take '(%expintegral_ei) (mul -1 arg)))) 1624 (take '(%log) (inv arg)) 1625 (mul -1 (take '(%log) (mul -1 (inv arg)))) 1626 (take '(%log) (mul -1 arg)) 1627 (mul -1 (take '(%log) arg))))) 1628 (%expintegral_li 1629 (add (mul (inv 2) 1630 (sub (take '(%expintegral_li) (power '$%e arg)) 1631 (take '(%expintegral_li) (power '$%e (mul -1 arg))))) 1632 (mul (div (mul '$%i '$%pi) -2) 1633 (take '(%signum) ($imagpart arg))))) 1634 ($expintegral_trig 1635 (mul -1 '$%i (take '(%expintegral_si) (mul '$%i arg)))))) 1636 1637 (t 1638 (eqtest (list '(%expintegral_shi) arg) expr))))) 1639 1640;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1641;;; Numerical evaluation of the Exponential Integral Shi 1642;;; 1643;;; We use the representation: 1644;;; 1645;;; Shi(z) = 1/2 * (E1(z) - E1(-z) - log(-z) + log(z)) 1646;;; 1647;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1648 1649(defun expintegral-shi (z) 1650 (* 1651 0.5 1652 (+ 1653 (expintegral-e 1 z) 1654 (- (expintegral-e 1 (- z))) 1655 (- (log (- z))) 1656 (log z)))) 1657 1658(defun bfloat-expintegral-shi (z) 1659 (let ((mz (mul -1 z))) 1660 (mul 1661 0.5 1662 (add 1663 (bfloat-expintegral-e 1 z) 1664 (mul -1 (bfloat-expintegral-e 1 mz)) 1665 (mul -1 ($log mz)) 1666 ($log z))))) 1667 1668;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1669;;; 1670;;; Part 7: The implementation of the Exponential Integral Ci 1671;;; 1672;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1673 1674(defmfun $expintegral_ci (z) 1675 (simplify (list '(%expintegral_ci) z))) 1676 1677;;; Set properties to give full support to the parser and display 1678 1679(defprop $expintegral_ci %expintegral_ci alias) 1680(defprop $expintegral_ci %expintegral_ci verb) 1681 1682(defprop %expintegral_ci $expintegral_ci reversealias) 1683(defprop %expintegral_ci $expintegral_ci noun) 1684 1685;;; Exponential Integral Ci is a simplifying function 1686 1687(defprop %expintegral_ci simp-expintegral-ci operators) 1688 1689;;; Exponential Integral Ci distributes over bags 1690 1691(defprop %expintegral_ci (mlist $matrix mequal) distribute_over) 1692 1693;;; Exponential Integral Ci has mirror symmetry, 1694;;; but not on the real negative axis. 1695 1696(defprop %expintegral_ci conjugate-expintegral-ci conjugate-function) 1697 1698(defun conjugate-expintegral-ci (args) 1699 (let ((z (first args))) 1700 (cond ((off-negative-real-axisp z) 1701 ;; Definitely not on the negative real axis for z. Mirror symmetry. 1702 (take '(%expintegral_ci) (take '($conjugate) z))) 1703 (t 1704 ;; On the negative real axis or no information. Unsimplified. 1705 (list '($conjugate simp) (take '(%expintegral_ci) z)))))) 1706 1707;;; Differentiation of Exponential Integral Ci 1708 1709(defprop %expintegral_ci 1710 ((x) 1711 ((mtimes) ((%cos) x) ((mexpt) x -1))) 1712 grad) 1713 1714;;; Integral of Exponential Ci 1715 1716(defprop %expintegral_ci 1717 ((x) 1718 ((mplus) 1719 ((mtimes) x ((%expintegral_ci) x)) 1720 ((mtimes) -1 ((%sin) x)))) 1721 integral) 1722 1723;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1724 1725;;; We support a simplim%function. The function is looked up in simplimit and 1726;;; handles specific values of the function. 1727 1728(defprop %expintegral_ci simplim%expintegral_ci simplim%function) 1729 1730(defun simplim%expintegral_ci (expr var val) 1731 ;; Look for the limit of the argument. 1732 (let ((z (limit (cadr expr) var val 'think))) 1733 (cond 1734 ;; Handle an argument 0 at this place 1735 ((or (zerop1 z) 1736 (eq z '$zeroa) 1737 (eq z '$zerob)) 1738 '$minf) 1739 (t 1740 ;; All other cases are handled by the simplifier of the function. 1741 (take '(%expintegral_ci) z))))) 1742 1743;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1744 1745(defun simp-expintegral-ci (expr ignored z) 1746 (declare (ignore ignored)) 1747 (oneargcheck expr) 1748 (let ((arg (simpcheck (cadr expr) z))) 1749 (cond 1750 ;; Check for special values 1751 ((zerop1 arg) 1752 (simp-domain-error 1753 (intl:gettext "expintegral_ci: expintegral_ci(~:M) is undefined.") arg)) 1754 ((or (eq arg '$inf) 1755 (alike1 arg '((mtimes) -1 $minf))) 1756 0) 1757 ((or (eq arg '$minf) 1758 (alike1 arg '((mtimes) -1 $inf))) 1759 (mul '$%i '$%pi)) 1760 1761 ;; Check for numerical evaluation 1762 ((complex-float-numerical-eval-p arg) 1763 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 1764 (complexify (expintegral-ci carg)))) 1765 1766 ((complex-bigfloat-numerical-eval-p arg) 1767 (let* (($ratprint nil) 1768 (carg (add ($bfloat ($realpart arg)) 1769 (mul '$%i ($bfloat ($imagpart arg))))) 1770 (result (bfloat-expintegral-ci carg))) 1771 (add (mul '$%i ($imagpart result)) 1772 ($realpart result)))) 1773 1774 ;; Check for argument simplifications and transformations 1775 ((taylorize (mop expr) (second expr))) 1776 1777 ((and $expintrep 1778 (member $expintrep *expintflag*) 1779 (not (eq $expintrep '$expintegral_trig))) 1780 (case $expintrep 1781 (%gamma_incomplete 1782 (sub (take '(%log) arg) 1783 (mul (inv 2) 1784 (add (take '(%gamma_incomplete) 0 (mul -1 '$%i arg)) 1785 (take '(%gamma_incomplete) 0 (mul '$%i arg)) 1786 (take '(%log) (mul -1 '$%i arg)) 1787 (take '(%log) (mul '$%i arg)))))) 1788 (%expintegral_e1 1789 (add (mul (inv -2) 1790 (add (take '(%expintegral_e1) (mul -1 '$%i arg)) 1791 (take '(%expintegral_e1) (mul '$%i arg))) 1792 (take '(%log) (mul -1 '$%i arg)) 1793 (take '(%log) (mul '$%i arg))) 1794 (take '(%log) arg))) 1795 (%expintegral_ei 1796 (add (mul (inv 4) 1797 (add (mul 2 1798 (add (take '(%expintegral_ei) (mul -1 '$%i arg)) 1799 (take '(%expintegral_ei) (mul '$%i arg)))) 1800 (take '(%log) (div '$%i arg)) 1801 (take '(%log) (mul -1 '$%i (inv arg))) 1802 (mul -1 (take '(%log) (mul -1 '$%i arg))) 1803 (mul -1 (take '(%log) (mul '$%i arg))))) 1804 (take '(%log) arg))) 1805 (%expintegral_li 1806 (add (mul (inv 2) 1807 (add (take '(%expintegral_li) 1808 (power '$%e (mul -1 '$%i arg))) 1809 (take '(%expintegral_li) 1810 (power '$%e (mul '$%i arg))))) 1811 (mul (div (mul '$%i '$%pi) 2) 1812 (take '(%signum) ($imagpart arg))) 1813 (sub 1 (take '(%signum) ($realpart arg))))) 1814 ($expintegral_hyp 1815 (add (take '(%expintegral_chi) (mul '$%i arg)) 1816 (mul -1 (take '(%log) (mul '$%i arg))) 1817 (take '(%log) arg))))) 1818 1819 (t 1820 (eqtest (list '(%expintegral_ci) arg) expr))))) 1821 1822;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1823;;; Numerical evaluation of the Exponential Integral Ci 1824;;; 1825;;; We use the representation: 1826;;; 1827;;; Ci(z) = -1/2 * (E1(-%i*z) + E1(%i*z) + log(-%i*z) + log(%i*z)) + log(z) 1828;;; 1829;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1830 1831(defun expintegral-ci (z) 1832 (let ((z (coerce z '(complex flonum)))) 1833 (+ (* -0.5 1834 (+ (expintegral-e 1 (* (complex 0 -1) z)) 1835 (expintegral-e 1 (* (complex 0 1) z)) 1836 (log (* (complex 0 -1) z)) 1837 (log (* (complex 0 1) z)))) 1838 (log z)))) 1839 1840(defun bfloat-expintegral-ci (z) 1841 (let ((z*%i (cmul '$%i z)) 1842 (mz*%i (cmul (mul -1 '$%i) z))) 1843 (add 1844 (cmul 1845 -0.5 1846 (add 1847 (bfloat-expintegral-e 1 mz*%i) 1848 (bfloat-expintegral-e 1 z*%i) 1849 ($log mz*%i) 1850 ($log z*%i))) 1851 ($log z)))) 1852 1853;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1854;;; 1855;;; Part 8: The implementation of the Exponential Integral Chi 1856;;; 1857;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1858 1859(defmfun $expintegral_chi (z) 1860 (simplify (list '(%expintegral_chi) z))) 1861 1862;;; Set properties to give full support to the parser and display 1863 1864(defprop $expintegral_chi %expintegral_chi alias) 1865(defprop $expintegral_chi %expintegral_chi verb) 1866 1867(defprop %expintegral_chi $expintegral_chi reversealias) 1868(defprop %expintegral_chi $expintegral_chi noun) 1869 1870;;; Exponential Integral Chi is a simplifying function 1871 1872(defprop %expintegral_chi simp-expintegral-chi operators) 1873 1874;;; Exponential Integral Chi distributes over bags 1875 1876(defprop %expintegral_chi (mlist $matrix mequal) distribute_over) 1877 1878;;; Exponential Integral Chi has mirror symmetry, 1879;;; but not on the real negative axis. 1880 1881(defprop %expintegral_chi conjugate-expintegral-chi conjugate-function) 1882 1883(defun conjugate-expintegral-chi (args) 1884 (let ((z (first args))) 1885 (cond ((off-negative-real-axisp z) 1886 ;; Definitely not on the negative real axis for z. Mirror symmetry. 1887 (take '(%expintegral_chi) (take '($conjugate) z))) 1888 (t 1889 ;; On the negative real axis or no information. Unsimplified. 1890 (list '($conjugate simp) (take '(%expintegral_chi) z)))))) 1891 1892;;; Differentiation of Exponential Integral Chi 1893 1894(defprop %expintegral_chi 1895 ((x) 1896 ((mtimes) ((%cosh) x) ((mexpt) x -1))) 1897 grad) 1898 1899;;; Integral of Exponential Chi 1900 1901(defprop %expintegral_chi 1902 ((x) 1903 ((mplus) 1904 ((mtimes) x ((%expintegral_chi) x)) 1905 ((mtimes) -1 ((%sinh) x)))) 1906 integral) 1907 1908;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1909 1910;;; We support a simplim%function. The function is looked up in simplimit and 1911;;; handles specific values of the function. 1912 1913(defprop %expintegral_chi simplim%expintegral_chi simplim%function) 1914 1915(defun simplim%expintegral_chi (expr var val) 1916 ;; Look for the limit of the argument. 1917 (let ((z (limit (cadr expr) var val 'think))) 1918 (cond 1919 ;; Handle an argument 0 at this place 1920 ((or (zerop1 z) 1921 (eq z '$zeroa) 1922 (eq z '$zerob)) 1923 '$minf) 1924 ((or (eq z '$inf) 1925 (eq z '$minf)) 1926 '$inf) 1927 (t 1928 ;; All other cases are handled by the simplifier of the function. 1929 (take '(%expintegral_chi) z))))) 1930 1931;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1932 1933(defun simp-expintegral-chi (expr ignored z) 1934 (declare (ignore ignored)) 1935 (oneargcheck expr) 1936 (let ((arg (simpcheck (cadr expr) z))) 1937 (cond 1938 ;; Check for special values 1939 ((zerop1 arg) 1940 ;; First check for zero argument. Throw Maxima error. 1941 (simp-domain-error 1942 (intl:gettext "expintegral_chi: expintegral_chi(~:M) is undefined.") 1943 arg)) 1944 ((or (alike1 arg '((mtimes) $%i $inf)) 1945 (alike1 arg '((mtimes) -1 $%i $minf))) 1946 (div (mul '$%pi '$%i) 2)) 1947 ((or (alike1 arg '((mtimes) $%i $minf)) 1948 (alike1 arg '((mtimes) -1 $%i $inf))) 1949 (div (mul -1 '$%pi '$%i) 2)) 1950 1951 ;; Check for numerical evaluation 1952 ((complex-float-numerical-eval-p arg) 1953 (let ((carg (complex ($float ($realpart arg)) ($float ($imagpart arg))))) 1954 (complexify (expintegral-chi carg)))) 1955 1956 ((complex-bigfloat-numerical-eval-p arg) 1957 (let* (($ratprint nil) 1958 (carg (add ($bfloat ($realpart arg)) 1959 (mul '$%i ($bfloat ($imagpart arg))))) 1960 (result (bfloat-expintegral-chi carg))) 1961 (add (mul '$%i ($imagpart result)) 1962 ($realpart result)))) 1963 1964 ;; Check for argument simplifications and transformations 1965 ((taylorize (mop expr) (second expr))) 1966 1967 ((and $expintrep 1968 (member $expintrep *expintflag*) 1969 (not (eq $expintrep '$expintegral_hyp))) 1970 (case $expintrep 1971 (%gamma_incomplete 1972 (mul (inv -2) 1973 (add (take '(%gamma_incomplete) 0 (mul -1 arg)) 1974 (take '(%gamma_incomplete) 0 arg) 1975 (take '(%log) (mul -1 arg)) 1976 (mul -1 (take '(%log) arg))))) 1977 (%expintegral_e1 1978 (mul (inv -2) 1979 (add (take '(%expintegral_e1) (mul -1 arg)) 1980 (take '(%expintegral_e1) arg) 1981 (take '(%log) (mul -1 arg)) 1982 (mul -1 (take '(%log) arg))))) 1983 (%expintegral_ei 1984 (mul (inv 4) 1985 (add (mul 2 1986 (add (take '(%expintegral_ei) (mul -1 arg)) 1987 (take '(%expintegral_ei) arg))) 1988 (take '(%log) (inv arg)) 1989 (take '(%log) (mul -1 (inv arg))) 1990 (mul -1 (take '(%log) (mul -1 arg))) 1991 (mul 3 (take '(%log) arg))))) 1992 (%expintegral_li 1993 (add (mul (inv 2) 1994 (add (take '(%expintegral_li) (power '$%e (mul -1 arg))) 1995 (take '(%expintegral_li) (power '$%e arg)))) 1996 (mul (div (mul '$%i '$%pi) 2) 1997 (take '(%signum) ($imagpart arg))) 1998 (mul (inv 2) 1999 (add (take '(%log) (inv arg)) 2000 (take '(%log) arg))))) 2001 ($expintegral_trig 2002 (add (take '(%expintegral_ci) (mul '$%i arg)) 2003 (take '(%log) arg) 2004 (mul -1 (take '(%log) (mul '$%i arg))))))) 2005 2006 (t 2007 (eqtest (list '(%expintegral_chi) arg) expr))))) 2008 2009;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2010;;; Numerical evaluation of the Exponential Integral Ci 2011;;; 2012;;; We use the representation: 2013;;; 2014;;; Chi(z) = -1/2 * (E1(-z) + E1(z) + log(-z) - log(z)) 2015;;; 2016;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2017 2018(defun expintegral-chi (z) 2019 (* 2020 -0.5 2021 (+ 2022 (expintegral-e 1 z) 2023 (expintegral-e 1 (- z)) 2024 (log (- z)) 2025 (- (log z))))) 2026 2027(defun bfloat-expintegral-chi (z) 2028 (let ((mz (mul -1 z))) 2029 (mul 2030 -0.5 2031 (add 2032 (bfloat-expintegral-e 1 z) 2033 (bfloat-expintegral-e 1 mz) 2034 ($log mz) 2035 (mul -1 ($log z)))))) 2036 2037;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2038;; Moved from bessel.lisp 2008-12-11. Consider deleting it. 2039;; 2040;; Exponential integral E1(x). The Cauchy principal value is used for 2041;; negative x. 2042(defmfun $expint (x) 2043 (cond ((numberp x) 2044 (values (slatec:de1 (float x)))) 2045 (t 2046 (list '($expint simp) x)))) 2047