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;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module csimp2) 14 15(load-macsyma-macros rzmac) 16 17(declare-top (special var %p%i varlist plogabs half%pi nn* dn* $factlim 18 $beta_expand)) 19 20(defmvar $gammalim 10000 21 "Controls simplification of gamma for rational number arguments.") 22 23(defmvar $gamma_expand nil 24 "Expand gamma(z+n) for n an integer when T.") 25 26;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 27 28;;; Implementation of the plog function 29 30(defun simpplog (x vestigial z) 31 (declare (ignore vestigial)) 32 (prog (varlist dd check y) 33 (oneargcheck x) 34 (setq check x) 35 (setq x (simpcheck (cadr x) z)) 36 (cond ((equal 0 x) (merror (intl:gettext "plog: plog(0) is undefined."))) 37 ((among var x) ;This is used in DEFINT. 1/19/81. -JIM 38 (return (eqtest (list '(%plog) x) check)))) 39 (newvar x) 40 (cond 41 ((and (member '$%i varlist) 42 (not (some #'(lambda (v) 43 (and (atom v) (not (eq v '$%i)))) 44 varlist))) 45 (setq dd (trisplit x)) 46 (cond ((setq z (patan (car dd) (cdr dd))) 47 (return (add2* (simpln (list '(%log) 48 (simpexpt (list '(mexpt) 49 ($expand (list '(mplus) 50 (list '(mexpt) (car dd) 2) 51 (list '(mexpt) (cdr dd) 2))) 52 '((rat) 1 2)) 1 nil)) 1 t) 53 (list '(mtimes) z '$%i)))))) 54 ((and (free x '$%i) (eq ($sign x) '$pnz)) 55 (return (eqtest (list '(%plog) x) check))) 56 ((and (equal ($imagpart x) 0) (setq y ($asksign x))) 57 (cond ((eq y '$pos) (return (simpln (list '(%log) x) 1 t))) 58 ((and plogabs (eq y '$neg)) 59 (return (simpln (list '(%log) (list '(mtimes) -1 x)) 1 nil))) 60 ((eq y '$neg) 61 (return (add2 %p%i 62 (simpln (list '(%log) (list '(mtimes) -1 x)) 1 nil)))) 63 (t (merror (intl:gettext "plog: plog(0) is undefined."))))) 64 ((and (equal ($imagpart (setq z (div* x '$%i))) 0) 65 (setq y ($asksign z))) 66 (cond 67 ((equal y '$zero) (merror (intl:gettext "plog: plog(0) is undefined."))) 68 (t (cond ((eq y '$pos) (setq y 1)) 69 ((eq y '$neg) (setq y -1))) 70 (return (add2* (simpln (list '(%log) 71 (list '(mtimes) y z)) 1 nil) 72 (list '(mtimes) y '((rat) 1 2) '$%i '$%pi))))))) 73 (return (eqtest (list '(%plog) x) check)))) 74 75(defun patan (r i) 76 (let (($numer $numer)) 77 (prog (a b var) 78 (setq i (simplifya i nil) r (simplifya r nil)) 79 (cond ((zerop1 r) 80 (if (floatp i) (setq $numer t)) 81 (setq i ($asksign i)) 82 (cond ((equal i '$pos) (return (simplify half%pi))) 83 ((equal i '$neg) 84 (return (mul2 -1 (simplify half%pi)))) 85 (t (merror (intl:gettext "plog: encountered atan(0/0)."))))) 86 ((zerop1 i) 87 (cond ((floatp r) (setq $numer t))) 88 (setq r ($asksign r)) 89 (cond ((equal r '$pos) (return 0)) 90 ((equal r '$neg) (return (simplify '$%pi))) 91 (t (merror (intl:gettext "plog: encountered atan(0/0)."))))) 92 ((and (among '%cos r) (among '%sin i)) 93 (setq var 'xz) 94 (numden (div* r i)) 95 (cond ((and (eq (caar nn*) '%cos) (eq (caar dn*) '%sin)) 96 (return (cadr nn*)))))) 97 (setq a ($sign r) b ($sign i)) 98 (cond ((eq a '$pos) (setq a 1)) 99 ((eq a '$neg) (setq a -1)) 100 ((eq a '$zero) (setq a 0))) 101 (cond ((eq b '$pos) (setq b 1)) 102 ((eq b '$neg) (setq b -1)) 103 ((eq a '$zero) (setq b 0))) 104 (cond ((equal i 0) 105 (return (if (equal a 1) 0 (simplify '$%pi)))) 106 ((equal r 0) 107 (return (cond ((equal b 1) (simplify half%pi)) 108 (t (mul2 '((rat simp) -1 2) 109 (simplify '$%pi))))))) 110 (setq r (simptimes (list '(mtimes) a b (div* i r)) 1 nil)) 111 (return (cond ((onep1 r) 112 (archk a b (list '(mtimes) '((rat) 1 4) '$%pi))) 113 ((alike1 r '((mexpt) 3 ((rat) 1 2))) 114 (archk a b (list '(mtimes) '((rat) 1 3) '$%pi))) 115 ((alike1 r '((mexpt) 3 ((rat) -1 2))) 116 (archk a b (list '(mtimes) '((rat) 1 6) '$%pi)))))))) 117 118;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 119 120;;; Implementation of the Binomial coefficient 121 122;; Verb function for the Binomial coefficient 123(defmfun $binomial (x y) 124 (simplify (list '(%binomial) x y))) 125 126;; Binomial has Mirror symmetry 127(defprop %binomial t commutes-with-conjugate) 128 129(defun simpbinocoef (x vestigial z) 130 (declare (ignore vestigial)) 131 (twoargcheck x) 132 (let ((u (simpcheck (cadr x) z)) 133 (v (simpcheck (caddr x) z)) 134 (y)) 135 (cond ((integerp v) 136 (cond ((minusp v) 137 (if (and (integerp u) (minusp u) (< v u)) 138 (bincomp u (- u v)) 139 0)) 140 ((or (zerop v) (equal u v)) 1) 141 ((and (integerp u) (not (minusp u))) 142 (bincomp u (min v (- u v)))) 143 (t (bincomp u v)))) 144 ((integerp (setq y (sub u v))) 145 (cond ((zerop1 y) 146 ;; u and v are equal, simplify not if argument can be negative 147 (if (member ($csign u) '($pnz $pn $neg $nz)) 148 (eqtest (list '(%binomial) u v) x) 149 (bincomp u y))) 150 (t (bincomp u y)))) 151 ((complex-float-numerical-eval-p u v) 152 ;; Numercial evaluation for real and complex floating point numbers. 153 (let (($numer t) ($float t)) 154 ($rectform 155 ($float 156 ($makegamma (list '(%binomial) ($float u) ($float v))))))) 157 ((complex-bigfloat-numerical-eval-p u v) 158 ;; Numerical evaluation for real and complex bigfloat numbers. 159 ($rectform 160 ($bfloat 161 ($makegamma (list '(%binomial) ($bfloat u) ($bfloat v)))))) 162 (t (eqtest (list '(%binomial) u v) x))))) 163 164(defun bincomp (u v) 165 (cond ((minusp v) 0) 166 ((zerop v) 1) 167 ((mnump u) (binocomp u v)) 168 (t (muln (bincomp1 u v) nil)))) 169 170(defun bincomp1 (u v) 171 (if (equal v 1) 172 (ncons u) 173 (list* u (list '(mexpt) v -1) (bincomp1 (add2 -1 u) (1- v))))) 174 175(defun binocomp (u v) 176 (prog (ans) 177 (setq ans 1) 178 loop (if (zerop v) (return ans)) 179 (setq ans (timesk (timesk u ans) (simplify (list '(rat) 1 v)))) 180 (setq u (addk -1 u) v (1- v)) 181 (go loop))) 182 183;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 184 185;;; Implementation of the Beta function 186 187(declare-top (special $numer $gammalim)) 188 189(defmvar $beta_args_sum_to_integer nil) 190 191;;; The Beta function has mirror symmetry 192(defprop $beta t commutes-with-conjugate) 193 194(defun simpbeta (x vestigial z &aux check) 195 (declare (ignore vestigial)) 196 (twoargcheck x) 197 (setq check x) 198 (let ((u (simpcheck (cadr x) z)) (v (simpcheck (caddr x) z))) 199 (cond ((or (zerop1 u) (zerop1 v)) 200 (if errorsw 201 (throw 'errorsw t) 202 (merror 203 (intl:gettext "beta: expected nonzero arguments; found ~M, ~M") 204 u v))) 205 206 ;; Check for numerical evaluation in float precision 207 ((complex-float-numerical-eval-p u v) 208 (cond 209 ;; We use gamma(u)*gamma(v)/gamma(u+v) for numerical evaluation. 210 ;; Therefore u, v or u+v can not be a negative integer or a 211 ;; floating point representation of a negative integer. 212 ((and (or (not (numberp u)) 213 (> u 0) 214 (not (= (nth-value 1 (truncate u)) 0))) 215 (and (or (not (numberp v)) 216 (> v 0) 217 (not (= (nth-value 1 (truncate v)) 0))) 218 (and (or (not (numberp (add u v))) 219 (> (add v u) 0) 220 (not (= (nth-value 1 ($truncate (add u v))) 0)))))) 221 ($rectform 222 (power ($float '$%e) 223 (add ($log_gamma ($float u)) 224 ($log_gamma ($float v)) 225 (mul -1 ($log_gamma ($float (add u v)))))))) 226 ((or (and (numberp u) 227 (> u 0) 228 (= (nth-value 1 (truncate u)) 0) 229 (not (and (mnump v) 230 (eq ($sign (sub ($truncate v) v)) '$zero) 231 (eq ($sign v) '$neg) 232 (eq ($sign (add u v)) '$pos))) 233 (setq u (truncate u))) 234 (and (numberp v) 235 (> v 0) 236 (= (nth-value 1 (truncate u)) 0) 237 (not (and (mnump u) 238 (eq ($sign (sub ($truncate u) u)) '$zero) 239 (eq ($sign u) '$neg) 240 (eq ($sign (add u v)) '$pos))) 241 (setq v (truncate v)))) 242 ;; One value is representing a negative integer, the other a 243 ;; positive integer and the sum is negative. Expand. 244 ($rectform ($float (beta-expand-integer u v)))) 245 (t 246 (eqtest (list '($beta) u v) check)))) 247 248 ;; Check for numerical evaluation in bigfloat precision 249 ((complex-bigfloat-numerical-eval-p u v) 250 (let (($ratprint nil)) 251 (cond 252 ((and (or (not (mnump u)) 253 (eq ($sign u) '$pos) 254 (not (eq ($sign (sub ($truncate u) u)) '$zero))) 255 (or (not (mnump v)) 256 (eq ($sign v) '$pos) 257 (not (eq ($sign (sub ($truncate v) v)) '$zero))) 258 (or (not (mnump (add u v))) 259 (eq ($sign (add u v)) '$pos) 260 (not (eq ($sign (sub ($truncate (add u v)) 261 (add u v))) 262 '$zero)))) 263 ($rectform 264 (power ($bfloat'$%e) 265 (add ($log_gamma ($bfloat u)) 266 ($log_gamma ($bfloat v)) 267 (mul -1 ($log_gamma ($bfloat (add u v)))))))) 268 ((or (and (mnump u) 269 (eq ($sign u) '$pos) 270 (eq ($sign (sub ($truncate u) u)) '$zero) 271 (not (and (mnump v) 272 (eq ($sign (sub ($truncate v) v)) '$zero) 273 (eq ($sign v) '$neg) 274 (eq ($sign (add u v)) '$pos))) 275 (setq u ($truncate u))) 276 (and (mnump v) 277 (eq ($sign v) '$pos) 278 (eq ($sign (sub ($truncate v) v)) '$zero) 279 (not (and (mnump u) 280 (eq ($sign (sub ($truncate u) u)) '$zero) 281 (eq ($sign u) '$neg) 282 (eq ($sign (add u v)) '$pos))) 283 (setq v ($truncate v)))) 284 ($rectform ($bfloat (beta-expand-integer u v)))) 285 (t 286 (eqtest (list '($beta) u v) check))))) 287 288 ((or (and (and (integerp u) 289 (plusp u)) 290 (not (and (mnump v) 291 (eq ($sign (sub ($truncate v) v)) '$zero) 292 (eq ($sign v) '$neg) 293 (eq ($sign (add u v)) '$pos)))) 294 (and (and (integerp v) 295 (plusp v)) 296 (not (and (mnump u) 297 (eq ($sign (sub ($truncate u) u)) '$zero) 298 (eq ($sign u) '$neg) 299 (eq ($sign (add u v)) '$pos))))) 300 ;; Expand for a positive integer. But not if the other argument is 301 ;; a negative integer and the sum of the integers is not negative. 302 (beta-expand-integer u v)) 303 304;;; At this point both integers are negative. This code does not work for 305;;; negative integers. The factorial function is not defined. 306; ((and (integerp u) (integerp v)) 307; (mul2* (div* (list '(mfactorial) (1- u)) 308; (list '(mfactorial) (+ u v -1))) 309; (list '(mfactorial) (1- v)))) 310 311 ((or (and (ratnump u) (ratnump v) (integerp (setq x (addk u v)))) 312 (and $beta_args_sum_to_integer 313 (integerp (setq x (expand1 (add2 u v) 1 1))))) 314 (let ((w (if (symbolp v) v u))) 315 (div* (mul2* '$%pi 316 (list '(%binomial) 317 (add2 (1- x) (neg w)) 318 (1- x))) 319 `((%sin) ((mtimes) ,w $%pi))))) 320 321 ((and $beta_expand (mplusp u) (integerp (cadr u))) 322 ;; Expand beta(a+n,b) where n is an integer. 323 (let ((n (cadr u)) 324 (u (simplify (cons '(mplus) (cddr u))))) 325 (beta-expand-add-integer n u v))) 326 327 ((and $beta_expand (mplusp v) (integerp (cadr v))) 328 ;; Expand beta(a,b+n) where n is an integer. 329 (let ((n (cadr v)) 330 (v (simplify (cons '(mplus) (cddr v))))) 331 (beta-expand-add-integer n v u))) 332 333 (t (eqtest (list '($beta) u v) check))))) 334 335(defun beta-expand-integer (u v) 336 ;; One of the arguments is a positive integer. Do an expansion. 337 ;; BUT for a negative integer as second argument the expansion is only 338 ;; possible when the sum of the integers is negative too. 339 ;; This routine expects that the calling routine has checked this. 340 (let ((x (add u v))) 341 (power 342 (mul (sub x 1) 343 (simplify 344 (list '(%binomial) 345 (sub x 2) 346 (sub (if (and (integerp u) (plusp u)) u v) 1)))) 347 -1))) 348 349(defun beta-expand-add-integer (n u v) 350 (if (plusp n) 351 (mul (simplify (list '($pochhammer) u n)) 352 (power (simplify (list '($pochhammer) (add u v) n)) -1) 353 (simplify (list '($beta) u v))) 354 (mul (simplify (list '($pochhammer) (add u v n) (- n))) 355 (power (simplify (list '($pochhammer) (add u n) (- n))) -1) 356 (simplify (list '($beta) u v))))) 357 358;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 359 360;;; Implementation of the Gamma function 361 362(defun simpgamma (x vestigial z) 363 (declare (ignore vestigial)) 364 (oneargcheck x) 365 (let ((j (simpcheck (cadr x) z))) 366 (cond ((and (floatp j) 367 (or (zerop j) 368 (and (< j 0) 369 (zerop (nth-value 1 (truncate j)))))) 370 (merror (intl:gettext "gamma: gamma(~:M) is undefined.") j)) 371 ((float-numerical-eval-p j) (gammafloat ($float j))) 372 ((and ($bfloatp j) 373 (or (zerop1 j) 374 (and (eq ($sign j) '$neg) 375 (zerop1 (sub j ($truncate j)))))) 376 (merror (intl:gettext "gamma: gamma(~:M) is undefined.") j)) 377 ((bigfloat-numerical-eval-p j) 378 ;; Adding 4 digits in the call to bffac. For $fpprec up to about 256 379 ;; and an argument up to about 500.0 the accuracy of the result is 380 ;; better than 10^(-$fpprec). 381 (let ((z (bigfloat:to ($bfloat j)))) 382 (cond 383 ((bigfloat:<= (bigfloat:abs z) (bigfloat:sqrt (bigfloat:epsilon z))) 384 ;; For small z, use gamma(z) = gamma(z+1)/z = z!/z 385 (div (mfuncall '$bffac 386 ($bfloat j) 387 (+ $fpprec 4)) 388 ($bfloat j))) 389 (t 390 (let ((result (mfuncall '$bffac (m+ ($bfloat j) -1) (+ $fpprec 4)))) 391 ;; bigfloatp will round the result to the correct fpprec 392 (bigfloatp result)))))) 393 ((complex-float-numerical-eval-p j) 394 (complexify (gamma-lanczos (complex ($float ($realpart j)) 395 ($float ($imagpart j)))))) 396 ((complex-bigfloat-numerical-eval-p j) 397 (let ((z (bigfloat:to ($bfloat j)))) 398 (cond 399 ((bigfloat:<= (bigfloat:abs z) 400 (bigfloat:sqrt (bigfloat:epsilon z))) 401 ;; For small z, use gamma(z) = gamma(z+1)/z = z!/z 402 (to (bigfloat:/ (bigfloat:to (mfuncall '$cbffac 403 (to z) 404 (+ $fpprec 4))) 405 z))) 406 (t 407 ;; Adding 4 digits in the call to cbffac. See comment above. 408 (let ((result 409 (mfuncall '$cbffac 410 (add -1 ($bfloat ($realpart j)) 411 (mul '$%i ($bfloat ($imagpart j)))) 412 (+ $fpprec 4)))) 413 (add (bigfloatp ($realpart result)) 414 (mul '$%i (bigfloatp ($imagpart result))))))))) 415 ((taylorize (mop x) (cadr x))) 416 ((eq j '$inf) '$inf) ; Simplify to $inf to be more consistent. 417 ((and $gamma_expand 418 (mplusp j) 419 (integerp (cadr j))) 420 ;; Expand gamma(z+n) for n an integer. 421 (let ((n (cadr j)) 422 (z (simplify (cons '(mplus) (cddr j))))) 423 (cond 424 ((> n 0) 425 (mul (simplify (list '($pochhammer) z n)) 426 (simplify (list '(%gamma) z)))) 427 ((< n 0) 428 (setq n (- n)) 429 (div (mul (power -1 n) (simplify (list '(%gamma) z))) 430 ;; We factor to get the order (z-1)*(z-2)*... 431 ;; and not (1-z)*(2-z)*... 432 ($factor 433 (simplify (list '($pochhammer) (sub 1 z) n)))))))) 434 ((integerp j) 435 (cond ((> j 0) 436 (cond ((<= j $factlim) 437 ;; Positive integer less than $factlim. Evaluate. 438 (simplify (list '(mfactorial) (1- j)))) 439 ;; Positive integer greater $factlim. Noun form. 440 (t (eqtest (list '(%gamma) j) x)))) 441 ;; Negative integer. Throw a Maxima error. 442 (errorsw (throw 'errorsw t)) 443 (t (merror (intl:gettext "gamma: gamma(~:M) is undefined.") j)))) 444 ((alike1 j '((rat) 1 2)) 445 (list '(mexpt simp) '$%pi j)) 446 ((and (mnump j) 447 (ratgreaterp $gammalim (simplify (list '(mabs) j))) 448 (or (ratgreaterp j 1) (ratgreaterp 0 j))) 449 ;; Expand for rational numbers less than $gammalim. 450 (gammared j)) 451 (t (eqtest (list '(%gamma) j) x))))) 452 453(defun gamma (y) ;;; numerical evaluation for 0 < y < 1 454 (prog (sum coefs) 455 (setq coefs (list 0.035868343 -0.193527817 0.48219939 456 -0.75670407 0.91820685 -0.89705693 457 0.98820588 -0.57719165)) 458 (unless (atom y) (setq y (fpcofrat y))) 459 (setq sum (car coefs) coefs (cdr coefs)) 460 loop (setq sum (+ (* sum y) (car coefs))) 461 (when (setq coefs (cdr coefs)) (go loop)) 462 (return (+ (/ y) sum)))) 463 464(defun gammared (a) ;A is assumed to 465 (prog (m q n) ;be '((RAT) M N) 466 (cond ((floatp a) (return (gammafloat a)))) 467 (setq m (cadr a) ;Numerator 468 n (caddr a) ;denominator 469 q (abs (truncate m n))) ;integer part 470 (cond ((minusp m) 471 (setq q (1+ q) m (+ m (* n q))) 472 (return 473 (simptimes (list '(mtimes) 474 (list '(mexpt) n q) 475 (simpgamma (list '(%gamma) 476 (list '(rat) m n)) 477 1 478 nil) 479 (list '(mexpt) (gammac m n q) -1)) 480 1 481 nil)))) 482 (return (m* (gammac m n q) 483 (simpgamma (list '(%gamma) 484 (list '(rat) (rem m n) n)) 485 1 nil) 486 (m^ n (- q)))))) 487 488(defun gammac (m n q) 489 (do ((ans 1)) 490 ((< q 1) ans) 491 (setq q (1- q) m (- m n) ans (* m ans)))) 492 493 494;; This implementation is based on Lanczos convergent formula for the 495;; gamma function for Re(z) > 0. We can use the reflection formula 496;; 497;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z) 498;; 499;; to handle the case of Re(z) <= 0. 500;; 501;; See http://my.fit.edu/~gabdo/gamma.m for some matlab code to 502;; compute this and http://my.fit.edu/~gabdo/gamma.txt for a nice 503;; discussion of Lanczos method and an improvement of Lanczos method. 504;; 505;; 506;; The document says this should give about 15 digits of accuracy for 507;; double-precision IEEE floats. The document also indicates how to 508;; compute a new set of coefficients if you need more range or 509;; accuracy. 510 511(defun gamma-lanczos (z) 512 (declare (type (complex flonum) z) 513 (optimize (safety 3))) 514 (let ((c (make-array 15 :element-type 'flonum 515 :initial-contents 516 '(0.99999999999999709182 517 57.156235665862923517 518 -59.597960355475491248 519 14.136097974741747174 520 -0.49191381609762019978 521 .33994649984811888699e-4 522 .46523628927048575665e-4 523 -.98374475304879564677e-4 524 .15808870322491248884e-3 525 -.21026444172410488319e-3 526 .21743961811521264320e-3 527 -.16431810653676389022e-3 528 .84418223983852743293e-4 529 -.26190838401581408670e-4 530 .36899182659531622704e-5)))) 531 (declare (type (simple-array flonum (15)) c)) 532 (cond 533 ((minusp (realpart z)) 534 ;; Use the reflection formula 535 ;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z) 536 ;; or 537 ;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z) 538 ;; 539 ;; If z is a negative integer, Gamma(z) is infinity. Should 540 ;; we test for this? Throw an error? 541 ;; The test must be done by the calling routine. 542 (/ (float pi) 543 (* (- z) (sin (* (float pi) z)) 544 (gamma-lanczos (- z))))) 545 ((<= (abs z) (sqrt flonum-epsilon)) 546 ;; For |z| small, use Gamma(z) = Gamma(z+1)/z 547 (/ (gamma-lanczos (+ 1 z)) 548 z)) 549 (t 550 (let* ((z (- z 1)) 551 (zh (+ z 1/2)) 552 (zgh (+ zh 607/128)) 553 (ss 554 (do ((sum 0.0) 555 (pp (1- (length c)) (1- pp))) 556 ((< pp 1) 557 sum) 558 (incf sum (/ (aref c pp) (+ z pp)))))) 559 (let ((result 560 ;; We check for an overflow. The last positive value in 561 ;; double-float precicsion for which Maxima can calculate 562 ;; gamma is ~171.6243 (CLISP 2.46 and GCL 2.6.8) 563 (ignore-errors 564 (let ((zp (expt zgh (/ zh 2)))) 565 (* (sqrt (float (* 2 pi))) 566 (+ ss (aref c 0)) 567 (* (/ zp (exp zgh)) zp)))))) 568 (cond ((null result) 569 ;; No result. Overflow. 570 (merror (intl:gettext "gamma: overflow in GAMMA-LANCZOS."))) 571 ((or (float-nan-p (realpart result)) 572 (float-inf-p (realpart result))) 573 ;; Result, but beyond extreme values. Overflow. 574 (merror (intl:gettext "gamma: overflow in GAMMA-LANCZOS."))) 575 (t result)))))))) 576 577(defun gammafloat (a) 578 (let ((a (float a))) 579 (cond ((minusp a) 580 ;; Reflection formula to make it positive: gamma(x) = 581 ;; %pi/sin(%pi*x)/x/gamma(-x) 582 (/ (float (- pi)) 583 (* a (sin (* (float pi) a))) 584 (gammafloat (- a)))) 585 ((<= a (sqrt flonum-epsilon)) 586 ;; Use gamma(x) = gamma(1+x)/x when x is very small 587 (/ (gammafloat (+ 1 a)) 588 a)) 589 ((< a 10) 590 (slatec::dgamma a)) 591 (t 592 (let ((result 593 (let ((c (* (sqrt (* 2 (float pi))) 594 (exp (slatec::d9lgmc a))))) 595 (let ((v (expt a (- (* .5e0 a) 0.25e0)))) 596 (* v 597 (/ v (exp a)) 598 c))))) 599 (if (or (float-nan-p result) 600 (float-inf-p result)) 601 (merror (intl:gettext "gamma: overflow in GAMMAFLOAT.")) 602 result)))))) 603 604(declare-top (special $numer $trigsign)) 605 606(defmfun $zeromatrix (m n) ($ematrix m n 0 1 1)) 607 608(defmfun $ematrix (m n var i j) 609 (prog (ans row) 610 (cond ((equal m 0) (return (ncons '($matrix simp)))) 611 ((and (equal n 0) (fixnump m) (> m 0)) 612 (return (cons '($matrix simp) (list-of-mlists m)))) 613 ((not (and (fixnump m) (fixnump n) 614 (fixnump i) (fixnump j) 615 (> m 0) (> n 0) (> i 0) (> j 0))) 616 (merror (intl:gettext "ematrix: arguments must be positive integers; found ~M") 617 (list '(mlist simp) m n i j) ))) 618 loop (cond ((= m i) (setq row (onen j n var 0)) (go on)) 619 ((zerop m) (return (cons '($matrix) (mxc ans))))) 620 (setq row nil) 621 (do ((n n (1- n))) ((zerop n)) (setq row (cons 0 row))) 622 on (setq ans (cons row ans) m (1- m)) 623 (go loop))) 624 625(defun list-of-mlists (n) 626 (do ((n n (1- n)) 627 (l nil (cons (ncons '(mlist simp)) l))) 628 ((= n 0) l))) 629 630(declare-top (special $ratmx)) 631 632(defmfun $coefmatrix (eql varl) (coefmatrix eql varl nil)) 633 634(defmfun $augcoefmatrix (eql varl) (coefmatrix eql varl t)) 635 636(defun coefmatrix (eql varl ind) 637 (prog (ans row a b elem) 638 (if (not ($listp eql)) (improper-arg-err eql '$coefmatrix)) 639 (if (not ($listp varl)) (improper-arg-err varl '$coefmatrix)) 640 (dolist (v (cdr varl)) 641 (if (and (not (atom v)) (member (caar v) '(mplus mtimes) :test #'eq)) 642 (merror (intl:gettext "coefmatrix: variables cannot be '+' or '*' expressions; found ~M") v))) 643 (setq eql (nreverse (mapcar #'meqhk (cdr eql))) 644 varl (reverse (cdr varl))) 645 loop1(if (null eql) (return (cons '($matrix) (mxc ans)))) 646 (setq a (car eql) eql (cdr eql) row nil) 647 (if ind (setq row (cons (const1 a varl) row))) 648 (setq b varl) 649 loop2(setq elem (ratcoef a (car b))) 650 (setq row (cons (if $ratmx elem (ratdisrep elem)) row)) 651 (if (setq b (cdr b)) (go loop2)) 652 (setq ans (cons row ans)) 653 (go loop1))) 654 655(defun const1 (e varl) 656 (dolist (v varl) (setq e (maxima-substitute 0 v e))) e) 657 658(defmfun $entermatrix (rows columns) 659 (prog (row column vector matrix sym symvector) 660 (cond ((or (not (fixnump rows)) 661 (not (fixnump columns))) 662 (merror (intl:gettext "entermatrix: arguments must be integers; found ~M, ~M") rows columns))) 663 (setq row 0) 664 (unless (= rows columns) (setq sym nil) (go oloop)) 665 quest (format t "~%Is the matrix 1. Diagonal 2. Symmetric 3. Antisymmetric 4. General~%") 666 (setq sym (retrieve "Answer 1, 2, 3 or 4 : " nil)) 667 (unless (member sym '(1 2 3 4)) (go quest)) 668 oloop (cond ((> (incf row) rows) 669 (format t "~%Matrix entered.~%") 670 (return (cons '($matrix) (mxc matrix))))) 671 (cond ((equal sym 1) 672 (setq column row) 673 (let ((prompt (format nil "Row ~a Column ~a: " row column))) 674 (setq matrix 675 (nconc matrix 676 (ncons (onen row columns 677 (meval (retrieve prompt nil)) 0))))) 678 (go oloop)) 679 ((equal sym 2) 680 (setq column (1- row)) 681 (cond ((equal row 1) (go iloop))) 682 (setq symvector 683 (cons (nthcdr column vector) symvector) 684 vector (nreverse (mapcar 'car symvector)) 685 symvector (mapcar 'cdr symvector)) 686 (go iloop)) 687 ((equal sym 3) 688 (setq column row) 689 (cond ((equal row 1) (setq vector (ncons 0)) (go iloop))) 690 (setq symvector 691 (cons (mapcar #'neg (nthcdr (1- column) vector)) 692 symvector) 693 vector (nreconc (mapcar 'car symvector) (ncons 0)) 694 symvector (mapcar 'cdr symvector)) 695 (go iloop))) 696 (setq column 0 vector nil) 697 iloop (cond ((> (incf column) columns) 698 (setq matrix (nconc matrix (ncons vector))) 699 (go oloop))) 700 (let ((prompt (format nil "Row ~a Column ~a: " row column))) 701 (setq vector (nconc vector (ncons (meval (retrieve prompt nil)))))) 702 (go iloop))) 703 704(declare-top (special sn* sd* rsn*)) 705 706(defmfun $xthru (e) 707 (cond ((atom e) e) 708 ((mtimesp e) (muln (mapcar '$xthru (cdr e)) nil)) 709 ((mplusp e) (simplify (comdenom (mapcar '$xthru (cdr e)) t))) 710 ((mexptp e) (power ($xthru (cadr e)) (caddr e))) 711 ((mbagp e) (cons (car e) (mapcar '$xthru (cdr e)))) 712 (t e))) 713 714(defun comdenom (l ind) 715 (prog (n d) 716 (prodnumden (car l)) 717 (setq n (m*l sn*) sn* nil) 718 (setq d (m*l sd*) sd* nil) 719 loop (setq l (cdr l)) 720 (cond ((null l) 721 (return (cond (ind (div* (cond (rsn* ($ratsimp n)) 722 (t n)) 723 d)) 724 (t (list n d)))))) 725 (prodnumden (car l)) 726 (setq d (comdenom1 n d (m*l sn*) (m*l sd*))) 727 (setq n (car d)) 728 (setq d (cadr d)) 729 (go loop))) 730 731(defun prodnumden (e) 732 (cond ((atom e) (prodnd (list e))) 733 ((eq (caar e) 'mtimes) (prodnd (cdr e))) 734 (t (prodnd (list e))))) 735 736(defun prodnd (l) 737 (prog (e) 738 (setq l (reverse l)) 739 (setq sn* nil sd* nil) 740 loop (cond ((null l) (return nil))) 741 (setq e (car l)) 742 (cond ((atom e) (setq sn* (cons e sn*))) 743 ((ratnump e) 744 (cond ((not (equal 1 (cadr e))) 745 (setq sn* (cons (cadr e) sn*)))) 746 (setq sd* (cons (caddr e) sd*))) 747 ((and (eq (caar e) 'mexpt) 748 (mnegp (caddr e))) 749 (setq sd* (cons (power (cadr e) 750 (timesk -1 (caddr e))) 751 sd*))) 752 (t (setq sn* (cons e sn*)))) 753 (setq l (cdr l)) 754 (go loop))) 755 756(defun comdenom1 (a b c d) 757 (prog (b1 c1) 758 (prodnumden (div* b d)) 759 (setq b1 (m*l sn*) sn* nil) 760 (setq c1 (m*l sd*) sd* nil) 761 (return 762 (list (add2 (m* a c1) (m* c b1)) 763 (mul2 d b1))))) 764 765(declare-top (special $globalsolve $backsubst $dispflag 766 $linsolve_params $%rnum_list ax *linelabel* $linechar 767 $linenum *mosesflag)) 768 769(defun xrutout (ax n m varl ind) 770 (let (($linsolve_params (and $backsubst $linsolve_params))) 771 (prog (ix imin ans zz m-1 sol tim chk zzz) 772 (setq ax (get-array-pointer ax) tim 0) 773 (if $linsolve_params (setq $%rnum_list (list '(mlist)))) 774 (setq imin (min (setq m-1 (1- m)) n)) 775 (setq ix (max imin (length varl))) 776 loop (if (zerop ix) (if ind (go out) (return (cons '(mlist) zz)))) 777 (when (or (> ix imin) (equal (car (aref ax ix ix)) 0)) 778 (setf (aref ax 0 ix) 779 (rform (if $linsolve_params (make-param) (ith varl ix)))) 780 (if $linsolve_params (go saval) (go next))) 781 (setq ans (aref ax ix m)) 782 (setf (aref ax ix m) nil) 783 (do ((j (1+ ix) (1+ j))) 784 ((> j m-1)) 785 (setq ans (ratdif ans (rattimes (aref ax ix j) (aref ax 0 j) t))) 786 (setf (aref ax ix j) nil)) 787 (setf (aref ax 0 ix) (ratquotient ans (aref ax ix ix))) 788 (setf (aref ax ix ix) nil) 789 (setq ans nil) 790 saval (push (cond (*mosesflag (aref ax 0 ix)) 791 (t (list (if $globalsolve '(msetq) '(mequal)) 792 (ith varl ix) 793 (simplify (rdis (aref ax 0 ix)))))) 794 zz) 795 (if (not $backsubst) 796 (setf (aref ax 0 ix) (rform (ith varl ix)))) 797 (and $globalsolve (meval (car zz))) 798 next (decf ix) 799 (go loop) 800 out 801 (cond ($dispflag (mtell (intl:gettext "Solution:~%")))) 802 (setq sol (list '(mlist)) chk (checklabel $linechar)) 803 (do ((ll zz (cdr ll))) 804 ((null ll)) 805 (setq zzz (car ll)) 806 (setq zzz (list '(mlabel) 807 (progn 808 (if chk 809 (setq chk nil) 810 (incf $linenum)) 811 (let (($nolabels nil)) 812 (makelabel $linechar)) 813 *linelabel*) 814 (setf (symbol-value *linelabel*) zzz))) 815 (nconc sol (ncons *linelabel*)) 816 (cond ($dispflag 817 (setq tim (get-internal-run-time)) 818 (mtell-open "~%~M" zzz) 819 (timeorg tim)) 820 (t 821 (putprop *linelabel* t 'nodisp)))) 822 (return sol)))) 823