src/bessel.lisp

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