code.delx.au - gnu-emacs/blob - lisp/calc/calc-funcs.el

   1 ;;; calc-funcs.el --- well-known functions for Calc
   2
   3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
   4
   5 ;; Author: David Gillespie <daveg@synaptics.com>
   6 ;; Maintainer: Colin Walters <walters@debian.org>
   7
   8 ;; This file is part of GNU Emacs.
   9
  10 ;; GNU Emacs is distributed in the hope that it will be useful,
  11 ;; but WITHOUT ANY WARRANTY.  No author or distributor
  12 ;; accepts responsibility to anyone for the consequences of using it
  13 ;; or for whether it serves any particular purpose or works at all,
  14 ;; unless he says so in writing.  Refer to the GNU Emacs General Public
  15 ;; License for full details.
  16
  17 ;; Everyone is granted permission to copy, modify and redistribute
  18 ;; GNU Emacs, but only under the conditions described in the
  19 ;; GNU Emacs General Public License.   A copy of this license is
  20 ;; supposed to have been given to you along with GNU Emacs so you
  21 ;; can know your rights and responsibilities.  It should be in a
  22 ;; file named COPYING.  Among other things, the copyright notice
  23 ;; and this notice must be preserved on all copies.
  24
  25 ;;; Commentary:
  26
  27 ;;; Code:
  28
  29 ;; This file is autoloaded from calc-ext.el.
  30 (require 'calc-ext)
  31
  32 (require 'calc-macs)
  33
  34 (defun calc-Need-calc-funcs () nil)
  35
  36
  37 (defun calc-inc-gamma (arg)
  38   (interactive "P")
  39   (calc-slow-wrapper
  40    (if (calc-is-inverse)
  41        (if (calc-is-hyperbolic)
  42            (calc-binary-op "gamG" 'calcFunc-gammaG arg)
  43          (calc-binary-op "gamQ" 'calcFunc-gammaQ arg))
  44        (if (calc-is-hyperbolic)
  45            (calc-binary-op "gamg" 'calcFunc-gammag arg)
  46          (calc-binary-op "gamP" 'calcFunc-gammaP arg)))))
  47
  48 (defun calc-erf (arg)
  49   (interactive "P")
  50   (calc-slow-wrapper
  51    (if (calc-is-inverse)
  52        (calc-unary-op "erfc" 'calcFunc-erfc arg)
  53      (calc-unary-op "erf" 'calcFunc-erf arg))))
  54
  55 (defun calc-erfc (arg)
  56   (interactive "P")
  57   (calc-invert-func)
  58   (calc-erf arg))
  59
  60 (defun calc-beta (arg)
  61   (interactive "P")
  62   (calc-slow-wrapper
  63    (calc-binary-op "beta" 'calcFunc-beta arg)))
  64
  65 (defun calc-inc-beta ()
  66   (interactive)
  67   (calc-slow-wrapper
  68    (if (calc-is-hyperbolic)
  69        (calc-enter-result 3 "betB" (cons 'calcFunc-betaB (calc-top-list-n 3)))
  70      (calc-enter-result 3 "betI" (cons 'calcFunc-betaI (calc-top-list-n 3))))))
  71
  72 (defun calc-bessel-J (arg)
  73   (interactive "P")
  74   (calc-slow-wrapper
  75    (calc-binary-op "besJ" 'calcFunc-besJ arg)))
  76
  77 (defun calc-bessel-Y (arg)
  78   (interactive "P")
  79   (calc-slow-wrapper
  80    (calc-binary-op "besY" 'calcFunc-besY arg)))
  81
  82 (defun calc-bernoulli-number (arg)
  83   (interactive "P")
  84   (calc-slow-wrapper
  85    (if (calc-is-hyperbolic)
  86        (calc-binary-op "bern" 'calcFunc-bern arg)
  87      (calc-unary-op "bern" 'calcFunc-bern arg))))
  88
  89 (defun calc-euler-number (arg)
  90   (interactive "P")
  91   (calc-slow-wrapper
  92    (if (calc-is-hyperbolic)
  93        (calc-binary-op "eulr" 'calcFunc-euler arg)
  94      (calc-unary-op "eulr" 'calcFunc-euler arg))))
  95
  96 (defun calc-stirling-number (arg)
  97   (interactive "P")
  98   (calc-slow-wrapper
  99    (if (calc-is-hyperbolic)
 100        (calc-binary-op "str2" 'calcFunc-stir2 arg)
 101      (calc-binary-op "str1" 'calcFunc-stir1 arg))))
 102
 103 (defun calc-utpb ()
 104   (interactive)
 105   (calc-prob-dist "b" 3))
 106
 107 (defun calc-utpc ()
 108   (interactive)
 109   (calc-prob-dist "c" 2))
 110
 111 (defun calc-utpf ()
 112   (interactive)
 113   (calc-prob-dist "f" 3))
 114
 115 (defun calc-utpn ()
 116   (interactive)
 117   (calc-prob-dist "n" 3))
 118
 119 (defun calc-utpp ()
 120   (interactive)
 121   (calc-prob-dist "p" 2))
 122
 123 (defun calc-utpt ()
 124   (interactive)
 125   (calc-prob-dist "t" 2))
 126
 127 (defun calc-prob-dist (letter nargs)
 128   (calc-slow-wrapper
 129    (if (calc-is-inverse)
 130        (calc-enter-result nargs (concat "ltp" letter)
 131                           (append (list (intern (concat "calcFunc-ltp" letter))
 132                                         (calc-top-n 1))
 133                                   (calc-top-list-n (1- nargs) 2)))
 134      (calc-enter-result nargs (concat "utp" letter)
 135                         (append (list (intern (concat "calcFunc-utp" letter))
 136                                       (calc-top-n 1))
 137                                 (calc-top-list-n (1- nargs) 2))))))
 138
 139
 140
 141
 142 ;;; Sources:  Numerical Recipes, Press et al;
 143 ;;;           Handbook of Mathematical Functions, Abramowitz & Stegun.
 144
 145
 146 ;;; Gamma function.
 147
 148 (defun calcFunc-gamma (x)
 149   (or (math-numberp x) (math-reject-arg x 'numberp))
 150   (calcFunc-fact (math-add x -1)))
 151
 152 (defun math-gammap1-raw (x &optional fprec nfprec)   ; compute gamma(1 + x)
 153   (or fprec
 154       (setq fprec (math-float calc-internal-prec)
 155             nfprec (math-float (- calc-internal-prec))))
 156   (cond ((math-lessp-float (calcFunc-re x) fprec)
 157          (if (math-lessp-float (calcFunc-re x) nfprec)
 158              (math-neg (math-div
 159                         (math-pi)
 160                         (math-mul (math-gammap1-raw
 161                                    (math-add (math-neg x)
 162                                              '(float -1 0))
 163                                    fprec nfprec)
 164                                   (math-sin-raw
 165                                    (math-mul (math-pi) x)))))
 166            (let ((xplus1 (math-add x '(float 1 0))))
 167              (math-div (math-gammap1-raw xplus1 fprec nfprec) xplus1))))
 168         ((and (math-realp x)
 169               (math-lessp-float '(float 736276 0) x))
 170          (math-overflow))
 171         (t   ; re(x) now >= 10.0
 172          (let ((xinv (math-div 1 x))
 173                (lnx (math-ln-raw x)))
 174            (math-mul (math-sqrt-two-pi)
 175                      (math-exp-raw
 176                       (math-gamma-series
 177                        (math-sub (math-mul (math-add x '(float 5 -1))
 178                                            lnx)
 179                                  x)
 180                        xinv
 181                        (math-sqr xinv)
 182                        '(float 0 0)
 183                        2)))))))
 184
 185 (defun math-gamma-series (sum x xinvsqr oterm n)
 186   (math-working "gamma" sum)
 187   (let* ((bn (math-bernoulli-number n))
 188          (term (math-mul (math-div-float (math-float (nth 1 bn))
 189                                          (math-float (* (nth 2 bn)
 190                                                         (* n (1- n)))))
 191                          x))
 192          (next (math-add sum term)))
 193     (if (math-nearly-equal sum next)
 194         next
 195       (if (> n (* 2 calc-internal-prec))
 196           (progn
 197             ;; Need this because series eventually diverges for large enough n.
 198             (calc-record-why
 199              "*Gamma computation stopped early, not all digits may be valid")
 200             next)
 201         (math-gamma-series next (math-mul x xinvsqr) xinvsqr term (+ n 2))))))
 202
 203
 204 ;;; Incomplete gamma function.
 205
 206 (defvar math-current-gamma-value nil)
 207 (defun calcFunc-gammaP (a x)
 208   (if (equal x '(var inf var-inf))
 209       '(float 1 0)
 210     (math-inexact-result)
 211     (or (Math-numberp a) (math-reject-arg a 'numberp))
 212     (or (math-numberp x) (math-reject-arg x 'numberp))
 213     (if (and (math-num-integerp a)
 214              (integerp (setq a (math-trunc a)))
 215              (> a 0) (< a 20))
 216         (math-sub 1 (calcFunc-gammaQ a x))
 217       (let ((math-current-gamma-value (calcFunc-gamma a)))
 218         (math-div (calcFunc-gammag a x) math-current-gamma-value)))))
 219
 220 (defun calcFunc-gammaQ (a x)
 221   (if (equal x '(var inf var-inf))
 222       '(float 0 0)
 223     (math-inexact-result)
 224     (or (Math-numberp a) (math-reject-arg a 'numberp))
 225     (or (math-numberp x) (math-reject-arg x 'numberp))
 226     (if (and (math-num-integerp a)
 227              (integerp (setq a (math-trunc a)))
 228              (> a 0) (< a 20))
 229         (let ((n 0)
 230               (sum '(float 1 0))
 231               (term '(float 1 0)))
 232           (math-with-extra-prec 1
 233             (while (< (setq n (1+ n)) a)
 234               (setq term (math-div (math-mul term x) n)
 235                     sum (math-add sum term))
 236               (math-working "gamma" sum))
 237             (math-mul sum (calcFunc-exp (math-neg x)))))
 238       (let ((math-current-gamma-value (calcFunc-gamma a)))
 239         (math-div (calcFunc-gammaG a x) math-current-gamma-value)))))
 240
 241 (defun calcFunc-gammag (a x)
 242   (if (equal x '(var inf var-inf))
 243       (calcFunc-gamma a)
 244     (math-inexact-result)
 245     (or (Math-numberp a) (math-reject-arg a 'numberp))
 246     (or (Math-numberp x) (math-reject-arg x 'numberp))
 247     (math-with-extra-prec 2
 248       (setq a (math-float a))
 249       (setq x (math-float x))
 250       (if (or (math-negp (calcFunc-re a))
 251               (math-lessp-float (calcFunc-re x)
 252                                 (math-add-float (calcFunc-re a)
 253                                                 '(float 1 0))))
 254           (math-inc-gamma-series a x)
 255         (math-sub (or math-current-gamma-value (calcFunc-gamma a))
 256                   (math-inc-gamma-cfrac a x))))))
 257
 258 (defun calcFunc-gammaG (a x)
 259   (if (equal x '(var inf var-inf))
 260       '(float 0 0)
 261     (math-inexact-result)
 262     (or (Math-numberp a) (math-reject-arg a 'numberp))
 263     (or (Math-numberp x) (math-reject-arg x 'numberp))
 264     (math-with-extra-prec 2
 265       (setq a (math-float a))
 266       (setq x (math-float x))
 267       (if (or (math-negp (calcFunc-re a))
 268               (math-lessp-float (calcFunc-re x)
 269                                 (math-add-float (math-abs-approx a)
 270                                                 '(float 1 0))))
 271           (math-sub (or math-current-gamma-value (calcFunc-gamma a))
 272                     (math-inc-gamma-series a x))
 273         (math-inc-gamma-cfrac a x)))))
 274
 275 (defun math-inc-gamma-series (a x)
 276   (if (Math-zerop x)
 277       '(float 0 0)
 278     (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x))
 279               (math-with-extra-prec 2
 280                 (let ((start (math-div '(float 1 0) a)))
 281                   (math-inc-gamma-series-step start start a x))))))
 282
 283 (defun math-inc-gamma-series-step (sum term a x)
 284   (math-working "gamma" sum)
 285   (setq a (math-add a '(float 1 0))
 286         term (math-div (math-mul term x) a))
 287   (let ((next (math-add sum term)))
 288     (if (math-nearly-equal sum next)
 289         next
 290       (math-inc-gamma-series-step next term a x))))
 291
 292 (defun math-inc-gamma-cfrac (a x)
 293   (if (Math-zerop x)
 294       (or math-current-gamma-value (calcFunc-gamma a))
 295     (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x))
 296               (math-inc-gamma-cfrac-step '(float 1 0) x
 297                                          '(float 0 0) '(float 1 0)
 298                                          '(float 1 0) '(float 1 0) '(float 0 0)
 299                                          a x))))
 300
 301 (defun math-inc-gamma-cfrac-step (a0 a1 b0 b1 n fac g a x)
 302   (let ((ana (math-sub n a))
 303         (anf (math-mul n fac)))
 304     (setq n (math-add n '(float 1 0))
 305           a0 (math-mul (math-add a1 (math-mul a0 ana)) fac)
 306           b0 (math-mul (math-add b1 (math-mul b0 ana)) fac)
 307           a1 (math-add (math-mul x a0) (math-mul anf a1))
 308           b1 (math-add (math-mul x b0) (math-mul anf b1)))
 309     (if (math-zerop a1)
 310         (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac g a x)
 311       (setq fac (math-div '(float 1 0) a1))
 312       (let ((next (math-mul b1 fac)))
 313         (math-working "gamma" next)
 314         (if (math-nearly-equal next g)
 315             next
 316           (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac next a x))))))
 317
 318
 319 ;;; Error function.
 320
 321 (defun calcFunc-erf (x)
 322   (if (equal x '(var inf var-inf))
 323       '(float 1 0)
 324     (if (equal x '(neg (var inf var-inf)))
 325         '(float -1 0)
 326       (if (Math-zerop x)
 327           x
 328         (let ((math-current-gamma-value (math-sqrt-pi)))
 329           (math-to-same-complex-quad
 330            (math-div (calcFunc-gammag '(float 5 -1)
 331                                       (math-sqr (math-to-complex-quad-one x)))
 332                      math-current-gamma-value)
 333            x))))))
 334
 335 (defun calcFunc-erfc (x)
 336   (if (equal x '(var inf var-inf))
 337       '(float 0 0)
 338     (if (math-posp x)
 339         (let ((math-current-gamma-value (math-sqrt-pi)))
 340           (math-div (calcFunc-gammaG '(float 5 -1) (math-sqr x))
 341                     math-current-gamma-value))
 342       (math-sub 1 (calcFunc-erf x)))))
 343
 344 (defun math-to-complex-quad-one (x)
 345   (if (eq (car-safe x) 'polar) (setq x (math-complex x)))
 346   (if (eq (car-safe x) 'cplx)
 347       (list 'cplx (math-abs (nth 1 x)) (math-abs (nth 2 x)))
 348     x))
 349
 350 (defun math-to-same-complex-quad (x y)
 351   (if (eq (car-safe y) 'cplx)
 352       (if (eq (car-safe x) 'cplx)
 353           (list 'cplx
 354                 (if (math-negp (nth 1 y)) (math-neg (nth 1 x)) (nth 1 x))
 355                 (if (math-negp (nth 2 y)) (math-neg (nth 2 x)) (nth 2 x)))
 356         (if (math-negp (nth 1 y)) (math-neg x) x))
 357     (if (math-negp y)
 358         (if (eq (car-safe x) 'cplx)
 359             (list 'cplx (math-neg (nth 1 x)) (nth 2 x))
 360           (math-neg x))
 361       x)))
 362
 363
 364 ;;; Beta function.
 365
 366 (defun calcFunc-beta (a b)
 367   (if (math-num-integerp a)
 368       (let ((am (math-add a -1)))
 369         (or (math-numberp b) (math-reject-arg b 'numberp))
 370         (math-div 1 (math-mul b (calcFunc-choose (math-add b am) am))))
 371     (if (math-num-integerp b)
 372         (calcFunc-beta b a)
 373       (math-div (math-mul (calcFunc-gamma a) (calcFunc-gamma b))
 374                 (calcFunc-gamma (math-add a b))))))
 375
 376
 377 ;;; Incomplete beta function.
 378
 379 (defvar math-current-beta-value nil)
 380 (defun calcFunc-betaI (x a b)
 381   (cond ((math-zerop x)
 382          '(float 0 0))
 383         ((math-equal-int x 1)
 384          '(float 1 0))
 385         ((or (math-zerop a)
 386              (and (math-num-integerp a)
 387                   (math-negp a)))
 388          (if (or (math-zerop b)
 389                  (and (math-num-integerp b)
 390                       (math-negp b)))
 391              (math-reject-arg b 'range)
 392            '(float 1 0)))
 393         ((or (math-zerop b)
 394              (and (math-num-integerp b)
 395                   (math-negp b)))
 396          '(float 0 0))
 397         ((not (math-numberp a)) (math-reject-arg a 'numberp))
 398         ((not (math-numberp b)) (math-reject-arg b 'numberp))
 399         ((math-inexact-result))
 400         (t (let ((math-current-beta-value (calcFunc-beta a b)))
 401              (math-div (calcFunc-betaB x a b) math-current-beta-value)))))
 402
 403 (defun calcFunc-betaB (x a b)
 404   (cond
 405    ((math-zerop x)
 406     '(float 0 0))
 407    ((math-equal-int x 1)
 408     (calcFunc-beta a b))
 409    ((not (math-numberp x)) (math-reject-arg x 'numberp))
 410    ((not (math-numberp a)) (math-reject-arg a 'numberp))
 411    ((not (math-numberp b)) (math-reject-arg b 'numberp))
 412    ((math-zerop a) (math-reject-arg a 'nonzerop))
 413    ((math-zerop b) (math-reject-arg b 'nonzerop))
 414    ((and (math-num-integerp b)
 415          (if (math-negp b)
 416              (math-reject-arg b 'range)
 417            (Math-natnum-lessp (setq b (math-trunc b)) 20)))
 418     (and calc-symbolic-mode (or (math-floatp a) (math-floatp b))
 419          (math-inexact-result))
 420     (math-mul
 421      (math-with-extra-prec 2
 422        (let* ((i 0)
 423               (term 1)
 424               (sum (math-div term a)))
 425          (while (< (setq i (1+ i)) b)
 426            (setq term (math-mul (math-div (math-mul term (- i b)) i) x)
 427                  sum (math-add sum (math-div term (math-add a i))))
 428            (math-working "beta" sum))
 429          sum))
 430      (math-pow x a)))
 431    ((and (math-num-integerp a)
 432          (if (math-negp a)
 433              (math-reject-arg a 'range)
 434            (Math-natnum-lessp (setq a (math-trunc a)) 20)))
 435     (math-sub (or math-current-beta-value (calcFunc-beta a b))
 436               (calcFunc-betaB (math-sub 1 x) b a)))
 437    (t
 438     (math-inexact-result)
 439     (math-with-extra-prec 2
 440       (setq x (math-float x))
 441       (setq a (math-float a))
 442       (setq b (math-float b))
 443       (let ((bt (math-exp-raw (math-add (math-mul a (math-ln-raw x))
 444                                         (math-mul b (math-ln-raw
 445                                                      (math-sub '(float 1 0)
 446                                                                x)))))))
 447         (if (Math-lessp x (math-div (math-add a '(float 1 0))
 448                                     (math-add (math-add a b) '(float 2 0))))
 449             (math-div (math-mul bt (math-beta-cfrac a b x)) a)
 450           (math-sub (or math-current-beta-value (calcFunc-beta a b))
 451                     (math-div (math-mul bt
 452                                         (math-beta-cfrac b a (math-sub 1 x)))
 453                               b))))))))
 454
 455 (defun math-beta-cfrac (a b x)
 456   (let ((qab (math-add a b))
 457         (qap (math-add a '(float 1 0)))
 458         (qam (math-add a '(float -1 0))))
 459     (math-beta-cfrac-step '(float 1 0)
 460                           (math-sub '(float 1 0)
 461                                     (math-div (math-mul qab x) qap))
 462                           '(float 1 0) '(float 1 0)
 463                           '(float 1 0)
 464                           qab qap qam a b x)))
 465
 466 (defun math-beta-cfrac-step (az bz am bm m qab qap qam a b x)
 467   (let* ((two-m (math-mul m '(float 2 0)))
 468          (d (math-div (math-mul (math-mul (math-sub b m) m) x)
 469                       (math-mul (math-add qam two-m) (math-add a two-m))))
 470          (ap (math-add az (math-mul d am)))
 471          (bp (math-add bz (math-mul d bm)))
 472          (d2 (math-neg
 473               (math-div (math-mul (math-mul (math-add a m) (math-add qab m)) x)
 474                         (math-mul (math-add qap two-m) (math-add a two-m)))))
 475          (app (math-add ap (math-mul d2 az)))
 476          (bpp (math-add bp (math-mul d2 bz)))
 477          (next (math-div app bpp)))
 478     (math-working "beta" next)
 479     (if (math-nearly-equal next az)
 480         next
 481       (math-beta-cfrac-step next '(float 1 0)
 482                             (math-div ap bpp) (math-div bp bpp)
 483                             (math-add m '(float 1 0))
 484                             qab qap qam a b x))))
 485
 486
 487 ;;; Bessel functions.
 488
 489 ;;; Should generalize this to handle arbitrary precision!
 490
 491 (defun calcFunc-besJ (v x)
 492   (or (math-numberp v) (math-reject-arg v 'numberp))
 493   (or (math-numberp x) (math-reject-arg x 'numberp))
 494   (let ((calc-internal-prec (min 8 calc-internal-prec)))
 495     (math-with-extra-prec 3
 496       (setq x (math-float (math-normalize x)))
 497       (setq v (math-float (math-normalize v)))
 498       (cond ((math-zerop x)
 499              (if (math-zerop v)
 500                  '(float 1 0)
 501                '(float 0 0)))
 502             ((math-inexact-result))
 503             ((not (math-num-integerp v))
 504              (let ((start (math-div 1 (calcFunc-fact v))))
 505                (math-mul (math-besJ-series start start
 506                                            0
 507                                            (math-mul '(float -25 -2)
 508                                                      (math-sqr x))
 509                                            v)
 510                          (math-pow (math-div x 2) v))))
 511             ((math-negp (setq v (math-trunc v)))
 512              (if (math-oddp v)
 513                  (math-neg (calcFunc-besJ (math-neg v) x))
 514                (calcFunc-besJ (math-neg v) x)))
 515             ((eq v 0)
 516              (math-besJ0 x))
 517             ((eq v 1)
 518              (math-besJ1 x))
 519             ((Math-lessp v (math-abs-approx x))
 520              (let ((j 0)
 521                    (bjm (math-besJ0 x))
 522                    (bj (math-besJ1 x))
 523                    (two-over-x (math-div 2 x))
 524                    bjp)
 525                (while (< (setq j (1+ j)) v)
 526                  (setq bjp (math-sub (math-mul (math-mul j two-over-x) bj)
 527                                      bjm)
 528                        bjm bj
 529                        bj bjp))
 530                bj))
 531             (t
 532              (if (Math-lessp 100 v) (math-reject-arg v 'range))
 533              (let* ((j (logior (+ v (math-isqrt-small (* 40 v))) 1))
 534                     (two-over-x (math-div 2 x))
 535                     (jsum nil)
 536                     (bjp '(float 0 0))
 537                     (sum '(float 0 0))
 538                     (bj '(float 1 0))
 539                     bjm ans)
 540                (while (> (setq j (1- j)) 0)
 541                  (setq bjm (math-sub (math-mul (math-mul j two-over-x) bj)
 542                                      bjp)
 543                        bjp bj
 544                        bj bjm)
 545                  (if (> (nth 2 (math-abs-approx bj)) 10)
 546                      (setq bj (math-mul bj '(float 1 -10))
 547                            bjp (math-mul bjp '(float 1 -10))
 548                            ans (and ans (math-mul ans '(float 1 -10)))
 549                            sum (math-mul sum '(float 1 -10))))
 550                  (or (setq jsum (not jsum))
 551                      (setq sum (math-add sum bj)))
 552                  (if (= j v)
 553                      (setq ans bjp)))
 554                (math-div ans (math-sub (math-mul 2 sum) bj))))))))
 555
 556 (defun math-besJ-series (sum term k zz vk)
 557   (math-working "besJ" sum)
 558   (setq k (1+ k)
 559         vk (math-add 1 vk)
 560         term (math-div (math-mul term zz) (math-mul k vk)))
 561   (let ((next (math-add sum term)))
 562     (if (math-nearly-equal next sum)
 563         next
 564       (math-besJ-series next term k zz vk))))
 565
 566 (defun math-besJ0 (x &optional yflag)
 567   (cond ((and (not yflag) (math-negp (calcFunc-re x)))
 568          (math-besJ0 (math-neg x)))
 569         ((Math-lessp '(float 8 0) (math-abs-approx x))
 570          (let* ((z (math-div '(float 8 0) x))
 571                 (y (math-sqr z))
 572                 (xx (math-add x '(float (bigneg 164 398 785) -9)))
 573                 (a1 (math-poly-eval y
 574                                     '((float (bigpos 211 887 093 2) -16)
 575                                       (float (bigneg 639 370 073 2) -15)
 576                                       (float (bigpos 407 510 734 2) -14)
 577                                       (float (bigneg 627 628 098 1) -12)
 578                                       (float 1 0))))
 579                 (a2 (math-poly-eval y
 580                                     '((float (bigneg 152 935 934) -16)
 581                                       (float (bigpos 161 095 621 7) -16)
 582                                       (float (bigneg 651 147 911 6) -15)
 583                                       (float (bigpos 765 488 430 1) -13)
 584                                       (float (bigneg 995 499 562 1) -11))))
 585                 (sc (math-sin-cos-raw xx)))
 586                (if yflag
 587                    (setq sc (cons (math-neg (cdr sc)) (car sc))))
 588                (math-mul (math-sqrt
 589                           (math-div '(float (bigpos 722 619 636) -9) x))
 590                          (math-sub (math-mul (cdr sc) a1)
 591                                    (math-mul (car sc) (math-mul z a2))))))
 592          (t
 593           (let ((y (math-sqr x)))
 594             (math-div (math-poly-eval y
 595                                       '((float (bigneg 456 052 849 1) -7)
 596                                         (float (bigpos 017 233 739 7) -5)
 597                                         (float (bigneg 418 442 121 1) -2)
 598                                         (float (bigpos 407 196 516 6) -1)
 599                                         (float (bigneg 354 590 362 13) 0)
 600                                         (float (bigpos 574 490 568 57) 0)))
 601                       (math-poly-eval y
 602                                       '((float 1 0)
 603                                         (float (bigpos 712 532 678 2) -7)
 604                                         (float (bigpos 853 264 927 5) -5)
 605                                         (float (bigpos 718 680 494 9) -3)
 606                                         (float (bigpos 985 532 029 1) 0)
 607                                         (float (bigpos 411 490 568 57) 0))))))))
 608
 609 (defun math-besJ1 (x &optional yflag)
 610   (cond ((and (math-negp (calcFunc-re x)) (not yflag))
 611          (math-neg (math-besJ1 (math-neg x))))
 612         ((Math-lessp '(float 8 0) (math-abs-approx x))
 613          (let* ((z (math-div '(float 8 0) x))
 614                 (y (math-sqr z))
 615                 (xx (math-add x '(float (bigneg 491 194 356 2) -9)))
 616                 (a1 (math-poly-eval y
 617                                     '((float (bigneg 019 337 240) -15)
 618                                       (float (bigpos 174 520 457 2) -15)
 619                                       (float (bigneg 496 396 516 3) -14)
 620                                       (float 183105 -8)
 621                                       (float 1 0))))
 622                 (a2 (math-poly-eval y
 623                                     '((float (bigpos 412 787 105) -15)
 624                                       (float (bigneg 987 228 88) -14)
 625                                       (float (bigpos 096 199 449 8) -15)
 626                                       (float (bigneg 873 690 002 2) -13)
 627                                       (float (bigpos 995 499 687 4) -11))))
 628                 (sc (math-sin-cos-raw xx)))
 629            (if yflag
 630                (setq sc (cons (math-neg (cdr sc)) (car sc)))
 631              (if (math-negp x)
 632                  (setq sc (cons (math-neg (car sc)) (math-neg (cdr sc))))))
 633            (math-mul (math-sqrt (math-div '(float (bigpos 722 619 636) -9) x))
 634                      (math-sub (math-mul (cdr sc) a1)
 635                                (math-mul (car sc) (math-mul z a2))))))
 636         (t
 637          (let ((y (math-sqr x)))
 638            (math-mul
 639             x
 640             (math-div (math-poly-eval y
 641                                       '((float (bigneg 606 036 016 3) -8)
 642                                         (float (bigpos 826 044 157) -4)
 643                                         (float (bigneg 439 611 972 2) -3)
 644                                         (float (bigpos 531 968 423 2) -1)
 645                                         (float (bigneg 235 059 895 7) 0)
 646                                         (float (bigpos 232 614 362 72) 0)))
 647                       (math-poly-eval y
 648                                       '((float 1 0)
 649                                         (float (bigpos 397 991 769 3) -7)
 650                                         (float (bigpos 394 743 944 9) -5)
 651                                         (float (bigpos 474 330 858 1) -2)
 652                                         (float (bigpos 178 535 300 2) 0)
 653                                         (float (bigpos 442 228 725 144)
 654                                                0)))))))))
 655
 656 (defun calcFunc-besY (v x)
 657   (math-inexact-result)
 658   (or (math-numberp v) (math-reject-arg v 'numberp))
 659   (or (math-numberp x) (math-reject-arg x 'numberp))
 660   (let ((calc-internal-prec (min 8 calc-internal-prec)))
 661     (math-with-extra-prec 3
 662       (setq x (math-float (math-normalize x)))
 663       (setq v (math-float (math-normalize v)))
 664       (cond ((not (math-num-integerp v))
 665              (let ((sc (math-sin-cos-raw (math-mul v (math-pi)))))
 666                (math-div (math-sub (math-mul (calcFunc-besJ v x) (cdr sc))
 667                                    (calcFunc-besJ (math-neg v) x))
 668                          (car sc))))
 669             ((math-negp (setq v (math-trunc v)))
 670              (if (math-oddp v)
 671                  (math-neg (calcFunc-besY (math-neg v) x))
 672                (calcFunc-besY (math-neg v) x)))
 673             ((eq v 0)
 674              (math-besY0 x))
 675             ((eq v 1)
 676              (math-besY1 x))
 677             (t
 678              (let ((j 0)
 679                    (bym (math-besY0 x))
 680                    (by (math-besY1 x))
 681                    (two-over-x (math-div 2 x))
 682                    byp)
 683                (while (< (setq j (1+ j)) v)
 684                  (setq byp (math-sub (math-mul (math-mul j two-over-x) by)
 685                                      bym)
 686                        bym by
 687                        by byp))
 688                by))))))
 689
 690 (defun math-besY0 (x)
 691   (cond ((Math-lessp (math-abs-approx x) '(float 8 0))
 692          (let ((y (math-sqr x)))
 693            (math-add
 694             (math-div (math-poly-eval y
 695                                       '((float (bigpos 733 622 284 2) -7)
 696                                         (float (bigneg 757 792 632 8) -5)
 697                                         (float (bigpos 129 988 087 1) -2)
 698                                         (float (bigneg 036 598 123 5) -1)
 699                                         (float (bigpos 065 834 062 7) 0)
 700                                         (float (bigneg 389 821 957 2) 0)))
 701                       (math-poly-eval y
 702                                       '((float 1 0)
 703                                         (float (bigpos 244 030 261 2) -7)
 704                                         (float (bigpos 647 472 474) -4)
 705                                         (float (bigpos 438 466 189 7) -3)
 706                                         (float (bigpos 648 499 452 7) -1)
 707                                         (float (bigpos 269 544 076 40) 0))))
 708             (math-mul '(float (bigpos 772 619 636) -9)
 709                       (math-mul (math-besJ0 x) (math-ln-raw x))))))
 710         ((math-negp (calcFunc-re x))
 711          (math-add (math-besJ0 (math-neg x) t)
 712                    (math-mul '(cplx 0 2)
 713                              (math-besJ0 (math-neg x)))))
 714         (t
 715          (math-besJ0 x t))))
 716
 717 (defun math-besY1 (x)
 718   (cond ((Math-lessp (math-abs-approx x) '(float 8 0))
 719          (let ((y (math-sqr x)))
 720            (math-add
 721             (math-mul
 722              x
 723              (math-div (math-poly-eval y
 724                                        '((float (bigpos 935 937 511 8) -6)
 725                                          (float (bigneg 726 922 237 4) -3)
 726                                          (float (bigpos 551 264 349 7) -1)
 727                                          (float (bigneg 139 438 153 5) 1)
 728                                          (float (bigpos 439 527 127) 4)
 729                                          (float (bigneg 943 604 900 4) 3)))
 730                        (math-poly-eval y
 731                                        '((float 1 0)
 732                                          (float (bigpos 885 632 549 3) -7)
 733                                          (float (bigpos 605 042 102) -3)
 734                                          (float (bigpos 002 904 245 2) -2)
 735                                          (float (bigpos 367 650 733 3) 0)
 736                                          (float (bigpos 664 419 244 4) 2)
 737                                          (float (bigpos 057 958 249) 5)))))
 738             (math-mul '(float (bigpos 772 619 636) -9)
 739                       (math-sub (math-mul (math-besJ1 x) (math-ln-raw x))
 740                                 (math-div 1 x))))))
 741         ((math-negp (calcFunc-re x))
 742          (math-neg
 743           (math-add (math-besJ1 (math-neg x) t)
 744                     (math-mul '(cplx 0 2)
 745                               (math-besJ1 (math-neg x))))))
 746         (t
 747          (math-besJ1 x t))))
 748
 749 (defun math-poly-eval (x coefs)
 750   (let ((accum (car coefs)))
 751     (while (setq coefs (cdr coefs))
 752       (setq accum (math-add (car coefs) (math-mul accum x))))
 753     accum))
 754
 755
 756 ;;;; Bernoulli and Euler polynomials and numbers.
 757
 758 (defun calcFunc-bern (n &optional x)
 759   (if (and x (not (math-zerop x)))
 760       (if (and calc-symbolic-mode (math-floatp x))
 761           (math-inexact-result)
 762         (math-build-polynomial-expr (math-bernoulli-coefs n) x))
 763     (or (math-num-natnump n) (math-reject-arg n 'natnump))
 764     (if (consp n)
 765         (progn
 766           (math-inexact-result)
 767           (math-float (math-bernoulli-number (math-trunc n))))
 768       (math-bernoulli-number n))))
 769
 770 (defun calcFunc-euler (n &optional x)
 771   (or (math-num-natnump n) (math-reject-arg n 'natnump))
 772   (if x
 773       (let* ((n1 (math-add n 1))
 774              (coefs (math-bernoulli-coefs n1))
 775              (fac (math-div (math-pow 2 n1) n1))
 776              (k -1)
 777              (x1 (math-div (math-add x 1) 2))
 778              (x2 (math-div x 2)))
 779         (if (math-numberp x)
 780             (if (and calc-symbolic-mode (math-floatp x))
 781                 (math-inexact-result)
 782               (math-mul fac
 783                         (math-sub (math-build-polynomial-expr coefs x1)
 784                                   (math-build-polynomial-expr coefs x2))))
 785           (calcFunc-collect
 786            (math-reduce-vec
 787             'math-add
 788             (cons 'vec
 789                   (mapcar (function
 790                            (lambda (c)
 791                              (setq k (1+ k))
 792                              (math-mul (math-mul fac c)
 793                                        (math-sub (math-pow x1 k)
 794                                                  (math-pow x2 k)))))
 795                           coefs)))
 796            x)))
 797     (math-mul (math-pow 2 n)
 798               (if (consp n)
 799                   (progn
 800                     (math-inexact-result)
 801                     (calcFunc-euler n '(float 5 -1)))
 802                 (calcFunc-euler n '(frac 1 2))))))
 803
 804 (defvar math-bernoulli-b-cache '((frac -174611
 805                                        (bigpos 0 200 291 698 662 857 802))
 806                                  (frac 43867 (bigpos 0 944 170 217 94 109 5))
 807                                  (frac -3617 (bigpos 0 880 842 622 670 10))
 808                                  (frac 1 (bigpos 600 249 724 74))
 809                                  (frac -691 (bigpos 0 368 674 307 1))
 810                                  (frac 1 (bigpos 160 900 47))
 811                                  (frac -1 (bigpos 600 209 1))
 812                                  (frac 1 30240) (frac -1 720)
 813                                  (frac 1 12) 1 ))
 814
 815 (defvar math-bernoulli-B-cache '((frac -174611 330) (frac 43867 798)
 816                                  (frac -3617 510) (frac 7 6) (frac -691 2730)
 817                                  (frac 5 66) (frac -1 30) (frac 1 42)
 818                                  (frac -1 30) (frac 1 6) 1 ))
 819
 820 (defvar math-bernoulli-cache-size 11)
 821 (defun math-bernoulli-coefs (n)
 822   (let* ((coefs (list (calcFunc-bern n)))
 823          (nn (math-trunc n))
 824          (k nn)
 825          (term nn)
 826          coef
 827          (calc-prefer-frac (or (integerp n) calc-prefer-frac)))
 828     (while (>= (setq k (1- k)) 0)
 829       (setq term (math-div term (- nn k))
 830             coef (math-mul term (math-bernoulli-number k))
 831             coefs (cons (if (consp n) (math-float coef) coef) coefs)
 832             term (math-mul term k)))
 833     (nreverse coefs)))
 834
 835 (defun math-bernoulli-number (n)
 836   (if (= (% n 2) 1)
 837       (if (= n 1)
 838           '(frac -1 2)
 839         0)
 840     (setq n (/ n 2))
 841     (while (>= n math-bernoulli-cache-size)
 842       (let* ((sum 0)
 843              (nk 1)     ; nk = n-k+1
 844              (fact 1)   ; fact = (n-k+1)!
 845              ofact
 846              (p math-bernoulli-b-cache)
 847              (calc-prefer-frac t))
 848         (math-working "bernoulli B" (* 2 math-bernoulli-cache-size))
 849         (while p
 850           (setq nk (+ nk 2)
 851                 ofact fact
 852                 fact (math-mul fact (* nk (1- nk)))
 853                 sum (math-add sum (math-div (car p) fact))
 854                 p (cdr p)))
 855         (setq ofact (math-mul ofact (1- nk))
 856               sum (math-sub (math-div '(frac 1 2) ofact) sum)
 857               math-bernoulli-b-cache (cons sum math-bernoulli-b-cache)
 858               math-bernoulli-B-cache (cons (math-mul sum ofact)
 859                                            math-bernoulli-B-cache)
 860               math-bernoulli-cache-size (1+ math-bernoulli-cache-size))))
 861     (nth (- math-bernoulli-cache-size n 1) math-bernoulli-B-cache)))
 862
 863 ;;;   Bn = n! bn
 864 ;;;   bn = - sum_k=0^n-1 bk / (n-k+1)!
 865
 866 ;;; A faster method would be to use "tangent numbers", c.f., Concrete
 867 ;;; Mathematics pg. 273.
 868
 869
 870 ;;; Probability distributions.
 871
 872 ;;; Binomial.
 873 (defun calcFunc-utpb (x n p)
 874   (if math-expand-formulas
 875       (math-normalize (list 'calcFunc-betaI p x (list '+ (list '- n x) 1)))
 876     (calcFunc-betaI p x (math-add (math-sub n x) 1))))
 877 (put 'calcFunc-utpb 'math-expandable t)
 878
 879 (defun calcFunc-ltpb (x n p)
 880   (math-sub 1 (calcFunc-utpb x n p)))
 881 (put 'calcFunc-ltpb 'math-expandable t)
 882
 883 ;;; Chi-square.
 884 (defun calcFunc-utpc (chisq v)
 885   (if math-expand-formulas
 886       (math-normalize (list 'calcFunc-gammaQ (list '/ v 2) (list '/ chisq 2)))
 887     (calcFunc-gammaQ (math-div v 2) (math-div chisq 2))))
 888 (put 'calcFunc-utpc 'math-expandable t)
 889
 890 (defun calcFunc-ltpc (chisq v)
 891   (if math-expand-formulas
 892       (math-normalize (list 'calcFunc-gammaP (list '/ v 2) (list '/ chisq 2)))
 893     (calcFunc-gammaP (math-div v 2) (math-div chisq 2))))
 894 (put 'calcFunc-ltpc 'math-expandable t)
 895
 896 ;;; F-distribution.
 897 (defun calcFunc-utpf (f v1 v2)
 898   (if math-expand-formulas
 899       (math-normalize (list 'calcFunc-betaI
 900                             (list '/ v2 (list '+ v2 (list '* v1 f)))
 901                             (list '/ v2 2)
 902                             (list '/ v1 2)))
 903     (calcFunc-betaI (math-div v2 (math-add v2 (math-mul v1 f)))
 904                     (math-div v2 2)
 905                     (math-div v1 2))))
 906 (put 'calcFunc-utpf 'math-expandable t)
 907
 908 (defun calcFunc-ltpf (f v1 v2)
 909   (math-sub 1 (calcFunc-utpf f v1 v2)))
 910 (put 'calcFunc-ltpf 'math-expandable t)
 911
 912 ;;; Normal.
 913 (defun calcFunc-utpn (x mean sdev)
 914   (if math-expand-formulas
 915       (math-normalize
 916        (list '/
 917              (list '+ 1
 918                    (list 'calcFunc-erf
 919                          (list '/ (list '- mean x)
 920                                (list '* sdev (list 'calcFunc-sqrt 2)))))
 921              2))
 922     (math-mul (math-add '(float 1 0)
 923                         (calcFunc-erf
 924                          (math-div (math-sub mean x)
 925                                    (math-mul sdev (math-sqrt-2)))))
 926               '(float 5 -1))))
 927 (put 'calcFunc-utpn 'math-expandable t)
 928
 929 (defun calcFunc-ltpn (x mean sdev)
 930   (if math-expand-formulas
 931       (math-normalize
 932        (list '/
 933              (list '+ 1
 934                    (list 'calcFunc-erf
 935                          (list '/ (list '- x mean)
 936                                (list '* sdev (list 'calcFunc-sqrt 2)))))
 937              2))
 938     (math-mul (math-add '(float 1 0)
 939                         (calcFunc-erf
 940                          (math-div (math-sub x mean)
 941                                    (math-mul sdev (math-sqrt-2)))))
 942               '(float 5 -1))))
 943 (put 'calcFunc-ltpn 'math-expandable t)
 944
 945 ;;; Poisson.
 946 (defun calcFunc-utpp (n x)
 947   (if math-expand-formulas
 948       (math-normalize (list 'calcFunc-gammaP x n))
 949     (calcFunc-gammaP x n)))
 950 (put 'calcFunc-utpp 'math-expandable t)
 951
 952 (defun calcFunc-ltpp (n x)
 953   (if math-expand-formulas
 954       (math-normalize (list 'calcFunc-gammaQ x n))
 955     (calcFunc-gammaQ x n)))
 956 (put 'calcFunc-ltpp 'math-expandable t)
 957
 958 ;;; Student's t.  (As defined in Abramowitz & Stegun and Numerical Recipes.)
 959 (defun calcFunc-utpt (tt v)
 960   (if math-expand-formulas
 961       (math-normalize (list 'calcFunc-betaI
 962                             (list '/ v (list '+ v (list '^ tt 2)))
 963                             (list '/ v 2)
 964                             '(float 5 -1)))
 965     (calcFunc-betaI (math-div v (math-add v (math-sqr tt)))
 966                     (math-div v 2)
 967                     '(float 5 -1))))
 968 (put 'calcFunc-utpt 'math-expandable t)
 969
 970 (defun calcFunc-ltpt (tt v)
 971   (math-sub 1 (calcFunc-utpt tt v)))
 972 (put 'calcFunc-ltpt 'math-expandable t)
 973
 974
 975 ;;; calc-funcs.el ends here