src/bessel.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8
   9 (in-package :maxima)
  10
  11 (declare-top (special $hypergeometric_representation))
  12
  13 ;; When non-NIL, the Bessel functions of half-integral order are
  14 ;; expanded in terms of elementary functions.
  15
  16 (defmvar $besselexpand nil)
  17
  18 ;; When T Bessel functions with an integer order are reduced to order 0 and 1
  19 (defmvar $bessel_reduce nil)
  20
  21 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  22
  23 ;;; Helper functions for this file
  24
  25 ;; If E is a maxima ratio with a denominator of DEN, return the ratio
  26 ;; as a Lisp rational.  Otherwise NIL.
  27 (defun max-numeric-ratio-p (e den)
  28   (if (and (listp e)
  29            (eq 'rat (caar e))
  30            (= den (third e))
  31            (integerp (second e)))
  32       (/ (second e) (third e))
  33       nil))
  34
  35 (defun bessel-numerical-eval-p (order arg)
  36   ;; Return non-NIL if we should numerically evaluate a bessel
  37   ;; function.  Basically, both args have to be numbers.  If both args
  38   ;; are integers, we don't evaluate unless $numer is true.
  39   (or (and (numberp order) (complex-number-p arg)
  40            (or (floatp order)
  41                (floatp ($realpart arg))
  42                (floatp ($imagpart arg))))
  43       (and $numer (numberp order)
  44            (complex-number-p arg))))
  45
  46 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  47 ;;;
  48 ;;; Implementation of the Bessel J function
  49 ;;;
  50 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  51
  52 (defmfun $bessel_j (v z)
  53   (simplify (list '(%bessel_j) v z)))
  54
  55 ;; Bessel J distributes over lists, matrices, and equations
  56
  57 (defprop %bessel_j (mlist $matrix mequal) distribute_over)
  58
  59 ;; Derivatives of the Bessel function.
  60 (defprop %bessel_j
  61     ((n x)
  62      ;; Derivative wrt to order n.  A&S 9.1.64.  Do we really want to
  63      ;; do this?  It's quite messy.
  64      ;;
  65      ;; J[n](x)*log(x/2)
  66      ;;       - (x/2)^n*sum((-1)^k*psi(n+k+1)/gamma(n+k+1)*(x^2/4)^k/k!,k,0,inf)
  67      ((mplus)
  68       ((mtimes)
  69        ((%bessel_j) n x)
  70        ((%log) ((mtimes) ((rat) 1 2) x)))
  71       ((mtimes) -1
  72        ((mexpt) ((mtimes) x ((rat) 1 2)) n)
  73        ((%sum)
  74         ((mtimes) ((mexpt) -1 $%k)
  75          ((mexpt) ((mfactorial) $%k) -1)
  76          ((mqapply) (($psi array) 0) ((mplus) 1 $%k n))
  77          ((mexpt) ((%gamma) ((mplus) 1 $%k n)) -1)
  78          ((mexpt) ((mtimes) x x ((rat) 1 4)) $%k))
  79         $%k 0 $inf)))
  80
  81      ;; Derivative wrt to arg x.
  82      ;; A&S 9.1.27; changed from 9.1.30 so that taylor works on Bessel functions
  83      ((mtimes)
  84       ((mplus)
  85        ((%bessel_j) ((mplus) -1 n) x)
  86        ((mtimes) -1 ((%bessel_j) ((mplus) 1 n) x)))
  87       ((rat) 1 2)))
  88   grad)
  89
  90 ;; Integral of the Bessel function wrt z
  91 (defun bessel-j-integral-2 (v z)
  92   (case v
  93     (0
  94      ;; integrate(bessel_j(0,z),z)
  95      ;; = (1/2)*z*(%pi*bessel_j(1,z)*struve_h(0,z)
  96      ;;            +bessel_j(0,z)*(2-%pi*struve_h(1,z)))
  97      `((mtimes) ((rat) 1 2) ,z
  98        ((mplus)
  99         ((mtimes) $%pi
 100          ((%bessel_j) 1 ,z)
 101          ((%struve_h) 0 ,z))
 102         ((mtimes)
 103          ((%bessel_j) 0 ,z)
 104          ((mplus) 2 ((mtimes) -1 $%pi ((%struve_h) 1 ,z)))))))
 105     (1
 106      ;; integrate(bessel_j(1,z),z) = -bessel_j(0,z)
 107      `((mtimes) -1 ((%bessel_j) 0 ,z)))
 108     (otherwise
 109      ;; http://functions.wolfram.com/03.01.21.0002.01
 110      ;; integrate(bessel_j(v,z),z)
 111      ;;  = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2)
 112      ;;   * hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],-z^2/4)
 113      ;;  = 2^(-v)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],-z^2/4)
 114      ;;   / gamma(v+2)
 115      `((mtimes)
 116        (($hypergeometric)
 117         ((mlist)
 118          ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) ,v)))
 119         ((mlist)
 120          ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) ,v))
 121          ((mplus) 1 ,v))
 122         ((mtimes) ((rat) -1 4) ((mexpt) ,z 2)))
 123        ((mexpt) 2 ((mtimes) -1 ,v))
 124        ((mexpt) ((%gamma) ((mplus) 2 ,v)) -1)
 125        ((mexpt) z ((mplus) 1 ,v))))))
 126
 127 (putprop '%bessel_j `((v z) nil ,#'bessel-j-integral-2) 'integral)
 128
 129 ;; Support a simplim%bessel_j function to handle specific limits
 130
 131 (defprop %bessel_j simplim%bessel_j simplim%function)
 132
 133 (defun simplim%bessel_j (expr var val)
 134   ;; Look for the limit of the arguments.
 135   (let ((v (limit (cadr expr) var val 'think))
 136         (z (limit (caddr expr) var val 'think)))
 137   (cond
 138     ;; Handle an argument 0 at this place.
 139     ((or (zerop1 z)
 140          (eq z '$zeroa)
 141          (eq z '$zerob))
 142      (let ((sgn ($sign ($realpart v))))
 143        (cond ((and (eq sgn '$neg)
 144                    (not (maxima-integerp v)))
 145               ;; bessel_j(v,0), Re(v)<0 and v not an integer
 146               '$infinity)
 147              ((and (eq sgn '$zero)
 148                    (not (zerop1 v)))
 149               ;; bessel_j(v,0), Re(v)=0 and v #0
 150               '$und)
 151              ;; Call the simplifier of the function.
 152              (t (simplify (list '(%bessel_j) v z))))))
 153     ((or (eq z '$inf)
 154          (eq z '$minf))
 155      ;; bessel_j(v,inf) or bessel_j(v,minf)
 156      0)
 157     (t
 158      ;; All other cases are handled by the simplifier of the function.
 159      (simplify (list '(%bessel_j) v z))))))
 160
 161 (def-simplifier bessel_j (order arg)
 162   (let ((rat-order nil))
 163     (cond
 164       ((zerop1 arg)
 165        ;; We handle the different case for zero arg carefully.
 166        (let ((sgn ($sign ($realpart order))))
 167          (cond ((and (eq sgn '$zero)
 168                      (zerop1 ($imagpart order)))
 169                 ;; bessel_j(0,0) = 1
 170                 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1))
 171                       ((or (floatp order) (floatp arg)) 1.0)
 172                       (t 1)))
 173                ((or (eq sgn '$pos)
 174                     (maxima-integerp order))
 175                 ;; bessel_j(v,0) and Re(v)>0 or v an integer
 176                 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0))
 177                       ((or (floatp order) (floatp arg)) 0.0)
 178                       (t 0)))
 179                ((and (eq sgn '$neg)
 180                      (not (maxima-integerp order)))
 181                 ;; bessel_j(v,0) and Re(v)<0 and v not an integer
 182                 (simp-domain-error
 183                   (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.")
 184                   order arg))
 185                ((and (eq sgn '$zero)
 186                      (not (zerop1 ($imagpart order))))
 187                 ;; bessel_j(v,0) and Re(v)=0 and v # 0
 188                 (simp-domain-error
 189                   (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.")
 190                   order arg))
 191                (t
 192                 ;; No information about the sign of the order
 193                 (give-up)))))
 194
 195       ((complex-float-numerical-eval-p order arg)
 196        ;; We have numeric order and arg and $numer is true, or we have either
 197        ;; the order or arg being floating-point, so let's evaluate it
 198        ;; numerically.
 199        ;; The numerical routine bessel-j returns a CL number, so we have
 200        ;; to add the conversion to a Maxima-complex-number.
 201        (cond ((= 0 ($imagpart order))
 202               ;; order is real, arg is real or complex
 203               (complexify (bessel-j ($float order)
 204                                     (complex ($float ($realpart arg))
 205                                              ($float ($imagpart arg))))))
 206              (t
 207               ;; order is complex, arg is real or complex
 208               (let (($numer t)
 209                     ($float t)
 210                     (order ($float order))
 211                     (arg ($float arg)))
 212                 ($float
 213                   ($rectform
 214                     (bessel-j-hypergeometric order arg)))))))
 215
 216       ((or (bigfloat-numerical-eval-p order arg)
 217            (complex-bigfloat-numerical-eval-p order arg))
 218        ;; We have the order or arg being a bigfloat, so evaluate it
 219        ;; numerically using the hypergeometric representation
 220        ;; bessel_j.
 221        (destructuring-bind (re . im)
 222            (risplit order)
 223          (cond ((and (zerop1 im)
 224                      (zerop1 (sub re (take '($floor) re)))
 225                      (mminusp re))
 226                 ;; bessel_j(-n,z) = (-1)^n*bessel_j(n,z)
 227                 ($rectform
 228                  (mul (power -1 re)
 229                       ($bfloat (bessel-j-hypergeometric (mul -1 re) arg)))))
 230                (t
 231                 ($rectform
 232                  ($bfloat (bessel-j-hypergeometric order arg)))))))
 233
 234       ((and (integerp order) (minusp order))
 235        ;; Some special cases when the order is an integer.
 236        ;; A&S 9.1.5: J[-n](x) = (-1)^n*J[n](x)
 237        (if (evenp order)
 238            (take '(%bessel_j) (- order) arg)
 239            (mul -1 (take '(%bessel_j) (- order) arg))))
 240
 241       ((and $besselexpand
 242             (setq rat-order (max-numeric-ratio-p order 2)))
 243        ;; When order is a fraction with a denominator of 2, we
 244        ;; can express the result in terms of elementary functions.
 245        (bessel-j-half-order rat-order arg))
 246
 247       ((and $bessel_reduce
 248             (and (integerp order)
 249                  (plusp order)
 250                  (> order 1)))
 251        ;; Reduce a bessel function of order > 2 to order 1 and 0.
 252        ;; A&S 9.1.27: bessel_j(v,z) = 2*(v-1)/z*bessel_j(v-1,z)-bessel_j(v-2,z)
 253        (sub (mul 2
 254                  (- order 1)
 255                  (inv arg)
 256                  (take '(%bessel_j) (- order 1) arg))
 257             (take '(%bessel_j) (- order 2) arg)))
 258
 259       ((and $%iargs (multiplep arg '$%i))
 260        ;; bessel_j(v, %i*x) = (%i*x)^v/(x^v) * bessel_i(v, x)
 261        ;; (From http://functions.wolfram.com/03.01.27.0002.01)
 262        (let ((x (coeff arg '$%i 1)))
 263          (mul (power (mul '$%i x) order)
 264               (inv (power x order))
 265               (take '(%bessel_i) order x))))
 266
 267       ($hypergeometric_representation
 268        (bessel-j-hypergeometric order arg))
 269
 270       (t
 271        (give-up)))))
 272
 273 ;; Returns the hypergeometric representation of bessel_j
 274 (defun bessel-j-hypergeometric (order arg)
 275   ;; A&S 9.1.69 / http://functions.wolfram.com/03.01.26.0002.01
 276   ;; hypergeometric([],[v+1],-z^2/4)*(z/2)^v/gamma(v+1)
 277   (mul (inv (take '(%gamma) (add order 1)))
 278        (power (div arg 2) order)
 279        (take '($hypergeometric)
 280              (list '(mlist))
 281              (list '(mlist) (add order 1))
 282              (neg (div (mul arg arg) 4)))))
 283
 284 ;; Compute value of Bessel function of the first kind of order ORDER.
 285 (defun bessel-j (order arg)
 286   (cond
 287     ((zerop (imagpart arg))
 288      ;; We have numeric args and the arg is purely real.
 289      ;; Call the real-valued Bessel function when possible.
 290      (let ((arg (realpart arg)))
 291        (cond
 292          ((= order 0)
 293           (slatec:dbesj0 (float arg)))
 294          ((= order 1)
 295           (slatec:dbesj1 (float arg)))
 296          ((minusp order)
 297           (cond ((zerop (nth-value 1 (truncate order)))
 298                  ;; The order is a negative integer.  We use A&S 9.1.5:
 299                  ;;   J[-n](z) = (-1)^n*J[n](z)
 300                  ;; and not the Hankel functions.
 301                  (if (evenp (floor order))
 302                      (bessel-j (- order) arg)
 303                      (- (bessel-j (- order) arg))))
 304                 (t
 305                  ;; Bessel function of negative order.  We use the Hankel
 306                  ;; functions to compute this:
 307                  ;;   J[v](z) = (H1[v](z) + H2[v](z)) / 2.
 308                  ;; This works for negative and positive arg and handles
 309                  ;; special cases correctly.
 310                  (let ((result (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg)))))
 311                    (cond ((= (nth-value 1 (floor order)) 1/2)
 312                           ;; ORDER is a half-integral-value or a float
 313                           ;; representation, thereof.
 314                           (if (minusp arg)
 315                               ;; arg is negative, the result is purely imaginary
 316                               (complex 0 (imagpart result))
 317                               ;; arg is positive, the result is purely real
 318                               (realpart result)))
 319                          ;; in all other cases general complex result
 320                          (t result))))))
 321          (t
 322           ;; We have a real arg and order > 0 and order not 0 or 1
 323           ;; for this case we can call the function dbesj
 324           (multiple-value-bind (n alpha) (floor (float order))
 325             (let ((jvals (make-array (1+ n) :element-type 'flonum)))
 326               (slatec:dbesj (abs (float arg)) alpha (1+ n) jvals 0)
 327               (cond ((>= arg 0)
 328                      (aref jvals n))
 329                     (t
 330                      ;; Use analytic continuation formula A&S 9.1.35:
 331                      ;; J[v](z*exp(m*%pi*%i)) = exp(m*%pi*%i*v)*J[v](z)
 332                      ;; for an integer m.  In particular, for m = 1:
 333                      ;; J[v](-x) = exp(v*%pi*%i)*J[v](x)
 334                      ;; and handle special cases
 335                      (cond ((zerop (nth-value 1 (truncate order)))
 336                             ;; order is an integer
 337                             (if (evenp (floor order))
 338                                 (aref jvals n)
 339                                 (- (aref jvals n))))
 340                            ((= (nth-value 1 (floor order)) 1/2)
 341                             ;; Order is a half-integral-value and we know that
 342                             ;; arg < 0, so the result is purely imginary.
 343                             (if (evenp (floor order))
 344                                 (complex 0 (aref jvals n))
 345                                 (complex 0 (- (aref jvals n)))))
 346                            ;; In all other cases a general complex result
 347                            (t
 348                             (* (cis (* order pi))
 349                                (aref jvals n))))))))))))
 350     (t
 351      ;; The arg is complex. Use the complex-valued Bessel function.
 352      (cond ((mminusp order)
 353             ;; Bessel function of negative order. We use the Hankel function to
 354             ;; compute this, because A&S 9.1.3 and 9.1.4 say
 355             ;;   H1[v](z) = J[v](z) + %i * Y[v](z)
 356             ;; and
 357             ;;   H2[v](z) = J[v](z) - %i * Y[v](z).
 358             ;; Thus
 359             ;;   J[v](z) = (H1[v](z) + H2[v](z)) / 2.
 360             ;; Not the most efficient way, but perhaps good enough for maxima.
 361             (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg))))
 362            (t
 363             (multiple-value-bind (n alpha) (floor (float order))
 364               (let ((cyr (make-array (1+ n) :element-type 'flonum))
 365                     (cyi (make-array (1+ n) :element-type 'flonum)))
 366                 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
 367                                            v-cyr v-cyi v-nz v-ierr)
 368                    (slatec:zbesj (float (realpart arg))
 369                                  (float (imagpart arg))
 370                                  alpha 1 (1+ n) cyr cyi 0 0)
 371                   (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz))
 372
 373                   ;; Should check the return status in v-ierr of this routine.
 374                   (when (plusp v-ierr)
 375                     (format t "zbesj ierr = ~A~%" v-ierr))
 376                   (complex (aref cyr n) (aref cyi n))))))))))
 377
 378 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 379 ;;;
 380 ;;; Implementation of the Bessel Y function
 381 ;;;
 382 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 383
 384 (defmfun $bessel_y (v z)
 385   (simplify (list '(%bessel_y) v z)))
 386
 387 ;; Bessel Y distributes over lists, matrices, and equations
 388
 389 (defprop %bessel_y (mlist $matrix mequal) distribute_over)
 390
 391 (defprop %bessel_y
 392     ((n x)
 393      ;; Derivative wrt order n.  A&S 9.1.65.
 394      ;;
 395      ;; cot(n*%pi)*[diff(bessel_j(n,x),n)-%pi*bessel_y(n,x)]
 396      ;;  - csc(n*%pi)*diff(bessel_j(-n,x),n)-%pi*bessel_j(n,x)
 397      ((mplus)
 398       ((mtimes) -1 $%pi ((%bessel_j) n x))
 399       ((mtimes)
 400        -1
 401        ((%csc) ((mtimes) $%pi n))
 402        ((%derivative) ((%bessel_j) ((mtimes) -1 n) x) n 1))
 403       ((mtimes)
 404        ((%cot) ((mtimes) $%pi n))
 405        ((mplus)
 406         ((mtimes) -1 $%pi ((%bessel_y) n x))
 407         ((%derivative) ((%bessel_j) n x) n 1))))
 408
 409      ;; Derivative wrt to arg x.  A&S 9.1.27; changed from A&S 9.1.30
 410      ;; to be consistent with bessel_j.
 411      ((mtimes)
 412       ((mplus)
 413        ((%bessel_y)((mplus) -1 n) x)
 414        ((mtimes) -1 ((%bessel_y) ((mplus) 1 n) x)))
 415       ((rat) 1 2)))
 416     grad)
 417
 418 ;; Integral of the Bessel Y function wrt z
 419 ;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/21/01/01/
 420 (defun bessel-y-integral-2 (n z)
 421   (cond
 422    ((and ($integerp n) (<= 0 n))
 423     (cond
 424      (($oddp n)
 425       ;; integrate(bessel_y(2*N+1,z),z) , N > 0
 426       ;; = -bessel_y(0,z) - 2 * sum(bessel_y(2*k,z),k,1,N)
 427       (let* ((k (gensym))
 428              (answer `((mplus) ((mtimes) -1 ((%bessel_y) 0 ,z))
 429                        ((mtimes) -2
 430                         ((%sum) ((%bessel_y) ((mtimes) 2 ,k) ,z) ,k 1
 431                          ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))))
 432         ;; Expand out the sum if n < 10.  Otherwise fix up the indices
 433         (if (< n 10)
 434             (meval `(($ev) ,answer $sum))   ; Is there a better way?
 435           (simplify ($niceindices answer)))))
 436      (($evenp n)
 437       ;; integrate(bessel_y(2*N,z),z) , N > 0
 438       ;; = (1/2)*%pi*z*(bessel_y(0,z)*struve_h(-1,z)
 439       ;;               +bessel_y(1,z)*struve_h(0,z))
 440       ;;    - 2 * sum(bessel_y(2*k+1,z),k,0,N-1)
 441       (let*
 442           ((k (gensym))
 443            (answer `((mplus)
 444                      ((mtimes) -2
 445                       ((%sum)
 446                        ((%bessel_y) ((mplus) 1 ((mtimes) 2 ,k)) ,z)
 447                        ,k 0
 448                        ((mplus)
 449                         -1
 450                         ((mtimes) ((rat) 1 2) ,n))))
 451                      ((mtimes) ((rat) 1 2) $%pi ,z
 452                       ((mplus)
 453                        ((mtimes)
 454                         ((%bessel_y) 0 ,z)
 455                         ((%struve_h) -1 ,z))
 456                        ((mtimes)
 457                         ((%bessel_y) 1 ,z)
 458                         ((%struve_h) 0 ,z)))))))
 459         ;; Expand out the sum if n < 10.  Otherwise fix up the indices
 460         (if (< n 10)
 461             (meval `(($ev) ,answer $sum))  ; Is there a better way?
 462           (simplify ($niceindices answer)))))))
 463    (t nil)))
 464
 465 (putprop '%bessel_y `((n z) nil ,#'bessel-y-integral-2) 'integral)
 466
 467 ;; Support a simplim%bessel_y function to handle specific limits
 468
 469 (defprop %bessel_y simplim%bessel_y simplim%function)
 470
 471 (defun simplim%bessel_y (expr var val)
 472   ;; Look for the limit of the arguments.
 473   (let ((v (limit (cadr expr) var val 'think))
 474         (z (limit (caddr expr) var val 'think)))
 475   (cond
 476     ;; Handle an argument 0 at this place.
 477     ((or (zerop1 z)
 478          (eq z '$zeroa)
 479          (eq z '$zerob))
 480      (cond ((zerop1 v)
 481             ;; bessel_y(0,0)
 482             '$minf)
 483            ((integerp v)
 484             ;; bessel_y(n,0), n an integer
 485             (cond ((evenp v) '$minf)
 486                   (t (cond ((eq z '$zeroa) '$minf)
 487                            ((eq z '$zerob) '$inf)
 488                            (t '$infinity)))))
 489            ((not (zerop1 ($realpart v)))
 490             ;; bessel_y(v,0), Re(v)#0
 491             '$infinity)
 492            ((and (zerop1 ($realpart v))
 493                  (not (zerop1 v)))
 494             ;; bessel_y(v,0), Re(v)=0 and v#0
 495             '$und)
 496            ;; Call the simplifier of the function.
 497            (t (simplify (list '(%bessel_y) v z)))))
 498     ((or (eq z '$inf)
 499          (eq z '$minf))
 500      ;; bessel_y(v,inf) or bessel_y(v,minf)
 501      0)
 502     (t
 503      ;; All other cases are handled by the simplifier of the function.
 504      (simplify (list '(%bessel_y) v z))))))
 505
 506 (def-simplifier bessel_y (order arg)
 507   (let ((rat-order nil))
 508     (cond
 509       ((zerop1 arg)
 510        ;; Domain error for a zero argument.
 511        (simp-domain-error
 512          (intl:gettext "bessel_y: bessel_y(~:M,~:M) is undefined.") order arg))
 513
 514       ((complex-float-numerical-eval-p order arg)
 515        ;; We have numeric order and arg and $numer is true, or
 516        ;; we have either the order or arg being floating-point,
 517        ;; so let's evaluate it numerically.
 518        (cond ((= 0 ($imagpart order))
 519               ;; order is real, arg is real or complex
 520               (complexify (bessel-y ($float order)
 521                                     (complex ($float ($realpart arg))
 522                                              ($float ($imagpart arg))))))
 523              (t
 524               ;; order is complex, arg is real or complex
 525               (let (($numer t)
 526                     ($float t)
 527                     (order ($float order))
 528                     (arg ($float arg)))
 529                 ($float
 530                   ($rectform
 531                     (bessel-y-hypergeometric order arg)))))))
 532
 533       ((or (bigfloat-numerical-eval-p order arg)
 534            (complex-bigfloat-numerical-eval-p order arg))
 535        ;; When the order is not an integer, we can use the
 536        ;; hypergeometric representation to evaluate it.
 537        (destructuring-bind (re . im)
 538            (risplit order)
 539          (cond ((and (zerop1 im)
 540                      (not (zerop1 (sub re (take '($floor) re)))))
 541                 ($rectform
 542                  ($bfloat (bessel-y-hypergeometric re arg))))
 543                (t
 544                 (give-up)))))
 545
 546       ((and (integerp order) (minusp order))
 547        ;; Special case when the order is an integer.
 548        ;; A&S 9.1.5: Y[-n](x) = (-1)^n*Y[n](x)
 549        (if (evenp order)
 550            (take '(%bessel_y) (- order) arg)
 551            (mul -1 (take '(%bessel_y) (- order) arg))))
 552
 553       ((and $besselexpand
 554             (setq rat-order (max-numeric-ratio-p order 2)))
 555        ;; When order is a fraction with a denominator of 2, we
 556        ;; can express the result in terms of elementary functions.
 557        ;; From A&S 10.1.1, 10.1.11-12 and 10.1.15:
 558        ;;
 559        ;; Y[1/2](z) = -J[-1/2](z) is a function of cos(z).
 560        ;; Y[-1/2](z) = J[1/2](z) is a function of sin(z).
 561        (bessel-y-half-order rat-order arg))
 562
 563       ((and $bessel_reduce
 564             (and (integerp order)
 565                  (plusp order)
 566                  (> order 1)))
 567        ;; Reduce a bessel function of order > 2 to order 1 and 0.
 568        ;; A&S 9.1.27: bessel_y(v,z) = 2*(v-1)/z*bessel_y(v-1,z)-bessel_y(v-2,z)
 569        (sub (mul 2
 570                  (- order 1)
 571                  (inv arg)
 572                  (take '(%bessel_y) (- order 1) arg))
 573             (take '(%bessel_y) (- order 2) arg)))
 574
 575       ($hypergeometric_representation
 576        (bessel-y-hypergeometric order arg))
 577
 578       (t
 579        (give-up)))))
 580
 581 ;; Returns the hypergeometric representation of bessel_y
 582 (defun bessel-y-hypergeometric (order arg)
 583   ;; http://functions.wolfram.com/03.03.26.0002.01
 584   ;; -hypergeometric([],[1-v],-z^2/4)*(2/z)^v*gamma(v)/%pi
 585   ;;   - hypergeometric([],[v+1],-z^2/4)*gamma(-v)*cos(%pi*v)*(z/2)^v/%pi
 586   (add (mul -1
 587             (power 2 order)
 588             (inv (power arg order))
 589             (take '(%gamma) order)
 590             (inv '$%pi)
 591             (take '($hypergeometric)
 592                   (list '(mlist))
 593                   (list '(mlist) (sub 1 order))
 594                   (neg (div (mul arg arg) 4))))
 595        (mul -1
 596             (inv (power 2 order))
 597             (power arg order)
 598             (take '(%cos) (mul order '$%pi))
 599             (take '(%gamma) (neg order))
 600             (inv '$%pi)
 601             (take '($hypergeometric)
 602                   (list '(mlist))
 603                   (list '(mlist) (add 1 order))
 604                   (neg (div (mul arg arg) 4))))))
 605
 606 ;; Bessel function of the second kind, Y[n](z), for real or complex z
 607 (defun bessel-y (order arg)
 608   (cond
 609     ((zerop (imagpart arg))
 610      ;; We have numeric args and the first arg is purely
 611      ;; real. Call the real-valued Bessel function.
 612      ;;
 613      ;; For negative values, use the analytic continuation formula
 614      ;; A&S 9.1.36:
 615      ;;
 616      ;; Y[v](z*exp(m*%pi*%i)) = exp(-v*m*%pi*%i) * Y[v](z)
 617      ;;       + 2*%i*sin(m*v*%pi)*cot(v*%pi)*J[v](z)
 618      ;;
 619      ;; In particular for m = 1:
 620      ;; Y[v](-z) = exp(-v*%pi*%i) * Y[v](z) + 2*%i*cos(v*%pi)*J[v](z)
 621      (let ((arg (realpart arg)))
 622        (cond
 623          ((zerop order)
 624           (cond ((>= arg 0)
 625                  (slatec:dbesy0 (float arg)))
 626                 (t
 627                  ;; For v = 0, this simplifies to
 628                  ;; Y[0](-z) = Y[0](z) + 2*%i*J[0](z)
 629                  ;; the return value has to be a CL number
 630                  (+ (slatec:dbesy0 (float (- arg)))
 631                     (complex 0 (* 2 (slatec:dbesj0 (float (- arg)))))))))
 632          ((= order 1)
 633           (cond ((>= arg 0)
 634                  (slatec:dbesy1 (float arg)))
 635                 (t
 636                  ;; For v = 1, this simplifies to
 637                  ;; Y[1](-z) = -Y[1](z) - 2*%i*J[1](v)
 638                  ;; the return value has to be a CL number
 639                  (+ (- (slatec:dbesy1 (float (- arg))))
 640                     (complex 0 (* -2 (slatec:dbesj1 (float (- arg)))))))))
 641          ((minusp order)
 642           (cond ((zerop (nth-value 1 (truncate order)))
 643                  ;; Order is a negative integer or float representation.
 644                  ;; Use A&S 9.1.5: Y[-n](z) = (-1)^n*Y[n](z).
 645                  (if (evenp (floor order))
 646                      (bessel-y (- order) arg)
 647                      (- (bessel-y (- order) arg))))
 648                 (t
 649                  ;; Bessel function of negative order. We use the Hankel
 650                  ;; functions to compute this:
 651                  ;;   Y[v](z) = (H1[v](z) - H2[v](z)) / 2 / %i
 652                  (let ((result (/ (- (hankel-1 order arg)
 653                                      (hankel-2 order arg))
 654                                   (complex 0 2))))
 655                    (cond ((= (nth-value 1 (floor order)) 1/2)
 656                           ;; ORDER is half-integral-value or a float
 657                           ;; representation thereof.
 658                           (if (minusp arg)
 659                               ;; arg is negative, the result is purely imaginary
 660                               (complex 0 (imagpart result))
 661                               ;; arg is positive, the result is purely real
 662                               (realpart result)))
 663                          ;; in all other cases general complex result
 664                          (t result))))))
 665          (t
 666           (multiple-value-bind (n alpha) (floor (float order))
 667             (let ((jvals (make-array (1+ n) :element-type 'flonum)))
 668               ;; First we do the calculation for an positive argument.
 669               (slatec:dbesy (abs (float arg)) alpha (1+ n) jvals)
 670
 671               ;; Now we look at the sign of the argument
 672               (cond ((>= arg 0)
 673                      (aref jvals n))
 674                     (t
 675                      (let* ((dpi (coerce pi 'flonum))
 676                             (s1 (cis (- (* order dpi))))
 677                             (s2 (* #c(0 2) (cos (* order dpi)))))
 678                        (let ((result (+ (* s1 (aref jvals n))
 679                                         (* s2 (bessel-j order (- arg))))))
 680                          (cond ((zerop (nth-value 1 (truncate order)))
 681                                 ;; ORDER is an integer or a float representation
 682                                 ;; of an integer, and the arg is positive the
 683                                 ;; result is general complex.
 684                                 result)
 685                                ((= (nth-value 1 (floor order)) 1/2)
 686                                 ;; ORDER is a half-integral-value or an float
 687                                 ;; representation and we have arg < 0. The
 688                                 ;; result is purely imaginary.
 689                                 (complex 0 (imagpart result)))
 690                                ;; in all other cases general complex result
 691                                (t result))))))))))))
 692     (t
 693      (cond ((minusp order)
 694             ;; Bessel function of negative order. We use the Hankel function to
 695             ;; compute this, because A&S 9.1.3 and 9.1.4 say
 696             ;;   H1[v](z) = J[v](z) + %i * Y[v](z)
 697             ;; and
 698             ;;   H2[v](z) = J[v](z) - %i * Y[v](z).
 699             ;; Thus
 700             ;;   Y[v](z) = (H1[v](z) - H2[v](z)) / 2 / %i
 701             (/ (- (hankel-1 order arg) (hankel-2 order arg))
 702                (complex 0 2)))
 703            (t
 704             (multiple-value-bind (n alpha) (floor (float order))
 705               (let ((cyr (make-array (1+ n) :element-type 'flonum))
 706                     (cyi (make-array (1+ n) :element-type 'flonum))
 707                     (cwrkr (make-array (1+ n) :element-type 'flonum))
 708                     (cwrki (make-array (1+ n) :element-type 'flonum)))
 709                 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi
 710                                            v-nz v-cwrkr v-cwrki v-ierr)
 711                    (slatec::zbesy (float (realpart arg))
 712                                   (float (imagpart arg))
 713                                   alpha 1 (1+ n) cyr cyi 0 cwrkr cwrki 0)
 714                   (declare (ignore v-zr v-zi v-fnu v-kode v-n
 715                                    v-cyr v-cyi v-cwrkr v-cwrki v-nz))
 716
 717                   ;; We should check for errors based on the value of v-ierr.
 718                   (when (plusp v-ierr)
 719                     (format t "zbesy ierr = ~A~%" v-ierr))
 720
 721                   (complex (aref cyr n) (aref cyi n))))))))))
 722
 723 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 724 ;;;
 725 ;;; Implementation of the Bessel I function
 726 ;;;
 727 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 728
 729 (defmfun $bessel_i (v z)
 730   (simplify (list '(%bessel_i) v z)))
 731
 732 ;; Bessel I distributes over lists, matrices, and equations
 733
 734 (defprop %bessel_i (mlist $matrix mequal) distribute_over)
 735
 736 (defprop %bessel_i
 737     ((n x)
 738      ;; Derivative wrt order n.  A&S 9.6.42.
 739      ;;
 740      ;; I[n](x)*log(x/2)
 741      ;;   - (x/2)^n*sum(psi(n+k+1)/gamma(n+k+1)*(x^2/4)^k/k!,k,0,inf)
 742      ((mplus)
 743       ((mtimes)
 744        ((%bessel_i) n x)
 745        ((%log) ((mtimes) ((rat) 1 2) x)))
 746       ((mtimes) -1
 747        ((mexpt) ((mtimes) x ((rat) 1 2)) n)
 748        ((%sum)
 749         ((mtimes)
 750          ((mexpt) ((mfactorial) $%k) -1)
 751          ((mqapply) (($psi array) 0) ((mplus) 1 $%k n))
 752          ((mexpt) ((%gamma) ((mplus) 1 $%k n)) -1)
 753          ((mexpt) ((mtimes) x x ((rat) 1 4)) $%k))
 754         $%k 0 $inf)))
 755      ;; Derivative wrt to x.  A&S 9.6.29.
 756      ((mtimes)
 757       ((mplus) ((%bessel_i) ((mplus) -1 n) x)
 758                ((%bessel_i) ((mplus) 1 n) x))
 759       ((rat) 1 2)))
 760   grad)
 761
 762 ;; Integral of the Bessel I function wrt z
 763 ;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselI/21/01/01/
 764 (defun bessel-i-integral-2 (v z)
 765   (case v
 766         (0
 767          ;; integrate(bessel_i(0,z),z)
 768          ;; = (1/2)*z*(bessel_i(0,z)*(%pi*struve_l(1,z)+2)
 769          ;;            -%pi*bessel_i(1,z)*struve_l(0,z))
 770          `((mtimes) ((rat) 1 2) ,z
 771            ((mplus)
 772             ((mtimes) -1 $%pi
 773              ((%bessel_i) 1 ,z)
 774              ((%struve_l) 0 ,z))
 775             ((mtimes)
 776              ((%bessel_i) 0 ,z)
 777              ((mplus) 2
 778               ((mtimes) $%pi ((%struve_l) 1 ,z)))))))
 779         (1
 780          ;; integrate(bessel_i(1,z),z) = bessel_i(0,z)
 781          `((%bessel_i) 0 ,z))
 782         (otherwise
 783          ;; http://functions.wolfram.com/03.02.21.0002.01
 784          ;; integrate(bessel_i(v,z),z)
 785          ;;  = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2)
 786          ;;   * hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],z^2/4)
 787          ;;  = 2^(-v)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],z^2/4)
 788          ;;   / gamma(v+2)
 789          `((mtimes)
 790            (($hypergeometric)
 791             ((mlist)
 792              ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) ,v)))
 793             ((mlist)
 794              ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) ,v))
 795              ((mplus) 1 ,v))
 796             ((mtimes) ((rat) 1 4) ((mexpt) ,z 2)))
 797            ((mexpt) 2 ((mtimes) -1 ,v))
 798            ((mexpt) ((%gamma) ((mplus) 2 ,v)) -1)
 799            ((mexpt) z ((mplus) 1 ,v))))))
 800
 801 (putprop '%bessel_i `((v z) nil ,#'bessel-i-integral-2) 'integral)
 802
 803 ;; Support a simplim%bessel_i function to handle specific limits
 804
 805 (defprop %bessel_i simplim%bessel_i simplim%function)
 806
 807 (defun simplim%bessel_i (expr var val)
 808   ;; Look for the limit of the arguments.
 809   (let ((v (limit (cadr expr) var val 'think))
 810         (z (limit (caddr expr) var val 'think)))
 811   (cond
 812     ;; Handle an argument 0 at this place.
 813     ((or (zerop1 z)
 814          (eq z '$zeroa)
 815          (eq z '$zerob))
 816      (let ((sgn ($sign ($realpart v))))
 817        (cond ((and (eq sgn '$neg)
 818                    (not (maxima-integerp v)))
 819               ;; bessel_i(v,0), Re(v)<0 and v not an integer
 820               '$infinity)
 821              ((and (eq sgn '$zero)
 822                    (not (zerop1 v)))
 823               ;; bessel_i(v,0), Re(v)=0 and v #0
 824               '$und)
 825              ;; Call the simplifier of the function.
 826              (t (simplify (list '(%bessel_i) v z))))))
 827     ((eq z '$inf)
 828      ;; bessel_i(v,inf)
 829      '$inf)
 830     ((eq z '$minf)
 831      ;; bessel_i(v,minf)
 832      '$infinity)
 833     (t
 834      ;; All other cases are handled by the simplifier of the function.
 835      (simplify (list '(%bessel_i) v z))))))
 836
 837 (def-simplifier bessel_i (order arg)
 838   (let ((rat-order nil))
 839     (cond
 840       ((zerop1 arg)
 841        ;; We handle the different case for zero arg carefully.
 842        (let ((sgn ($sign ($realpart order))))
 843          (cond ((and (eq sgn '$zero)
 844                      (zerop1 ($imagpart order)))
 845                 ;; bessel_i(0,0) = 1
 846                 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1))
 847                       ((or (floatp order) (floatp arg)) 1.0)
 848                       (t 1)))
 849                ((or (eq sgn '$pos)
 850                     (maxima-integerp order))
 851                 ;; bessel_i(v,0) and Re(v)>0 or v an integer
 852                 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0))
 853                       ((or (floatp order) (floatp arg)) 0.0)
 854                       (t 0)))
 855                ((and (eq sgn '$neg)
 856                      (not (maxima-integerp order)))
 857                 ;; bessel_i(v,0) and Re(v)<0 and v not an integer
 858                 (simp-domain-error
 859                   (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.")
 860                   order arg))
 861                ((and (eq sgn '$zero)
 862                      (not (zerop1 ($imagpart order))))
 863                 ;; bessel_i(v,0) and Re(v)=0 and v # 0
 864                 (simp-domain-error
 865                   (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.")
 866                   order arg))
 867                (t
 868                 ;; No information about the sign of the order
 869                 (give-up)))))
 870
 871       ((complex-float-numerical-eval-p order arg)
 872        (cond ((= 0 ($imagpart order))
 873               ;; order is real, arg is real or complex
 874               (complexify (bessel-i ($float order)
 875                                     (complex ($float ($realpart arg))
 876                                              ($float ($imagpart arg))))))
 877              (t
 878               ;; order is complex, arg is real or complex
 879               (let (($numer t)
 880                     ($float t)
 881                     (order ($float order))
 882                     (arg ($float arg)))
 883                 ($float
 884                   ($rectform
 885                     (bessel-i-hypergeometric order arg)))))))
 886
 887       ((or (bigfloat-numerical-eval-p order arg)
 888            (complex-bigfloat-numerical-eval-p order arg))
 889        ;; We have the order or arg being a bigfloat, so evaluate it
 890        ;; numerically using the hypergeometric representation
 891        ;; bessel_j.
 892        ;;
 893        ;; When the order is a negative integer, the hypergeometric
 894        ;; representation is invalid.  However,
 895        ;; https://dlmf.nist.gov/10.27.E1 says bessel_i(-n,z) =
 896        ;; bessel_i(n,z).
 897        (destructuring-bind (re . im)
 898            (risplit order)
 899          (cond ((and (zerop im)
 900                      (zerop1 (sub re (take '($floor) re)))
 901                      (mminusp re))
 902                 ($rectform
 903                  ($bfloat (bessel-i-hypergeometric (mul -1 order) arg))))
 904                (t
 905                 ($rectform
 906                  ($bfloat (bessel-i-hypergeometric order arg)))))))
 907
 908       ((and (integerp order) (minusp order))
 909        ;; Some special cases when the order is an integer
 910        ;; A&S 9.6.6: I[-n](x) = I[n](x)
 911        (take '(%bessel_i) (- order) arg))
 912
 913       ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
 914        ;; When order is a fraction with a denominator of 2, we
 915        ;; can express the result in terms of elementary functions.
 916        ;; From A&S 10.2.13 and 10.2.14:
 917        ;;
 918        ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
 919        ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
 920        (bessel-i-half-order rat-order arg))
 921
 922       ((and $bessel_reduce
 923             (and (integerp order)
 924                  (plusp order)
 925                  (> order 1)))
 926        ;; Reduce a bessel function of order > 2 to order 1 and 0.
 927        ;; A&S 9.6.26: bessel_i(v,z) = -2*(v-1)/z*bessel_i(v-1,z)+bessel_i(v-2,z)
 928        (add (mul -2
 929                  (- order 1)
 930                  (inv arg)
 931                  (take '(%bessel_i) (- order 1) arg))
 932             (take '(%bessel_i) (- order 2) arg)))
 933
 934       ((and $%iargs (multiplep arg '$%i))
 935        ;; bessel_i(v, %i*x) = (%i*x)^v/(x^v) * bessel_j(v, x)
 936        ;; (From http://functions.wolfram.com/03.02.27.0002.01)
 937        (let ((x (coeff arg '$%i 1)))
 938          (mul (power (mul '$%i x) order)
 939               (inv (power x order))
 940               (take '(%bessel_j) order x))))
 941
 942       ($hypergeometric_representation
 943        (bessel-i-hypergeometric order arg))
 944
 945       (t
 946        (give-up)))))
 947
 948 ;; Returns the hypergeometric representation of bessel_i
 949 (defun bessel-i-hypergeometric (order arg)
 950   ;; A&S 9.6.47 / http://functions.wolfram.com/03.02.26.0002.01
 951   ;; hypergeometric([],[v+1],z^2/4)*(z/2)^v/gamma(v+1)
 952   (mul (inv (take '(%gamma) (add order 1)))
 953        (power (div arg 2) order)
 954        (take '($hypergeometric)
 955              (list '(mlist))
 956              (list '(mlist) (add order 1))
 957              (div (mul arg arg) 4))))
 958
 959 ;; Compute value of Modified Bessel function of the first kind of order ORDER
 960 (defun bessel-i (order arg)
 961   (cond
 962     ((zerop (imagpart arg))
 963      ;; We have numeric args and the first arg is purely
 964      ;; real. Call the real-valued Bessel function.  Use special
 965      ;; routines for order 0 and 1, when possible
 966      (let ((arg (realpart arg)))
 967        (cond
 968          ((zerop order)
 969           (slatec:dbesi0 (float arg)))
 970          ((= order 1)
 971           (slatec:dbesi1 (float arg)))
 972          ((or (minusp order) (< arg 0))
 973           (multiple-value-bind (order-int order-frac) (floor order)
 974             (cond ((zerop order-frac)
 975                    ;; Order is an integer. From A&S 9.6.6 and 9.6.30 we have
 976                    ;;   I[-n](z) = I[n](z)
 977                    ;; and
 978                    ;;   I[n](-z) = (-1)^n*I[n](z)
 979                    (if (< arg 0)
 980                        (if (evenp order-int)
 981                            (bessel-i (abs order) (abs arg))
 982                            (- (bessel-i (abs order) (abs arg))))
 983                        (bessel-i (abs order) arg)))
 984                   (t
 985                    ;; Order or arg is negative and order is not an integer, use
 986                    ;; the bessel-j function for calculation.  We know from
 987                    ;; the definition I[v](x) = x^v/(%i*x)^v * J[v](%i*x).
 988                    (let* ((arg (float arg))
 989                           (result (* (expt arg order)
 990                                      (expt (complex 0 arg) (- order))
 991                                      (bessel-j order (complex 0 arg)))))
 992                      ;; Try to clean up result if we know the result is
 993                      ;; purely real or purely imaginary.
 994                      (cond ((>= arg 0)
 995                             ;; Result is purely real for arg >= 0
 996                             (realpart result))
 997                            ((zerop order-frac)
 998                             ;; Order is an integer or a float representation of
 999                             ;; an integer, the result is purely real.
1000                             (realpart result))
1001                            ((= order-frac 1/2)
1002                             ;; Order is half-integral-value or a float
1003                             ;; representation and arg < 0, the result
1004                             ;; is purely imaginary.
1005                             (complex 0 (imagpart result)))
1006                            (t result)))))))
1007          (t
1008           ;; Now the case order > 0 and arg >= 0
1009           (multiple-value-bind (n alpha) (floor (float order))
1010             (let ((jvals (make-array (1+ n) :element-type 'flonum)))
1011               (slatec:dbesi (float (realpart arg)) alpha 1 (1+ n) jvals 0)
1012               (aref jvals n)))))))
1013
1014     ((and (zerop (realpart arg))
1015           (zerop (rem order 1)))
1016      ;; Handle the case for a pure imaginary arg and integer order.
1017      ;; In this case, the result is purely real or purely imaginary,
1018      ;; so we want to make sure that happens.
1019      ;;
1020      ;; bessel_i(n, %i*x) = (%i*x)^n/x^n * bessel_j(n,x)
1021      ;;    = %i^n * bessel_j(n,x)
1022      (let* ((n (floor order))
1023             (result (bessel-j (float n) (imagpart arg))))
1024        (cond ((evenp n)
1025               ;; %i^(2*m) = (-1)^m, where n=2*m
1026               (if (evenp (/ n 2))
1027                   result
1028                   (- result)))
1029              ((oddp n)
1030               ;; %i^(2*m+1) = %i*(-1)^m, where n = 2*m+1
1031               (if (evenp (floor n 2))
1032                   (complex 0 result)
1033                   (complex 0 (- result)))))))
1034
1035     (t
1036      ;; The arg is complex.  Use the complex-valued Bessel function.
1037      (multiple-value-bind (n alpha) (floor (abs (float order)))
1038        ;; Evaluate the function for positive order and fixup the result later.
1039        (let ((cyr (make-array (1+ n) :element-type 'flonum))
1040              (cyi (make-array (1+ n) :element-type 'flonum)))
1041          (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
1042                                     v-cyr v-cyi v-nz v-ierr)
1043             (slatec::zbesi (float (realpart arg))
1044                            (float (imagpart arg))
1045                            alpha 1 (1+ n) cyr cyi 0 0)
1046            (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz))
1047            ;; We should check for errors here based on the value of v-ierr.
1048            (when (plusp v-ierr)
1049              (format t "zbesi ierr = ~A~%" v-ierr))
1050
1051            ;; We have evaluated I[|order|](arg), now we look at
1052            ;; the the sign of the order.
1053            (cond ((minusp order)
1054                   ;; From A&S 9.6.2:
1055                   ;; I[-a](z) = I[a](z) + (2/%pi)*sin(%pi*a)*K[a](z)
1056                   (+ (complex (aref cyr n) (aref cyi n))
1057                      (let ((dpi (coerce pi 'flonum)))
1058                        (* (/ 2.0 dpi)
1059                           (sin (* dpi (- order)))
1060                           (bessel-k (- order) arg)))))
1061                  (t
1062                   (complex (aref cyr n) (aref cyi n))))))))))
1063
1064 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1065 ;;;
1066 ;;; Implementation of the Bessel K function
1067 ;;;
1068 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1069
1070 (defmfun $bessel_k (v z)
1071   (simplify (list '(%bessel_k) v z)))
1072
1073 ;; Bessel K distributes over lists, matrices, and equations
1074
1075 (defprop %bessel_k (mlist $matrix mequal) distribute_over)
1076
1077 (defprop %bessel_k
1078     ((n x)
1079      ;; Derivativer wrt order n.  A&S 9.6.43.
1080      ;;
1081      ;; %pi/2*csc(n*%pi)*['diff(bessel_i(-n,x),n)-'diff(bessel_i(n,x),n)]
1082      ;;    - %pi*cot(n*%pi)*bessel_k(n,x)
1083      ((mplus)
1084       ((mtimes) -1 $%pi
1085        ((%bessel_k) n x)
1086        ((%cot) ((mtimes) $%pi n)))
1087       ((mtimes)
1088        ((rat) 1 2)
1089        $%pi
1090        ((%csc) ((mtimes) $%pi n))
1091        ((mplus)
1092         ((%derivative) ((%bessel_i) ((mtimes) -1 n) x) n 1)
1093         ((mtimes) -1
1094          ((%derivative) ((%bessel_i) n x) n 1)))))
1095      ;; Derivative wrt to x.  A&S 9.6.29.
1096      ((mtimes)
1097       -1
1098       ((mplus) ((%bessel_k) ((mplus) -1 n) x)
1099                ((%bessel_k) ((mplus) 1 n) x))
1100       ((rat) 1 2)))
1101   grad)
1102
1103 ;; Integral of the Bessel K function wrt z
1104 ;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/21/01/01/
1105 (defun bessel-k-integral-2 (n z)
1106   (cond
1107    ((and ($integerp n) (<= 0 n))
1108     (cond
1109      (($oddp n)
1110       ;; integrate(bessel_k(2*N+1,z),z) , N > 0
1111       ;; = -(-1)^((n-1)/2)*bessel_k(0,z)
1112       ;;   + 2*sum((-1)^(k+(n-1)/2-1)*bessel_k(2*k,z),k,1,(n-1)/2)
1113       (let* ((k (gensym))
1114              (answer `((mplus)
1115                        ((mtimes) -1 ((%bessel_k) 0 ,z)
1116                         ((mexpt) -1
1117                          ((mtimes) ((rat) 1 2) ((mplus) -1 ,n))))
1118                        ((mtimes) 2
1119                         ((%sum)
1120                          ((mtimes) ((%bessel_k) ((mtimes) 2 ,k) ,z)
1121                           ((mexpt) -1
1122                            ((mplus) -1 ,k
1123                             ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))
1124                          ,k 1 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))))
1125         ;; Expand out the sum if n < 10.  Otherwise fix up the indices
1126         (if (< n 10)
1127             (meval `(($ev) ,answer $sum))   ; Is there a better way?
1128           (simplify ($niceindices answer)))))
1129      (($evenp n)
1130       ;; integrate(bessel_k(2*N,z),z) , N > 0
1131       ;; = (1/2)*(-1)^(n/2)*%pi*z*(bessel_k(0,z)*struve_l(-1,z)
1132       ;;               +bessel_k(1,z)*struve_l(0,z))
1133       ;;    + 2 * sum((-1)^(k+n/2)*bessel_k(2*k+1,z),k,0,n/2-1)
1134       (let*
1135           ((k (gensym))
1136            (answer `((mplus)
1137                      ((mtimes) 2
1138                       ((%sum)
1139                        ((mtimes)
1140                         ((%bessel_k) ((mplus) 1 ((mtimes) 2 ,k)) ,z)
1141                         ((mexpt) -1
1142                          ((mplus) ,k ((mtimes) ((rat) 1 2) ,n))))
1143                        ,k 0 ((mplus) -1 ((mtimes) ((rat) 1 2) ,n))))
1144                      ((mtimes)
1145                       ((rat) 1 2)
1146                       $%pi
1147                       ((mexpt) -1 ((mtimes) ((rat) 1 2) ,n))
1148                       z
1149                       ((mplus)
1150                        ((mtimes)
1151                         ((%bessel_k) 0 ,z)
1152                         ((%struve_l) -1 ,z))
1153                        ((mtimes)
1154                         ((%bessel_k) 1 ,z)
1155                         ((%struve_l) 0 ,z)))))))
1156         ;; expand out the sum if n < 10.  Otherwise fix up the indices
1157         (if (< n 10)
1158             (meval `(($ev) ,answer $sum))  ; Is there a better way?
1159           (simplify ($niceindices answer)))))))
1160    (t nil)))
1161
1162 (putprop '%bessel_k `((n z) nil ,#'bessel-k-integral-2) 'integral)
1163
1164 ;; Support a simplim%bessel_k function to handle specific limits
1165
1166 (defprop %bessel_k simplim%bessel_k simplim%function)
1167
1168 (defun simplim%bessel_k (expr var val)
1169   ;; Look for the limit of the arguments.
1170   (let ((v (limit (cadr expr) var val 'think))
1171         (z (limit (caddr expr) var val 'think)))
1172   (cond
1173     ;; Handle an argument 0 at this place.
1174     ((or (zerop1 z)
1175          (eq z '$zeroa)
1176          (eq z '$zerob))
1177      (cond ((zerop1 v)
1178             ;; bessel_k(0,0)
1179             '$inf)
1180            ((integerp v)
1181             ;; bessel_k(n,0), n an integer
1182             (cond ((evenp v) '$inf)
1183                   (t (cond ((eq z '$zeroa) '$inf)
1184                            ((eq z '$zerob) '$minf)
1185                            (t '$infinity)))))
1186            ((not (zerop1 ($realpart v)))
1187             ;; bessel_k(v,0), Re(v)#0
1188             '$infinity)
1189            ((and (zerop1 ($realpart v))
1190                  (not (zerop1 v)))
1191             ;; bessel_k(v,0), Re(v)=0 and v#0
1192             '$und)
1193            ;; Call the simplifier of the function.
1194            (t (simplify (list '(%bessel_k) v z)))))
1195     ((eq z '$inf)
1196      ;; bessel_k(v,inf)
1197      0)
1198     ((eq z '$minf)
1199      ;; bessel_k(v,minf)
1200      '$infinity)
1201     (t
1202      ;; All other cases are handled by the simplifier of the function.
1203      (simplify (list '(%bessel_k) v z))))))
1204
1205 (def-simplifier bessel_k (order arg)
1206   (let ((rat-order nil))
1207     (cond
1208       ((zerop1 arg)
1209        ;; Domain error for a zero argument.
1210        (simp-domain-error
1211          (intl:gettext "bessel_k: bessel_k(~:M,~:M) is undefined.") order arg))
1212
1213       ((complex-float-numerical-eval-p order arg)
1214        (cond ((= 0 ($imagpart order))
1215               ;; order is real, arg is real or complex
1216               (complexify (bessel-k ($float order)
1217                                     (complex ($float ($realpart arg))
1218                                              ($float ($imagpart arg))))))
1219              (t
1220               ;; order is complex, arg is real or complex
1221               (let (($numer t)
1222                     ($float t)
1223                     (order ($float order))
1224                     (arg ($float arg)))
1225                 ($float
1226                   ($rectform
1227                     (bessel-k-hypergeometric order arg)))))))
1228
1229       ((or (bigfloat-numerical-eval-p order arg)
1230            (complex-bigfloat-numerical-eval-p order arg))
1231        ;; When the order is not an integer, we can use the
1232        ;; hypergeometric representation to evaluate it.
1233        (destructuring-bind (re . im)
1234            (risplit order)
1235          (cond ((and (zerop1 im)
1236                      (not (zerop1 (sub re (take '($floor) re)))))
1237                 ($rectform
1238                  ($bfloat (bessel-k-hypergeometric re arg))))
1239                (t
1240                 (give-up)))))
1241
1242       ((mminusp order)
1243        ;; A&S 9.6.6: K[-v](x) = K[v](x)
1244        (take '(%bessel_k) (mul -1 order) arg))
1245
1246       ((and $besselexpand
1247             (setq rat-order (max-numeric-ratio-p order 2)))
1248        ;; When order is a fraction with a denominator of 2, we
1249        ;; can express the result in terms of elementary
1250        ;; functions.  From A&S 10.2.16 and 10.2.17:
1251        ;;
1252        ;; K[1/2](z) = sqrt(%pi/2/z)*exp(-z)
1253        ;;           = K[-1/2](z)
1254        (bessel-k-half-order rat-order arg))
1255
1256       ((and $bessel_reduce
1257             (and (integerp order)
1258                  (plusp order)
1259                  (> order 1)))
1260        ;; Reduce a bessel function of order > 2 to order 1 and 0.
1261        ;; A&S 9.6.26: bessel_k(v,z) = 2*(v-1)/z*bessel_k(v-1,z)+bessel_k(v-2,z)
1262        (add (mul 2
1263                  (- order 1)
1264                  (inv arg)
1265                  (take '(%bessel_k) (- order 1) arg))
1266             (take '(%bessel_k) (- order 2) arg)))
1267
1268       ($hypergeometric_representation
1269        (bessel-k-hypergeometric order arg))
1270
1271       (t
1272        (give-up)))))
1273
1274 ;; Returns the hypergeometric representation of bessel_k
1275 (defun bessel-k-hypergeometric (order arg)
1276   ;; http://functions.wolfram.com/03.04.26.0002.01
1277   ;; hypergeometric([],[1-v],z^2/4)*2^(v-1)*gamma(v)/z^v
1278   ;;   + hypergeometric([],[v+1],z^2/4)*2^(-v-1)*gamma(-v)*z^v
1279   (add (mul (power 2 (sub order 1))
1280             (take '(%gamma) order)
1281             (inv (power arg order))
1282             (take '($hypergeometric)
1283                   (list '(mlist))
1284                   (list '(mlist) (sub 1 order))
1285                   (div (mul arg arg) 4)))
1286        (mul (inv (power 2 (add order 1)))
1287             (power arg order)
1288             (take '(%gamma) (neg order))
1289             (take '($hypergeometric)
1290                   (list '(mlist))
1291                   (list '(mlist) (add order 1))
1292                   (div (mul arg arg) 4)))))
1293
1294 ;; Compute value of Modified Bessel function of the second kind of order ORDER
1295 (defun bessel-k (order arg)
1296   (cond
1297     ((zerop (imagpart arg))
1298      ;; We have numeric args and the first arg is purely real. Call the
1299      ;; real-valued Bessel function. Handle orders 0 and 1 specially, when
1300      ;; possible.
1301      (let ((arg (realpart arg)))
1302        (cond
1303          ((< arg 0)
1304           ;; This is the extension for negative arg.
1305           ;; We use A&S 9.6.6 and 9.6.31 for evaluation:
1306           ;; K[v](-z) = K[|v|](-z)
1307           ;;          = exp(-%i*%pi*|v|)*K[|v|](z) - %i*%pi*I[|v|](z)
1308           (let* ((dpi (coerce pi 'flonum))
1309                  (s1 (cis (* dpi (- (abs order)))))
1310                  (s2 (* (complex 0 -1) dpi))
1311                  (result (+ (* s1 (bessel-k (abs order) (- arg)))
1312                             (* s2 (bessel-i (abs order) (- arg)))))
1313                  (rem (nth-value 1 (floor order))))
1314             (cond ((zerop rem)
1315                    ;; order is an integer or a float representation of an
1316                    ;; integer, the result is a general complex
1317                    result)
1318                   ((= rem 1/2)
1319                    ;; order is half-integral-value or an float representation
1320                    ;; and arg  < 0, the result is pure imaginary
1321                    (complex 0 (imagpart result)))
1322                   ;; in all other cases general complex result
1323                   (t result))))
1324          ((= order 0)
1325           (slatec:dbesk0 (float arg)))
1326          ((= order 1)
1327           (slatec:dbesk1 (float arg)))
1328          (t
1329           ;; From A&S 9.6.6, K[-v](z) = K[v](z), so take the
1330           ;; absolute value of the order.
1331           (multiple-value-bind (n alpha) (floor (abs (float order)))
1332             (let ((jvals (make-array (1+ n) :element-type 'flonum)))
1333               (slatec:dbesk (float arg) alpha 1 (1+ n) jvals 0)
1334               (aref jvals n)))))))
1335     (t
1336      ;; The arg is complex.  Use the complex-valued Bessel function. From
1337      ;; A&S 9.6.6, K[-v](z) = K[v](z), so take the absolute value of the order.
1338      (multiple-value-bind (n alpha) (floor (abs (float order)))
1339        (let ((cyr (make-array (1+ n) :element-type 'flonum))
1340              (cyi (make-array (1+ n) :element-type 'flonum)))
1341          (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
1342                                     v-cyr v-cyi v-nz v-ierr)
1343             (slatec::zbesk (float (realpart arg))
1344                            (float (imagpart arg))
1345                            alpha 1 (1+ n) cyr cyi 0 0)
1346            (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz))
1347
1348            ;; We should check for errors here based on the value of v-ierr.
1349            (when (plusp v-ierr)
1350              (format t "zbesk ierr = ~A~%" v-ierr))
1351            (complex (aref cyr n) (aref cyi n))))))))
1352
1353 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1354 ;;;
1355 ;;; Expansion of Bessel J and Y functions of half-integral order.
1356 ;;;
1357 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1358 ;;
1359 ;; From A&S 10.1.1, we have
1360 ;;
1361 ;; J[n+1/2](z) = sqrt(2*z/%pi)*j[n](z)
1362 ;; Y[n+1/2](z) = sqrt(2*z/%pi)*y[n](z)
1363 ;;
1364 ;; where j[n](z) is the spherical bessel function of the first kind
1365 ;; and y[n](z) is the spherical bessel function of the second kind.
1366 ;;
1367 ;; A&S 10.1.8 and 10.1.9 give
1368 ;;
1369 ;; j[n](z) = 1/z*[P(n+1/2,z)*sin(z-n*%pi/2) + Q(n+1/2,z)*cos(z-n*%pi/2)]
1370 ;;
1371 ;; y[n](z) = (-1)^(n+1)*1/z*[P(n+1/2,z)*cos(z+n*%pi/2) - Q(n+1/2,z)*sin(z+n*%pi/2)]
1372 ;;
1373
1374 ;; A&S 10.1.10
1375 ;;
1376 ;; j[n](z) = f[n](z)*sin(z) + (-1)^(n+1)*f[-n-1](z)*cos(z)
1377 ;;
1378 ;; f[0](z) = 1/z, f[1](z) = 1/z^2
1379 ;;
1380 ;; f[n-1](z) + f[n+1](z) = (2*n+1)/z*f[n](z)
1381 ;;
1382 (defun f-fun (n z)
1383   (cond ((= n 0)
1384          (div 1 z))
1385         ((= n 1)
1386          (div 1 (mul z z)))
1387         ((= n -1)
1388          0)
1389         ((>= n 2)
1390          ;; f[n+1](z) = (2*n+1)/z*f[n](z) - f[n-1](z) or
1391          ;; f[n](z) = (2*n-1)/z*f[n-1](z) - f[n-2](z)
1392          (sub (mul (div (+ n n -1) z)
1393                    (f-fun (1- n) z))
1394               (f-fun (- n 2) z)))
1395         (t
1396          ;; Negative n
1397          ;;
1398          ;; f[n-1](z) = (2*n+1)/z*f[n](z) - f[n+1](z) or
1399          ;; f[n](z) = (2*n+3)/z*f[n+1](z) - f[n+2](z)
1400          (sub (mul (div (+ n n 3) z)
1401                    (f-fun (1+ n) z))
1402               (f-fun (+ n 2) z)))))
1403
1404 ;; Compute sqrt(2*z/%pi)
1405 (defun root-2z/pi (z)
1406   (power (mul 2 z (inv '$%pi)) '((rat simp) 1 2)))
1407
1408 (defun bessel-j-half-order (order arg)
1409   "Compute J[n+1/2](z)"
1410   (let* ((n (floor order))
1411          (sign (if (oddp n) -1 1))
1412          (jn (sub (mul ($expand (f-fun n arg))
1413                        (take '(%sin) arg))
1414                   (mul sign
1415                        ($expand (f-fun (- (- n) 1) arg))
1416                        (take '(%cos) arg)))))
1417     (mul (root-2z/pi arg)
1418          jn)))
1419
1420 (defun bessel-y-half-order (order arg)
1421   "Compute Y[n+1/2](z)"
1422   ;; A&S 10.1.1:
1423   ;; Y[n+1/2](z) = sqrt(2*z/%pi)*y[n](z)
1424   ;;
1425   ;; A&S 10.1.15:
1426   ;; y[n](z) = (-1)^(n+1)*j[-n-1](z)
1427   ;;
1428   ;; So
1429   ;; Y[n+1/2](z) = sqrt(2*z/%pi)*(-1)^(n+1)*j[-n-1](z)
1430   ;;             = (-1)^(n+1)*sqrt(2*z/%pi)*j[-n-1](z)
1431   ;;             = (-1)^(n+1)*J[-n-1/2](z)
1432   (let* ((n (floor order))
1433          (jn (bessel-j-half-order (- (- order) 1/2) arg)))
1434     (if (evenp n)
1435         (mul -1 jn)
1436         jn)))
1437
1438 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1439 ;;;
1440 ;;; Expansion of Bessel I and K functions of half-integral order.
1441 ;;;
1442 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1443
1444 ;; See A&S 10.2.12
1445 ;;
1446 ;; I[n+1/2](z) = sqrt(2*z/%pi)*[g[n](z)*sinh(z) + g[-n-1](z)*cosh(z)]
1447 ;;
1448 ;; g[0](z) = 1/z, g[1](z) = -1/z^2
1449 ;;
1450 ;; g[n-1](z) - g[n+1](z) = (2*n+1)/z*g[n](z)
1451 ;;
1452 ;;
1453 (defun g-fun (n z)
1454   (declare (type integer n))
1455   (cond ((= n 0)
1456          (div 1 z))
1457         ((= n 1)
1458          (div -1 (mul z z)))
1459         ((>= n 2)
1460          ;; g[n](z) = g[n-2](z) - (2*n-1)/z*g[n-1](z)
1461          (sub (g-fun (- n 2) z)
1462               (mul (div (+ n n -1) z)
1463                    (g-fun (- n 1) z))))
1464         (t
1465          ;; n is negative
1466          ;;
1467          ;; g[n](z) = (2*n+3)/z*g[n+1](z) + g[n+2](z)
1468          (add (mul (div (+ n n 3) z)
1469                    (g-fun (+ n 1) z))
1470               (g-fun (+ n 2) z)))))
1471
1472 ;; See A&S 10.2.12
1473 ;;
1474 ;; I[n+1/2](z) = sqrt(2*z/%pi)*[g[n](z)*sinh(z) + g[-n-1](z)*cosh(z)]
1475
1476 (defun bessel-i-half-order (order arg)
1477   (let ((order (floor order)))
1478     (mul (root-2z/pi arg)
1479          (add (mul ($expand (g-fun order arg))
1480                    (take '(%sinh) arg))
1481               (mul ($expand (g-fun (- (+ order 1)) arg))
1482                    (take '(%cosh) arg))))))
1483
1484 ;; See A&S 10.2.15
1485 ;;
1486 ;; sqrt(%pi/2/z)*K[n+1/2](z) = (%pi/2/z)*exp(-z)*sum((n+1/2,k)/(2*z)^k,k,0,n)
1487 ;;
1488 ;; or
1489 ;;
1490 ;; K[n+1/2](z) = sqrt(%pi/2)/sqrt(z)*exp(-z)*sum((n+1/2,k)/(2*z)^k,k,0,n)
1491 ;;
1492 ;; where (A&S 10.1.9)
1493 ;;
1494 ;; (n+1/2,k) = (n+k)!/k!/(n-k)!
1495
1496 (defun k-fun (n z)
1497   (declare (type unsigned-byte n))
1498   ;; Computes the sum above
1499   (let ((sum 1)
1500         (term 1))
1501     (loop for k from 0 upto n do
1502           (setf term (mul term
1503                           (div (div (* (- n k) (+ n k 1))
1504                                     (+ k 1))
1505                                (mul 2 z))))
1506           (setf sum (add sum term)))
1507     sum))
1508
1509 (defun bessel-k-half-order (order arg)
1510   (let ((order (truncate (abs order))))
1511     (mul (mul (power (div '$%pi 2) '((rat simp) 1 2))
1512               (power arg '((rat simp) -1 2))
1513               (power '$%e (neg arg)))
1514          (k-fun (abs order) arg))))
1515
1516 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1517 ;;;
1518 ;;; Realpart und imagpart of Bessel functions
1519 ;;;
1520 ;;; We handle the special cases when we know that the Bessel function
1521 ;;; is pure real or pure imaginary. In all other cases Maxima generates
1522 ;;; a general noun form as result.
1523 ;;;
1524 ;;; To get the complex sign of the order an argument of the Bessel function
1525 ;;; the function $csign is used which calls $sign in a complex mode.
1526 ;;;
1527 ;;; This is an extension of of the algorithm of the function risplit in
1528 ;;; the file rpart.lisp. risplit looks for a risplit-function on the
1529 ;;; property list and call it if available.
1530 ;;;
1531 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1532
1533 (defvar *debug-bessel* nil)
1534
1535 ;;; Put the risplit-function for Bessel J and Bessel I on the property list
1536
1537 (defprop %bessel_j risplit-bessel-j-or-i risplit-function)
1538 (defprop %bessel_i risplit-bessel-j-or-i risplit-function)
1539
1540 ;;; realpart und imagpart for Bessel J and Bessel I function
1541
1542 (defun risplit-bessel-j-or-i (expr)
1543   (when *debug-bessel*
1544     (format t "~&RISPLIT-BESSEL-J with ~A~%" expr))
1545   (let ((order (cadr  expr))
1546         (arg   (caddr expr))
1547         (sign-order ($csign (cadr  expr)))
1548         (sign-arg   ($csign (caddr expr))))
1549
1550     (when *debug-bessel*
1551       (format t "~&   : order      = ~A~%" order)
1552       (format t "~&   : arg        = ~A~%" arg)
1553       (format t "~&   : sign-order = ~A~%" sign-order)
1554       (format t "~&   : sign-arg   = ~A~%" sign-arg))
1555
1556     (cond
1557       ((or (member sign-order '($complex $imaginary))
1558            (eq sign-arg '$complex))
1559        ;; order or arg are complex, return general noun-form
1560        (risplit-noun expr))
1561       ((eq sign-arg '$imaginary)
1562        ;; arg is pure imaginary
1563        (cond
1564          ((or ($oddp  order)
1565               ($featurep order '$odd))
1566           ;; order is an odd integer, pure imaginary noun-form
1567           (cons 0 (mul -1 '$%i expr)))
1568          ((or ($evenp order)
1569               ($featurep order '$even))
1570            ;; order is an even integer, real noun-form
1571            (cons expr 0))
1572          (t
1573           ;; order is not an odd or even integer, or Maxima can not
1574           ;; determine it, return general noun-form
1575           (risplit-noun expr))))
1576       ;; At this point order and arg are real.
1577       ;; We have to look for some special cases involing a negative arg
1578       ((or (maxima-integerp order)
1579            (member sign-arg '($pos $pz)))
1580        ;; arg is positive or order an integer, real noun-form
1581        (cons expr 0))
1582       ;; At this point we know that arg is negative or the sign is not known
1583       ((zerop1 (sub ($truncate ($multthru 2 order)) ($multthru 2 order)))
1584        ;; order is half integral
1585        (cond
1586           ((eq sign-arg '$neg)
1587            ;; order is half integral and arg negative, imaginary noun-form
1588            (cons 0 (mul -1 '$%i expr)))
1589           (t
1590             ;; the sign of arg or order is not known
1591             (risplit-noun expr))))
1592       (t
1593         ;; the sign of arg or order is not known
1594         (risplit-noun expr)))))
1595
1596 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1597
1598 ;;; Put the risplit-function for Bessel K and Bessel Y on the property list
1599
1600 (defprop %bessel_k risplit-bessel-k-or-y risplit-function)
1601 (defprop %bessel_y risplit-bessel-k-or-y risplit-function)
1602
1603 ;;; realpart und imagpart for Bessel K and Bessel Y function
1604
1605 (defun risplit-bessel-k-or-y (expr)
1606   (when *debug-bessel*
1607     (format t "~&RISPLIT-BESSEL-K with ~A~%" expr))
1608   (let ((order (cadr  expr))
1609         (arg   (caddr expr))
1610         (sign-order ($csign (cadr  expr)))
1611         (sign-arg   ($csign (caddr expr))))
1612
1613     (when *debug-bessel*
1614       (format t "~&   : order      = ~A~%" order)
1615       (format t "~&   : arg        = ~A~%" arg)
1616       (format t "~&   : sign-order = ~A~%" sign-order)
1617       (format t "~&   : sign-arg   = ~A~%" sign-arg))
1618
1619     (cond
1620       ((or (member sign-order '($complex $imaginary))
1621            (member sign-arg '($complex '$imaginary)))
1622        (risplit-noun expr))
1623       ;; At this point order and arg are real valued.
1624       ;; We have to look for some special cases involing a negative arg
1625       ((member sign-arg '($pos $pz))
1626        ;; arg is positive
1627        (cons expr 0))
1628       ;; At this point we know that arg is negative or the sign is not known
1629       ((and (not (maxima-integerp order))
1630             (zerop1 (sub ($truncate ($multthru 2 order)) ($multthru 2 order))))
1631        ;; order is half integral
1632        (cond
1633           ((eq sign-arg '$neg)
1634            (cons 0 (mul -1 '$%i expr)))
1635           (t
1636             ;; the sign of arg is not known
1637             (risplit-noun expr))))
1638       (t
1639         ;; the sign of arg is not known
1640         (risplit-noun expr)))))
1641
1642 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1643 ;;;
1644 ;;; Implementation of scaled Bessel I functions
1645 ;;;
1646 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1647
1648 ;; FIXME: The following scaled functions need work.  They should be
1649 ;; extended to match bessel_i, but still carefully compute the scaled
1650 ;; value.
1651
1652 ;; I think g0(x) = exp(-x)*I[0](x), g1(x) = exp(-x)*I[1](x), and
1653 ;; gn(x,n) = exp(-x)*I[n](x), based on some simple numerical
1654 ;; evaluations.
1655
1656 (defmfun $scaled_bessel_i0 ($x)
1657   (cond ((mnump $x)
1658          ;; XXX Should we return noun forms if $x is rational?
1659          (slatec:dbsi0e ($float $x)))
1660         (t
1661          (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil)))
1662               `((%bessel_i) 0 ,$x)))))
1663
1664 (defmfun $scaled_bessel_i1 ($x)
1665   (cond ((mnump $x)
1666          ;; XXX Should we return noun forms if $x is rational?
1667          (slatec:dbsi1e ($float $x)))
1668         (t
1669          (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil)))
1670               `((%bessel_i) 1 ,$x)))))
1671
1672 (defmfun $scaled_bessel_i ($n $x)
1673   (cond ((and (mnump $x) (mnump $n))
1674          ;; XXX Should we return noun forms if $n and $x are rational?
1675          (multiple-value-bind (n alpha) (floor ($float $n))
1676            (let ((iarray (make-array (1+ n) :element-type 'flonum)))
1677            (slatec:dbesi ($float $x) alpha 2 (1+ n) iarray 0)
1678            (aref iarray n))))
1679         (t
1680          (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil)))
1681               ($bessel_i $n $x)))))
1682
1683 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1684 ;;;
1685 ;;; Implementation of the Hankel 1 function
1686 ;;;
1687 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1688
1689 (defmfun $hankel_1 (v z)
1690   (simplify (list '(%hankel_1) v z)))
1691
1692 ;; hankel_1 distributes over lists, matrices, and equations
1693
1694 (defprop %hankel_1 (mlist $matrix mequal) distribute_over)
1695
1696 #+nil
1697 (defprop %hankel_1 simp-hankel-1 operators)
1698
1699 ; Derivatives of the Hankel 1 function
1700 (defprop %hankel_1
1701     ((n x)
1702      nil
1703
1704      ;; Derivative wrt arg x.  A&S 9.1.27.
1705      ((mtimes)
1706        ((mplus)
1707          ((%hankel_1) ((mplus) -1 n) x)
1708          ((mtimes) -1 ((%hankel_1) ((mplus) 1 n) x)))
1709        ((rat) 1 2)))
1710     grad)
1711
1712 (def-simplifier hankel_1 (order arg)
1713   (let (rat-order)
1714     (cond
1715       ((zerop1 arg)
1716        (simp-domain-error
1717          (intl:gettext "hankel_1: hankel_1(~:M,~:M) is undefined.")
1718          order arg))
1719       ((complex-float-numerical-eval-p order arg)
1720        (cond ((= 0 ($imagpart order))
1721               ;; order is real, arg is real or complex
1722               (complexify (hankel-1 ($float order)
1723                                     (complex ($float ($realpart arg))
1724                                              ($float ($imagpart arg))))))
1725              (t
1726               ;; The order is complex.  Use
1727               ;;   hankel_1(v,z) = bessel_j(v,z) + %i*bessel_y(v,z)
1728               ;; and evaluate using the hypergeometric function
1729               (let (($numer t)
1730                     ($float t)
1731                     (order ($float order))
1732                     (arg ($float arg)))
1733                 ($float
1734                   ($rectform
1735                     (add (bessel-j-hypergeometric order arg)
1736                          (mul '$%i
1737                               (bessel-y-hypergeometric order arg)))))))))
1738       ((and $besselexpand
1739             (setq rat-order (max-numeric-ratio-p order 2)))
1740        ;; When order is a fraction with a denominator of 2, we can express
1741        ;; the result in terms of elementary functions.
1742        ;; Use the definition hankel_1(v,z) = bessel_j(v,z)+%i*bessel_y(v,z)
1743        (sratsimp
1744          (add (bessel-j-half-order rat-order arg)
1745               (mul '$%i
1746                    (bessel-y-half-order rat-order arg)))))
1747       (t (give-up)))))
1748
1749 ;; Numerically compute H1[v](z).
1750 ;;
1751 ;; A&S 9.1.3 says H1[v](z) = J[v](z) + %i * Y[v](z)
1752 ;;
1753 (defun hankel-1 (v z)
1754   (let ((v (float v))
1755         (z (coerce z '(complex flonum))))
1756     (cond ((minusp v)
1757            ;; A&S 9.1.6:
1758            ;;
1759            ;; H1[-v](z) = exp(v*%pi*%i)*H1[v](z)
1760            ;;
1761            ;; or
1762            ;;
1763            ;; H1[v](z) = exp(-v*%pi*%i)*H1[-v](z)
1764
1765            (* (cis (* (float pi) (- v))) (hankel-1 (- v) z)))
1766           (t
1767            (multiple-value-bind (n fnu) (floor v)
1768              (let ((zr (realpart z))
1769                    (zi (imagpart z))
1770                    (cyr (make-array (1+ n) :element-type 'flonum))
1771                    (cyi (make-array (1+ n) :element-type 'flonum)))
1772                (multiple-value-bind (dzr dzi df dk dm dn dcyr dcyi nz ierr)
1773                   (slatec::zbesh zr zi fnu 1 1 (1+ n) cyr cyi 0 0)
1774                  (declare (ignore dzr dzi df dk dm dn dcyr dcyi nz ierr))
1775                  (complex (aref cyr n) (aref cyi n)))))))))
1776
1777 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1778 ;;;
1779 ;;; Implementation of the Hankel 2 function
1780 ;;;
1781 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1782
1783 (defmfun $hankel_2 (v z)
1784   (simplify (list '(%hankel_2) v z)))
1785
1786 ;; hankel_2 distributes over lists, matrices, and equations
1787
1788 (defprop %hankel_2 (mlist $matrix mequal) distribute_over)
1789
1790 #+nil
1791 (defprop %hankel_2 simp-hankel-2 operators)
1792
1793 ; Derivatives of the Hankel 2 function
1794 (defprop %hankel_2
1795     ((n x)
1796      nil
1797
1798      ;; Derivative wrt arg x.  A&S 9.1.27.
1799      ((mtimes)
1800        ((mplus)
1801          ((%hankel_2) ((mplus) -1 n) x)
1802          ((mtimes) -1 ((%hankel_2) ((mplus) 1 n) x)))
1803        ((rat) 1 2)))
1804     grad)
1805
1806 (def-simplifier hankel_2 (order arg)
1807   (let (rat-order)
1808     (cond
1809       ((zerop1 arg)
1810        (simp-domain-error
1811          (intl:gettext "hankel_2: hankel_2(~:M,~:M) is undefined.")
1812          order arg))
1813       ((complex-float-numerical-eval-p order arg)
1814        (cond ((= 0 ($imagpart order))
1815               ;; order is real, arg is real or complex
1816               (complexify (hankel-2 ($float order)
1817                                     (complex ($float ($realpart arg))
1818                                              ($float ($imagpart arg))))))
1819              (t
1820               ;; The order is complex.  Use
1821               ;;   hankel_2(v,z) = bessel_j(v,z) - %i*bessel_y(v,z)
1822               ;; and evaluate using the hypergeometric function
1823               (let (($numer t)
1824                     ($float t)
1825                     (order ($float order))
1826                     (arg ($float arg)))
1827                 ($float
1828                   ($rectform
1829                     (sub (bessel-j-hypergeometric order arg)
1830                          (mul '$%i
1831                               (bessel-y-hypergeometric order arg)))))))))
1832       ((and $besselexpand
1833             (setq rat-order (max-numeric-ratio-p order 2)))
1834        ;; When order is a fraction with a denominator of 2, we can express
1835        ;; the result in terms of elementary functions.
1836        ;; Use the definition hankel_2(v,z) = bessel_j(v,z)-%i*bessel_y(v,z)
1837        (sratsimp
1838          (sub (bessel-j-half-order rat-order arg)
1839               (mul '$%i
1840                    (bessel-y-half-order rat-order arg)))))
1841       (t (give-up)))))
1842
1843 ;; Numerically compute H2[v](z).
1844 ;;
1845 ;; A&S 9.1.4 says H2[v](z) = J[v](z) - %i * Y[v](z)
1846 ;;
1847 (defun hankel-2 (v z)
1848   (let ((v (float v))
1849         (z (coerce z '(complex flonum))))
1850     (cond ((minusp v)
1851            ;; A&S 9.1.6:
1852            ;;
1853            ;; H2[-v](z) = exp(-v*%pi*%i)*H2[v](z)
1854            ;;
1855            ;; or
1856            ;;
1857            ;; H2[v](z) = exp(v*%pi*%i)*H2[-v](z)
1858
1859            (* (cis (* (float pi) v)) (hankel-2 (- v) z)))
1860           (t
1861            (multiple-value-bind (n fnu) (floor v)
1862              (let ((zr (realpart z))
1863                    (zi (imagpart z))
1864                    (cyr (make-array (1+ n) :element-type 'flonum))
1865                    (cyi (make-array (1+ n) :element-type 'flonum)))
1866                (multiple-value-bind (dzr dzi df dk dm dn dcyr dcyi nz ierr)
1867                   (slatec::zbesh zr zi fnu 1 2 (1+ n) cyr cyi 0 0)
1868                  (declare (ignore dzr dzi df dk dm dn dcyr dcyi nz ierr))
1869                  (complex (aref cyr n) (aref cyi n)))))))))
1870
1871 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1872 ;;;
1873 ;;; Implementation of Struve H function
1874 ;;;
1875 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1876
1877 (defmfun $struve_h (v z)
1878   (simplify (list '(%struve_h) v z)))
1879
1880 ;; struve_h distributes over lists, matrices, and equations
1881
1882 (defprop %struve_h (mlist $matrix mequal) distribute_over)
1883
1884 #+nil
1885 (defprop %struve_h simp-struve-h operators)
1886
1887 ; Derivatives of the Struve H function
1888 (defprop %struve_h
1889   ((v z)
1890    nil
1891
1892    ;; Derivative wrt arg z.  A&S 12.1.10.
1893    ((mtimes)
1894     ((rat) 1 2)
1895     ((mplus)
1896      ((%struve_h) ((mplus) -1 v) z)
1897      ((mtimes) -1 ((%struve_h) ((mplus) 1 v) z))
1898      ((mtimes)
1899       ((mexpt) $%pi ((rat) -1 2))
1900       ((mexpt) 2 ((mtimes) -1 v))
1901       ((mexpt) ((%gamma) ((mplus) ((rat) 3 2) v)) -1)
1902       ((mexpt) z v)))))
1903   grad)
1904
1905 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1906
1907 (def-simplifier struve_h (order arg)
1908   (cond
1909
1910     ;; Check for special values
1911
1912     ((zerop1 arg)
1913      (cond ((eq ($csign (add order 1)) '$zero)
1914                                         ; http://functions.wolfram.com/03.09.03.0018.01
1915             (cond ((or ($bfloatp order)
1916                        ($bfloatp arg))
1917                    ($bfloat (div 2 '$%pi)))
1918                   ((or (floatp order)
1919                        (floatp arg))
1920                    ($float (div 2 '$%pi)))
1921                   (t
1922                    (div 2 '$%pi))))
1923            ((eq ($sign (add ($realpart order) 1)) '$pos)
1924                                         ; http://functions.wolfram.com/03.09.03.0001.01
1925             arg)
1926            ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz))
1927             (simp-domain-error
1928              (intl:gettext "struve_h: struve_h(~:M,~:M) is undefined.")
1929              order arg))
1930            (t
1931             (give-up))))
1932
1933     ;; Check for numerical evaluation
1934
1935     ((complex-float-numerical-eval-p order arg)
1936                                         ; A&S 12.1.21
1937      (let (($numer t) ($float t))
1938        ($rectform
1939         ($float
1940          (mul
1941           ($rectform (power arg (add order 1.0)))
1942           ($rectform (inv (power 2.0 order)))
1943           (inv (power ($float '$%pi) 0.5))
1944           (inv (simplify (list '(%gamma) (add order 1.5))))
1945           (simplify (list '($hypergeometric)
1946                           (list '(mlist) 1)
1947                           (list '(mlist) '((rat simp) 3 2)
1948                                 (add order '((rat simp) 3 2)))
1949                           (div (mul arg arg) -4.0))))))))
1950
1951     ((complex-bigfloat-numerical-eval-p order arg)
1952                                         ; A&S 12.1.21
1953      (let (($ratprint nil)
1954            (arg ($bfloat arg))
1955            (order ($bfloat order)))
1956        ($rectform
1957         ($bfloat
1958          (mul
1959           ($rectform (power arg (add order 1)))
1960           ($rectform (inv (power 2 order)))
1961           (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2))))
1962           (inv (simplify (list '(%gamma)
1963                                (add order ($bfloat '((rat simp) 3 2))))))
1964           (simplify (list '($hypergeometric)
1965                           (list '(mlist) 1)
1966                           (list '(mlist) '((rat simp) 3 2)
1967                                 (add order '((rat simp) 3 2)))
1968                           (div (mul arg arg) ($bfloat -4)))))))))
1969
1970     ;; Transformations and argument simplifications
1971
1972     ((and $besselexpand
1973           (ratnump order)
1974           (integerp (mul 2 order)))
1975      (cond
1976        ((eq ($sign order) '$pos)
1977         ;; Expansion of Struve H for a positive half integral order.
1978         (sratsimp
1979          (add
1980           (mul
1981            (inv (simplify (list '(mfactorial) (sub order
1982                                                    '((rat simp) 1 2)))))
1983            (inv (power '$%pi '((rat simp) 1 2 )))
1984            (power (div arg 2) (add order -1))
1985            (let ((index (gensumindex)))
1986              (dosum
1987               (mul
1988                (simplify (list '($pochhammer) '((rat simp) 1 2) index))
1989                (simplify (list '($pochhammer)
1990                                (sub '((rat simp) 1 2) order)
1991                                index))
1992                (power (mul -1 arg arg (inv 4)) (mul -1 index)))
1993               index 0 (sub order '((rat simp) 1 2)) t)))
1994           (mul
1995            (power (div 2 '$%pi) '((rat simp) 1 2))
1996            (power -1 (add order '((rat simp) 1 2)))
1997            (inv (power arg '((rat simp) 1 2)))
1998            (add
1999             (mul
2000              (simplify
2001               (list '(%sin)
2002                     (add (mul '((rat simp) 1 2)
2003                               '$%pi
2004                               (add order '((rat simp) 1 2)))
2005                          arg)))
2006              (let ((index (gensumindex)))
2007                (dosum
2008                 (mul
2009                  (power -1 index)
2010                  (simplify (list '(mfactorial)
2011                                  (add (mul 2 index)
2012                                       order
2013                                       '((rat simp) -1 2))))
2014                  (inv (simplify (list '(mfactorial) (mul 2 index))))
2015                  (inv (simplify (list '(mfactorial)
2016                                       (add (mul -2 index)
2017                                            order
2018                                            '((rat simp) -1 2)))))
2019                  (inv (power (mul 2 arg) (mul 2 index))))
2020                 index 0
2021                 (simplify (list '($floor)
2022                                 (div (sub (mul 2 order) 1) 4)))
2023                 t)))
2024             (mul
2025              (simplify (list '(%cos)
2026                              (add (mul '((rat simp) 1 2)
2027                                        '$%pi
2028                                        (add order '((rat simp) 1 2)))
2029                                   arg)))
2030              (let ((index (gensumindex)))
2031                (dosum
2032                 (mul
2033                  (power -1 index)
2034                  (simplify (list '(mfactorial)
2035                                  (add (mul 2 index)
2036                                       order
2037                                       '((rat simp) 1 2))))
2038                  (power (mul 2 arg) (mul -1 (add (mul 2 index) 1)))
2039                  (inv (simplify (list '(mfactorial)
2040                                       (add (mul 2 index) 1))))
2041                  (inv (simplify (list '(mfactorial)
2042                                       (add (mul -2 index)
2043                                            order
2044                                            '((rat simp) -3 2))))))
2045                 index 0
2046                 (simplify (list '($floor)
2047                                 (div (sub (mul 2 order) 3) 4)))
2048                 t))))))))
2049
2050        ((eq ($sign order) '$neg)
2051         ;; Expansion of Struve H for a negative half integral order.
2052         (sratsimp
2053          (add
2054           (mul
2055            (power (div 2 '$%pi) '((rat simp) 1 2))
2056            (power -1 (add order '((rat simp) 1 2)))
2057            (inv (power arg '((rat simp) 1 2)))
2058            (add
2059             (mul
2060              (simplify (list '(%sin)
2061                              (add
2062                               (mul
2063                                '((rat simp) 1 2)
2064                                '$%pi
2065                                (add order '((rat simp) 1 2)))
2066                               arg)))
2067              (let ((index (gensumindex)))
2068                (dosum
2069                 (mul
2070                  (power -1 index)
2071                  (simplify (list '(mfactorial)
2072                                  (add (mul 2 index)
2073                                       (neg order)
2074                                       '((rat simp) -1 2))))
2075                  (inv (simplify (list '(mfactorial) (mul 2 index))))
2076                  (inv (simplify (list '(mfactorial)
2077                                       (add (mul -2 index)
2078                                            (neg order)
2079                                            '((rat simp) -1 2)))))
2080                  (inv (power (mul 2 arg) (mul 2 index))))
2081                 index 0
2082                 (simplify (list '($floor)
2083                                 (div (add (mul 2 order) 1) -4)))
2084                 t)))
2085             (mul
2086              (simplify (list '(%cos)
2087                              (add
2088                               (mul
2089                                '((rat simp) 1 2)
2090                                '$%pi
2091                                (add order '((rat simp) 1 2)))
2092                               arg)))
2093              (let ((index (gensumindex)))
2094                (dosum
2095                 (mul
2096                  (power -1 index)
2097                  (simplify (list '(mfactorial)
2098                                  (add (mul 2 index)
2099                                       (neg order)
2100                                       '((rat simp) 1 2))))
2101                  (power (mul 2 arg) (mul -1 (add (mul 2 index) 1)))
2102                  (inv (simplify (list '(mfactorial)
2103                                       (add (mul 2 index) 1))))
2104                  (inv (simplify (list '(mfactorial)
2105                                       (add (mul -2 index)
2106                                            (neg order)
2107                                            '((rat simp) -3 2))))))
2108                 index 0
2109                 (simplify (list '($floor)
2110                                 (div (add (mul 2 order) 3) -4)))
2111                 t))))))))))
2112     (t
2113      (give-up))))
2114
2115 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2116 ;;;
2117 ;;; Implementation of Struve L function
2118 ;;;
2119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2120
2121 (defmfun $struve_l (v z)
2122   (simplify (list '(%struve_l) v z)))
2123
2124 ;; struve_l distributes over lists, matrices, and equations
2125
2126 (defprop %struve_l (mlist $matrix mequal) distribute_over)
2127
2128 ; Derivatives of the Struve L function
2129 (defprop %struve_l
2130   ((v z)
2131    nil
2132
2133    ;; Derivative wrt arg z.  A&S 12.2.5.
2134    ((mtimes)
2135     ((rat) 1 2)
2136     ((mplus)
2137      ((%struve_l) ((mplus) -1 v) z)
2138      ((%struve_l) ((mplus) 1 v) z)
2139      ((mtimes)
2140       ((mexpt) $%pi ((rat) -1 2))
2141       ((mexpt) 2 ((mtimes) -1 v))
2142       ((mexpt) ((%gamma) ((mplus) ((rat) 3 2) v)) -1)
2143       ((mexpt) z v)))))
2144   grad)
2145
2146 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2147
2148 (def-simplifier struve_l (order arg)
2149   (cond
2150
2151     ;; Check for special values
2152
2153     ((zerop1 arg)
2154      (cond ((eq ($csign (add order 1)) '$zero)
2155                                         ; http://functions.wolfram.com/03.10.03.0018.01
2156             (cond ((or ($bfloatp order)
2157                        ($bfloatp arg))
2158                    ($bfloat (div 2 '$%pi)))
2159                   ((or (floatp order)
2160                        (floatp arg))
2161                    ($float (div 2 '$%pi)))
2162                   (t
2163                    (div 2 '$%pi))))
2164            ((eq ($sign (add ($realpart order) 1)) '$pos)
2165                                         ; http://functions.wolfram.com/03.10.03.0001.01
2166             arg)
2167            ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz))
2168             (simp-domain-error
2169              (intl:gettext "struve_l: struve_l(~:M,~:M) is undefined.")
2170              order arg))
2171            (t
2172             (give-up))))
2173
2174     ;; Check for numerical evaluation
2175
2176     ((complex-float-numerical-eval-p order arg)
2177                                         ; http://functions.wolfram.com/03.10.26.0002.01
2178      (let (($numer t) ($float t))
2179        ($rectform
2180         ($float
2181          (mul
2182           ($rectform (power arg (add order 1.0)))
2183           ($rectform (inv (power 2.0 order)))
2184           (inv (power ($float '$%pi) 0.5))
2185           (inv (simplify (list '(%gamma) (add order 1.5))))
2186           (simplify (list '($hypergeometric)
2187                           (list '(mlist) 1)
2188                           (list '(mlist) '((rat simp) 3 2)
2189                                 (add order '((rat simp) 3 2)))
2190                           (div (mul arg arg) 4.0))))))))
2191
2192     ((complex-bigfloat-numerical-eval-p order arg)
2193                                         ; http://functions.wolfram.com/03.10.26.0002.01
2194      (let (($ratprint nil))
2195        ($rectform
2196         ($bfloat
2197          (mul
2198           ($rectform (power arg (add order 1)))
2199           ($rectform (inv (power 2 order)))
2200           (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2))))
2201           (inv (simplify (list '(%gamma) (add order '((rat simp) 3 2)))))
2202           (simplify (list '($hypergeometric)
2203                           (list '(mlist) 1)
2204                           (list '(mlist) '((rat simp) 3 2)
2205                                 (add order '((rat simp) 3 2)))
2206                           (div (mul arg arg) 4))))))))
2207
2208     ;; Transformations and argument simplifications
2209
2210     ((and $besselexpand
2211           (ratnump order)
2212           (integerp (mul 2 order)))
2213      (cond
2214        ((eq ($sign order) '$pos)
2215         ;; Expansion of Struve L for a positive half integral order.
2216         (sratsimp
2217          (add
2218           (mul -1
2219                (power 2 (sub 1 order))
2220                (power arg (sub order 1))
2221                (inv (power '$%pi '((rat simp) 1 2)))
2222                (inv (simplify (list '(mfactorial) (sub order
2223                                                        '((rat simp) 1 2)))))
2224                (let ((index (gensumindex)))
2225                  (dosum
2226                   (mul
2227                    (simplify (list '($pochhammer) '((rat simp) 1 2) index))
2228                    (simplify (list '($pochhammer)
2229                                    (sub '((rat simp) 1 2) order)
2230                                    index))
2231                    (power (mul arg arg (inv 4)) (mul -1 index)))
2232                   index 0 (sub order '((rat simp) 1 2)) t)))
2233           (mul -1
2234                (power (div 2 '$%pi) '((rat simp) 1 2))
2235                (inv (power arg '((rat simp) 1 2)))
2236                ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2))
2237                (add
2238                 (mul
2239                  (let (($trigexpand t))
2240                    (simplify (list '(%sinh)
2241                                    (sub (mul '((rat simp) 1 2)
2242                                              '$%pi
2243                                              '$%i
2244                                              (add order '((rat simp) 1 2)))
2245                                         arg))))
2246                  (let ((index (gensumindex)))
2247                    (dosum
2248                     (mul
2249                      (simplify (list '(mfactorial)
2250                                      (add (mul 2 index)
2251                                           (simplify (list '(mabs) order))
2252                                           '((rat simp) -1 2))))
2253                      (inv (simplify (list '(mfactorial) (mul 2 index))))
2254                      (inv (simplify (list '(mfactorial)
2255                                           (add (simplify (list '(mabs)
2256                                                                order))
2257                                                (mul -2 index)
2258                                                '((rat simp) -1 2)))))
2259                      (inv (power (mul 2 arg) (mul 2 index))))
2260                     index 0
2261                     (simplify (list '($floor)
2262                                     (div (sub (mul 2
2263                                                    (simplify (list '(mabs)
2264                                                                    order)))
2265                                               1)
2266                                          4)))
2267                     t)))
2268                 (mul
2269                  (let (($trigexpand t))
2270                    (simplify (list '(%cosh)
2271                                    (sub (mul '((rat simp) 1 2)
2272                                              '$%pi
2273                                              '$%i
2274                                              (add order '((rat simp) 1 2)))
2275                                         arg))))
2276                  (let ((index (gensumindex)))
2277                    (dosum
2278                     (mul
2279                      (simplify (list '(mfactorial)
2280                                      (add (mul 2 index)
2281                                           (simplify (list '(mabs) order))
2282                                           '((rat simp) 1 2))))
2283                      (power (mul 2 arg) (neg (add (mul 2 index) 1)))
2284                      (inv (simplify (list '(mfactorial)
2285                                           (add (mul 2 index) 1))))
2286                      (inv (simplify (list '(mfactorial)
2287                                           (add (simplify (list '(mabs)
2288                                                                order))
2289                                                (mul -2 index)
2290                                                '((rat simp) -3 2))))))
2291                     index 0
2292                     (simplify (list '($floor)
2293                                     (div (sub (mul 2
2294                                                    (simplify (list '(mabs)
2295                                                                    order)))
2296                                               3)
2297                                          4)))
2298                     t))))))))
2299        ((eq ($sign order) '$neg)
2300         ;; Expansion of Struve L for a negative half integral order.
2301         (sratsimp
2302          (add
2303           (mul -1
2304                (power (div 2 '$%pi) '((rat simp) 1 2))
2305                (inv (power arg '((rat simp) 1 2)))
2306                ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2))
2307                (add
2308                 (mul
2309                  (let (($trigexpand t))
2310                    (simplify (list '(%sinh)
2311                                    (sub (mul '((rat simp) 1 2)
2312                                              '$%pi
2313                                              '$%i
2314                                              (add order '((rat simp) 1 2)))
2315                                         arg))))
2316                  (let ((index (gensumindex)))
2317                    (dosum
2318                     (mul
2319                      (simplify (list '(mfactorial)
2320                                      (add (mul 2 index)
2321                                           (neg order)
2322                                           '((rat simp) -1 2))))
2323                      (inv (simplify (list '(mfactorial) (mul 2 index))))
2324                      (inv (simplify (list '(mfactorial)
2325                                           (add (neg order)
2326                                                (mul -2 index)
2327                                                '((rat simp) -1 2)))))
2328                      (inv (power (mul 2 arg) (mul 2 index))))
2329                     index 0
2330                     (simplify (list '($floor)
2331                                     (div (add (mul 2 order) 1) -4)))
2332                     t)))
2333                 (mul
2334                  (let (($trigexpand t))
2335                    (simplify (list '(%cosh)
2336                                    (sub (mul '((rat simp) 1 2)
2337                                              '$%pi
2338                                              '$%i
2339                                              (add order '((rat simp) 1 2)))
2340                                         arg))))
2341                  (let ((index (gensumindex)))
2342                    (dosum
2343                     (mul
2344                      (simplify (list '(mfactorial)
2345                                      (add (mul 2 index)
2346                                           (neg order)
2347                                           '((rat simp) 1 2))))
2348                      (power (mul 2 arg) (neg (add (mul 2 index) 1)))
2349                      (inv (simplify (list '(mfactorial)
2350                                           (add (mul 2 index) 1))))
2351                      (inv (simplify (list '(mfactorial)
2352                                           (add (neg order)
2353                                                (mul -2 index)
2354                                                '((rat simp) -3 2))))))
2355                     index 0
2356                     (simplify (list '($floor)
2357                                     (div (add (mul 2 order) 3) -4)))
2358                     t))))))))))
2359     (t
2360      (give-up))))
2361
2362 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2363
2364 (defmspec $gauss (form)
2365   (format t
2366 "NOTE: The gauss function is superseded by random_normal in the `distrib' package.
2367 Perhaps you meant to enter `~a'.~%"
2368     (print-invert-case (implode (mstring `(($random_normal) ,@ (cdr form))))))
2369   '$done)