src/bessel.lisp

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