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

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