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