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