1;;; calcalg2.el --- more algebraic functions for Calc 2 3;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004, 4;; 2005, 2006, 2007 Free Software Foundation, Inc. 5 6;; Author: David Gillespie <daveg@synaptics.com> 7;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com> 8 9;; This file is part of GNU Emacs. 10 11;; GNU Emacs is free software; you can redistribute it and/or modify 12;; it under the terms of the GNU General Public License as published by 13;; the Free Software Foundation; either version 2, or (at your option) 14;; any later version. 15 16;; GNU Emacs is distributed in the hope that it will be useful, 17;; but WITHOUT ANY WARRANTY; without even the implied warranty of 18;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19;; GNU General Public License for more details. 20 21;; You should have received a copy of the GNU General Public License 22;; along with GNU Emacs; see the file COPYING. If not, write to the 23;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 24;; Boston, MA 02110-1301, USA. 25 26;;; Commentary: 27 28;;; Code: 29 30;; This file is autoloaded from calc-ext.el. 31 32(require 'calc-ext) 33(require 'calc-macs) 34 35(defun calc-derivative (var num) 36 (interactive "sDifferentiate with respect to: \np") 37 (calc-slow-wrapper 38 (when (< num 0) 39 (error "Order of derivative must be positive")) 40 (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv)) 41 n expr) 42 (if (or (equal var "") (equal var "$")) 43 (setq n 2 44 expr (calc-top-n 2) 45 var (calc-top-n 1)) 46 (setq var (math-read-expr var)) 47 (when (eq (car-safe var) 'error) 48 (error "Bad format in expression: %s" (nth 1 var))) 49 (setq n 1 50 expr (calc-top-n 1))) 51 (while (>= (setq num (1- num)) 0) 52 (setq expr (list func expr var))) 53 (calc-enter-result n "derv" expr)))) 54 55(defun calc-integral (var &optional arg) 56 (interactive "sIntegration variable: \nP") 57 (if arg 58 (calc-tabular-command 'calcFunc-integ "Integration" "intg" nil var nil nil) 59 (calc-slow-wrapper 60 (if (or (equal var "") (equal var "$")) 61 (calc-enter-result 2 "intg" (list 'calcFunc-integ 62 (calc-top-n 2) 63 (calc-top-n 1))) 64 (let ((var (math-read-expr var))) 65 (if (eq (car-safe var) 'error) 66 (error "Bad format in expression: %s" (nth 1 var))) 67 (calc-enter-result 1 "intg" (list 'calcFunc-integ 68 (calc-top-n 1) 69 var))))))) 70 71(defun calc-num-integral (&optional varname lowname highname) 72 (interactive "sIntegration variable: ") 73 (calc-tabular-command 'calcFunc-ninteg "Integration" "nint" 74 nil varname lowname highname)) 75 76(defun calc-summation (arg &optional varname lowname highname) 77 (interactive "P\nsSummation variable: ") 78 (calc-tabular-command 'calcFunc-sum "Summation" "sum" 79 arg varname lowname highname)) 80 81(defun calc-alt-summation (arg &optional varname lowname highname) 82 (interactive "P\nsSummation variable: ") 83 (calc-tabular-command 'calcFunc-asum "Summation" "asum" 84 arg varname lowname highname)) 85 86(defun calc-product (arg &optional varname lowname highname) 87 (interactive "P\nsIndex variable: ") 88 (calc-tabular-command 'calcFunc-prod "Index" "prod" 89 arg varname lowname highname)) 90 91(defun calc-tabulate (arg &optional varname lowname highname) 92 (interactive "P\nsIndex variable: ") 93 (calc-tabular-command 'calcFunc-table "Index" "tabl" 94 arg varname lowname highname)) 95 96(defun calc-tabular-command (func prompt prefix arg varname lowname highname) 97 (calc-slow-wrapper 98 (let (var (low nil) (high nil) (step nil) stepname stepnum (num 1) expr) 99 (if (consp arg) 100 (setq stepnum 1) 101 (setq stepnum 0)) 102 (if (or (equal varname "") (equal varname "$") (null varname)) 103 (setq high (calc-top-n (+ stepnum 1)) 104 low (calc-top-n (+ stepnum 2)) 105 var (calc-top-n (+ stepnum 3)) 106 num (+ stepnum 4)) 107 (setq var (if (stringp varname) (math-read-expr varname) varname)) 108 (if (eq (car-safe var) 'error) 109 (error "Bad format in expression: %s" (nth 1 var))) 110 (or lowname 111 (setq lowname (read-string (concat prompt " variable: " varname 112 ", from: ")))) 113 (if (or (equal lowname "") (equal lowname "$")) 114 (setq high (calc-top-n (+ stepnum 1)) 115 low (calc-top-n (+ stepnum 2)) 116 num (+ stepnum 3)) 117 (setq low (if (stringp lowname) (math-read-expr lowname) lowname)) 118 (if (eq (car-safe low) 'error) 119 (error "Bad format in expression: %s" (nth 1 low))) 120 (or highname 121 (setq highname (read-string (concat prompt " variable: " varname 122 ", from: " lowname 123 ", to: ")))) 124 (if (or (equal highname "") (equal highname "$")) 125 (setq high (calc-top-n (+ stepnum 1)) 126 num (+ stepnum 2)) 127 (setq high (if (stringp highname) (math-read-expr highname) 128 highname)) 129 (if (eq (car-safe high) 'error) 130 (error "Bad format in expression: %s" (nth 1 high))) 131 (if (consp arg) 132 (progn 133 (setq stepname (read-string (concat prompt " variable: " 134 varname 135 ", from: " lowname 136 ", to: " highname 137 ", step: "))) 138 (if (or (equal stepname "") (equal stepname "$")) 139 (setq step (calc-top-n 1) 140 num 2) 141 (setq step (math-read-expr stepname)) 142 (if (eq (car-safe step) 'error) 143 (error "Bad format in expression: %s" 144 (nth 1 step))))))))) 145 (or step 146 (if (consp arg) 147 (setq step (calc-top-n 1)) 148 (if arg 149 (setq step (prefix-numeric-value arg))))) 150 (setq expr (calc-top-n num)) 151 (calc-enter-result num prefix (append (list func expr var low high) 152 (and step (list step))))))) 153 154(defun calc-solve-for (var) 155 (interactive "sVariable(s) to solve for: ") 156 (calc-slow-wrapper 157 (let ((func (if (calc-is-inverse) 158 (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv) 159 (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve)))) 160 (if (or (equal var "") (equal var "$")) 161 (calc-enter-result 2 "solv" (list func 162 (calc-top-n 2) 163 (calc-top-n 1))) 164 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var) 165 (not (string-match "\\[" var))) 166 (math-read-expr (concat "[" var "]")) 167 (math-read-expr var)))) 168 (if (eq (car-safe var) 'error) 169 (error "Bad format in expression: %s" (nth 1 var))) 170 (calc-enter-result 1 "solv" (list func 171 (calc-top-n 1) 172 var))))))) 173 174(defun calc-poly-roots (var) 175 (interactive "sVariable to solve for: ") 176 (calc-slow-wrapper 177 (if (or (equal var "") (equal var "$")) 178 (calc-enter-result 2 "prts" (list 'calcFunc-roots 179 (calc-top-n 2) 180 (calc-top-n 1))) 181 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var) 182 (not (string-match "\\[" var))) 183 (math-read-expr (concat "[" var "]")) 184 (math-read-expr var)))) 185 (if (eq (car-safe var) 'error) 186 (error "Bad format in expression: %s" (nth 1 var))) 187 (calc-enter-result 1 "prts" (list 'calcFunc-roots 188 (calc-top-n 1) 189 var)))))) 190 191(defun calc-taylor (var nterms) 192 (interactive "sTaylor expansion variable: \nNNumber of terms: ") 193 (calc-slow-wrapper 194 (let ((var (math-read-expr var))) 195 (if (eq (car-safe var) 'error) 196 (error "Bad format in expression: %s" (nth 1 var))) 197 (calc-enter-result 1 "tylr" (list 'calcFunc-taylor 198 (calc-top-n 1) 199 var 200 (prefix-numeric-value nterms)))))) 201 202 203;; The following are global variables used by math-derivative and some 204;; related functions 205(defvar math-deriv-var) 206(defvar math-deriv-total) 207(defvar math-deriv-symb) 208(defvar math-decls-cache) 209(defvar math-decls-all) 210 211(defun math-derivative (expr) 212 (cond ((equal expr math-deriv-var) 213 1) 214 ((or (Math-scalarp expr) 215 (eq (car expr) 'sdev) 216 (and (eq (car expr) 'var) 217 (or (not math-deriv-total) 218 (math-const-var expr) 219 (progn 220 (math-setup-declarations) 221 (memq 'const (nth 1 (or (assq (nth 2 expr) 222 math-decls-cache) 223 math-decls-all))))))) 224 0) 225 ((eq (car expr) '+) 226 (math-add (math-derivative (nth 1 expr)) 227 (math-derivative (nth 2 expr)))) 228 ((eq (car expr) '-) 229 (math-sub (math-derivative (nth 1 expr)) 230 (math-derivative (nth 2 expr)))) 231 ((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt 232 calcFunc-gt calcFunc-leq calcFunc-geq)) 233 (list (car expr) 234 (math-derivative (nth 1 expr)) 235 (math-derivative (nth 2 expr)))) 236 ((eq (car expr) 'neg) 237 (math-neg (math-derivative (nth 1 expr)))) 238 ((eq (car expr) '*) 239 (math-add (math-mul (nth 2 expr) 240 (math-derivative (nth 1 expr))) 241 (math-mul (nth 1 expr) 242 (math-derivative (nth 2 expr))))) 243 ((eq (car expr) '/) 244 (math-sub (math-div (math-derivative (nth 1 expr)) 245 (nth 2 expr)) 246 (math-div (math-mul (nth 1 expr) 247 (math-derivative (nth 2 expr))) 248 (math-sqr (nth 2 expr))))) 249 ((eq (car expr) '^) 250 (let ((du (math-derivative (nth 1 expr))) 251 (dv (math-derivative (nth 2 expr)))) 252 (or (Math-zerop du) 253 (setq du (math-mul (nth 2 expr) 254 (math-mul (math-normalize 255 (list '^ 256 (nth 1 expr) 257 (math-add (nth 2 expr) -1))) 258 du)))) 259 (or (Math-zerop dv) 260 (setq dv (math-mul (math-normalize 261 (list 'calcFunc-ln (nth 1 expr))) 262 (math-mul expr dv)))) 263 (math-add du dv))) 264 ((eq (car expr) '%) 265 (math-derivative (nth 1 expr))) ; a reasonable definition 266 ((eq (car expr) 'vec) 267 (math-map-vec 'math-derivative expr)) 268 ((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im)) 269 (= (length expr) 2)) 270 (list (car expr) (math-derivative (nth 1 expr)))) 271 ((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol)) 272 (= (length expr) 3)) 273 (let ((d (math-derivative (nth 1 expr)))) 274 (if (math-numberp d) 275 0 ; assume x and x_1 are independent vars 276 (list (car expr) d (nth 2 expr))))) 277 (t (or (and (symbolp (car expr)) 278 (if (= (length expr) 2) 279 (let ((handler (get (car expr) 'math-derivative))) 280 (and handler 281 (let ((deriv (math-derivative (nth 1 expr)))) 282 (if (Math-zerop deriv) 283 deriv 284 (math-mul (funcall handler (nth 1 expr)) 285 deriv))))) 286 (let ((handler (get (car expr) 'math-derivative-n))) 287 (and handler 288 (funcall handler expr))))) 289 (and (not (eq math-deriv-symb 'pre-expand)) 290 (let ((exp (math-expand-formula expr))) 291 (and exp 292 (or (let ((math-deriv-symb 'pre-expand)) 293 (catch 'math-deriv (math-derivative expr))) 294 (math-derivative exp))))) 295 (if (or (Math-objvecp expr) 296 (eq (car expr) 'var) 297 (not (symbolp (car expr)))) 298 (if math-deriv-symb 299 (throw 'math-deriv nil) 300 (list (if math-deriv-total 'calcFunc-tderiv 'calcFunc-deriv) 301 expr 302 math-deriv-var)) 303 (let ((accum 0) 304 (arg expr) 305 (n 1) 306 derv) 307 (while (setq arg (cdr arg)) 308 (or (Math-zerop (setq derv (math-derivative (car arg)))) 309 (let ((func (intern (concat (symbol-name (car expr)) 310 "'" 311 (if (> n 1) 312 (int-to-string n) 313 "")))) 314 (prop (cond ((= (length expr) 2) 315 'math-derivative-1) 316 ((= (length expr) 3) 317 'math-derivative-2) 318 ((= (length expr) 4) 319 'math-derivative-3) 320 ((= (length expr) 5) 321 'math-derivative-4) 322 ((= (length expr) 6) 323 'math-derivative-5)))) 324 (setq accum 325 (math-add 326 accum 327 (math-mul 328 derv 329 (let ((handler (get func prop))) 330 (or (and prop handler 331 (apply handler (cdr expr))) 332 (if (and math-deriv-symb 333 (not (get func 334 'calc-user-defn))) 335 (throw 'math-deriv nil) 336 (cons func (cdr expr)))))))))) 337 (setq n (1+ n))) 338 accum)))))) 339 340(defun calcFunc-deriv (expr math-deriv-var &optional deriv-value math-deriv-symb) 341 (let* ((math-deriv-total nil) 342 (res (catch 'math-deriv (math-derivative expr)))) 343 (or (eq (car-safe res) 'calcFunc-deriv) 344 (null res) 345 (setq res (math-normalize res))) 346 (and res 347 (if deriv-value 348 (math-expr-subst res math-deriv-var deriv-value) 349 res)))) 350 351(defun calcFunc-tderiv (expr math-deriv-var &optional deriv-value math-deriv-symb) 352 (math-setup-declarations) 353 (let* ((math-deriv-total t) 354 (res (catch 'math-deriv (math-derivative expr)))) 355 (or (eq (car-safe res) 'calcFunc-tderiv) 356 (null res) 357 (setq res (math-normalize res))) 358 (and res 359 (if deriv-value 360 (math-expr-subst res math-deriv-var deriv-value) 361 res)))) 362 363(put 'calcFunc-inv\' 'math-derivative-1 364 (function (lambda (u) (math-neg (math-div 1 (math-sqr u)))))) 365 366(put 'calcFunc-sqrt\' 'math-derivative-1 367 (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u)))))) 368 369(put 'calcFunc-deg\' 'math-derivative-1 370 (function (lambda (u) (math-div-float '(float 18 1) (math-pi))))) 371 372(put 'calcFunc-rad\' 'math-derivative-1 373 (function (lambda (u) (math-pi-over-180)))) 374 375(put 'calcFunc-ln\' 'math-derivative-1 376 (function (lambda (u) (math-div 1 u)))) 377 378(put 'calcFunc-log10\' 'math-derivative-1 379 (function (lambda (u) 380 (math-div (math-div 1 (math-normalize '(calcFunc-ln 10))) 381 u)))) 382 383(put 'calcFunc-lnp1\' 'math-derivative-1 384 (function (lambda (u) (math-div 1 (math-add u 1))))) 385 386(put 'calcFunc-log\' 'math-derivative-2 387 (function (lambda (x b) 388 (and (not (Math-zerop b)) 389 (let ((lnv (math-normalize 390 (list 'calcFunc-ln b)))) 391 (math-div 1 (math-mul lnv x))))))) 392 393(put 'calcFunc-log\'2 'math-derivative-2 394 (function (lambda (x b) 395 (let ((lnv (list 'calcFunc-ln b))) 396 (math-neg (math-div (list 'calcFunc-log x b) 397 (math-mul lnv b))))))) 398 399(put 'calcFunc-exp\' 'math-derivative-1 400 (function (lambda (u) (math-normalize (list 'calcFunc-exp u))))) 401 402(put 'calcFunc-expm1\' 'math-derivative-1 403 (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u))))) 404 405(put 'calcFunc-sin\' 'math-derivative-1 406 (function (lambda (u) (math-to-radians-2 (math-normalize 407 (list 'calcFunc-cos u)))))) 408 409(put 'calcFunc-cos\' 'math-derivative-1 410 (function (lambda (u) (math-neg (math-to-radians-2 411 (math-normalize 412 (list 'calcFunc-sin u))))))) 413 414(put 'calcFunc-tan\' 'math-derivative-1 415 (function (lambda (u) (math-to-radians-2 416 (math-sqr 417 (math-normalize 418 (list 'calcFunc-sec u))))))) 419 420(put 'calcFunc-sec\' 'math-derivative-1 421 (function (lambda (u) (math-to-radians-2 422 (math-mul 423 (math-normalize 424 (list 'calcFunc-sec u)) 425 (math-normalize 426 (list 'calcFunc-tan u))))))) 427 428(put 'calcFunc-csc\' 'math-derivative-1 429 (function (lambda (u) (math-neg 430 (math-to-radians-2 431 (math-mul 432 (math-normalize 433 (list 'calcFunc-csc u)) 434 (math-normalize 435 (list 'calcFunc-cot u)))))))) 436 437(put 'calcFunc-cot\' 'math-derivative-1 438 (function (lambda (u) (math-neg 439 (math-to-radians-2 440 (math-sqr 441 (math-normalize 442 (list 'calcFunc-csc u)))))))) 443 444(put 'calcFunc-arcsin\' 'math-derivative-1 445 (function (lambda (u) 446 (math-from-radians-2 447 (math-div 1 (math-normalize 448 (list 'calcFunc-sqrt 449 (math-sub 1 (math-sqr u))))))))) 450 451(put 'calcFunc-arccos\' 'math-derivative-1 452 (function (lambda (u) 453 (math-from-radians-2 454 (math-div -1 (math-normalize 455 (list 'calcFunc-sqrt 456 (math-sub 1 (math-sqr u))))))))) 457 458(put 'calcFunc-arctan\' 'math-derivative-1 459 (function (lambda (u) (math-from-radians-2 460 (math-div 1 (math-add 1 (math-sqr u))))))) 461 462(put 'calcFunc-sinh\' 'math-derivative-1 463 (function (lambda (u) (math-normalize (list 'calcFunc-cosh u))))) 464 465(put 'calcFunc-cosh\' 'math-derivative-1 466 (function (lambda (u) (math-normalize (list 'calcFunc-sinh u))))) 467 468(put 'calcFunc-tanh\' 'math-derivative-1 469 (function (lambda (u) (math-sqr 470 (math-normalize 471 (list 'calcFunc-sech u)))))) 472 473(put 'calcFunc-sech\' 'math-derivative-1 474 (function (lambda (u) (math-neg 475 (math-mul 476 (math-normalize (list 'calcFunc-sech u)) 477 (math-normalize (list 'calcFunc-tanh u))))))) 478 479(put 'calcFunc-csch\' 'math-derivative-1 480 (function (lambda (u) (math-neg 481 (math-mul 482 (math-normalize (list 'calcFunc-csch u)) 483 (math-normalize (list 'calcFunc-coth u))))))) 484 485(put 'calcFunc-coth\' 'math-derivative-1 486 (function (lambda (u) (math-neg 487 (math-sqr 488 (math-normalize 489 (list 'calcFunc-csch u))))))) 490 491(put 'calcFunc-arcsinh\' 'math-derivative-1 492 (function (lambda (u) 493 (math-div 1 (math-normalize 494 (list 'calcFunc-sqrt 495 (math-add (math-sqr u) 1))))))) 496 497(put 'calcFunc-arccosh\' 'math-derivative-1 498 (function (lambda (u) 499 (math-div 1 (math-normalize 500 (list 'calcFunc-sqrt 501 (math-add (math-sqr u) -1))))))) 502 503(put 'calcFunc-arctanh\' 'math-derivative-1 504 (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u)))))) 505 506(put 'calcFunc-bern\'2 'math-derivative-2 507 (function (lambda (n x) 508 (math-mul n (list 'calcFunc-bern (math-add n -1) x))))) 509 510(put 'calcFunc-euler\'2 'math-derivative-2 511 (function (lambda (n x) 512 (math-mul n (list 'calcFunc-euler (math-add n -1) x))))) 513 514(put 'calcFunc-gammag\'2 'math-derivative-2 515 (function (lambda (a x) (math-deriv-gamma a x 1)))) 516 517(put 'calcFunc-gammaG\'2 'math-derivative-2 518 (function (lambda (a x) (math-deriv-gamma a x -1)))) 519 520(put 'calcFunc-gammaP\'2 'math-derivative-2 521 (function (lambda (a x) (math-deriv-gamma a x 522 (math-div 523 1 (math-normalize 524 (list 'calcFunc-gamma 525 a))))))) 526 527(put 'calcFunc-gammaQ\'2 'math-derivative-2 528 (function (lambda (a x) (math-deriv-gamma a x 529 (math-div 530 -1 (math-normalize 531 (list 'calcFunc-gamma 532 a))))))) 533 534(defun math-deriv-gamma (a x scale) 535 (math-mul scale 536 (math-mul (math-pow x (math-add a -1)) 537 (list 'calcFunc-exp (math-neg x))))) 538 539(put 'calcFunc-betaB\' 'math-derivative-3 540 (function (lambda (x a b) (math-deriv-beta x a b 1)))) 541 542(put 'calcFunc-betaI\' 'math-derivative-3 543 (function (lambda (x a b) (math-deriv-beta x a b 544 (math-div 545 1 (list 'calcFunc-beta 546 a b)))))) 547 548(defun math-deriv-beta (x a b scale) 549 (math-mul (math-mul (math-pow x (math-add a -1)) 550 (math-pow (math-sub 1 x) (math-add b -1))) 551 scale)) 552 553(put 'calcFunc-erf\' 'math-derivative-1 554 (function (lambda (x) (math-div 2 555 (math-mul (list 'calcFunc-exp 556 (math-sqr x)) 557 (if calc-symbolic-mode 558 '(calcFunc-sqrt 559 (var pi var-pi)) 560 (math-sqrt-pi))))))) 561 562(put 'calcFunc-erfc\' 'math-derivative-1 563 (function (lambda (x) (math-div -2 564 (math-mul (list 'calcFunc-exp 565 (math-sqr x)) 566 (if calc-symbolic-mode 567 '(calcFunc-sqrt 568 (var pi var-pi)) 569 (math-sqrt-pi))))))) 570 571(put 'calcFunc-besJ\'2 'math-derivative-2 572 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ 573 (math-add v -1) 574 z) 575 (list 'calcFunc-besJ 576 (math-add v 1) 577 z)) 578 2)))) 579 580(put 'calcFunc-besY\'2 'math-derivative-2 581 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY 582 (math-add v -1) 583 z) 584 (list 'calcFunc-besY 585 (math-add v 1) 586 z)) 587 2)))) 588 589(put 'calcFunc-sum 'math-derivative-n 590 (function 591 (lambda (expr) 592 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var) 593 (throw 'math-deriv nil) 594 (cons 'calcFunc-sum 595 (cons (math-derivative (nth 1 expr)) 596 (cdr (cdr expr)))))))) 597 598(put 'calcFunc-prod 'math-derivative-n 599 (function 600 (lambda (expr) 601 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var) 602 (throw 'math-deriv nil) 603 (math-mul expr 604 (cons 'calcFunc-sum 605 (cons (math-div (math-derivative (nth 1 expr)) 606 (nth 1 expr)) 607 (cdr (cdr expr))))))))) 608 609(put 'calcFunc-integ 'math-derivative-n 610 (function 611 (lambda (expr) 612 (if (= (length expr) 3) 613 (if (equal (nth 2 expr) math-deriv-var) 614 (nth 1 expr) 615 (math-normalize 616 (list 'calcFunc-integ 617 (math-derivative (nth 1 expr)) 618 (nth 2 expr)))) 619 (if (= (length expr) 5) 620 (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr) 621 (nth 3 expr))) 622 (upper (math-expr-subst (nth 1 expr) (nth 2 expr) 623 (nth 4 expr)))) 624 (math-add (math-sub (math-mul upper 625 (math-derivative (nth 4 expr))) 626 (math-mul lower 627 (math-derivative (nth 3 expr)))) 628 (if (equal (nth 2 expr) math-deriv-var) 629 0 630 (math-normalize 631 (list 'calcFunc-integ 632 (math-derivative (nth 1 expr)) (nth 2 expr) 633 (nth 3 expr) (nth 4 expr))))))))))) 634 635(put 'calcFunc-if 'math-derivative-n 636 (function 637 (lambda (expr) 638 (and (= (length expr) 4) 639 (list 'calcFunc-if (nth 1 expr) 640 (math-derivative (nth 2 expr)) 641 (math-derivative (nth 3 expr))))))) 642 643(put 'calcFunc-subscr 'math-derivative-n 644 (function 645 (lambda (expr) 646 (and (= (length expr) 3) 647 (list 'calcFunc-subscr (nth 1 expr) 648 (math-derivative (nth 2 expr))))))) 649 650 651(defvar math-integ-var '(var X ---)) 652(defvar math-integ-var-2 '(var Y ---)) 653(defvar math-integ-vars (list 'f math-integ-var math-integ-var-2)) 654(defvar math-integ-var-list (list math-integ-var)) 655(defvar math-integ-var-list-list (list math-integ-var-list)) 656 657;; math-integ-depth is a local variable for math-try-integral, but is used 658;; by math-integral and math-tracing-integral 659;; which are called (directly or indirectly) by math-try-integral. 660(defvar math-integ-depth) 661;; math-integ-level is a local variable for math-try-integral, but is used 662;; by math-integral, math-do-integral, math-tracing-integral, 663;; math-sub-integration, math-integrate-by-parts and 664;; math-integrate-by-substitution, which are called (directly or 665;; indirectly) by math-try-integral. 666(defvar math-integ-level) 667;; math-integral-limit is a local variable for calcFunc-integ, but is 668;; used by math-tracing-integral, math-sub-integration and 669;; math-try-integration. 670(defvar math-integral-limit) 671 672(defmacro math-tracing-integral (&rest parts) 673 (list 'and 674 'trace-buffer 675 (list 'save-excursion 676 '(set-buffer trace-buffer) 677 '(goto-char (point-max)) 678 (list 'and 679 '(bolp) 680 '(insert (make-string (- math-integral-limit 681 math-integ-level) 32) 682 (format "%2d " math-integ-depth) 683 (make-string math-integ-level 32))) 684 ;;(list 'condition-case 'err 685 (cons 'insert parts) 686 ;; '(error (insert (prin1-to-string err)))) 687 '(sit-for 0)))) 688 689;;; The following wrapper caches results and avoids infinite recursion. 690;;; Each cache entry is: ( A B ) Integral of A is B; 691;;; ( A N ) Integral of A failed at level N; 692;;; ( A busy ) Currently working on integral of A; 693;;; ( A parts ) Currently working, integ-by-parts; 694;;; ( A parts2 ) Currently working, integ-by-parts; 695;;; ( A cancelled ) Ignore this cache entry; 696;;; ( A [B] ) Same result as for math-cur-record = B. 697 698;; math-cur-record is a local variable for math-try-integral, but is used 699;; by math-integral, math-replace-integral-parts and math-integrate-by-parts 700;; which are called (directly or indirectly) by math-try-integral, as well as 701;; by calc-dump-integral-cache 702(defvar math-cur-record) 703;; math-enable-subst and math-any-substs are local variables for 704;; calcFunc-integ, but are used by math-integral and math-try-integral. 705(defvar math-enable-subst) 706(defvar math-any-substs) 707 708;; math-integ-msg is a local variable for math-try-integral, but is 709;; used (both locally and non-locally) by math-integral. 710(defvar math-integ-msg) 711 712(defvar math-integral-cache nil) 713(defvar math-integral-cache-state nil) 714 715(defun math-integral (expr &optional simplify same-as-above) 716 (let* ((simp math-cur-record) 717 (math-cur-record (assoc expr math-integral-cache)) 718 (math-integ-depth (1+ math-integ-depth)) 719 (val 'cancelled)) 720 (math-tracing-integral "Integrating " 721 (math-format-value expr 1000) 722 "...\n") 723 (and math-cur-record 724 (progn 725 (math-tracing-integral "Found " 726 (math-format-value (nth 1 math-cur-record) 1000)) 727 (and (consp (nth 1 math-cur-record)) 728 (math-replace-integral-parts math-cur-record)) 729 (math-tracing-integral " => " 730 (math-format-value (nth 1 math-cur-record) 1000) 731 "\n"))) 732 (or (and math-cur-record 733 (not (eq (nth 1 math-cur-record) 'cancelled)) 734 (or (not (integerp (nth 1 math-cur-record))) 735 (>= (nth 1 math-cur-record) math-integ-level))) 736 (and (math-integral-contains-parts expr) 737 (progn 738 (setq val nil) 739 t)) 740 (unwind-protect 741 (progn 742 (let (math-integ-msg) 743 (if (eq calc-display-working-message 'lots) 744 (progn 745 (calc-set-command-flag 'clear-message) 746 (setq math-integ-msg (format 747 "Working... Integrating %s" 748 (math-format-flat-expr expr 0))) 749 (message math-integ-msg))) 750 (if math-cur-record 751 (setcar (cdr math-cur-record) 752 (if same-as-above (vector simp) 'busy)) 753 (setq math-cur-record 754 (list expr (if same-as-above (vector simp) 'busy)) 755 math-integral-cache (cons math-cur-record 756 math-integral-cache))) 757 (if (eq simplify 'yes) 758 (progn 759 (math-tracing-integral "Simplifying...") 760 (setq simp (math-simplify expr)) 761 (setq val (if (equal simp expr) 762 (progn 763 (math-tracing-integral " no change\n") 764 (math-do-integral expr)) 765 (math-tracing-integral " simplified\n") 766 (math-integral simp 'no t)))) 767 (or (setq val (math-do-integral expr)) 768 (eq simplify 'no) 769 (let ((simp (math-simplify expr))) 770 (or (equal simp expr) 771 (progn 772 (math-tracing-integral "Trying again after " 773 "simplification...\n") 774 (setq val (math-integral simp 'no t)))))))) 775 (if (eq calc-display-working-message 'lots) 776 (message math-integ-msg))) 777 (setcar (cdr math-cur-record) (or val 778 (if (or math-enable-subst 779 (not math-any-substs)) 780 math-integ-level 781 'cancelled))))) 782 (setq val math-cur-record) 783 (while (vectorp (nth 1 val)) 784 (setq val (aref (nth 1 val) 0))) 785 (setq val (if (memq (nth 1 val) '(parts parts2)) 786 (progn 787 (setcar (cdr val) 'parts2) 788 (list 'var 'PARTS val)) 789 (and (consp (nth 1 val)) 790 (nth 1 val)))) 791 (math-tracing-integral "Integral of " 792 (math-format-value expr 1000) 793 " is " 794 (math-format-value val 1000) 795 "\n") 796 val)) 797 798(defun math-integral-contains-parts (expr) 799 (if (Math-primp expr) 800 (and (eq (car-safe expr) 'var) 801 (eq (nth 1 expr) 'PARTS) 802 (listp (nth 2 expr))) 803 (while (and (setq expr (cdr expr)) 804 (not (math-integral-contains-parts (car expr))))) 805 expr)) 806 807(defun math-replace-integral-parts (expr) 808 (or (Math-primp expr) 809 (while (setq expr (cdr expr)) 810 (and (consp (car expr)) 811 (if (eq (car (car expr)) 'var) 812 (and (eq (nth 1 (car expr)) 'PARTS) 813 (consp (nth 2 (car expr))) 814 (if (listp (nth 1 (nth 2 (car expr)))) 815 (progn 816 (setcar expr (nth 1 (nth 2 (car expr)))) 817 (math-replace-integral-parts (cons 'foo expr))) 818 (setcar (cdr math-cur-record) 'cancelled))) 819 (math-replace-integral-parts (car expr))))))) 820 821(defvar math-linear-subst-tried t 822 "Non-nil means that a linear substitution has been tried.") 823 824;; The variable math-has-rules is a local variable for math-try-integral, 825;; but is used by math-do-integral, which is called (non-directly) by 826;; math-try-integral. 827(defvar math-has-rules) 828 829;; math-old-integ is a local variable for math-do-integral, but is 830;; used by math-sub-integration. 831(defvar math-old-integ) 832 833;; The variables math-t1, math-t2 and math-t3 are local to 834;; math-do-integral, math-try-solve-for and math-decompose-poly, but 835;; are used by functions they call (directly or indirectly); 836;; math-do-integral calls math-do-integral-methods; 837;; math-try-solve-for calls math-try-solve-prod, 838;; math-solve-find-root-term and math-solve-find-root-in-prod; 839;; math-decompose-poly calls math-solve-poly-funny-powers and 840;; math-solve-crunch-poly. 841(defvar math-t1) 842(defvar math-t2) 843(defvar math-t3) 844 845(defun math-do-integral (expr) 846 (let ((math-linear-subst-tried nil) 847 math-t1 math-t2) 848 (or (cond ((not (math-expr-contains expr math-integ-var)) 849 (math-mul expr math-integ-var)) 850 ((equal expr math-integ-var) 851 (math-div (math-sqr expr) 2)) 852 ((eq (car expr) '+) 853 (and (setq math-t1 (math-integral (nth 1 expr))) 854 (setq math-t2 (math-integral (nth 2 expr))) 855 (math-add math-t1 math-t2))) 856 ((eq (car expr) '-) 857 (and (setq math-t1 (math-integral (nth 1 expr))) 858 (setq math-t2 (math-integral (nth 2 expr))) 859 (math-sub math-t1 math-t2))) 860 ((eq (car expr) 'neg) 861 (and (setq math-t1 (math-integral (nth 1 expr))) 862 (math-neg math-t1))) 863 ((eq (car expr) '*) 864 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var)) 865 (and (setq math-t1 (math-integral (nth 2 expr))) 866 (math-mul (nth 1 expr) math-t1))) 867 ((not (math-expr-contains (nth 2 expr) math-integ-var)) 868 (and (setq math-t1 (math-integral (nth 1 expr))) 869 (math-mul math-t1 (nth 2 expr)))) 870 ((memq (car-safe (nth 1 expr)) '(+ -)) 871 (math-integral (list (car (nth 1 expr)) 872 (math-mul (nth 1 (nth 1 expr)) 873 (nth 2 expr)) 874 (math-mul (nth 2 (nth 1 expr)) 875 (nth 2 expr))) 876 'yes t)) 877 ((memq (car-safe (nth 2 expr)) '(+ -)) 878 (math-integral (list (car (nth 2 expr)) 879 (math-mul (nth 1 (nth 2 expr)) 880 (nth 1 expr)) 881 (math-mul (nth 2 (nth 2 expr)) 882 (nth 1 expr))) 883 'yes t)))) 884 ((eq (car expr) '/) 885 (cond ((and (not (math-expr-contains (nth 1 expr) 886 math-integ-var)) 887 (not (math-equal-int (nth 1 expr) 1))) 888 (and (setq math-t1 (math-integral (math-div 1 (nth 2 expr)))) 889 (math-mul (nth 1 expr) math-t1))) 890 ((not (math-expr-contains (nth 2 expr) math-integ-var)) 891 (and (setq math-t1 (math-integral (nth 1 expr))) 892 (math-div math-t1 (nth 2 expr)))) 893 ((and (eq (car-safe (nth 1 expr)) '*) 894 (not (math-expr-contains (nth 1 (nth 1 expr)) 895 math-integ-var))) 896 (and (setq math-t1 (math-integral 897 (math-div (nth 2 (nth 1 expr)) 898 (nth 2 expr)))) 899 (math-mul math-t1 (nth 1 (nth 1 expr))))) 900 ((and (eq (car-safe (nth 1 expr)) '*) 901 (not (math-expr-contains (nth 2 (nth 1 expr)) 902 math-integ-var))) 903 (and (setq math-t1 (math-integral 904 (math-div (nth 1 (nth 1 expr)) 905 (nth 2 expr)))) 906 (math-mul math-t1 (nth 2 (nth 1 expr))))) 907 ((and (eq (car-safe (nth 2 expr)) '*) 908 (not (math-expr-contains (nth 1 (nth 2 expr)) 909 math-integ-var))) 910 (and (setq math-t1 (math-integral 911 (math-div (nth 1 expr) 912 (nth 2 (nth 2 expr))))) 913 (math-div math-t1 (nth 1 (nth 2 expr))))) 914 ((and (eq (car-safe (nth 2 expr)) '*) 915 (not (math-expr-contains (nth 2 (nth 2 expr)) 916 math-integ-var))) 917 (and (setq math-t1 (math-integral 918 (math-div (nth 1 expr) 919 (nth 1 (nth 2 expr))))) 920 (math-div math-t1 (nth 2 (nth 2 expr))))) 921 ((eq (car-safe (nth 2 expr)) 'calcFunc-exp) 922 (math-integral 923 (math-mul (nth 1 expr) 924 (list 'calcFunc-exp 925 (math-neg (nth 1 (nth 2 expr))))))))) 926 ((eq (car expr) '^) 927 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var)) 928 (or (and (setq math-t1 (math-is-polynomial (nth 2 expr) 929 math-integ-var 1)) 930 (math-div expr 931 (math-mul (nth 1 math-t1) 932 (math-normalize 933 (list 'calcFunc-ln 934 (nth 1 expr)))))) 935 (math-integral 936 (list 'calcFunc-exp 937 (math-mul (nth 2 expr) 938 (math-normalize 939 (list 'calcFunc-ln 940 (nth 1 expr))))) 941 'yes t))) 942 ((not (math-expr-contains (nth 2 expr) math-integ-var)) 943 (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0)) 944 (math-integral 945 (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr)))) 946 nil t) 947 (or (and (setq math-t1 (math-is-polynomial (nth 1 expr) 948 math-integ-var 949 1)) 950 (setq math-t2 (math-add (nth 2 expr) 1)) 951 (math-div (math-pow (nth 1 expr) math-t2) 952 (math-mul math-t2 (nth 1 math-t1)))) 953 (and (Math-negp (nth 2 expr)) 954 (math-integral 955 (math-div 1 956 (math-pow (nth 1 expr) 957 (math-neg 958 (nth 2 expr)))) 959 nil t)) 960 nil)))))) 961 962 ;; Integral of a polynomial. 963 (and (setq math-t1 (math-is-polynomial expr math-integ-var 20)) 964 (let ((accum 0) 965 (n 1)) 966 (while math-t1 967 (if (setq accum (math-add accum 968 (math-div (math-mul (car math-t1) 969 (math-pow 970 math-integ-var 971 n)) 972 n)) 973 math-t1 (cdr math-t1)) 974 (setq n (1+ n)))) 975 accum)) 976 977 ;; Try looking it up! 978 (cond ((= (length expr) 2) 979 (and (symbolp (car expr)) 980 (setq math-t1 (get (car expr) 'math-integral)) 981 (progn 982 (while (and math-t1 983 (not (setq math-t2 (funcall (car math-t1) 984 (nth 1 expr))))) 985 (setq math-t1 (cdr math-t1))) 986 (and math-t2 (math-normalize math-t2))))) 987 ((= (length expr) 3) 988 (and (symbolp (car expr)) 989 (setq math-t1 (get (car expr) 'math-integral-2)) 990 (progn 991 (while (and math-t1 992 (not (setq math-t2 (funcall (car math-t1) 993 (nth 1 expr) 994 (nth 2 expr))))) 995 (setq math-t1 (cdr math-t1))) 996 (and math-t2 (math-normalize math-t2)))))) 997 998 ;; Integral of a rational function. 999 (and (math-ratpoly-p expr math-integ-var) 1000 (setq math-t1 (calcFunc-apart expr math-integ-var)) 1001 (not (equal math-t1 expr)) 1002 (math-integral math-t1)) 1003 1004 ;; Try user-defined integration rules. 1005 (and math-has-rules 1006 (let ((math-old-integ (symbol-function 'calcFunc-integ)) 1007 (input (list 'calcFunc-integtry expr math-integ-var)) 1008 res part) 1009 (unwind-protect 1010 (progn 1011 (fset 'calcFunc-integ 'math-sub-integration) 1012 (setq res (math-rewrite input 1013 '(var IntegRules var-IntegRules) 1014 1)) 1015 (fset 'calcFunc-integ math-old-integ) 1016 (and (not (equal res input)) 1017 (if (setq part (math-expr-calls 1018 res '(calcFunc-integsubst))) 1019 (and (memq (length part) '(3 4 5)) 1020 (let ((parts (mapcar 1021 (function 1022 (lambda (x) 1023 (math-expr-subst 1024 x (nth 2 part) 1025 math-integ-var))) 1026 (cdr part)))) 1027 (math-integrate-by-substitution 1028 expr (car parts) t 1029 (or (nth 2 parts) 1030 (list 'calcFunc-integfailed 1031 math-integ-var)) 1032 (nth 3 parts)))) 1033 (if (not (math-expr-calls res 1034 '(calcFunc-integtry 1035 calcFunc-integfailed))) 1036 res)))) 1037 (fset 'calcFunc-integ math-old-integ)))) 1038 1039 ;; See if the function is a symbolic derivative. 1040 (and (string-match "'" (symbol-name (car expr))) 1041 (let ((name (symbol-name (car expr))) 1042 (p expr) (n 0) (which nil) (bad nil)) 1043 (while (setq n (1+ n) p (cdr p)) 1044 (if (equal (car p) math-integ-var) 1045 (if which (setq bad t) (setq which n)) 1046 (if (math-expr-contains (car p) math-integ-var) 1047 (setq bad t)))) 1048 (and which (not bad) 1049 (let ((prime (if (= which 1) "'" (format "'%d" which)))) 1050 (and (string-match (concat prime "\\('['0-9]*\\|$\\)") 1051 name) 1052 (cons (intern 1053 (concat 1054 (substring name 0 (match-beginning 0)) 1055 (substring name (+ (match-beginning 0) 1056 (length prime))))) 1057 (cdr expr))))))) 1058 1059 ;; Try transformation methods (parts, substitutions). 1060 (and (> math-integ-level 0) 1061 (math-do-integral-methods expr)) 1062 1063 ;; Try expanding the function's definition. 1064 (let ((res (math-expand-formula expr))) 1065 (and res 1066 (math-integral res)))))) 1067 1068(defun math-sub-integration (expr &rest rest) 1069 (or (if (or (not rest) 1070 (and (< math-integ-level math-integral-limit) 1071 (eq (car rest) math-integ-var))) 1072 (math-integral expr) 1073 (let ((res (apply math-old-integ expr rest))) 1074 (and (or (= math-integ-level math-integral-limit) 1075 (not (math-expr-calls res 'calcFunc-integ))) 1076 res))) 1077 (list 'calcFunc-integfailed expr))) 1078 1079;; math-so-far is a local variable for math-do-integral-methods, but 1080;; is used by math-integ-try-linear-substitutions and 1081;; math-integ-try-substitutions. 1082(defvar math-so-far) 1083 1084;; math-integ-expr is a local variable for math-do-integral-methods, 1085;; but is used by math-integ-try-linear-substitutions and 1086;; math-integ-try-substitutions. 1087(defvar math-integ-expr) 1088 1089(defun math-do-integral-methods (math-integ-expr) 1090 (let ((math-so-far math-integ-var-list-list) 1091 rat-in) 1092 1093 ;; Integration by substitution, for various likely sub-expressions. 1094 ;; (In first pass, we look only for sub-exprs that are linear in X.) 1095 (or (math-integ-try-linear-substitutions math-integ-expr) 1096 (math-integ-try-substitutions math-integ-expr) 1097 1098 ;; If function has sines and cosines, try tan(x/2) substitution. 1099 (and (let ((p (setq rat-in (math-expr-rational-in math-integ-expr)))) 1100 (while (and p 1101 (memq (car (car p)) '(calcFunc-sin 1102 calcFunc-cos 1103 calcFunc-tan 1104 calcFunc-sec 1105 calcFunc-csc 1106 calcFunc-cot)) 1107 (equal (nth 1 (car p)) math-integ-var)) 1108 (setq p (cdr p))) 1109 (null p)) 1110 (or (and (math-integ-parts-easy math-integ-expr) 1111 (math-integ-try-parts math-integ-expr t)) 1112 (math-integrate-by-good-substitution 1113 math-integ-expr (list 'calcFunc-tan (math-div math-integ-var 2))))) 1114 1115 ;; If function has sinh and cosh, try tanh(x/2) substitution. 1116 (and (let ((p rat-in)) 1117 (while (and p 1118 (memq (car (car p)) '(calcFunc-sinh 1119 calcFunc-cosh 1120 calcFunc-tanh 1121 calcFunc-sech 1122 calcFunc-csch 1123 calcFunc-coth 1124 calcFunc-exp)) 1125 (equal (nth 1 (car p)) math-integ-var)) 1126 (setq p (cdr p))) 1127 (null p)) 1128 (or (and (math-integ-parts-easy math-integ-expr) 1129 (math-integ-try-parts math-integ-expr t)) 1130 (math-integrate-by-good-substitution 1131 math-integ-expr (list 'calcFunc-tanh (math-div math-integ-var 2))))) 1132 1133 ;; If function has square roots, try sin, tan, or sec substitution. 1134 (and (let ((p rat-in)) 1135 (setq math-t1 nil) 1136 (while (and p 1137 (or (equal (car p) math-integ-var) 1138 (and (eq (car (car p)) 'calcFunc-sqrt) 1139 (setq math-t1 (math-is-polynomial 1140 (nth 1 (setq math-t2 (car p))) 1141 math-integ-var 2))))) 1142 (setq p (cdr p))) 1143 (and (null p) math-t1)) 1144 (if (cdr (cdr math-t1)) 1145 (if (math-guess-if-neg (nth 2 math-t1)) 1146 (let* ((c (math-sqrt (math-neg (nth 2 math-t1)))) 1147 (d (math-div (nth 1 math-t1) (math-mul -2 c))) 1148 (a (math-sqrt (math-add (car math-t1) (math-sqr d))))) 1149 (math-integrate-by-good-substitution 1150 math-integ-expr (list 'calcFunc-arcsin 1151 (math-div-thru 1152 (math-add (math-mul c math-integ-var) d) 1153 a)))) 1154 (let* ((c (math-sqrt (nth 2 math-t1))) 1155 (d (math-div (nth 1 math-t1) (math-mul 2 c))) 1156 (aa (math-sub (car math-t1) (math-sqr d)))) 1157 (if (and nil (not (and (eq d 0) (eq c 1)))) 1158 (math-integrate-by-good-substitution 1159 math-integ-expr (math-add (math-mul c math-integ-var) d)) 1160 (if (math-guess-if-neg aa) 1161 (math-integrate-by-good-substitution 1162 math-integ-expr (list 'calcFunc-arccosh 1163 (math-div-thru 1164 (math-add (math-mul c math-integ-var) 1165 d) 1166 (math-sqrt (math-neg aa))))) 1167 (math-integrate-by-good-substitution 1168 math-integ-expr (list 'calcFunc-arcsinh 1169 (math-div-thru 1170 (math-add (math-mul c math-integ-var) 1171 d) 1172 (math-sqrt aa)))))))) 1173 (math-integrate-by-good-substitution math-integ-expr math-t2)) ) 1174 1175 ;; Try integration by parts. 1176 (math-integ-try-parts math-integ-expr) 1177 1178 ;; Give up. 1179 nil))) 1180 1181(defun math-integ-parts-easy (expr) 1182 (cond ((Math-primp expr) t) 1183 ((memq (car expr) '(+ - *)) 1184 (and (math-integ-parts-easy (nth 1 expr)) 1185 (math-integ-parts-easy (nth 2 expr)))) 1186 ((eq (car expr) '/) 1187 (and (math-integ-parts-easy (nth 1 expr)) 1188 (math-atomic-factorp (nth 2 expr)))) 1189 ((eq (car expr) '^) 1190 (and (natnump (nth 2 expr)) 1191 (math-integ-parts-easy (nth 1 expr)))) 1192 ((eq (car expr) 'neg) 1193 (math-integ-parts-easy (nth 1 expr))) 1194 (t t))) 1195 1196;; math-prev-parts-v is local to calcFunc-integ (as well as 1197;; math-integrate-by-parts), but is used by math-integ-try-parts. 1198(defvar math-prev-parts-v) 1199 1200;; math-good-parts is local to calcFunc-integ (as well as 1201;; math-integ-try-parts), but is used by math-integrate-by-parts. 1202(defvar math-good-parts) 1203 1204 1205(defun math-integ-try-parts (expr &optional math-good-parts) 1206 ;; Integration by parts: 1207 ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x) 1208 ;; where h(x) = integ(g(x),x). 1209 (or (let ((exp (calcFunc-expand expr))) 1210 (and (not (equal exp expr)) 1211 (math-integral exp))) 1212 (and (eq (car expr) '*) 1213 (let ((first-bad (or (math-polynomial-p (nth 1 expr) 1214 math-integ-var) 1215 (equal (nth 2 expr) math-prev-parts-v)))) 1216 (or (and first-bad ; so try this one first 1217 (math-integrate-by-parts (nth 1 expr) (nth 2 expr))) 1218 (math-integrate-by-parts (nth 2 expr) (nth 1 expr)) 1219 (and (not first-bad) 1220 (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))))) 1221 (and (eq (car expr) '/) 1222 (math-expr-contains (nth 1 expr) math-integ-var) 1223 (let ((recip (math-div 1 (nth 2 expr)))) 1224 (or (math-integrate-by-parts (nth 1 expr) recip) 1225 (math-integrate-by-parts recip (nth 1 expr))))) 1226 (and (eq (car expr) '^) 1227 (math-integrate-by-parts (math-pow (nth 1 expr) 1228 (math-sub (nth 2 expr) 1)) 1229 (nth 1 expr))))) 1230 1231(defun math-integrate-by-parts (u vprime) 1232 (let ((math-integ-level (if (or math-good-parts 1233 (math-polynomial-p u math-integ-var)) 1234 math-integ-level 1235 (1- math-integ-level))) 1236 (math-doing-parts t) 1237 v temp) 1238 (and (>= math-integ-level 0) 1239 (unwind-protect 1240 (progn 1241 (setcar (cdr math-cur-record) 'parts) 1242 (math-tracing-integral "Integrating by parts, u = " 1243 (math-format-value u 1000) 1244 ", v' = " 1245 (math-format-value vprime 1000) 1246 "\n") 1247 (and (setq v (math-integral vprime)) 1248 (setq temp (calcFunc-deriv u math-integ-var nil t)) 1249 (setq temp (let ((math-prev-parts-v v)) 1250 (math-integral (math-mul v temp) 'yes))) 1251 (setq temp (math-sub (math-mul u v) temp)) 1252 (if (eq (nth 1 math-cur-record) 'parts) 1253 (calcFunc-expand temp) 1254 (setq v (list 'var 'PARTS math-cur-record) 1255 temp (let (calc-next-why) 1256 (math-simplify-extended 1257 (math-solve-for (math-sub v temp) 0 v nil))) 1258 temp (if (and (eq (car-safe temp) '/) 1259 (math-zerop (nth 2 temp))) 1260 nil temp))))) 1261 (setcar (cdr math-cur-record) 'busy))))) 1262 1263;;; This tries two different formulations, hoping the algebraic simplifier 1264;;; will be strong enough to handle at least one. 1265(defun math-integrate-by-substitution (expr u &optional user uinv uinvprime) 1266 (and (> math-integ-level 0) 1267 (let ((math-integ-level (max (- math-integ-level 2) 0))) 1268 (math-integrate-by-good-substitution expr u user uinv uinvprime)))) 1269 1270(defun math-integrate-by-good-substitution (expr u &optional user 1271 uinv uinvprime) 1272 (let ((math-living-dangerously t) 1273 deriv temp) 1274 (and (setq uinv (if uinv 1275 (math-expr-subst uinv math-integ-var 1276 math-integ-var-2) 1277 (let (calc-next-why) 1278 (math-solve-for u 1279 math-integ-var-2 1280 math-integ-var nil)))) 1281 (progn 1282 (math-tracing-integral "Integrating by substitution, u = " 1283 (math-format-value u 1000) 1284 "\n") 1285 (or (and (setq deriv (calcFunc-deriv u 1286 math-integ-var nil 1287 (not user))) 1288 (setq temp (math-integral (math-expr-subst 1289 (math-expr-subst 1290 (math-expr-subst 1291 (math-div expr deriv) 1292 u 1293 math-integ-var-2) 1294 math-integ-var 1295 uinv) 1296 math-integ-var-2 1297 math-integ-var) 1298 'yes))) 1299 (and (setq deriv (or uinvprime 1300 (calcFunc-deriv uinv 1301 math-integ-var-2 1302 math-integ-var 1303 (not user)))) 1304 (setq temp (math-integral (math-mul 1305 (math-expr-subst 1306 (math-expr-subst 1307 (math-expr-subst 1308 expr 1309 u 1310 math-integ-var-2) 1311 math-integ-var 1312 uinv) 1313 math-integ-var-2 1314 math-integ-var) 1315 deriv) 1316 'yes))))) 1317 (math-simplify-extended 1318 (math-expr-subst temp math-integ-var u))))) 1319 1320;;; Look for substitutions of the form u = a x + b. 1321(defun math-integ-try-linear-substitutions (sub-expr) 1322 (setq math-linear-subst-tried t) 1323 (and (not (Math-primp sub-expr)) 1324 (or (and (not (memq (car sub-expr) '(+ - * / neg))) 1325 (not (and (eq (car sub-expr) '^) 1326 (integerp (nth 2 sub-expr)))) 1327 (math-expr-contains sub-expr math-integ-var) 1328 (let ((res nil)) 1329 (while (and (setq sub-expr (cdr sub-expr)) 1330 (or (not (math-linear-in (car sub-expr) 1331 math-integ-var)) 1332 (assoc (car sub-expr) math-so-far) 1333 (progn 1334 (setq math-so-far (cons (list (car sub-expr)) 1335 math-so-far)) 1336 (not (setq res 1337 (math-integrate-by-substitution 1338 math-integ-expr (car sub-expr)))))))) 1339 res)) 1340 (let ((res nil)) 1341 (while (and (setq sub-expr (cdr sub-expr)) 1342 (not (setq res (math-integ-try-linear-substitutions 1343 (car sub-expr)))))) 1344 res)))) 1345 1346;;; Recursively try different substitutions based on various sub-expressions. 1347(defun math-integ-try-substitutions (sub-expr &optional allow-rat) 1348 (and (not (Math-primp sub-expr)) 1349 (not (assoc sub-expr math-so-far)) 1350 (math-expr-contains sub-expr math-integ-var) 1351 (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg))) 1352 (not (and (eq (car sub-expr) '^) 1353 (integerp (nth 2 sub-expr))))) 1354 (setq allow-rat t) 1355 (prog1 allow-rat (setq allow-rat nil))) 1356 (not (eq sub-expr math-integ-expr)) 1357 (or (math-integrate-by-substitution math-integ-expr sub-expr) 1358 (and (eq (car sub-expr) '^) 1359 (integerp (nth 2 sub-expr)) 1360 (< (nth 2 sub-expr) 0) 1361 (math-integ-try-substitutions 1362 (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr))) 1363 t)))) 1364 (let ((res nil)) 1365 (setq math-so-far (cons (list sub-expr) math-so-far)) 1366 (while (and (setq sub-expr (cdr sub-expr)) 1367 (not (setq res (math-integ-try-substitutions 1368 (car sub-expr) allow-rat))))) 1369 res)))) 1370 1371;; The variable math-expr-parts is local to math-expr-rational-in, 1372;; but is used by math-expr-rational-in-rec 1373(defvar math-expr-parts) 1374 1375(defun math-expr-rational-in (expr) 1376 (let ((math-expr-parts nil)) 1377 (math-expr-rational-in-rec expr) 1378 (mapcar 'car math-expr-parts))) 1379 1380(defun math-expr-rational-in-rec (expr) 1381 (cond ((Math-primp expr) 1382 (and (equal expr math-integ-var) 1383 (not (assoc expr math-expr-parts)) 1384 (setq math-expr-parts (cons (list expr) math-expr-parts)))) 1385 ((or (memq (car expr) '(+ - * / neg)) 1386 (and (eq (car expr) '^) (integerp (nth 2 expr)))) 1387 (math-expr-rational-in-rec (nth 1 expr)) 1388 (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr)))) 1389 ((and (eq (car expr) '^) 1390 (eq (math-quarter-integer (nth 2 expr)) 2)) 1391 (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr)))) 1392 (t 1393 (and (not (assoc expr math-expr-parts)) 1394 (math-expr-contains expr math-integ-var) 1395 (setq math-expr-parts (cons (list expr) math-expr-parts)))))) 1396 1397(defun math-expr-calls (expr funcs &optional arg-contains) 1398 (if (consp expr) 1399 (if (or (memq (car expr) funcs) 1400 (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt) 1401 (eq (math-quarter-integer (nth 2 expr)) 2))) 1402 (and (or (not arg-contains) 1403 (math-expr-contains expr arg-contains)) 1404 expr) 1405 (and (not (Math-primp expr)) 1406 (let ((res nil)) 1407 (while (and (setq expr (cdr expr)) 1408 (not (setq res (math-expr-calls 1409 (car expr) funcs arg-contains))))) 1410 res))))) 1411 1412(defun math-fix-const-terms (expr except-vars) 1413 (cond ((not (math-expr-depends expr except-vars)) 0) 1414 ((Math-primp expr) expr) 1415 ((eq (car expr) '+) 1416 (math-add (math-fix-const-terms (nth 1 expr) except-vars) 1417 (math-fix-const-terms (nth 2 expr) except-vars))) 1418 ((eq (car expr) '-) 1419 (math-sub (math-fix-const-terms (nth 1 expr) except-vars) 1420 (math-fix-const-terms (nth 2 expr) except-vars))) 1421 (t expr))) 1422 1423;; Command for debugging the Calculator's symbolic integrator. 1424(defun calc-dump-integral-cache (&optional arg) 1425 (interactive "P") 1426 (let ((buf (current-buffer))) 1427 (unwind-protect 1428 (let ((p math-integral-cache) 1429 math-cur-record) 1430 (display-buffer (get-buffer-create "*Integral Cache*")) 1431 (set-buffer (get-buffer "*Integral Cache*")) 1432 (erase-buffer) 1433 (while p 1434 (setq math-cur-record (car p)) 1435 (or arg (math-replace-integral-parts math-cur-record)) 1436 (insert (math-format-flat-expr (car math-cur-record) 0) 1437 " --> " 1438 (if (symbolp (nth 1 math-cur-record)) 1439 (concat "(" (symbol-name (nth 1 math-cur-record)) ")") 1440 (math-format-flat-expr (nth 1 math-cur-record) 0)) 1441 "\n") 1442 (setq p (cdr p))) 1443 (goto-char (point-min))) 1444 (set-buffer buf)))) 1445 1446;; The variable math-max-integral-limit is local to calcFunc-integ, 1447;; but is used by math-try-integral. 1448(defvar math-max-integral-limit) 1449 1450(defun math-try-integral (expr) 1451 (let ((math-integ-level math-integral-limit) 1452 (math-integ-depth 0) 1453 (math-integ-msg "Working...done") 1454 (math-cur-record nil) ; a technicality 1455 (math-integrating t) 1456 (calc-prefer-frac t) 1457 (calc-symbolic-mode t) 1458 (math-has-rules (calc-has-rules 'var-IntegRules))) 1459 (or (math-integral expr 'yes) 1460 (and math-any-substs 1461 (setq math-enable-subst t) 1462 (math-integral expr 'yes)) 1463 (and (> math-max-integral-limit math-integral-limit) 1464 (setq math-integral-limit math-max-integral-limit 1465 math-integ-level math-integral-limit) 1466 (math-integral expr 'yes))))) 1467 1468(defvar var-IntegLimit nil) 1469 1470(defun calcFunc-integ (expr var &optional low high) 1471 (cond 1472 ;; Do these even if the parts turn out not to be integrable. 1473 ((eq (car-safe expr) '+) 1474 (math-add (calcFunc-integ (nth 1 expr) var low high) 1475 (calcFunc-integ (nth 2 expr) var low high))) 1476 ((eq (car-safe expr) '-) 1477 (math-sub (calcFunc-integ (nth 1 expr) var low high) 1478 (calcFunc-integ (nth 2 expr) var low high))) 1479 ((eq (car-safe expr) 'neg) 1480 (math-neg (calcFunc-integ (nth 1 expr) var low high))) 1481 ((and (eq (car-safe expr) '*) 1482 (not (math-expr-contains (nth 1 expr) var))) 1483 (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high))) 1484 ((and (eq (car-safe expr) '*) 1485 (not (math-expr-contains (nth 2 expr) var))) 1486 (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr))) 1487 ((and (eq (car-safe expr) '/) 1488 (not (math-expr-contains (nth 1 expr) var)) 1489 (not (math-equal-int (nth 1 expr) 1))) 1490 (math-mul (nth 1 expr) 1491 (calcFunc-integ (math-div 1 (nth 2 expr)) var low high))) 1492 ((and (eq (car-safe expr) '/) 1493 (not (math-expr-contains (nth 2 expr) var))) 1494 (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr))) 1495 ((and (eq (car-safe expr) '/) 1496 (eq (car-safe (nth 1 expr)) '*) 1497 (not (math-expr-contains (nth 1 (nth 1 expr)) var))) 1498 (math-mul (nth 1 (nth 1 expr)) 1499 (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr)) 1500 var low high))) 1501 ((and (eq (car-safe expr) '/) 1502 (eq (car-safe (nth 1 expr)) '*) 1503 (not (math-expr-contains (nth 2 (nth 1 expr)) var))) 1504 (math-mul (nth 2 (nth 1 expr)) 1505 (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr)) 1506 var low high))) 1507 ((and (eq (car-safe expr) '/) 1508 (eq (car-safe (nth 2 expr)) '*) 1509 (not (math-expr-contains (nth 1 (nth 2 expr)) var))) 1510 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr))) 1511 var low high) 1512 (nth 1 (nth 2 expr)))) 1513 ((and (eq (car-safe expr) '/) 1514 (eq (car-safe (nth 2 expr)) '*) 1515 (not (math-expr-contains (nth 2 (nth 2 expr)) var))) 1516 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr))) 1517 var low high) 1518 (nth 2 (nth 2 expr)))) 1519 ((eq (car-safe expr) 'vec) 1520 (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high))) 1521 (cdr expr)))) 1522 (t 1523 (let ((state (list calc-angle-mode 1524 ;;calc-symbolic-mode 1525 ;;calc-prefer-frac 1526 calc-internal-prec 1527 (calc-var-value 'var-IntegRules) 1528 (calc-var-value 'var-IntegSimpRules)))) 1529 (or (equal state math-integral-cache-state) 1530 (setq math-integral-cache-state state 1531 math-integral-cache nil))) 1532 (let* ((math-max-integral-limit (or (and (natnump var-IntegLimit) 1533 var-IntegLimit) 1534 3)) 1535 (math-integral-limit 1) 1536 (sexpr (math-expr-subst expr var math-integ-var)) 1537 (trace-buffer (get-buffer "*Trace*")) 1538 (calc-language (if (eq calc-language 'big) nil calc-language)) 1539 (math-any-substs t) 1540 (math-enable-subst nil) 1541 (math-prev-parts-v nil) 1542 (math-doing-parts nil) 1543 (math-good-parts nil) 1544 (res 1545 (if trace-buffer 1546 (let ((calcbuf (current-buffer)) 1547 (calcwin (selected-window))) 1548 (unwind-protect 1549 (progn 1550 (if (get-buffer-window trace-buffer) 1551 (select-window (get-buffer-window trace-buffer))) 1552 (set-buffer trace-buffer) 1553 (goto-char (point-max)) 1554 (or (assq 'scroll-stop (buffer-local-variables)) 1555 (progn 1556 (make-local-variable 'scroll-step) 1557 (setq scroll-step 3))) 1558 (insert "\n\n\n") 1559 (set-buffer calcbuf) 1560 (math-try-integral sexpr)) 1561 (select-window calcwin) 1562 (set-buffer calcbuf))) 1563 (math-try-integral sexpr)))) 1564 (if res 1565 (progn 1566 (if (calc-has-rules 'var-IntegAfterRules) 1567 (setq res (math-rewrite res '(var IntegAfterRules 1568 var-IntegAfterRules)))) 1569 (math-simplify 1570 (if (and low high) 1571 (math-sub (math-expr-subst res math-integ-var high) 1572 (math-expr-subst res math-integ-var low)) 1573 (setq res (math-fix-const-terms res math-integ-vars)) 1574 (if low 1575 (math-expr-subst res math-integ-var low) 1576 (math-expr-subst res math-integ-var var))))) 1577 (append (list 'calcFunc-integ expr var) 1578 (and low (list low)) 1579 (and high (list high)))))))) 1580 1581 1582(math-defintegral calcFunc-inv 1583 (math-integral (math-div 1 u))) 1584 1585(math-defintegral calcFunc-conj 1586 (let ((int (math-integral u))) 1587 (and int 1588 (list 'calcFunc-conj int)))) 1589 1590(math-defintegral calcFunc-deg 1591 (let ((int (math-integral u))) 1592 (and int 1593 (list 'calcFunc-deg int)))) 1594 1595(math-defintegral calcFunc-rad 1596 (let ((int (math-integral u))) 1597 (and int 1598 (list 'calcFunc-rad int)))) 1599 1600(math-defintegral calcFunc-re 1601 (let ((int (math-integral u))) 1602 (and int 1603 (list 'calcFunc-re int)))) 1604 1605(math-defintegral calcFunc-im 1606 (let ((int (math-integral u))) 1607 (and int 1608 (list 'calcFunc-im int)))) 1609 1610(math-defintegral calcFunc-sqrt 1611 (and (equal u math-integ-var) 1612 (math-mul '(frac 2 3) 1613 (list 'calcFunc-sqrt (math-pow u 3))))) 1614 1615(math-defintegral calcFunc-exp 1616 (or (and (equal u math-integ-var) 1617 (list 'calcFunc-exp u)) 1618 (let ((p (math-is-polynomial u math-integ-var 2))) 1619 (and (nth 2 p) 1620 (let ((sqa (math-sqrt (math-neg (nth 2 p))))) 1621 (math-div 1622 (math-mul 1623 (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi)) 1624 sqa) 1625 (math-normalize 1626 (list 'calcFunc-exp 1627 (math-div (math-sub (math-mul (car p) 1628 (nth 2 p)) 1629 (math-div 1630 (math-sqr (nth 1 p)) 1631 4)) 1632 (nth 2 p))))) 1633 (list 'calcFunc-erf 1634 (math-sub (math-mul sqa math-integ-var) 1635 (math-div (nth 1 p) (math-mul 2 sqa))))) 1636 2)))))) 1637 1638(math-defintegral calcFunc-ln 1639 (or (and (equal u math-integ-var) 1640 (math-sub (math-mul u (list 'calcFunc-ln u)) u)) 1641 (and (eq (car u) '*) 1642 (math-integral (math-add (list 'calcFunc-ln (nth 1 u)) 1643 (list 'calcFunc-ln (nth 2 u))))) 1644 (and (eq (car u) '/) 1645 (math-integral (math-sub (list 'calcFunc-ln (nth 1 u)) 1646 (list 'calcFunc-ln (nth 2 u))))) 1647 (and (eq (car u) '^) 1648 (math-integral (math-mul (nth 2 u) 1649 (list 'calcFunc-ln (nth 1 u))))))) 1650 1651(math-defintegral calcFunc-log10 1652 (and (equal u math-integ-var) 1653 (math-sub (math-mul u (list 'calcFunc-ln u)) 1654 (math-div u (list 'calcFunc-ln 10))))) 1655 1656(math-defintegral-2 calcFunc-log 1657 (math-integral (math-div (list 'calcFunc-ln u) 1658 (list 'calcFunc-ln v)))) 1659 1660(math-defintegral calcFunc-sin 1661 (or (and (equal u math-integ-var) 1662 (math-neg (math-from-radians-2 (list 'calcFunc-cos u)))) 1663 (and (nth 2 (math-is-polynomial u math-integ-var 2)) 1664 (math-integral (math-to-exponentials (list 'calcFunc-sin u)))))) 1665 1666(math-defintegral calcFunc-cos 1667 (or (and (equal u math-integ-var) 1668 (math-from-radians-2 (list 'calcFunc-sin u))) 1669 (and (nth 2 (math-is-polynomial u math-integ-var 2)) 1670 (math-integral (math-to-exponentials (list 'calcFunc-cos u)))))) 1671 1672(math-defintegral calcFunc-tan 1673 (and (equal u math-integ-var) 1674 (math-from-radians-2 1675 (list 'calcFunc-ln (list 'calcFunc-sec u))))) 1676 1677(math-defintegral calcFunc-sec 1678 (and (equal u math-integ-var) 1679 (math-from-radians-2 1680 (list 'calcFunc-ln 1681 (math-add 1682 (list 'calcFunc-sec u) 1683 (list 'calcFunc-tan u)))))) 1684 1685(math-defintegral calcFunc-csc 1686 (and (equal u math-integ-var) 1687 (math-from-radians-2 1688 (list 'calcFunc-ln 1689 (math-sub 1690 (list 'calcFunc-csc u) 1691 (list 'calcFunc-cot u)))))) 1692 1693(math-defintegral calcFunc-cot 1694 (and (equal u math-integ-var) 1695 (math-from-radians-2 1696 (list 'calcFunc-ln (list 'calcFunc-sin u))))) 1697 1698(math-defintegral calcFunc-arcsin 1699 (and (equal u math-integ-var) 1700 (math-add (math-mul u (list 'calcFunc-arcsin u)) 1701 (math-from-radians-2 1702 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))) 1703 1704(math-defintegral calcFunc-arccos 1705 (and (equal u math-integ-var) 1706 (math-sub (math-mul u (list 'calcFunc-arccos u)) 1707 (math-from-radians-2 1708 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))) 1709 1710(math-defintegral calcFunc-arctan 1711 (and (equal u math-integ-var) 1712 (math-sub (math-mul u (list 'calcFunc-arctan u)) 1713 (math-from-radians-2 1714 (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u))) 1715 2))))) 1716 1717(math-defintegral calcFunc-sinh 1718 (and (equal u math-integ-var) 1719 (list 'calcFunc-cosh u))) 1720 1721(math-defintegral calcFunc-cosh 1722 (and (equal u math-integ-var) 1723 (list 'calcFunc-sinh u))) 1724 1725(math-defintegral calcFunc-tanh 1726 (and (equal u math-integ-var) 1727 (list 'calcFunc-ln (list 'calcFunc-cosh u)))) 1728 1729(math-defintegral calcFunc-sech 1730 (and (equal u math-integ-var) 1731 (list 'calcFunc-arctan (list 'calcFunc-sinh u)))) 1732 1733(math-defintegral calcFunc-csch 1734 (and (equal u math-integ-var) 1735 (list 'calcFunc-ln (list 'calcFunc-tanh (math-div u 2))))) 1736 1737(math-defintegral calcFunc-coth 1738 (and (equal u math-integ-var) 1739 (list 'calcFunc-ln (list 'calcFunc-sinh u)))) 1740 1741(math-defintegral calcFunc-arcsinh 1742 (and (equal u math-integ-var) 1743 (math-sub (math-mul u (list 'calcFunc-arcsinh u)) 1744 (list 'calcFunc-sqrt (math-add (math-sqr u) 1))))) 1745 1746(math-defintegral calcFunc-arccosh 1747 (and (equal u math-integ-var) 1748 (math-sub (math-mul u (list 'calcFunc-arccosh u)) 1749 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))) 1750 1751(math-defintegral calcFunc-arctanh 1752 (and (equal u math-integ-var) 1753 (math-sub (math-mul u (list 'calcFunc-arctan u)) 1754 (math-div (list 'calcFunc-ln 1755 (math-add 1 (math-sqr u))) 1756 2)))) 1757 1758;;; (Ax + B) / (ax^2 + bx + c)^n forms. 1759(math-defintegral-2 / 1760 (math-integral-rational-funcs u v)) 1761 1762(defun math-integral-rational-funcs (u v) 1763 (let ((pu (math-is-polynomial u math-integ-var 1)) 1764 (vpow 1) pv) 1765 (and pu 1766 (catch 'int-rat 1767 (if (and (eq (car-safe v) '^) (natnump (nth 2 v))) 1768 (setq vpow (nth 2 v) 1769 v (nth 1 v))) 1770 (and (setq pv (math-is-polynomial v math-integ-var 2)) 1771 (let ((int (math-mul-thru 1772 (car pu) 1773 (math-integral-q02 (car pv) (nth 1 pv) 1774 (nth 2 pv) v vpow)))) 1775 (if (cdr pu) 1776 (setq int (math-add int 1777 (math-mul-thru 1778 (nth 1 pu) 1779 (math-integral-q12 1780 (car pv) (nth 1 pv) 1781 (nth 2 pv) v vpow))))) 1782 int)))))) 1783 1784(defun math-integral-q12 (a b c v vpow) 1785 (let (q) 1786 (cond ((not c) 1787 (cond ((= vpow 1) 1788 (math-sub (math-div math-integ-var b) 1789 (math-mul (math-div a (math-sqr b)) 1790 (list 'calcFunc-ln v)))) 1791 ((= vpow 2) 1792 (math-div (math-add (list 'calcFunc-ln v) 1793 (math-div a v)) 1794 (math-sqr b))) 1795 (t 1796 (let ((nm1 (math-sub vpow 1)) 1797 (nm2 (math-sub vpow 2))) 1798 (math-div (math-sub 1799 (math-div a (math-mul nm1 (math-pow v nm1))) 1800 (math-div 1 (math-mul nm2 (math-pow v nm2)))) 1801 (math-sqr b)))))) 1802 ((math-zerop 1803 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b)))) 1804 (let ((part (math-div b (math-mul 2 c)))) 1805 (math-mul-thru (math-pow c vpow) 1806 (math-integral-q12 part 1 nil 1807 (math-add math-integ-var part) 1808 (* vpow 2))))) 1809 ((= vpow 1) 1810 (and (math-ratp q) (math-negp q) 1811 (let ((calc-symbolic-mode t)) 1812 (math-ratp (math-sqrt (math-neg q)))) 1813 (throw 'int-rat nil)) ; should have used calcFunc-apart first 1814 (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c)) 1815 (math-mul-thru (math-div b (math-mul 2 c)) 1816 (math-integral-q02 a b c v 1)))) 1817 (t 1818 (let ((n (1- vpow))) 1819 (math-sub (math-neg (math-div 1820 (math-add (math-mul b math-integ-var) 1821 (math-mul 2 a)) 1822 (math-mul n (math-mul q (math-pow v n))))) 1823 (math-mul-thru (math-div (math-mul b (1- (* 2 n))) 1824 (math-mul n q)) 1825 (math-integral-q02 a b c v n)))))))) 1826 1827(defun math-integral-q02 (a b c v vpow) 1828 (let (q rq part) 1829 (cond ((not c) 1830 (cond ((= vpow 1) 1831 (math-div (list 'calcFunc-ln v) b)) 1832 (t 1833 (math-div (math-pow v (- 1 vpow)) 1834 (math-mul (- 1 vpow) b))))) 1835 ((math-zerop 1836 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b)))) 1837 (let ((part (math-div b (math-mul 2 c)))) 1838 (math-mul-thru (math-pow c vpow) 1839 (math-integral-q02 part 1 nil 1840 (math-add math-integ-var part) 1841 (* vpow 2))))) 1842 ((progn 1843 (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b)) 1844 (> vpow 1)) 1845 (let ((n (1- vpow))) 1846 (math-add (math-div part (math-mul n (math-mul q (math-pow v n)))) 1847 (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c) 1848 (math-mul n q)) 1849 (math-integral-q02 a b c v n))))) 1850 ((math-guess-if-neg q) 1851 (setq rq (list 'calcFunc-sqrt (math-neg q))) 1852 ;;(math-div-thru (list 'calcFunc-ln 1853 ;; (math-div (math-sub part rq) 1854 ;; (math-add part rq))) 1855 ;; rq) 1856 (math-div (math-mul -2 (list 'calcFunc-arctanh 1857 (math-div part rq))) 1858 rq)) 1859 (t 1860 (setq rq (list 'calcFunc-sqrt q)) 1861 (math-div (math-mul 2 (math-to-radians-2 1862 (list 'calcFunc-arctan 1863 (math-div part rq)))) 1864 rq))))) 1865 1866 1867(math-defintegral calcFunc-erf 1868 (and (equal u math-integ-var) 1869 (math-add (math-mul u (list 'calcFunc-erf u)) 1870 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u)) 1871 (list 'calcFunc-sqrt 1872 '(var pi var-pi))))))) 1873 1874(math-defintegral calcFunc-erfc 1875 (and (equal u math-integ-var) 1876 (math-sub (math-mul u (list 'calcFunc-erfc u)) 1877 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u)) 1878 (list 'calcFunc-sqrt 1879 '(var pi var-pi))))))) 1880 1881 1882 1883 1884(defvar math-tabulate-initial nil) 1885(defvar math-tabulate-function nil) 1886 1887;; The variables calc-low and calc-high are local to calcFunc-table, 1888;; but are used by math-scan-for-limits. 1889(defvar calc-low) 1890(defvar calc-high) 1891 1892(defun calcFunc-table (expr var &optional calc-low calc-high step) 1893 (or calc-low 1894 (setq calc-low '(neg (var inf var-inf)) calc-high '(var inf var-inf))) 1895 (or calc-high (setq calc-high calc-low calc-low 1)) 1896 (and (or (math-infinitep calc-low) (math-infinitep calc-high)) 1897 (not step) 1898 (math-scan-for-limits expr)) 1899 (and step (math-zerop step) (math-reject-arg step 'nonzerop)) 1900 (let ((known (+ (if (Math-objectp calc-low) 1 0) 1901 (if (Math-objectp calc-high) 1 0) 1902 (if (or (null step) (Math-objectp step)) 1 0))) 1903 (count '(var inf var-inf)) 1904 vec) 1905 (or (= known 2) ; handy optimization 1906 (equal calc-high '(var inf var-inf)) 1907 (progn 1908 (setq count (math-div (math-sub calc-high calc-low) (or step 1))) 1909 (or (Math-objectp count) 1910 (setq count (math-simplify count))) 1911 (if (Math-messy-integerp count) 1912 (setq count (math-trunc count))))) 1913 (if (Math-negp count) 1914 (setq count -1)) 1915 (if (integerp count) 1916 (let ((var-DUMMY nil) 1917 (vec math-tabulate-initial) 1918 (math-working-step-2 (1+ count)) 1919 (math-working-step 0)) 1920 (setq expr (math-evaluate-expr 1921 (math-expr-subst expr var '(var DUMMY var-DUMMY)))) 1922 (while (>= count 0) 1923 (setq math-working-step (1+ math-working-step) 1924 var-DUMMY calc-low 1925 vec (cond ((eq math-tabulate-function 'calcFunc-sum) 1926 (math-add vec (math-evaluate-expr expr))) 1927 ((eq math-tabulate-function 'calcFunc-prod) 1928 (math-mul vec (math-evaluate-expr expr))) 1929 (t 1930 (cons (math-evaluate-expr expr) vec))) 1931 calc-low (math-add calc-low (or step 1)) 1932 count (1- count))) 1933 (if math-tabulate-function 1934 vec 1935 (cons 'vec (nreverse vec)))) 1936 (if (Math-integerp count) 1937 (calc-record-why 'fixnump calc-high) 1938 (if (Math-num-integerp calc-low) 1939 (if (Math-num-integerp calc-high) 1940 (calc-record-why 'integerp step) 1941 (calc-record-why 'integerp calc-high)) 1942 (calc-record-why 'integerp calc-low))) 1943 (append (list (or math-tabulate-function 'calcFunc-table) 1944 expr var) 1945 (and (not (and (equal calc-low '(neg (var inf var-inf))) 1946 (equal calc-high '(var inf var-inf)))) 1947 (list calc-low calc-high)) 1948 (and step (list step)))))) 1949 1950(defun math-scan-for-limits (x) 1951 (cond ((Math-primp x)) 1952 ((and (eq (car x) 'calcFunc-subscr) 1953 (Math-vectorp (nth 1 x)) 1954 (math-expr-contains (nth 2 x) var)) 1955 (let* ((calc-next-why nil) 1956 (low-val (math-solve-for (nth 2 x) 1 var nil)) 1957 (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x))) 1958 var nil)) 1959 temp) 1960 (and low-val (math-realp low-val) 1961 high-val (math-realp high-val)) 1962 (and (Math-lessp high-val low-val) 1963 (setq temp low-val low-val high-val high-val temp)) 1964 (setq calc-low (math-max calc-low (math-ceiling low-val)) 1965 calc-high (math-min calc-high (math-floor high-val))))) 1966 (t 1967 (while (setq x (cdr x)) 1968 (math-scan-for-limits (car x)))))) 1969 1970 1971(defvar math-disable-sums nil) 1972(defun calcFunc-sum (expr var &optional low high step) 1973 (if math-disable-sums (math-reject-arg)) 1974 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2))) 1975 (math-sum-rec expr var low high step))) 1976 (math-disable-sums t)) 1977 (math-normalize res))) 1978 1979(defun math-sum-rec (expr var &optional low high step) 1980 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) 1981 (and low (not high) (setq high low low 1)) 1982 (let (t1 t2 val) 1983 (setq val 1984 (cond 1985 ((not (math-expr-contains expr var)) 1986 (math-mul expr (math-add (math-div (math-sub high low) (or step 1)) 1987 1))) 1988 ((and step (not (math-equal-int step 1))) 1989 (if (math-negp step) 1990 (math-sum-rec expr var high low (math-neg step)) 1991 (let ((lo (math-simplify (math-div low step)))) 1992 (if (math-known-num-integerp lo) 1993 (math-sum-rec (math-normalize 1994 (math-expr-subst expr var 1995 (math-mul step var))) 1996 var lo (math-simplify (math-div high step))) 1997 (math-sum-rec (math-normalize 1998 (math-expr-subst expr var 1999 (math-add (math-mul step var) 2000 low))) 2001 var 0 2002 (math-simplify (math-div (math-sub high low) 2003 step))))))) 2004 ((memq (setq t1 (math-compare low high)) '(0 1)) 2005 (if (eq t1 0) 2006 (math-expr-subst expr var low) 2007 0)) 2008 ((setq t1 (math-is-polynomial expr var 20)) 2009 (let ((poly nil) 2010 (n 0)) 2011 (while t1 2012 (setq poly (math-poly-mix poly 1 2013 (math-sum-integer-power n) (car t1)) 2014 n (1+ n) 2015 t1 (cdr t1))) 2016 (setq n (math-build-polynomial-expr poly high)) 2017 (if (= low 1) 2018 n 2019 (math-sub n (math-build-polynomial-expr poly 2020 (math-sub low 1)))))) 2021 ((and (memq (car expr) '(+ -)) 2022 (setq t1 (math-sum-rec (nth 1 expr) var low high) 2023 t2 (math-sum-rec (nth 2 expr) var low high)) 2024 (not (and (math-expr-calls t1 '(calcFunc-sum)) 2025 (math-expr-calls t2 '(calcFunc-sum))))) 2026 (list (car expr) t1 t2)) 2027 ((and (eq (car expr) '*) 2028 (setq t1 (math-sum-const-factors expr var))) 2029 (math-mul (car t1) (math-sum-rec (cdr t1) var low high))) 2030 ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -))) 2031 (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr)) 2032 (nth 2 expr)) 2033 (math-mul (nth 2 (nth 1 expr)) 2034 (nth 2 expr)) 2035 nil (eq (car (nth 1 expr)) '-)) 2036 var low high)) 2037 ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -))) 2038 (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr) 2039 (nth 1 (nth 2 expr))) 2040 (math-mul (nth 1 expr) 2041 (nth 2 (nth 2 expr))) 2042 nil (eq (car (nth 2 expr)) '-)) 2043 var low high)) 2044 ((and (eq (car expr) '/) 2045 (not (math-primp (nth 1 expr))) 2046 (setq t1 (math-sum-const-factors (nth 1 expr) var))) 2047 (math-mul (car t1) 2048 (math-sum-rec (math-div (cdr t1) (nth 2 expr)) 2049 var low high))) 2050 ((and (eq (car expr) '/) 2051 (setq t1 (math-sum-const-factors (nth 2 expr) var))) 2052 (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1)) 2053 var low high) 2054 (car t1))) 2055 ((eq (car expr) 'neg) 2056 (math-neg (math-sum-rec (nth 1 expr) var low high))) 2057 ((and (eq (car expr) '^) 2058 (not (math-expr-contains (nth 1 expr) var)) 2059 (setq t1 (math-is-polynomial (nth 2 expr) var 1))) 2060 (let ((x (math-pow (nth 1 expr) (nth 1 t1)))) 2061 (math-div (math-mul (math-sub (math-pow x (math-add 1 high)) 2062 (math-pow x low)) 2063 (math-pow (nth 1 expr) (car t1))) 2064 (math-sub x 1)))) 2065 ((and (setq t1 (math-to-exponentials expr)) 2066 (setq t1 (math-sum-rec t1 var low high)) 2067 (not (math-expr-calls t1 '(calcFunc-sum)))) 2068 (math-to-exps t1)) 2069 ((memq (car expr) '(calcFunc-ln calcFunc-log10)) 2070 (list (car expr) (calcFunc-prod (nth 1 expr) var low high))) 2071 ((and (eq (car expr) 'calcFunc-log) 2072 (= (length expr) 3) 2073 (not (math-expr-contains (nth 2 expr) var))) 2074 (list 'calcFunc-log 2075 (calcFunc-prod (nth 1 expr) var low high) 2076 (nth 2 expr))))) 2077 (if (equal val '(var nan var-nan)) (setq val nil)) 2078 (or val 2079 (let* ((math-tabulate-initial 0) 2080 (math-tabulate-function 'calcFunc-sum)) 2081 (calcFunc-table expr var low high))))) 2082 2083(defun calcFunc-asum (expr var low &optional high step no-mul-flag) 2084 (or high (setq high low low 1)) 2085 (if (and step (not (math-equal-int step 1))) 2086 (if (math-negp step) 2087 (math-mul (math-pow -1 low) 2088 (calcFunc-asum expr var high low (math-neg step) t)) 2089 (let ((lo (math-simplify (math-div low step)))) 2090 (if (math-num-integerp lo) 2091 (calcFunc-asum (math-normalize 2092 (math-expr-subst expr var 2093 (math-mul step var))) 2094 var lo (math-simplify (math-div high step))) 2095 (calcFunc-asum (math-normalize 2096 (math-expr-subst expr var 2097 (math-add (math-mul step var) 2098 low))) 2099 var 0 2100 (math-simplify (math-div (math-sub high low) 2101 step)))))) 2102 (math-mul (if no-mul-flag 1 (math-pow -1 low)) 2103 (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high)))) 2104 2105(defun math-sum-const-factors (expr var) 2106 (let ((const nil) 2107 (not-const nil) 2108 (p expr)) 2109 (while (eq (car-safe p) '*) 2110 (if (math-expr-contains (nth 1 p) var) 2111 (setq not-const (cons (nth 1 p) not-const)) 2112 (setq const (cons (nth 1 p) const))) 2113 (setq p (nth 2 p))) 2114 (if (math-expr-contains p var) 2115 (setq not-const (cons p not-const)) 2116 (setq const (cons p const))) 2117 (and const 2118 (cons (let ((temp (car const))) 2119 (while (setq const (cdr const)) 2120 (setq temp (list '* (car const) temp))) 2121 temp) 2122 (let ((temp (or (car not-const) 1))) 2123 (while (setq not-const (cdr not-const)) 2124 (setq temp (list '* (car not-const) temp))) 2125 temp))))) 2126 2127(defvar math-sum-int-pow-cache (list '(0 1))) 2128;; Following is from CRC Math Tables, 27th ed, pp. 52-53. 2129(defun math-sum-integer-power (pow) 2130 (let ((calc-prefer-frac t) 2131 (n (length math-sum-int-pow-cache))) 2132 (while (<= n pow) 2133 (let* ((new (list 0 0)) 2134 (lin new) 2135 (pp (cdr (nth (1- n) math-sum-int-pow-cache))) 2136 (p 2) 2137 (sum 0) 2138 q) 2139 (while pp 2140 (setq q (math-div (car pp) p) 2141 new (cons (math-mul q n) new) 2142 sum (math-add sum q) 2143 p (1+ p) 2144 pp (cdr pp))) 2145 (setcar lin (math-sub 1 (math-mul n sum))) 2146 (setq math-sum-int-pow-cache 2147 (nconc math-sum-int-pow-cache (list (nreverse new))) 2148 n (1+ n)))) 2149 (nth pow math-sum-int-pow-cache))) 2150 2151(defun math-to-exponentials (expr) 2152 (and (consp expr) 2153 (= (length expr) 2) 2154 (let ((x (nth 1 expr)) 2155 (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi))) 2156 (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1)))) 2157 (cond ((eq (car expr) 'calcFunc-exp) 2158 (list '^ '(var e var-e) x)) 2159 ((eq (car expr) 'calcFunc-sin) 2160 (or (eq calc-angle-mode 'rad) 2161 (setq x (list '/ (list '* x pi) 180))) 2162 (list '/ (list '- 2163 (list '^ '(var e var-e) (list '* x i)) 2164 (list '^ '(var e var-e) 2165 (list 'neg (list '* x i)))) 2166 (list '* 2 i))) 2167 ((eq (car expr) 'calcFunc-cos) 2168 (or (eq calc-angle-mode 'rad) 2169 (setq x (list '/ (list '* x pi) 180))) 2170 (list '/ (list '+ 2171 (list '^ '(var e var-e) 2172 (list '* x i)) 2173 (list '^ '(var e var-e) 2174 (list 'neg (list '* x i)))) 2175 2)) 2176 ((eq (car expr) 'calcFunc-sinh) 2177 (list '/ (list '- 2178 (list '^ '(var e var-e) x) 2179 (list '^ '(var e var-e) (list 'neg x))) 2180 2)) 2181 ((eq (car expr) 'calcFunc-cosh) 2182 (list '/ (list '+ 2183 (list '^ '(var e var-e) x) 2184 (list '^ '(var e var-e) (list 'neg x))) 2185 2)) 2186 (t nil))))) 2187 2188(defun math-to-exps (expr) 2189 (cond (calc-symbolic-mode expr) 2190 ((Math-primp expr) 2191 (if (equal expr '(var e var-e)) (math-e) expr)) 2192 ((and (eq (car expr) '^) 2193 (equal (nth 1 expr) '(var e var-e))) 2194 (list 'calcFunc-exp (nth 2 expr))) 2195 (t 2196 (cons (car expr) (mapcar 'math-to-exps (cdr expr)))))) 2197 2198 2199(defvar math-disable-prods nil) 2200(defun calcFunc-prod (expr var &optional low high step) 2201 (if math-disable-prods (math-reject-arg)) 2202 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2))) 2203 (math-prod-rec expr var low high step))) 2204 (math-disable-prods t)) 2205 (math-normalize res))) 2206 2207(defun math-prod-rec (expr var &optional low high step) 2208 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) 2209 (and low (not high) (setq high '(var inf var-inf))) 2210 (let (t1 t2 t3 val) 2211 (setq val 2212 (cond 2213 ((not (math-expr-contains expr var)) 2214 (math-pow expr (math-add (math-div (math-sub high low) (or step 1)) 2215 1))) 2216 ((and step (not (math-equal-int step 1))) 2217 (if (math-negp step) 2218 (math-prod-rec expr var high low (math-neg step)) 2219 (let ((lo (math-simplify (math-div low step)))) 2220 (if (math-known-num-integerp lo) 2221 (math-prod-rec (math-normalize 2222 (math-expr-subst expr var 2223 (math-mul step var))) 2224 var lo (math-simplify (math-div high step))) 2225 (math-prod-rec (math-normalize 2226 (math-expr-subst expr var 2227 (math-add (math-mul step 2228 var) 2229 low))) 2230 var 0 2231 (math-simplify (math-div (math-sub high low) 2232 step))))))) 2233 ((and (memq (car expr) '(* /)) 2234 (setq t1 (math-prod-rec (nth 1 expr) var low high) 2235 t2 (math-prod-rec (nth 2 expr) var low high)) 2236 (not (and (math-expr-calls t1 '(calcFunc-prod)) 2237 (math-expr-calls t2 '(calcFunc-prod))))) 2238 (list (car expr) t1 t2)) 2239 ((and (eq (car expr) '^) 2240 (not (math-expr-contains (nth 2 expr) var))) 2241 (math-pow (math-prod-rec (nth 1 expr) var low high) 2242 (nth 2 expr))) 2243 ((and (eq (car expr) '^) 2244 (not (math-expr-contains (nth 1 expr) var))) 2245 (math-pow (nth 1 expr) 2246 (calcFunc-sum (nth 2 expr) var low high))) 2247 ((eq (car expr) 'sqrt) 2248 (math-normalize (list 'calcFunc-sqrt 2249 (list 'calcFunc-prod (nth 1 expr) 2250 var low high)))) 2251 ((eq (car expr) 'neg) 2252 (math-mul (math-pow -1 (math-add (math-sub high low) 1)) 2253 (math-prod-rec (nth 1 expr) var low high))) 2254 ((eq (car expr) 'calcFunc-exp) 2255 (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high))) 2256 ((and (setq t1 (math-is-polynomial expr var 1)) 2257 (setq t2 2258 (cond 2259 ((or (and (math-equal-int (nth 1 t1) 1) 2260 (setq low (math-simplify 2261 (math-add low (car t1))) 2262 high (math-simplify 2263 (math-add high (car t1))))) 2264 (and (math-equal-int (nth 1 t1) -1) 2265 (setq t2 low 2266 low (math-simplify 2267 (math-sub (car t1) high)) 2268 high (math-simplify 2269 (math-sub (car t1) t2))))) 2270 (if (or (math-zerop low) (math-zerop high)) 2271 0 2272 (if (and (or (math-negp low) (math-negp high)) 2273 (or (math-num-integerp low) 2274 (math-num-integerp high))) 2275 (if (math-posp high) 2276 0 2277 (math-mul (math-pow -1 2278 (math-add 2279 (math-add low high) 1)) 2280 (list '/ 2281 (list 'calcFunc-fact 2282 (math-neg low)) 2283 (list 'calcFunc-fact 2284 (math-sub -1 high))))) 2285 (list '/ 2286 (list 'calcFunc-fact high) 2287 (list 'calcFunc-fact (math-sub low 1)))))) 2288 ((and (or (and (math-equal-int (nth 1 t1) 2) 2289 (setq t2 (math-simplify 2290 (math-add (math-mul low 2) 2291 (car t1))) 2292 t3 (math-simplify 2293 (math-add (math-mul high 2) 2294 (car t1))))) 2295 (and (math-equal-int (nth 1 t1) -2) 2296 (setq t2 (math-simplify 2297 (math-sub (car t1) 2298 (math-mul high 2))) 2299 t3 (math-simplify 2300 (math-sub (car t1) 2301 (math-mul low 2302 2)))))) 2303 (or (math-integerp t2) 2304 (and (math-messy-integerp t2) 2305 (setq t2 (math-trunc t2))) 2306 (math-integerp t3) 2307 (and (math-messy-integerp t3) 2308 (setq t3 (math-trunc t3))))) 2309 (if (or (math-zerop t2) (math-zerop t3)) 2310 0 2311 (if (or (math-evenp t2) (math-evenp t3)) 2312 (if (or (math-negp t2) (math-negp t3)) 2313 (if (math-posp high) 2314 0 2315 (list '/ 2316 (list 'calcFunc-dfact 2317 (math-neg t2)) 2318 (list 'calcFunc-dfact 2319 (math-sub -2 t3)))) 2320 (list '/ 2321 (list 'calcFunc-dfact t3) 2322 (list 'calcFunc-dfact 2323 (math-sub t2 2)))) 2324 (if (math-negp t3) 2325 (list '* 2326 (list '^ -1 2327 (list '/ (list '- (list '- t2 t3) 2328 2) 2329 2)) 2330 (list '/ 2331 (list 'calcFunc-dfact 2332 (math-neg t2)) 2333 (list 'calcFunc-dfact 2334 (math-sub -2 t3)))) 2335 (if (math-posp t2) 2336 (list '/ 2337 (list 'calcFunc-dfact t3) 2338 (list 'calcFunc-dfact 2339 (math-sub t2 2))) 2340 nil)))))))) 2341 t2))) 2342 (if (equal val '(var nan var-nan)) (setq val nil)) 2343 (or val 2344 (let* ((math-tabulate-initial 1) 2345 (math-tabulate-function 'calcFunc-prod)) 2346 (calcFunc-table expr var low high))))) 2347 2348 2349 2350 2351(defvar math-solve-ranges nil) 2352(defvar math-solve-sign) 2353;;; Attempt to reduce math-solve-lhs = math-solve-rhs to 2354;;; math-solve-var = math-solve-rhs', where math-solve-var appears 2355;;; in math-solve-lhs but not in math-solve-rhs or math-solve-rhs'; 2356;;; return math-solve-rhs'. 2357;;; Uses global values: math-solve-var, math-solve-full. 2358(defvar math-solve-var) 2359(defvar math-solve-full) 2360 2361;; The variables math-solve-lhs, math-solve-rhs and math-try-solve-sign 2362;; are local to math-try-solve-for, but are used by math-try-solve-prod. 2363;; (math-solve-lhs and math-solve-rhs are is also local to 2364;; math-decompose-poly, but used by math-solve-poly-funny-powers.) 2365(defvar math-solve-lhs) 2366(defvar math-solve-rhs) 2367(defvar math-try-solve-sign) 2368 2369(defun math-try-solve-for 2370 (math-solve-lhs math-solve-rhs &optional math-try-solve-sign no-poly) 2371 (let (math-t1 math-t2 math-t3) 2372 (cond ((equal math-solve-lhs math-solve-var) 2373 (setq math-solve-sign math-try-solve-sign) 2374 (if (eq math-solve-full 'all) 2375 (let ((vec (list 'vec (math-evaluate-expr math-solve-rhs))) 2376 newvec var p) 2377 (while math-solve-ranges 2378 (setq p (car math-solve-ranges) 2379 var (car p) 2380 newvec (list 'vec)) 2381 (while (setq p (cdr p)) 2382 (setq newvec (nconc newvec 2383 (cdr (math-expr-subst 2384 vec var (car p)))))) 2385 (setq vec newvec 2386 math-solve-ranges (cdr math-solve-ranges))) 2387 (math-normalize vec)) 2388 math-solve-rhs)) 2389 ((Math-primp math-solve-lhs) 2390 nil) 2391 ((and (eq (car math-solve-lhs) '-) 2392 (eq (car-safe (nth 1 math-solve-lhs)) (car-safe (nth 2 math-solve-lhs))) 2393 (Math-zerop math-solve-rhs) 2394 (= (length (nth 1 math-solve-lhs)) 2) 2395 (= (length (nth 2 math-solve-lhs)) 2) 2396 (setq math-t1 (get (car (nth 1 math-solve-lhs)) 'math-inverse)) 2397 (setq math-t2 (funcall math-t1 '(var SOLVEDUM SOLVEDUM))) 2398 (eq (math-expr-contains-count math-t2 '(var SOLVEDUM SOLVEDUM)) 1) 2399 (setq math-t3 (math-solve-above-dummy math-t2)) 2400 (setq math-t1 (math-try-solve-for 2401 (math-sub (nth 1 (nth 1 math-solve-lhs)) 2402 (math-expr-subst 2403 math-t2 math-t3 2404 (nth 1 (nth 2 math-solve-lhs)))) 2405 0))) 2406 math-t1) 2407 ((eq (car math-solve-lhs) 'neg) 2408 (math-try-solve-for (nth 1 math-solve-lhs) (math-neg math-solve-rhs) 2409 (and math-try-solve-sign (- math-try-solve-sign)))) 2410 ((and (not (eq math-solve-full 't)) (math-try-solve-prod))) 2411 ((and (not no-poly) 2412 (setq math-t2 2413 (math-decompose-poly math-solve-lhs 2414 math-solve-var 15 math-solve-rhs))) 2415 (setq math-t1 (cdr (nth 1 math-t2)) 2416 math-t1 (let ((math-solve-ranges math-solve-ranges)) 2417 (cond ((= (length math-t1) 5) 2418 (apply 'math-solve-quartic (car math-t2) math-t1)) 2419 ((= (length math-t1) 4) 2420 (apply 'math-solve-cubic (car math-t2) math-t1)) 2421 ((= (length math-t1) 3) 2422 (apply 'math-solve-quadratic (car math-t2) math-t1)) 2423 ((= (length math-t1) 2) 2424 (apply 'math-solve-linear 2425 (car math-t2) math-try-solve-sign math-t1)) 2426 (math-solve-full 2427 (math-poly-all-roots (car math-t2) math-t1)) 2428 (calc-symbolic-mode nil) 2429 (t 2430 (math-try-solve-for 2431 (car math-t2) 2432 (math-poly-any-root (reverse math-t1) 0 t) 2433 nil t))))) 2434 (if math-t1 2435 (if (eq (nth 2 math-t2) 1) 2436 math-t1 2437 (math-solve-prod math-t1 (math-try-solve-for (nth 2 math-t2) 0 nil t))) 2438 (calc-record-why "*Unable to find a symbolic solution") 2439 nil)) 2440 ((and (math-solve-find-root-term math-solve-lhs nil) 2441 (eq (math-expr-contains-count math-solve-lhs math-t1) 1)) ; just in case 2442 (math-try-solve-for (math-simplify 2443 (math-sub (if (or math-t3 (math-evenp math-t2)) 2444 (math-pow math-t1 math-t2) 2445 (math-neg (math-pow math-t1 math-t2))) 2446 (math-expand-power 2447 (math-sub (math-normalize 2448 (math-expr-subst 2449 math-solve-lhs math-t1 0)) 2450 math-solve-rhs) 2451 math-t2 math-solve-var))) 2452 0)) 2453 ((eq (car math-solve-lhs) '+) 2454 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) 2455 (math-try-solve-for (nth 2 math-solve-lhs) 2456 (math-sub math-solve-rhs (nth 1 math-solve-lhs)) 2457 math-try-solve-sign)) 2458 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) 2459 (math-try-solve-for (nth 1 math-solve-lhs) 2460 (math-sub math-solve-rhs (nth 2 math-solve-lhs)) 2461 math-try-solve-sign)))) 2462 ((eq (car math-solve-lhs) 'calcFunc-eq) 2463 (math-try-solve-for (math-sub (nth 1 math-solve-lhs) (nth 2 math-solve-lhs)) 2464 math-solve-rhs math-try-solve-sign no-poly)) 2465 ((eq (car math-solve-lhs) '-) 2466 (cond ((or (and (eq (car-safe (nth 1 math-solve-lhs)) 'calcFunc-sin) 2467 (eq (car-safe (nth 2 math-solve-lhs)) 'calcFunc-cos)) 2468 (and (eq (car-safe (nth 1 math-solve-lhs)) 'calcFunc-cos) 2469 (eq (car-safe (nth 2 math-solve-lhs)) 'calcFunc-sin))) 2470 (math-try-solve-for (math-sub (nth 1 math-solve-lhs) 2471 (list (car (nth 1 math-solve-lhs)) 2472 (math-sub 2473 (math-quarter-circle t) 2474 (nth 1 (nth 2 math-solve-lhs))))) 2475 math-solve-rhs)) 2476 ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) 2477 (math-try-solve-for (nth 2 math-solve-lhs) 2478 (math-sub (nth 1 math-solve-lhs) math-solve-rhs) 2479 (and math-try-solve-sign 2480 (- math-try-solve-sign)))) 2481 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) 2482 (math-try-solve-for (nth 1 math-solve-lhs) 2483 (math-add math-solve-rhs (nth 2 math-solve-lhs)) 2484 math-try-solve-sign)))) 2485 ((and (eq math-solve-full 't) (math-try-solve-prod))) 2486 ((and (eq (car math-solve-lhs) '%) 2487 (not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))) 2488 (math-try-solve-for (nth 1 math-solve-lhs) (math-add math-solve-rhs 2489 (math-solve-get-int 2490 (nth 2 math-solve-lhs))))) 2491 ((eq (car math-solve-lhs) 'calcFunc-log) 2492 (cond ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) 2493 (math-try-solve-for (nth 1 math-solve-lhs) 2494 (math-pow (nth 2 math-solve-lhs) math-solve-rhs))) 2495 ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) 2496 (math-try-solve-for (nth 2 math-solve-lhs) (math-pow 2497 (nth 1 math-solve-lhs) 2498 (math-div 1 math-solve-rhs)))))) 2499 ((and (= (length math-solve-lhs) 2) 2500 (symbolp (car math-solve-lhs)) 2501 (setq math-t1 (get (car math-solve-lhs) 'math-inverse)) 2502 (setq math-t2 (funcall math-t1 math-solve-rhs))) 2503 (setq math-t1 (get (car math-solve-lhs) 'math-inverse-sign)) 2504 (math-try-solve-for (nth 1 math-solve-lhs) (math-normalize math-t2) 2505 (and math-try-solve-sign math-t1 2506 (if (integerp math-t1) 2507 (* math-t1 math-try-solve-sign) 2508 (funcall math-t1 math-solve-lhs 2509 math-try-solve-sign))))) 2510 ((and (symbolp (car math-solve-lhs)) 2511 (setq math-t1 (get (car math-solve-lhs) 'math-inverse-n)) 2512 (setq math-t2 (funcall math-t1 math-solve-lhs math-solve-rhs))) 2513 math-t2) 2514 ((setq math-t1 (math-expand-formula math-solve-lhs)) 2515 (math-try-solve-for math-t1 math-solve-rhs math-try-solve-sign)) 2516 (t 2517 (calc-record-why "*No inverse known" math-solve-lhs) 2518 nil)))) 2519 2520 2521(defun math-try-solve-prod () 2522 (cond ((eq (car math-solve-lhs) '*) 2523 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) 2524 (math-try-solve-for (nth 2 math-solve-lhs) 2525 (math-div math-solve-rhs (nth 1 math-solve-lhs)) 2526 (math-solve-sign math-try-solve-sign 2527 (nth 1 math-solve-lhs)))) 2528 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) 2529 (math-try-solve-for (nth 1 math-solve-lhs) 2530 (math-div math-solve-rhs (nth 2 math-solve-lhs)) 2531 (math-solve-sign math-try-solve-sign 2532 (nth 2 math-solve-lhs)))) 2533 ((Math-zerop math-solve-rhs) 2534 (math-solve-prod (let ((math-solve-ranges math-solve-ranges)) 2535 (math-try-solve-for (nth 2 math-solve-lhs) 0)) 2536 (math-try-solve-for (nth 1 math-solve-lhs) 0))))) 2537 ((eq (car math-solve-lhs) '/) 2538 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) 2539 (math-try-solve-for (nth 2 math-solve-lhs) 2540 (math-div (nth 1 math-solve-lhs) math-solve-rhs) 2541 (math-solve-sign math-try-solve-sign 2542 (nth 1 math-solve-lhs)))) 2543 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) 2544 (math-try-solve-for (nth 1 math-solve-lhs) 2545 (math-mul math-solve-rhs (nth 2 math-solve-lhs)) 2546 (math-solve-sign math-try-solve-sign 2547 (nth 2 math-solve-lhs)))) 2548 ((setq math-t1 (math-try-solve-for (math-sub (nth 1 math-solve-lhs) 2549 (math-mul (nth 2 math-solve-lhs) 2550 math-solve-rhs)) 2551 0)) 2552 math-t1))) 2553 ((eq (car math-solve-lhs) '^) 2554 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) 2555 (math-try-solve-for 2556 (nth 2 math-solve-lhs) 2557 (math-add (math-normalize 2558 (list 'calcFunc-log math-solve-rhs (nth 1 math-solve-lhs))) 2559 (math-div 2560 (math-mul 2 2561 (math-mul '(var pi var-pi) 2562 (math-solve-get-int 2563 '(var i var-i)))) 2564 (math-normalize 2565 (list 'calcFunc-ln (nth 1 math-solve-lhs))))))) 2566 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) 2567 (cond ((and (integerp (nth 2 math-solve-lhs)) 2568 (>= (nth 2 math-solve-lhs) 2) 2569 (setq math-t1 (math-integer-log2 (nth 2 math-solve-lhs)))) 2570 (setq math-t2 math-solve-rhs) 2571 (if (and (eq math-solve-full t) 2572 (math-known-realp (nth 1 math-solve-lhs))) 2573 (progn 2574 (while (>= (setq math-t1 (1- math-t1)) 0) 2575 (setq math-t2 (list 'calcFunc-sqrt math-t2))) 2576 (setq math-t2 (math-solve-get-sign math-t2))) 2577 (while (>= (setq math-t1 (1- math-t1)) 0) 2578 (setq math-t2 (math-solve-get-sign 2579 (math-normalize 2580 (list 'calcFunc-sqrt math-t2)))))) 2581 (math-try-solve-for 2582 (nth 1 math-solve-lhs) 2583 (math-normalize math-t2))) 2584 ((math-looks-negp (nth 2 math-solve-lhs)) 2585 (math-try-solve-for 2586 (list '^ (nth 1 math-solve-lhs) 2587 (math-neg (nth 2 math-solve-lhs))) 2588 (math-div 1 math-solve-rhs))) 2589 ((and (eq math-solve-full t) 2590 (Math-integerp (nth 2 math-solve-lhs)) 2591 (math-known-realp (nth 1 math-solve-lhs))) 2592 (setq math-t1 (math-normalize 2593 (list 'calcFunc-nroot math-solve-rhs 2594 (nth 2 math-solve-lhs)))) 2595 (if (math-evenp (nth 2 math-solve-lhs)) 2596 (setq math-t1 (math-solve-get-sign math-t1))) 2597 (math-try-solve-for 2598 (nth 1 math-solve-lhs) math-t1 2599 (and math-try-solve-sign 2600 (math-oddp (nth 2 math-solve-lhs)) 2601 (math-solve-sign math-try-solve-sign 2602 (nth 2 math-solve-lhs))))) 2603 (t (math-try-solve-for 2604 (nth 1 math-solve-lhs) 2605 (math-mul 2606 (math-normalize 2607 (list 'calcFunc-exp 2608 (if (Math-realp (nth 2 math-solve-lhs)) 2609 (math-div (math-mul 2610 '(var pi var-pi) 2611 (math-solve-get-int 2612 '(var i var-i) 2613 (and (integerp (nth 2 math-solve-lhs)) 2614 (math-abs 2615 (nth 2 math-solve-lhs))))) 2616 (math-div (nth 2 math-solve-lhs) 2)) 2617 (math-div (math-mul 2618 2 2619 (math-mul 2620 '(var pi var-pi) 2621 (math-solve-get-int 2622 '(var i var-i) 2623 (and (integerp (nth 2 math-solve-lhs)) 2624 (math-abs 2625 (nth 2 math-solve-lhs)))))) 2626 (nth 2 math-solve-lhs))))) 2627 (math-normalize 2628 (list 'calcFunc-nroot 2629 math-solve-rhs 2630 (nth 2 math-solve-lhs)))) 2631 (and math-try-solve-sign 2632 (math-oddp (nth 2 math-solve-lhs)) 2633 (math-solve-sign math-try-solve-sign 2634 (nth 2 math-solve-lhs))))))))) 2635 (t nil))) 2636 2637(defun math-solve-prod (lsoln rsoln) 2638 (cond ((null lsoln) 2639 rsoln) 2640 ((null rsoln) 2641 lsoln) 2642 ((eq math-solve-full 'all) 2643 (cons 'vec (append (cdr lsoln) (cdr rsoln)))) 2644 (math-solve-full 2645 (list 'calcFunc-if 2646 (list 'calcFunc-gt (math-solve-get-sign 1) 0) 2647 lsoln 2648 rsoln)) 2649 (t lsoln))) 2650 2651;;; This deals with negative, fractional, and symbolic powers of "x". 2652;; The variable math-solve-b is local to math-decompose-poly, 2653;; but is used by math-solve-poly-funny-powers. 2654(defvar math-solve-b) 2655 2656(defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2" 2657 (setq math-t1 math-solve-lhs) 2658 (let ((pp math-poly-neg-powers) 2659 fac) 2660 (while pp 2661 (setq fac (math-pow (car pp) (or math-poly-mult-powers 1)) 2662 math-t1 (math-mul math-t1 fac) 2663 math-solve-rhs (math-mul math-solve-rhs fac) 2664 pp (cdr pp)))) 2665 (if sub-rhs (setq math-t1 (math-sub math-t1 math-solve-rhs))) 2666 (let ((math-poly-neg-powers nil)) 2667 (setq math-t2 (math-mul (or math-poly-mult-powers 1) 2668 (let ((calc-prefer-frac t)) 2669 (math-div 1 math-poly-frac-powers))) 2670 math-t1 (math-is-polynomial 2671 (math-simplify (calcFunc-expand math-t1)) math-solve-b 50)))) 2672 2673;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2". 2674(defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3" 2675 (let ((count 0)) 2676 (while (and math-t1 (Math-zerop (car math-t1))) 2677 (setq math-t1 (cdr math-t1) 2678 count (1+ count))) 2679 (and math-t1 2680 (let* ((degree (1- (length math-t1))) 2681 (scale degree)) 2682 (while (and (> scale 1) (= (car math-t3) 1)) 2683 (and (= (% degree scale) 0) 2684 (let ((p math-t1) 2685 (n 0) 2686 (new-t1 nil) 2687 (okay t)) 2688 (while (and p okay) 2689 (if (= (% n scale) 0) 2690 (setq new-t1 (nconc new-t1 (list (car p)))) 2691 (or (Math-zerop (car p)) 2692 (setq okay nil))) 2693 (setq p (cdr p) 2694 n (1+ n))) 2695 (if okay 2696 (setq math-t3 (cons scale (cdr math-t3)) 2697 math-t1 new-t1)))) 2698 (setq scale (1- scale))) 2699 (setq math-t3 (list (math-mul (car math-t3) math-t2) 2700 (math-mul count math-t2))) 2701 (<= (1- (length math-t1)) max-degree))))) 2702 2703(defun calcFunc-poly (expr var &optional degree) 2704 (if degree 2705 (or (natnump degree) (math-reject-arg degree 'fixnatnump)) 2706 (setq degree 50)) 2707 (let ((p (math-is-polynomial expr var degree 'gen))) 2708 (if p 2709 (if (equal p '(0)) 2710 (list 'vec) 2711 (cons 'vec p)) 2712 (math-reject-arg expr "Expected a polynomial")))) 2713 2714(defun calcFunc-gpoly (expr var &optional degree) 2715 (if degree 2716 (or (natnump degree) (math-reject-arg degree 'fixnatnump)) 2717 (setq degree 50)) 2718 (let* ((math-poly-base-variable var) 2719 (d (math-decompose-poly expr var degree nil))) 2720 (if d 2721 (cons 'vec d) 2722 (math-reject-arg expr "Expected a polynomial")))) 2723 2724(defun math-decompose-poly (math-solve-lhs math-solve-var degree sub-rhs) 2725 (let ((math-solve-rhs (or sub-rhs 1)) 2726 math-t1 math-t2 math-t3) 2727 (setq math-t2 (math-polynomial-base 2728 math-solve-lhs 2729 (function 2730 (lambda (math-solve-b) 2731 (let ((math-poly-neg-powers '(1)) 2732 (math-poly-mult-powers nil) 2733 (math-poly-frac-powers 1) 2734 (math-poly-exp-base t)) 2735 (and (not (equal math-solve-b math-solve-lhs)) 2736 (or (not (memq (car-safe math-solve-b) '(+ -))) sub-rhs) 2737 (setq math-t3 '(1 0) math-t2 1 2738 math-t1 (math-is-polynomial math-solve-lhs 2739 math-solve-b 50)) 2740 (if (and (equal math-poly-neg-powers '(1)) 2741 (memq math-poly-mult-powers '(nil 1)) 2742 (eq math-poly-frac-powers 1) 2743 sub-rhs) 2744 (setq math-t1 (cons (math-sub (car math-t1) math-solve-rhs) 2745 (cdr math-t1))) 2746 (math-solve-poly-funny-powers sub-rhs)) 2747 (math-solve-crunch-poly degree) 2748 (or (math-expr-contains math-solve-b math-solve-var) 2749 (math-expr-contains (car math-t3) math-solve-var)))))))) 2750 (if math-t2 2751 (list (math-pow math-t2 (car math-t3)) 2752 (cons 'vec math-t1) 2753 (if sub-rhs 2754 (math-pow math-t2 (nth 1 math-t3)) 2755 (math-div (math-pow math-t2 (nth 1 math-t3)) math-solve-rhs)))))) 2756 2757(defun math-solve-linear (var sign b a) 2758 (math-try-solve-for var 2759 (math-div (math-neg b) a) 2760 (math-solve-sign sign a) 2761 t)) 2762 2763(defun math-solve-quadratic (var c b a) 2764 (math-try-solve-for 2765 var 2766 (if (math-looks-evenp b) 2767 (let ((halfb (math-div b 2))) 2768 (math-div 2769 (math-add 2770 (math-neg halfb) 2771 (math-solve-get-sign 2772 (math-normalize 2773 (list 'calcFunc-sqrt 2774 (math-add (math-sqr halfb) 2775 (math-mul (math-neg c) a)))))) 2776 a)) 2777 (math-div 2778 (math-add 2779 (math-neg b) 2780 (math-solve-get-sign 2781 (math-normalize 2782 (list 'calcFunc-sqrt 2783 (math-add (math-sqr b) 2784 (math-mul 4 (math-mul (math-neg c) a))))))) 2785 (math-mul 2 a))) 2786 nil t)) 2787 2788(defun math-solve-cubic (var d c b a) 2789 (let* ((p (math-div b a)) 2790 (q (math-div c a)) 2791 (r (math-div d a)) 2792 (psqr (math-sqr p)) 2793 (aa (math-sub q (math-div psqr 3))) 2794 (bb (math-add r 2795 (math-div (math-sub (math-mul 2 (math-mul psqr p)) 2796 (math-mul 9 (math-mul p q))) 2797 27))) 2798 m) 2799 (if (Math-zerop aa) 2800 (math-try-solve-for (math-pow (math-add var (math-div p 3)) 3) 2801 (math-neg bb) nil t) 2802 (if (Math-zerop bb) 2803 (math-try-solve-for 2804 (math-mul (math-add var (math-div p 3)) 2805 (math-add (math-sqr (math-add var (math-div p 3))) 2806 aa)) 2807 0 nil t) 2808 (setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3)))) 2809 (math-try-solve-for 2810 var 2811 (math-sub 2812 (math-normalize 2813 (math-mul 2814 m 2815 (list 'calcFunc-cos 2816 (math-div 2817 (math-sub (list 'calcFunc-arccos 2818 (math-div (math-mul 3 bb) 2819 (math-mul aa m))) 2820 (math-mul 2 2821 (math-mul 2822 (math-add 1 (math-solve-get-int 2823 1 3)) 2824 (math-half-circle 2825 calc-symbolic-mode)))) 2826 3)))) 2827 (math-div p 3)) 2828 nil t))))) 2829 2830(defun math-solve-quartic (var d c b a aa) 2831 (setq a (math-div a aa)) 2832 (setq b (math-div b aa)) 2833 (setq c (math-div c aa)) 2834 (setq d (math-div d aa)) 2835 (math-try-solve-for 2836 var 2837 (let* ((asqr (math-sqr a)) 2838 (asqr4 (math-div asqr 4)) 2839 (y (let ((math-solve-full nil) 2840 calc-next-why) 2841 (math-solve-cubic math-solve-var 2842 (math-sub (math-sub 2843 (math-mul 4 (math-mul b d)) 2844 (math-mul asqr d)) 2845 (math-sqr c)) 2846 (math-sub (math-mul a c) 2847 (math-mul 4 d)) 2848 (math-neg b) 2849 1))) 2850 (rsqr (math-add (math-sub asqr4 b) y)) 2851 (r (list 'calcFunc-sqrt rsqr)) 2852 (sign1 (math-solve-get-sign 1)) 2853 (de (list 'calcFunc-sqrt 2854 (math-add 2855 (math-sub (math-mul 3 asqr4) 2856 (math-mul 2 b)) 2857 (if (Math-zerop rsqr) 2858 (math-mul 2859 2 2860 (math-mul sign1 2861 (list 'calcFunc-sqrt 2862 (math-sub (math-sqr y) 2863 (math-mul 4 d))))) 2864 (math-sub 2865 (math-mul sign1 2866 (math-div 2867 (math-sub (math-sub 2868 (math-mul 4 (math-mul a b)) 2869 (math-mul 8 c)) 2870 (math-mul asqr a)) 2871 (math-mul 4 r))) 2872 rsqr)))))) 2873 (math-normalize 2874 (math-sub (math-add (math-mul sign1 (math-div r 2)) 2875 (math-solve-get-sign (math-div de 2))) 2876 (math-div a 4)))) 2877 nil t)) 2878 2879(defvar math-symbolic-solve nil) 2880(defvar math-int-coefs nil) 2881 2882;; The variable math-int-threshold is local to math-poly-all-roots, 2883;; but is used by math-poly-newton-root. 2884(defvar math-int-threshold) 2885;; The variables math-int-scale, math-int-factors and math-double-roots 2886;; are local to math-poly-all-roots, but are used by math-poly-integer-root. 2887(defvar math-int-scale) 2888(defvar math-int-factors) 2889(defvar math-double-roots) 2890 2891(defun math-poly-all-roots (var p &optional math-factoring) 2892 (catch 'ouch 2893 (let* ((math-symbolic-solve calc-symbolic-mode) 2894 (roots nil) 2895 (deg (1- (length p))) 2896 (orig-p (reverse p)) 2897 (math-int-coefs nil) 2898 (math-int-scale nil) 2899 (math-double-roots nil) 2900 (math-int-factors nil) 2901 (math-int-threshold nil) 2902 (pp p)) 2903 ;; If rational coefficients, look for exact rational factors. 2904 (while (and pp (Math-ratp (car pp))) 2905 (setq pp (cdr pp))) 2906 (if pp 2907 (if (or math-factoring math-symbolic-solve) 2908 (throw 'ouch nil)) 2909 (let ((lead (car orig-p)) 2910 (calc-prefer-frac t) 2911 (scale (apply 'math-lcm-denoms p))) 2912 (setq math-int-scale (math-abs (math-mul scale lead)) 2913 math-int-threshold (math-div '(float 5 -2) math-int-scale) 2914 math-int-coefs (cdr (math-div (cons 'vec orig-p) lead))))) 2915 (if (> deg 4) 2916 (let ((calc-prefer-frac nil) 2917 (calc-symbolic-mode nil) 2918 (pp p) 2919 (def-p (copy-sequence orig-p))) 2920 (while pp 2921 (if (Math-numberp (car pp)) 2922 (setq pp (cdr pp)) 2923 (throw 'ouch nil))) 2924 (while (> deg (if math-symbolic-solve 2 4)) 2925 (let* ((x (math-poly-any-root def-p '(float 0 0) nil)) 2926 b c pp) 2927 (if (and (eq (car-safe x) 'cplx) 2928 (math-nearly-zerop (nth 2 x) (nth 1 x))) 2929 (setq x (calcFunc-re x))) 2930 (or math-factoring 2931 (setq roots (cons x roots))) 2932 (or (math-numberp x) 2933 (setq x (math-evaluate-expr x))) 2934 (setq pp def-p 2935 b (car def-p)) 2936 (while (setq pp (cdr pp)) 2937 (setq c (car pp)) 2938 (setcar pp b) 2939 (setq b (math-add (math-mul x b) c))) 2940 (setq def-p (cdr def-p) 2941 deg (1- deg)))) 2942 (setq p (reverse def-p)))) 2943 (if (> deg 1) 2944 (let ((math-solve-var '(var DUMMY var-DUMMY)) 2945 (math-solve-sign nil) 2946 (math-solve-ranges nil) 2947 (math-solve-full 'all)) 2948 (if (= (length p) (length math-int-coefs)) 2949 (setq p (reverse math-int-coefs))) 2950 (setq roots (append (cdr (apply (cond ((= deg 2) 2951 'math-solve-quadratic) 2952 ((= deg 3) 2953 'math-solve-cubic) 2954 (t 2955 'math-solve-quartic)) 2956 math-solve-var p)) 2957 roots))) 2958 (if (> deg 0) 2959 (setq roots (cons (math-div (math-neg (car p)) (nth 1 p)) 2960 roots)))) 2961 (if math-factoring 2962 (progn 2963 (while roots 2964 (math-poly-integer-root (car roots)) 2965 (setq roots (cdr roots))) 2966 (list math-int-factors (nreverse math-int-coefs) math-int-scale)) 2967 (let ((vec nil) res) 2968 (while roots 2969 (let ((root (car roots)) 2970 (math-solve-full (and math-solve-full 'all))) 2971 (if (math-floatp root) 2972 (setq root (math-poly-any-root orig-p root t))) 2973 (setq vec (append vec 2974 (cdr (or (math-try-solve-for var root nil t) 2975 (throw 'ouch nil)))))) 2976 (setq roots (cdr roots))) 2977 (setq vec (cons 'vec (nreverse vec))) 2978 (if math-symbolic-solve 2979 (setq vec (math-normalize vec))) 2980 (if (eq math-solve-full t) 2981 (list 'calcFunc-subscr 2982 vec 2983 (math-solve-get-int 1 (1- (length orig-p)) 1)) 2984 vec)))))) 2985 2986(defun math-lcm-denoms (&rest fracs) 2987 (let ((den 1)) 2988 (while fracs 2989 (if (eq (car-safe (car fracs)) 'frac) 2990 (setq den (calcFunc-lcm den (nth 2 (car fracs))))) 2991 (setq fracs (cdr fracs))) 2992 den)) 2993 2994(defun math-poly-any-root (p x polish) ; p is a reverse poly coeff list 2995 (let* ((newt (if (math-zerop x) 2996 (math-poly-newton-root 2997 p '(cplx (float 123 -6) (float 1 -4)) 4) 2998 (math-poly-newton-root p x 4))) 2999 (res (if (math-zerop (cdr newt)) 3000 (car newt) 3001 (if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish)) 3002 (setq newt (math-poly-newton-root p (car newt) 30))) 3003 (if (math-zerop (cdr newt)) 3004 (car newt) 3005 (math-poly-laguerre-root p x polish))))) 3006 (and math-symbolic-solve (math-floatp res) 3007 (throw 'ouch nil)) 3008 res)) 3009 3010(defun math-poly-newton-root (p x iters) 3011 (let* ((calc-prefer-frac nil) 3012 (calc-symbolic-mode nil) 3013 (try-integer math-int-coefs) 3014 (dx x) b d) 3015 (while (and (> (setq iters (1- iters)) 0) 3016 (let ((pp p)) 3017 (math-working "newton" x) 3018 (setq b (car p) 3019 d 0) 3020 (while (setq pp (cdr pp)) 3021 (setq d (math-add (math-mul x d) b) 3022 b (math-add (math-mul x b) (car pp)))) 3023 (not (math-zerop d))) 3024 (progn 3025 (setq dx (math-div b d) 3026 x (math-sub x dx)) 3027 (if try-integer 3028 (let ((adx (math-abs-approx dx))) 3029 (and (math-lessp adx math-int-threshold) 3030 (let ((iroot (math-poly-integer-root x))) 3031 (if iroot 3032 (setq x iroot dx 0) 3033 (setq try-integer nil)))))) 3034 (or (not (or (eq dx 0) 3035 (math-nearly-zerop dx (math-abs-approx x)))) 3036 (progn (setq dx 0) nil))))) 3037 (cons x (if (math-zerop x) 3038 1 (math-div (math-abs-approx dx) (math-abs-approx x)))))) 3039 3040(defun math-poly-integer-root (x) 3041 (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec) 3042 math-int-coefs 3043 (let* ((calc-prefer-frac t) 3044 (xre (calcFunc-re x)) 3045 (xim (calcFunc-im x)) 3046 (xresq (math-sqr xre)) 3047 (ximsq (math-sqr xim))) 3048 (if (math-lessp ximsq (calcFunc-scf xresq -1)) 3049 ;; Look for linear factor 3050 (let* ((rnd (math-div (math-round (math-mul xre math-int-scale)) 3051 math-int-scale)) 3052 (icp math-int-coefs) 3053 (rem (car icp)) 3054 (newcoef nil)) 3055 (while (setq icp (cdr icp)) 3056 (setq newcoef (cons rem newcoef) 3057 rem (math-add (car icp) 3058 (math-mul rem rnd)))) 3059 (and (math-zerop rem) 3060 (progn 3061 (setq math-int-coefs (nreverse newcoef) 3062 math-int-factors (cons (list (math-neg rnd)) 3063 math-int-factors)) 3064 rnd))) 3065 ;; Look for irreducible quadratic factor 3066 (let* ((rnd1 (math-div (math-round 3067 (math-mul xre (math-mul -2 math-int-scale))) 3068 math-int-scale)) 3069 (sqscale (math-sqr math-int-scale)) 3070 (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq) 3071 sqscale)) 3072 sqscale)) 3073 (rem1 (car math-int-coefs)) 3074 (icp (cdr math-int-coefs)) 3075 (rem0 (car icp)) 3076 (newcoef nil) 3077 (found (assoc (list rnd0 rnd1 (math-posp xim)) 3078 math-double-roots)) 3079 this) 3080 (if found 3081 (setq math-double-roots (delq found math-double-roots) 3082 rem0 0 rem1 0) 3083 (while (setq icp (cdr icp)) 3084 (setq this rem1 3085 newcoef (cons rem1 newcoef) 3086 rem1 (math-sub rem0 (math-mul this rnd1)) 3087 rem0 (math-sub (car icp) (math-mul this rnd0))))) 3088 (and (math-zerop rem0) 3089 (math-zerop rem1) 3090 (let ((aa (math-div rnd1 -2))) 3091 (or found (setq math-int-coefs (reverse newcoef) 3092 math-double-roots (cons (list 3093 (list 3094 rnd0 rnd1 3095 (math-negp xim))) 3096 math-double-roots) 3097 math-int-factors (cons (cons rnd0 rnd1) 3098 math-int-factors))) 3099 (math-add aa 3100 (let ((calc-symbolic-mode math-symbolic-solve)) 3101 (math-mul (math-sqrt (math-sub (math-sqr aa) 3102 rnd0)) 3103 (if (math-negp xim) -1 1))))))))))) 3104 3105;;; The following routine is from Numerical Recipes, section 9.5. 3106(defun math-poly-laguerre-root (p x polish) 3107 (let* ((calc-prefer-frac nil) 3108 (calc-symbolic-mode nil) 3109 (iters 0) 3110 (m (1- (length p))) 3111 (try-newt (not polish)) 3112 (tried-newt nil) 3113 b d f x1 dx dxold) 3114 (while 3115 (and (or (< (setq iters (1+ iters)) 50) 3116 (math-reject-arg x "*Laguerre's method failed to converge")) 3117 (let ((err (math-abs-approx (car p))) 3118 (abx (math-abs-approx x)) 3119 (pp p)) 3120 (setq b (car p) 3121 d 0 f 0) 3122 (while (setq pp (cdr pp)) 3123 (setq f (math-add (math-mul x f) d) 3124 d (math-add (math-mul x d) b) 3125 b (math-add (math-mul x b) (car pp)) 3126 err (math-add (math-abs-approx b) (math-mul abx err)))) 3127 (math-lessp (calcFunc-scf err (- -2 calc-internal-prec)) 3128 (math-abs-approx b))) 3129 (or (not (math-zerop d)) 3130 (not (math-zerop f)) 3131 (progn 3132 (setq x (math-pow (math-neg b) (list 'frac 1 m))) 3133 nil)) 3134 (let* ((g (math-div d b)) 3135 (g2 (math-sqr g)) 3136 (h (math-sub g2 (math-mul 2 (math-div f b)))) 3137 (sq (math-sqrt 3138 (math-mul (1- m) (math-sub (math-mul m h) g2)))) 3139 (gp (math-add g sq)) 3140 (gm (math-sub g sq))) 3141 (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm)) 3142 (setq gp gm)) 3143 (setq dx (math-div m gp) 3144 x1 (math-sub x dx)) 3145 (if (and try-newt 3146 (math-lessp (math-abs-approx dx) 3147 (calcFunc-scf (math-abs-approx x) -3))) 3148 (let ((newt (math-poly-newton-root p x1 7))) 3149 (setq tried-newt t 3150 try-newt nil) 3151 (if (math-zerop (cdr newt)) 3152 (setq x (car newt) x1 x) 3153 (if (math-lessp (cdr newt) '(float 1 -6)) 3154 (let ((newt2 (math-poly-newton-root 3155 p (car newt) 20))) 3156 (if (math-zerop (cdr newt2)) 3157 (setq x (car newt2) x1 x) 3158 (setq x (car newt)))))))) 3159 (not (or (eq x x1) 3160 (math-nearly-equal x x1)))) 3161 (let ((cdx (math-abs-approx dx))) 3162 (setq x x1 3163 tried-newt nil) 3164 (prog1 3165 (or (<= iters 6) 3166 (math-lessp cdx dxold) 3167 (progn 3168 (if polish 3169 (let ((digs (calcFunc-xpon 3170 (math-div (math-abs-approx x) cdx)))) 3171 (calc-record-why 3172 "*Could not attain full precision") 3173 (if (natnump digs) 3174 (let ((calc-internal-prec (max 3 digs))) 3175 (setq x (math-normalize x)))))) 3176 nil)) 3177 (setq dxold cdx))) 3178 (or polish 3179 (math-lessp (calcFunc-scf (math-abs-approx x) 3180 (- calc-internal-prec)) 3181 dxold)))) 3182 (or (and (math-floatp x) 3183 (math-poly-integer-root x)) 3184 x))) 3185 3186(defun math-solve-above-dummy (x) 3187 (and (not (Math-primp x)) 3188 (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM)) 3189 (= (length x) 2)) 3190 x 3191 (let ((res nil)) 3192 (while (and (setq x (cdr x)) 3193 (not (setq res (math-solve-above-dummy (car x)))))) 3194 res)))) 3195 3196(defun math-solve-find-root-term (x neg) ; sets "t2", "t3" 3197 (if (math-solve-find-root-in-prod x) 3198 (setq math-t3 neg 3199 math-t1 x) 3200 (and (memq (car-safe x) '(+ -)) 3201 (or (math-solve-find-root-term (nth 1 x) neg) 3202 (math-solve-find-root-term (nth 2 x) 3203 (if (eq (car x) '-) (not neg) neg)))))) 3204 3205(defun math-solve-find-root-in-prod (x) 3206 (and (consp x) 3207 (math-expr-contains x math-solve-var) 3208 (or (and (eq (car x) 'calcFunc-sqrt) 3209 (setq math-t2 2)) 3210 (and (eq (car x) '^) 3211 (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3)) 3212 (setq math-t2 2)) 3213 (and (eq (car-safe (nth 2 x)) 'frac) 3214 (eq (nth 2 (nth 2 x)) 3) 3215 (setq math-t2 3)))) 3216 (and (memq (car x) '(* /)) 3217 (or (and (not (math-expr-contains (nth 1 x) math-solve-var)) 3218 (math-solve-find-root-in-prod (nth 2 x))) 3219 (and (not (math-expr-contains (nth 2 x) math-solve-var)) 3220 (math-solve-find-root-in-prod (nth 1 x)))))))) 3221 3222;; The variable math-solve-vars is local to math-solve-system, 3223;; but is used by math-solve-system-rec. 3224(defvar math-solve-vars) 3225 3226;; The variable math-solve-simplifying is local to math-solve-system 3227;; and math-solve-system-rec, but is used by math-solve-system-subst. 3228(defvar math-solve-simplifying) 3229 3230(defun math-solve-system (exprs math-solve-vars math-solve-full) 3231 (setq exprs (mapcar 'list (if (Math-vectorp exprs) 3232 (cdr exprs) 3233 (list exprs))) 3234 math-solve-vars (if (Math-vectorp math-solve-vars) 3235 (cdr math-solve-vars) 3236 (list math-solve-vars))) 3237 (or (let ((math-solve-simplifying nil)) 3238 (math-solve-system-rec exprs math-solve-vars nil)) 3239 (let ((math-solve-simplifying t)) 3240 (math-solve-system-rec exprs math-solve-vars nil)))) 3241 3242;;; The following backtracking solver works by choosing a variable 3243;;; and equation, and trying to solve the equation for the variable. 3244;;; If it succeeds it calls itself recursively with that variable and 3245;;; equation removed from their respective lists, and with the solution 3246;;; added to solns as well as being substituted into all existing 3247;;; equations. The algorithm terminates when any solution path 3248;;; manages to remove all the variables from var-list. 3249 3250;;; To support calcFunc-roots, entries in eqn-list and solns are 3251;;; actually lists of equations. 3252 3253;; The variables math-solve-system-res and math-solve-system-vv are 3254;; local to math-solve-system-rec, but are used by math-solve-system-subst. 3255(defvar math-solve-system-vv) 3256(defvar math-solve-system-res) 3257 3258 3259(defun math-solve-system-rec (eqn-list var-list solns) 3260 (if var-list 3261 (let ((v var-list) 3262 (math-solve-system-res nil)) 3263 3264 ;; Try each variable in turn. 3265 (while 3266 (and 3267 v 3268 (let* ((math-solve-system-vv (car v)) 3269 (e eqn-list) 3270 (elim (eq (car-safe math-solve-system-vv) 'calcFunc-elim))) 3271 (if elim 3272 (setq math-solve-system-vv (nth 1 math-solve-system-vv))) 3273 3274 ;; Try each equation in turn. 3275 (while 3276 (and 3277 e 3278 (let ((e2 (car e)) 3279 (eprev nil) 3280 res2) 3281 (setq math-solve-system-res nil) 3282 3283 ;; Try to solve for math-solve-system-vv the list of equations e2. 3284 (while (and e2 3285 (setq res2 (or (and (eq (car e2) eprev) 3286 res2) 3287 (math-solve-for (car e2) 0 3288 math-solve-system-vv 3289 math-solve-full)))) 3290 (setq eprev (car e2) 3291 math-solve-system-res (cons (if (eq math-solve-full 'all) 3292 (cdr res2) 3293 (list res2)) 3294 math-solve-system-res) 3295 e2 (cdr e2))) 3296 (if e2 3297 (setq math-solve-system-res nil) 3298 3299 ;; Found a solution. Now try other variables. 3300 (setq math-solve-system-res (nreverse math-solve-system-res) 3301 math-solve-system-res (math-solve-system-rec 3302 (mapcar 3303 'math-solve-system-subst 3304 (delq (car e) 3305 (copy-sequence eqn-list))) 3306 (delq (car v) (copy-sequence var-list)) 3307 (let ((math-solve-simplifying nil) 3308 (s (mapcar 3309 (function 3310 (lambda (x) 3311 (cons 3312 (car x) 3313 (math-solve-system-subst 3314 (cdr x))))) 3315 solns))) 3316 (if elim 3317 s 3318 (cons (cons 3319 math-solve-system-vv 3320 (apply 'append math-solve-system-res)) 3321 s))))) 3322 (not math-solve-system-res)))) 3323 (setq e (cdr e))) 3324 (not math-solve-system-res))) 3325 (setq v (cdr v))) 3326 math-solve-system-res) 3327 3328 ;; Eliminated all variables, so now put solution into the proper format. 3329 (setq solns (sort solns 3330 (function 3331 (lambda (x y) 3332 (not (memq (car x) (memq (car y) math-solve-vars))))))) 3333 (if (eq math-solve-full 'all) 3334 (math-transpose 3335 (math-normalize 3336 (cons 'vec 3337 (if solns 3338 (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns) 3339 (mapcar (function (lambda (x) (cons 'vec x))) eqn-list))))) 3340 (math-normalize 3341 (cons 'vec 3342 (if solns 3343 (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns) 3344 (mapcar 'car eqn-list))))))) 3345 3346(defun math-solve-system-subst (x) ; uses "res" and "v" 3347 (let ((accum nil) 3348 (res2 math-solve-system-res)) 3349 (while x 3350 (setq accum (nconc accum 3351 (mapcar (function 3352 (lambda (r) 3353 (if math-solve-simplifying 3354 (math-simplify 3355 (math-expr-subst 3356 (car x) math-solve-system-vv r)) 3357 (math-expr-subst 3358 (car x) math-solve-system-vv r)))) 3359 (car res2))) 3360 x (cdr x) 3361 res2 (cdr res2))) 3362 accum)) 3363 3364 3365;; calc-command-flags is declared in calc.el 3366(defvar calc-command-flags) 3367 3368(defun math-get-from-counter (name) 3369 (let ((ctr (assq name calc-command-flags))) 3370 (if ctr 3371 (setcdr ctr (1+ (cdr ctr))) 3372 (setq ctr (cons name 1) 3373 calc-command-flags (cons ctr calc-command-flags))) 3374 (cdr ctr))) 3375 3376(defvar var-GenCount) 3377 3378(defun math-solve-get-sign (val) 3379 (setq val (math-simplify val)) 3380 (if (and (eq (car-safe val) '*) 3381 (Math-numberp (nth 1 val))) 3382 (list '* (nth 1 val) (math-solve-get-sign (nth 2 val))) 3383 (and (eq (car-safe val) 'calcFunc-sqrt) 3384 (eq (car-safe (nth 1 val)) '^) 3385 (setq val (math-normalize (list '^ 3386 (nth 1 (nth 1 val)) 3387 (math-div (nth 2 (nth 1 val)) 2))))) 3388 (if math-solve-full 3389 (if (and (calc-var-value 'var-GenCount) 3390 (Math-natnump var-GenCount) 3391 (not (eq math-solve-full 'all))) 3392 (prog1 3393 (math-mul (list 'calcFunc-as var-GenCount) val) 3394 (setq var-GenCount (math-add var-GenCount 1)) 3395 (calc-refresh-evaltos 'var-GenCount)) 3396 (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign)))) 3397 (var2 (list 'var (intern var) (intern (concat "var-" var))))) 3398 (if (eq math-solve-full 'all) 3399 (setq math-solve-ranges (cons (list var2 1 -1) 3400 math-solve-ranges))) 3401 (math-mul var2 val))) 3402 (calc-record-why "*Choosing positive solution") 3403 val))) 3404 3405(defun math-solve-get-int (val &optional range first) 3406 (if math-solve-full 3407 (if (and (calc-var-value 'var-GenCount) 3408 (Math-natnump var-GenCount) 3409 (not (eq math-solve-full 'all))) 3410 (prog1 3411 (math-mul val (list 'calcFunc-an var-GenCount)) 3412 (setq var-GenCount (math-add var-GenCount 1)) 3413 (calc-refresh-evaltos 'var-GenCount)) 3414 (let* ((var (concat "n" (int-to-string 3415 (math-get-from-counter 'solve-int)))) 3416 (var2 (list 'var (intern var) (intern (concat "var-" var))))) 3417 (if (and range (eq math-solve-full 'all)) 3418 (setq math-solve-ranges (cons (cons var2 3419 (cdr (calcFunc-index 3420 range (or first 0)))) 3421 math-solve-ranges))) 3422 (math-mul val var2))) 3423 (calc-record-why "*Choosing 0 for arbitrary integer in solution") 3424 0)) 3425 3426(defun math-solve-sign (sign expr) 3427 (and sign 3428 (let ((s1 (math-possible-signs expr))) 3429 (cond ((memq s1 '(4 6)) 3430 sign) 3431 ((memq s1 '(1 3)) 3432 (- sign)))))) 3433 3434(defun math-looks-evenp (expr) 3435 (if (Math-integerp expr) 3436 (math-evenp expr) 3437 (if (memq (car expr) '(* /)) 3438 (math-looks-evenp (nth 1 expr))))) 3439 3440(defun math-solve-for (lhs rhs math-solve-var math-solve-full &optional sign) 3441 (if (math-expr-contains rhs math-solve-var) 3442 (math-solve-for (math-sub lhs rhs) 0 math-solve-var math-solve-full) 3443 (and (math-expr-contains lhs math-solve-var) 3444 (math-with-extra-prec 1 3445 (let* ((math-poly-base-variable math-solve-var) 3446 (res (math-try-solve-for lhs rhs sign))) 3447 (if (and (eq math-solve-full 'all) 3448 (math-known-realp math-solve-var)) 3449 (let ((old-len (length res)) 3450 new-len) 3451 (setq res (delq nil 3452 (mapcar (function 3453 (lambda (x) 3454 (and (not (memq (car-safe x) 3455 '(cplx polar))) 3456 x))) 3457 res)) 3458 new-len (length res)) 3459 (if (< new-len old-len) 3460 (calc-record-why (if (= new-len 1) 3461 "*All solutions were complex" 3462 (format 3463 "*Omitted %d complex solutions" 3464 (- old-len new-len))))))) 3465 res))))) 3466 3467(defun math-solve-eqn (expr var full) 3468 (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt 3469 calcFunc-leq calcFunc-geq)) 3470 (let ((res (math-solve-for (cons '- (cdr expr)) 3471 0 var full 3472 (if (eq (car expr) 'calcFunc-neq) nil 1)))) 3473 (and res 3474 (if (eq math-solve-sign 1) 3475 (list (car expr) var res) 3476 (if (eq math-solve-sign -1) 3477 (list (car expr) res var) 3478 (or (eq (car expr) 'calcFunc-neq) 3479 (calc-record-why 3480 "*Can't determine direction of inequality")) 3481 (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt)) 3482 (list 'calcFunc-neq var res)))))) 3483 (let ((res (math-solve-for expr 0 var full))) 3484 (and res 3485 (list 'calcFunc-eq var res))))) 3486 3487(defun math-reject-solution (expr var func) 3488 (if (math-expr-contains expr var) 3489 (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution")) 3490 (calc-record-why "*Unable to find a solution"))) 3491 (list func expr var)) 3492 3493(defun calcFunc-solve (expr var) 3494 (or (if (or (Math-vectorp expr) (Math-vectorp var)) 3495 (math-solve-system expr var nil) 3496 (math-solve-eqn expr var nil)) 3497 (math-reject-solution expr var 'calcFunc-solve))) 3498 3499(defun calcFunc-fsolve (expr var) 3500 (or (if (or (Math-vectorp expr) (Math-vectorp var)) 3501 (math-solve-system expr var t) 3502 (math-solve-eqn expr var t)) 3503 (math-reject-solution expr var 'calcFunc-fsolve))) 3504 3505(defun calcFunc-roots (expr var) 3506 (let ((math-solve-ranges nil)) 3507 (or (if (or (Math-vectorp expr) (Math-vectorp var)) 3508 (math-solve-system expr var 'all) 3509 (math-solve-for expr 0 var 'all)) 3510 (math-reject-solution expr var 'calcFunc-roots)))) 3511 3512(defun calcFunc-finv (expr var) 3513 (let ((res (math-solve-for expr math-integ-var var nil))) 3514 (if res 3515 (math-normalize (math-expr-subst res math-integ-var var)) 3516 (math-reject-solution expr var 'calcFunc-finv)))) 3517 3518(defun calcFunc-ffinv (expr var) 3519 (let ((res (math-solve-for expr math-integ-var var t))) 3520 (if res 3521 (math-normalize (math-expr-subst res math-integ-var var)) 3522 (math-reject-solution expr var 'calcFunc-finv)))) 3523 3524 3525(put 'calcFunc-inv 'math-inverse 3526 (function (lambda (x) (math-div 1 x)))) 3527(put 'calcFunc-inv 'math-inverse-sign -1) 3528 3529(put 'calcFunc-sqrt 'math-inverse 3530 (function (lambda (x) (math-sqr x)))) 3531 3532(put 'calcFunc-conj 'math-inverse 3533 (function (lambda (x) (list 'calcFunc-conj x)))) 3534 3535(put 'calcFunc-abs 'math-inverse 3536 (function (lambda (x) (math-solve-get-sign x)))) 3537 3538(put 'calcFunc-deg 'math-inverse 3539 (function (lambda (x) (list 'calcFunc-rad x)))) 3540(put 'calcFunc-deg 'math-inverse-sign 1) 3541 3542(put 'calcFunc-rad 'math-inverse 3543 (function (lambda (x) (list 'calcFunc-deg x)))) 3544(put 'calcFunc-rad 'math-inverse-sign 1) 3545 3546(put 'calcFunc-ln 'math-inverse 3547 (function (lambda (x) (list 'calcFunc-exp x)))) 3548(put 'calcFunc-ln 'math-inverse-sign 1) 3549 3550(put 'calcFunc-log10 'math-inverse 3551 (function (lambda (x) (list 'calcFunc-exp10 x)))) 3552(put 'calcFunc-log10 'math-inverse-sign 1) 3553 3554(put 'calcFunc-lnp1 'math-inverse 3555 (function (lambda (x) (list 'calcFunc-expm1 x)))) 3556(put 'calcFunc-lnp1 'math-inverse-sign 1) 3557 3558(put 'calcFunc-exp 'math-inverse 3559 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x)) 3560 (math-mul 2 3561 (math-mul '(var pi var-pi) 3562 (math-solve-get-int 3563 '(var i var-i)))))))) 3564(put 'calcFunc-exp 'math-inverse-sign 1) 3565 3566(put 'calcFunc-expm1 'math-inverse 3567 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x)) 3568 (math-mul 2 3569 (math-mul '(var pi var-pi) 3570 (math-solve-get-int 3571 '(var i var-i)))))))) 3572(put 'calcFunc-expm1 'math-inverse-sign 1) 3573 3574(put 'calcFunc-sin 'math-inverse 3575 (function (lambda (x) (let ((n (math-solve-get-int 1))) 3576 (math-add (math-mul (math-normalize 3577 (list 'calcFunc-arcsin x)) 3578 (math-pow -1 n)) 3579 (math-mul (math-half-circle t) 3580 n)))))) 3581 3582(put 'calcFunc-cos 'math-inverse 3583 (function (lambda (x) (math-add (math-solve-get-sign 3584 (math-normalize 3585 (list 'calcFunc-arccos x))) 3586 (math-solve-get-int 3587 (math-full-circle t)))))) 3588 3589(put 'calcFunc-tan 'math-inverse 3590 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x)) 3591 (math-solve-get-int 3592 (math-half-circle t)))))) 3593 3594(put 'calcFunc-arcsin 'math-inverse 3595 (function (lambda (x) (math-normalize (list 'calcFunc-sin x))))) 3596 3597(put 'calcFunc-arccos 'math-inverse 3598 (function (lambda (x) (math-normalize (list 'calcFunc-cos x))))) 3599 3600(put 'calcFunc-arctan 'math-inverse 3601 (function (lambda (x) (math-normalize (list 'calcFunc-tan x))))) 3602 3603(put 'calcFunc-sinh 'math-inverse 3604 (function (lambda (x) (let ((n (math-solve-get-int 1))) 3605 (math-add (math-mul (math-normalize 3606 (list 'calcFunc-arcsinh x)) 3607 (math-pow -1 n)) 3608 (math-mul (math-half-circle t) 3609 (math-mul 3610 '(var i var-i) 3611 n))))))) 3612(put 'calcFunc-sinh 'math-inverse-sign 1) 3613 3614(put 'calcFunc-cosh 'math-inverse 3615 (function (lambda (x) (math-add (math-solve-get-sign 3616 (math-normalize 3617 (list 'calcFunc-arccosh x))) 3618 (math-mul (math-full-circle t) 3619 (math-solve-get-int 3620 '(var i var-i))))))) 3621 3622(put 'calcFunc-tanh 'math-inverse 3623 (function (lambda (x) (math-add (math-normalize 3624 (list 'calcFunc-arctanh x)) 3625 (math-mul (math-half-circle t) 3626 (math-solve-get-int 3627 '(var i var-i))))))) 3628(put 'calcFunc-tanh 'math-inverse-sign 1) 3629 3630(put 'calcFunc-arcsinh 'math-inverse 3631 (function (lambda (x) (math-normalize (list 'calcFunc-sinh x))))) 3632(put 'calcFunc-arcsinh 'math-inverse-sign 1) 3633 3634(put 'calcFunc-arccosh 'math-inverse 3635 (function (lambda (x) (math-normalize (list 'calcFunc-cosh x))))) 3636 3637(put 'calcFunc-arctanh 'math-inverse 3638 (function (lambda (x) (math-normalize (list 'calcFunc-tanh x))))) 3639(put 'calcFunc-arctanh 'math-inverse-sign 1) 3640 3641 3642 3643(defun calcFunc-taylor (expr var num) 3644 (let ((x0 0) (v var)) 3645 (if (memq (car-safe var) '(+ - calcFunc-eq)) 3646 (setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var)) 3647 v (nth 1 var))) 3648 (or (and (eq (car-safe v) 'var) 3649 (math-expr-contains expr v) 3650 (natnump num) 3651 (let ((accum (math-expr-subst expr v x0)) 3652 (var2 (if (eq (car var) 'calcFunc-eq) 3653 (cons '- (cdr var)) 3654 var)) 3655 (n 0) 3656 (nfac 1) 3657 (fprime expr)) 3658 (while (and (<= (setq n (1+ n)) num) 3659 (setq fprime (calcFunc-deriv fprime v nil t))) 3660 (setq fprime (math-simplify fprime) 3661 nfac (math-mul nfac n) 3662 accum (math-add accum 3663 (math-div (math-mul (math-pow var2 n) 3664 (math-expr-subst 3665 fprime v x0)) 3666 nfac)))) 3667 (and fprime 3668 (math-normalize accum)))) 3669 (list 'calcFunc-taylor expr var num)))) 3670 3671(provide 'calcalg2) 3672 3673;;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4 3674;;; calcalg2.el ends here 3675