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

   1 ;;; calcalg2.el --- more algebraic functions for Calc
   2
   3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
   4 ;;   2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   5
   6 ;; Author: David Gillespie <daveg@synaptics.com>
   7 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
   8
   9 ;; This file is part of GNU Emacs.
  10
  11 ;; GNU Emacs is free software: you can redistribute it and/or modify
  12 ;; it under the terms of the GNU General Public License as published by
  13 ;; the Free Software Foundation, either version 3 of the License, or
  14 ;; (at your option) any later version.
  15
  16 ;; GNU Emacs is distributed in the hope that it will be useful,
  17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 ;; GNU General Public License for more details.
  20
  21 ;; You should have received a copy of the GNU General Public License
  22 ;; along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 ;;; Commentary:
  25
  26 ;;; Code:
  27
  28 ;; This file is autoloaded from calc-ext.el.
  29
  30 (require 'calc-ext)
  31 (require 'calc-macs)
  32
  33 (defun calc-derivative (var num)
  34   (interactive "sDifferentiate with respect to: \np")
  35   (calc-slow-wrapper
  36    (when (< num 0)
  37      (error "Order of derivative must be positive"))
  38    (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv))
  39          n expr)
  40      (if (or (equal var "") (equal var "$"))
  41          (setq n 2
  42                expr (calc-top-n 2)
  43                var (calc-top-n 1))
  44        (setq var (math-read-expr var))
  45        (when (eq (car-safe var) 'error)
  46          (error "Bad format in expression: %s" (nth 1 var)))
  47        (setq n 1
  48              expr (calc-top-n 1)))
  49      (while (>= (setq num (1- num)) 0)
  50        (setq expr (list func expr var)))
  51      (calc-enter-result n "derv" expr))))
  52
  53 (defun calc-integral (var &optional arg)
  54   (interactive "sIntegration variable: \nP")
  55   (if arg
  56       (calc-tabular-command 'calcFunc-integ "Integration" "intg" nil var nil nil)
  57     (calc-slow-wrapper
  58      (if (or (equal var "") (equal var "$"))
  59          (calc-enter-result 2 "intg" (list 'calcFunc-integ
  60                                            (calc-top-n 2)
  61                                            (calc-top-n 1)))
  62        (let ((var (math-read-expr var)))
  63          (if (eq (car-safe var) 'error)
  64              (error "Bad format in expression: %s" (nth 1 var)))
  65          (calc-enter-result 1 "intg" (list 'calcFunc-integ
  66                                            (calc-top-n 1)
  67                                            var)))))))
  68
  69 (defun calc-num-integral (&optional varname lowname highname)
  70   (interactive "sIntegration variable: ")
  71   (calc-tabular-command 'calcFunc-ninteg "Integration" "nint"
  72                         nil varname lowname highname))
  73
  74 (defun calc-summation (arg &optional varname lowname highname)
  75   (interactive "P\nsSummation variable: ")
  76   (calc-tabular-command 'calcFunc-sum "Summation" "sum"
  77                         arg varname lowname highname))
  78
  79 (defun calc-alt-summation (arg &optional varname lowname highname)
  80   (interactive "P\nsSummation variable: ")
  81   (calc-tabular-command 'calcFunc-asum "Summation" "asum"
  82                         arg varname lowname highname))
  83
  84 (defun calc-product (arg &optional varname lowname highname)
  85   (interactive "P\nsIndex variable: ")
  86   (calc-tabular-command 'calcFunc-prod "Index" "prod"
  87                         arg varname lowname highname))
  88
  89 (defun calc-tabulate (arg &optional varname lowname highname)
  90   (interactive "P\nsIndex variable: ")
  91   (calc-tabular-command 'calcFunc-table "Index" "tabl"
  92                         arg varname lowname highname))
  93
  94 (defun calc-tabular-command (func prompt prefix arg varname lowname highname)
  95   (calc-slow-wrapper
  96    (let (var (low nil) (high nil) (step nil) stepname stepnum (num 1) expr)
  97      (if (consp arg)
  98          (setq stepnum 1)
  99        (setq stepnum 0))
 100      (if (or (equal varname "") (equal varname "$") (null varname))
 101          (setq high (calc-top-n (+ stepnum 1))
 102                low (calc-top-n (+ stepnum 2))
 103                var (calc-top-n (+ stepnum 3))
 104                num (+ stepnum 4))
 105        (setq var (if (stringp varname) (math-read-expr varname) varname))
 106        (if (eq (car-safe var) 'error)
 107            (error "Bad format in expression: %s" (nth 1 var)))
 108        (or lowname
 109            (setq lowname (read-string (concat prompt " variable: " varname
 110                                               ", from: "))))
 111        (if (or (equal lowname "") (equal lowname "$"))
 112            (setq high (calc-top-n (+ stepnum 1))
 113                  low (calc-top-n (+ stepnum 2))
 114                  num (+ stepnum 3))
 115          (setq low (if (stringp lowname) (math-read-expr lowname) lowname))
 116          (if (eq (car-safe low) 'error)
 117              (error "Bad format in expression: %s" (nth 1 low)))
 118          (or highname
 119              (setq highname (read-string (concat prompt " variable: " varname
 120                                                  ", from: " lowname
 121                                                  ", to: "))))
 122          (if (or (equal highname "") (equal highname "$"))
 123              (setq high (calc-top-n (+ stepnum 1))
 124                    num (+ stepnum 2))
 125            (setq high (if (stringp highname) (math-read-expr highname)
 126                         highname))
 127            (if (eq (car-safe high) 'error)
 128                (error "Bad format in expression: %s" (nth 1 high)))
 129            (if (consp arg)
 130                (progn
 131                  (setq stepname (read-string (concat prompt " variable: "
 132                                                      varname
 133                                                      ", from: " lowname
 134                                                      ", to: " highname
 135                                                      ", step: ")))
 136                  (if (or (equal stepname "") (equal stepname "$"))
 137                      (setq step (calc-top-n 1)
 138                            num 2)
 139                    (setq step (math-read-expr stepname))
 140                    (if (eq (car-safe step) 'error)
 141                        (error "Bad format in expression: %s"
 142                               (nth 1 step)))))))))
 143      (or step
 144          (if (consp arg)
 145              (setq step (calc-top-n 1))
 146            (if arg
 147                (setq step (prefix-numeric-value arg)))))
 148      (setq expr (calc-top-n num))
 149      (calc-enter-result num prefix (append (list func expr var low high)
 150                                            (and step (list step)))))))
 151
 152 (defun calc-solve-for (var)
 153   (interactive "sVariable(s) to solve for: ")
 154   (calc-slow-wrapper
 155    (let ((func (if (calc-is-inverse)
 156                    (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv)
 157                  (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve))))
 158      (if (or (equal var "") (equal var "$"))
 159          (calc-enter-result 2 "solv" (list func
 160                                            (calc-top-n 2)
 161                                            (calc-top-n 1)))
 162        (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
 163                            (not (string-match "\\[" var)))
 164                       (math-read-expr (concat "[" var "]"))
 165                     (math-read-expr var))))
 166          (if (eq (car-safe var) 'error)
 167              (error "Bad format in expression: %s" (nth 1 var)))
 168          (calc-enter-result 1 "solv" (list func
 169                                            (calc-top-n 1)
 170                                            var)))))))
 171
 172 (defun calc-poly-roots (var)
 173   (interactive "sVariable to solve for: ")
 174   (calc-slow-wrapper
 175    (if (or (equal var "") (equal var "$"))
 176        (calc-enter-result 2 "prts" (list 'calcFunc-roots
 177                                          (calc-top-n 2)
 178                                          (calc-top-n 1)))
 179      (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
 180                          (not (string-match "\\[" var)))
 181                     (math-read-expr (concat "[" var "]"))
 182                   (math-read-expr var))))
 183        (if (eq (car-safe var) 'error)
 184            (error "Bad format in expression: %s" (nth 1 var)))
 185        (calc-enter-result 1 "prts" (list 'calcFunc-roots
 186                                          (calc-top-n 1)
 187                                          var))))))
 188
 189 (defun calc-taylor (var nterms)
 190   (interactive "sTaylor expansion variable: \nNNumber of terms: ")
 191   (calc-slow-wrapper
 192    (let ((var (math-read-expr var)))
 193      (if (eq (car-safe var) 'error)
 194          (error "Bad format in expression: %s" (nth 1 var)))
 195      (calc-enter-result 1 "tylr" (list 'calcFunc-taylor
 196                                        (calc-top-n 1)
 197                                        var
 198                                        (prefix-numeric-value nterms))))))
 199
 200
 201 ;; The following are global variables used by math-derivative and some
 202 ;; related functions
 203 (defvar math-deriv-var)
 204 (defvar math-deriv-total)
 205 (defvar math-deriv-symb)
 206 (defvar math-decls-cache)
 207 (defvar math-decls-all)
 208
 209 (defun math-derivative (expr)
 210   (cond ((equal expr math-deriv-var)
 211          1)
 212         ((or (Math-scalarp expr)
 213              (eq (car expr) 'sdev)
 214              (and (eq (car expr) 'var)
 215                   (or (not math-deriv-total)
 216                       (math-const-var expr)
 217                       (progn
 218                         (math-setup-declarations)
 219                         (memq 'const (nth 1 (or (assq (nth 2 expr)
 220                                                       math-decls-cache)
 221                                                 math-decls-all)))))))
 222          0)
 223         ((eq (car expr) '+)
 224          (math-add (math-derivative (nth 1 expr))
 225                    (math-derivative (nth 2 expr))))
 226         ((eq (car expr) '-)
 227          (math-sub (math-derivative (nth 1 expr))
 228                    (math-derivative (nth 2 expr))))
 229         ((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt
 230                                         calcFunc-gt calcFunc-leq calcFunc-geq))
 231          (list (car expr)
 232                (math-derivative (nth 1 expr))
 233                (math-derivative (nth 2 expr))))
 234         ((eq (car expr) 'neg)
 235          (math-neg (math-derivative (nth 1 expr))))
 236         ((eq (car expr) '*)
 237          (math-add (math-mul (nth 2 expr)
 238                              (math-derivative (nth 1 expr)))
 239                    (math-mul (nth 1 expr)
 240                              (math-derivative (nth 2 expr)))))
 241         ((eq (car expr) '/)
 242          (math-sub (math-div (math-derivative (nth 1 expr))
 243                              (nth 2 expr))
 244                    (math-div (math-mul (nth 1 expr)
 245                                        (math-derivative (nth 2 expr)))
 246                              (math-sqr (nth 2 expr)))))
 247         ((eq (car expr) '^)
 248          (let ((du (math-derivative (nth 1 expr)))
 249                (dv (math-derivative (nth 2 expr))))
 250            (or (Math-zerop du)
 251                (setq du (math-mul (nth 2 expr)
 252                                   (math-mul (math-normalize
 253                                              (list '^
 254                                                    (nth 1 expr)
 255                                                    (math-add (nth 2 expr) -1)))
 256                                             du))))
 257            (or (Math-zerop dv)
 258                (setq dv (math-mul (math-normalize
 259                                    (list 'calcFunc-ln (nth 1 expr)))
 260                                   (math-mul expr dv))))
 261            (math-add du dv)))
 262         ((eq (car expr) '%)
 263          (math-derivative (nth 1 expr)))   ; a reasonable definition
 264         ((eq (car expr) 'vec)
 265          (math-map-vec 'math-derivative expr))
 266         ((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im))
 267               (= (length expr) 2))
 268          (list (car expr) (math-derivative (nth 1 expr))))
 269         ((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol))
 270               (= (length expr) 3))
 271          (let ((d (math-derivative (nth 1 expr))))
 272            (if (math-numberp d)
 273                0    ; assume x and x_1 are independent vars
 274              (list (car expr) d (nth 2 expr)))))
 275         (t (or (and (symbolp (car expr))
 276                     (if (= (length expr) 2)
 277                         (let ((handler (get (car expr) 'math-derivative)))
 278                           (and handler
 279                                (let ((deriv (math-derivative (nth 1 expr))))
 280                                  (if (Math-zerop deriv)
 281                                      deriv
 282                                    (math-mul (funcall handler (nth 1 expr))
 283                                              deriv)))))
 284                       (let ((handler (get (car expr) 'math-derivative-n)))
 285                         (and handler
 286                              (funcall handler expr)))))
 287                (and (not (eq math-deriv-symb 'pre-expand))
 288                     (let ((exp (math-expand-formula expr)))
 289                       (and exp
 290                            (or (let ((math-deriv-symb 'pre-expand))
 291                                  (catch 'math-deriv (math-derivative expr)))
 292                                (math-derivative exp)))))
 293                (if (or (Math-objvecp expr)
 294                        (eq (car expr) 'var)
 295                        (not (symbolp (car expr))))
 296                    (if math-deriv-symb
 297                        (throw 'math-deriv nil)
 298                      (list (if math-deriv-total 'calcFunc-tderiv 'calcFunc-deriv)
 299                            expr
 300                            math-deriv-var))
 301                  (let ((accum 0)
 302                        (arg expr)
 303                        (n 1)
 304                        derv)
 305                    (while (setq arg (cdr arg))
 306                      (or (Math-zerop (setq derv (math-derivative (car arg))))
 307                          (let ((func (intern (concat (symbol-name (car expr))
 308                                                      "'"
 309                                                      (if (> n 1)
 310                                                          (int-to-string n)
 311                                                        ""))))
 312                                (prop (cond ((= (length expr) 2)
 313                                             'math-derivative-1)
 314                                            ((= (length expr) 3)
 315                                             'math-derivative-2)
 316                                            ((= (length expr) 4)
 317                                             'math-derivative-3)
 318                                            ((= (length expr) 5)
 319                                             'math-derivative-4)
 320                                            ((= (length expr) 6)
 321                                             'math-derivative-5))))
 322                            (setq accum
 323                                  (math-add
 324                                   accum
 325                                   (math-mul
 326                                    derv
 327                                    (let ((handler (get func prop)))
 328                                      (or (and prop handler
 329                                               (apply handler (cdr expr)))
 330                                          (if (and math-deriv-symb
 331                                                   (not (get func
 332                                                             'calc-user-defn)))
 333                                              (throw 'math-deriv nil)
 334                                            (cons func (cdr expr))))))))))
 335                      (setq n (1+ n)))
 336                    accum))))))
 337
 338 (defun calcFunc-deriv (expr math-deriv-var &optional deriv-value math-deriv-symb)
 339   (let* ((math-deriv-total nil)
 340          (res (catch 'math-deriv (math-derivative expr))))
 341     (or (eq (car-safe res) 'calcFunc-deriv)
 342         (null res)
 343         (setq res (math-normalize res)))
 344     (and res
 345          (if deriv-value
 346              (math-expr-subst res math-deriv-var deriv-value)
 347            res))))
 348
 349 (defun calcFunc-tderiv (expr math-deriv-var &optional deriv-value math-deriv-symb)
 350   (math-setup-declarations)
 351   (let* ((math-deriv-total t)
 352          (res (catch 'math-deriv (math-derivative expr))))
 353     (or (eq (car-safe res) 'calcFunc-tderiv)
 354         (null res)
 355         (setq res (math-normalize res)))
 356     (and res
 357          (if deriv-value
 358              (math-expr-subst res math-deriv-var deriv-value)
 359            res))))
 360
 361 (put 'calcFunc-inv\' 'math-derivative-1
 362      (function (lambda (u) (math-neg (math-div 1 (math-sqr u))))))
 363
 364 (put 'calcFunc-sqrt\' 'math-derivative-1
 365      (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u))))))
 366
 367 (put 'calcFunc-deg\' 'math-derivative-1
 368      (function (lambda (u) (math-div-float '(float 18 1) (math-pi)))))
 369
 370 (put 'calcFunc-rad\' 'math-derivative-1
 371      (function (lambda (u) (math-pi-over-180))))
 372
 373 (put 'calcFunc-ln\' 'math-derivative-1
 374      (function (lambda (u) (math-div 1 u))))
 375
 376 (put 'calcFunc-log10\' 'math-derivative-1
 377      (function (lambda (u)
 378                  (math-div (math-div 1 (math-normalize '(calcFunc-ln 10)))
 379                            u))))
 380
 381 (put 'calcFunc-lnp1\' 'math-derivative-1
 382      (function (lambda (u) (math-div 1 (math-add u 1)))))
 383
 384 (put 'calcFunc-log\' 'math-derivative-2
 385      (function (lambda (x b)
 386                  (and (not (Math-zerop b))
 387                       (let ((lnv (math-normalize
 388                                   (list 'calcFunc-ln b))))
 389                         (math-div 1 (math-mul lnv x)))))))
 390
 391 (put 'calcFunc-log\'2 'math-derivative-2
 392      (function (lambda (x b)
 393                  (let ((lnv (list 'calcFunc-ln b)))
 394                    (math-neg (math-div (list 'calcFunc-log x b)
 395                                        (math-mul lnv b)))))))
 396
 397 (put 'calcFunc-exp\' 'math-derivative-1
 398      (function (lambda (u) (math-normalize (list 'calcFunc-exp u)))))
 399
 400 (put 'calcFunc-expm1\' 'math-derivative-1
 401      (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u)))))
 402
 403 (put 'calcFunc-sin\' 'math-derivative-1
 404      (function (lambda (u) (math-to-radians-2 (math-normalize
 405                                                (list 'calcFunc-cos u))))))
 406
 407 (put 'calcFunc-cos\' 'math-derivative-1
 408      (function (lambda (u) (math-neg (math-to-radians-2
 409                                       (math-normalize
 410                                        (list 'calcFunc-sin u)))))))
 411
 412 (put 'calcFunc-tan\' 'math-derivative-1
 413      (function (lambda (u) (math-to-radians-2
 414                             (math-sqr
 415                              (math-normalize
 416                               (list 'calcFunc-sec u)))))))
 417
 418 (put 'calcFunc-sec\' 'math-derivative-1
 419      (function (lambda (u) (math-to-radians-2
 420                             (math-mul
 421                              (math-normalize
 422                               (list 'calcFunc-sec u))
 423                              (math-normalize
 424                               (list 'calcFunc-tan u)))))))
 425
 426 (put 'calcFunc-csc\' 'math-derivative-1
 427      (function (lambda (u) (math-neg
 428                             (math-to-radians-2
 429                              (math-mul
 430                               (math-normalize
 431                                (list 'calcFunc-csc u))
 432                               (math-normalize
 433                                (list 'calcFunc-cot u))))))))
 434
 435 (put 'calcFunc-cot\' 'math-derivative-1
 436      (function (lambda (u) (math-neg
 437                             (math-to-radians-2
 438                              (math-sqr
 439                               (math-normalize
 440                                (list 'calcFunc-csc u))))))))
 441
 442 (put 'calcFunc-arcsin\' 'math-derivative-1
 443      (function (lambda (u)
 444                  (math-from-radians-2
 445                   (math-div 1 (math-normalize
 446                                (list 'calcFunc-sqrt
 447                                      (math-sub 1 (math-sqr u)))))))))
 448
 449 (put 'calcFunc-arccos\' 'math-derivative-1
 450      (function (lambda (u)
 451                  (math-from-radians-2
 452                   (math-div -1 (math-normalize
 453                                 (list 'calcFunc-sqrt
 454                                       (math-sub 1 (math-sqr u)))))))))
 455
 456 (put 'calcFunc-arctan\' 'math-derivative-1
 457      (function (lambda (u) (math-from-radians-2
 458                             (math-div 1 (math-add 1 (math-sqr u)))))))
 459
 460 (put 'calcFunc-sinh\' 'math-derivative-1
 461      (function (lambda (u) (math-normalize (list 'calcFunc-cosh u)))))
 462
 463 (put 'calcFunc-cosh\' 'math-derivative-1
 464      (function (lambda (u) (math-normalize (list 'calcFunc-sinh u)))))
 465
 466 (put 'calcFunc-tanh\' 'math-derivative-1
 467      (function (lambda (u) (math-sqr
 468                             (math-normalize
 469                              (list 'calcFunc-sech u))))))
 470
 471 (put 'calcFunc-sech\' 'math-derivative-1
 472      (function (lambda (u) (math-neg
 473                             (math-mul
 474                              (math-normalize (list 'calcFunc-sech u))
 475                              (math-normalize (list 'calcFunc-tanh u)))))))
 476
 477 (put 'calcFunc-csch\' 'math-derivative-1
 478      (function (lambda (u) (math-neg
 479                             (math-mul
 480                              (math-normalize (list 'calcFunc-csch u))
 481                              (math-normalize (list 'calcFunc-coth u)))))))
 482
 483 (put 'calcFunc-coth\' 'math-derivative-1
 484      (function (lambda (u) (math-neg
 485                             (math-sqr
 486                              (math-normalize
 487                               (list 'calcFunc-csch u)))))))
 488
 489 (put 'calcFunc-arcsinh\' 'math-derivative-1
 490      (function (lambda (u)
 491                  (math-div 1 (math-normalize
 492                               (list 'calcFunc-sqrt
 493                                     (math-add (math-sqr u) 1)))))))
 494
 495 (put 'calcFunc-arccosh\' 'math-derivative-1
 496      (function (lambda (u)
 497                   (math-div 1 (math-normalize
 498                                (list 'calcFunc-sqrt
 499                                      (math-add (math-sqr u) -1)))))))
 500
 501 (put 'calcFunc-arctanh\' 'math-derivative-1
 502      (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u))))))
 503
 504 (put 'calcFunc-bern\'2 'math-derivative-2
 505      (function (lambda (n x)
 506                  (math-mul n (list 'calcFunc-bern (math-add n -1) x)))))
 507
 508 (put 'calcFunc-euler\'2 'math-derivative-2
 509      (function (lambda (n x)
 510                  (math-mul n (list 'calcFunc-euler (math-add n -1) x)))))
 511
 512 (put 'calcFunc-gammag\'2 'math-derivative-2
 513      (function (lambda (a x) (math-deriv-gamma a x 1))))
 514
 515 (put 'calcFunc-gammaG\'2 'math-derivative-2
 516      (function (lambda (a x) (math-deriv-gamma a x -1))))
 517
 518 (put 'calcFunc-gammaP\'2 'math-derivative-2
 519      (function (lambda (a x) (math-deriv-gamma a x
 520                                                (math-div
 521                                                 1 (math-normalize
 522                                                    (list 'calcFunc-gamma
 523                                                          a)))))))
 524
 525 (put 'calcFunc-gammaQ\'2 'math-derivative-2
 526      (function (lambda (a x) (math-deriv-gamma a x
 527                                                (math-div
 528                                                 -1 (math-normalize
 529                                                     (list 'calcFunc-gamma
 530                                                           a)))))))
 531
 532 (defun math-deriv-gamma (a x scale)
 533   (math-mul scale
 534             (math-mul (math-pow x (math-add a -1))
 535                       (list 'calcFunc-exp (math-neg x)))))
 536
 537 (put 'calcFunc-betaB\' 'math-derivative-3
 538      (function (lambda (x a b) (math-deriv-beta x a b 1))))
 539
 540 (put 'calcFunc-betaI\' 'math-derivative-3
 541      (function (lambda (x a b) (math-deriv-beta x a b
 542                                                 (math-div
 543                                                  1 (list 'calcFunc-beta
 544                                                          a b))))))
 545
 546 (defun math-deriv-beta (x a b scale)
 547   (math-mul (math-mul (math-pow x (math-add a -1))
 548                       (math-pow (math-sub 1 x) (math-add b -1)))
 549             scale))
 550
 551 (put 'calcFunc-erf\' 'math-derivative-1
 552      (function (lambda (x) (math-div 2
 553                                      (math-mul (list 'calcFunc-exp
 554                                                      (math-sqr x))
 555                                                (if calc-symbolic-mode
 556                                                    '(calcFunc-sqrt
 557                                                      (var pi var-pi))
 558                                                  (math-sqrt-pi)))))))
 559
 560 (put 'calcFunc-erfc\' 'math-derivative-1
 561      (function (lambda (x) (math-div -2
 562                                      (math-mul (list 'calcFunc-exp
 563                                                      (math-sqr x))
 564                                                (if calc-symbolic-mode
 565                                                    '(calcFunc-sqrt
 566                                                      (var pi var-pi))
 567                                                  (math-sqrt-pi)))))))
 568
 569 (put 'calcFunc-besJ\'2 'math-derivative-2
 570      (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ
 571                                                        (math-add v -1)
 572                                                        z)
 573                                                  (list 'calcFunc-besJ
 574                                                        (math-add v 1)
 575                                                        z))
 576                                        2))))
 577
 578 (put 'calcFunc-besY\'2 'math-derivative-2
 579      (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY
 580                                                        (math-add v -1)
 581                                                        z)
 582                                                  (list 'calcFunc-besY
 583                                                        (math-add v 1)
 584                                                        z))
 585                                        2))))
 586
 587 (put 'calcFunc-sum 'math-derivative-n
 588      (function
 589       (lambda (expr)
 590         (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var)
 591             (throw 'math-deriv nil)
 592           (cons 'calcFunc-sum
 593                 (cons (math-derivative (nth 1 expr))
 594                       (cdr (cdr expr))))))))
 595
 596 (put 'calcFunc-prod 'math-derivative-n
 597      (function
 598       (lambda (expr)
 599         (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var)
 600             (throw 'math-deriv nil)
 601           (math-mul expr
 602                     (cons 'calcFunc-sum
 603                           (cons (math-div (math-derivative (nth 1 expr))
 604                                           (nth 1 expr))
 605                                 (cdr (cdr expr)))))))))
 606
 607 (put 'calcFunc-integ 'math-derivative-n
 608      (function
 609       (lambda (expr)
 610         (if (= (length expr) 3)
 611             (if (equal (nth 2 expr) math-deriv-var)
 612                 (nth 1 expr)
 613               (math-normalize
 614                (list 'calcFunc-integ
 615                      (math-derivative (nth 1 expr))
 616                      (nth 2 expr))))
 617           (if (= (length expr) 5)
 618               (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr)
 619                                             (nth 3 expr)))
 620                     (upper (math-expr-subst (nth 1 expr) (nth 2 expr)
 621                                             (nth 4 expr))))
 622                 (math-add (math-sub (math-mul upper
 623                                               (math-derivative (nth 4 expr)))
 624                                     (math-mul lower
 625                                               (math-derivative (nth 3 expr))))
 626                           (if (equal (nth 2 expr) math-deriv-var)
 627                               0
 628                             (math-normalize
 629                              (list 'calcFunc-integ
 630                                    (math-derivative (nth 1 expr)) (nth 2 expr)
 631                                    (nth 3 expr) (nth 4 expr)))))))))))
 632
 633 (put 'calcFunc-if 'math-derivative-n
 634      (function
 635       (lambda (expr)
 636         (and (= (length expr) 4)
 637              (list 'calcFunc-if (nth 1 expr)
 638                    (math-derivative (nth 2 expr))
 639                    (math-derivative (nth 3 expr)))))))
 640
 641 (put 'calcFunc-subscr 'math-derivative-n
 642      (function
 643       (lambda (expr)
 644         (and (= (length expr) 3)
 645              (list 'calcFunc-subscr (nth 1 expr)
 646                    (math-derivative (nth 2 expr)))))))
 647
 648
 649 (defvar math-integ-var '(var X ---))
 650 (defvar math-integ-var-2 '(var Y ---))
 651 (defvar math-integ-vars (list 'f math-integ-var math-integ-var-2))
 652 (defvar math-integ-var-list (list math-integ-var))
 653 (defvar math-integ-var-list-list (list math-integ-var-list))
 654
 655 ;; math-integ-depth is a local variable for math-try-integral, but is used
 656 ;; by math-integral and math-tracing-integral
 657 ;; which are called (directly or indirectly) by math-try-integral.
 658 (defvar math-integ-depth)
 659 ;; math-integ-level is a local variable for math-try-integral, but is used
 660 ;; by math-integral, math-do-integral, math-tracing-integral,
 661 ;; math-sub-integration, math-integrate-by-parts and
 662 ;; math-integrate-by-substitution, which are called (directly or
 663 ;; indirectly) by math-try-integral.
 664 (defvar math-integ-level)
 665 ;; math-integral-limit is a local variable for calcFunc-integ, but is
 666 ;; used by math-tracing-integral, math-sub-integration and
 667 ;; math-try-integration.
 668 (defvar math-integral-limit)
 669
 670 (defmacro math-tracing-integral (&rest parts)
 671   (list 'and
 672         'trace-buffer
 673         (list 'save-excursion
 674               '(set-buffer trace-buffer)
 675               '(goto-char (point-max))
 676               (list 'and
 677                     '(bolp)
 678                     '(insert (make-string (- math-integral-limit
 679                                              math-integ-level) 32)
 680                              (format "%2d " math-integ-depth)
 681                              (make-string math-integ-level 32)))
 682               ;;(list 'condition-case 'err
 683                     (cons 'insert parts)
 684                 ;;    '(error (insert (prin1-to-string err))))
 685               '(sit-for 0))))
 686
 687 ;;; The following wrapper caches results and avoids infinite recursion.
 688 ;;; Each cache entry is: ( A B )          Integral of A is B;
 689 ;;;                      ( A N )          Integral of A failed at level N;
 690 ;;;                      ( A busy )       Currently working on integral of A;
 691 ;;;                      ( A parts )      Currently working, integ-by-parts;
 692 ;;;                      ( A parts2 )     Currently working, integ-by-parts;
 693 ;;;                      ( A cancelled )  Ignore this cache entry;
 694 ;;;                      ( A [B] )        Same result as for math-cur-record = B.
 695
 696 ;; math-cur-record is a local variable for math-try-integral, but is used
 697 ;; by math-integral, math-replace-integral-parts and math-integrate-by-parts
 698 ;; which are called (directly or indirectly) by math-try-integral, as well as
 699 ;; by calc-dump-integral-cache
 700 (defvar math-cur-record)
 701 ;; math-enable-subst and math-any-substs are local variables for
 702 ;; calcFunc-integ, but are used by math-integral and math-try-integral.
 703 (defvar math-enable-subst)
 704 (defvar math-any-substs)
 705
 706 ;; math-integ-msg is a local variable for math-try-integral, but is
 707 ;; used (both locally and non-locally) by math-integral.
 708 (defvar math-integ-msg)
 709
 710 (defvar math-integral-cache nil)
 711 (defvar math-integral-cache-state nil)
 712
 713 (defun math-integral (expr &optional simplify same-as-above)
 714   (let* ((simp math-cur-record)
 715          (math-cur-record (assoc expr math-integral-cache))
 716          (math-integ-depth (1+ math-integ-depth))
 717          (val 'cancelled))
 718     (math-tracing-integral "Integrating "
 719                            (math-format-value expr 1000)
 720                            "...\n")
 721     (and math-cur-record
 722          (progn
 723            (math-tracing-integral "Found "
 724                                   (math-format-value (nth 1 math-cur-record) 1000))
 725            (and (consp (nth 1 math-cur-record))
 726                 (math-replace-integral-parts math-cur-record))
 727            (math-tracing-integral " => "
 728                                   (math-format-value (nth 1 math-cur-record) 1000)
 729                                   "\n")))
 730     (or (and math-cur-record
 731              (not (eq (nth 1 math-cur-record) 'cancelled))
 732              (or (not (integerp (nth 1 math-cur-record)))
 733                  (>= (nth 1 math-cur-record) math-integ-level)))
 734         (and (math-integral-contains-parts expr)
 735              (progn
 736                (setq val nil)
 737                t))
 738         (unwind-protect
 739             (progn
 740               (let (math-integ-msg)
 741                 (if (eq calc-display-working-message 'lots)
 742                     (progn
 743                       (calc-set-command-flag 'clear-message)
 744                       (setq math-integ-msg (format
 745                                             "Working... Integrating %s"
 746                                             (math-format-flat-expr expr 0)))
 747                       (message "%s" math-integ-msg)))
 748                 (if math-cur-record
 749                     (setcar (cdr math-cur-record)
 750                             (if same-as-above (vector simp) 'busy))
 751                   (setq math-cur-record
 752                         (list expr (if same-as-above (vector simp) 'busy))
 753                         math-integral-cache (cons math-cur-record
 754                                                   math-integral-cache)))
 755                 (if (eq simplify 'yes)
 756                     (progn
 757                       (math-tracing-integral "Simplifying...")
 758                       (setq simp (math-simplify expr))
 759                       (setq val (if (equal simp expr)
 760                                     (progn
 761                                       (math-tracing-integral " no change\n")
 762                                       (math-do-integral expr))
 763                                   (math-tracing-integral " simplified\n")
 764                                   (math-integral simp 'no t))))
 765                   (or (setq val (math-do-integral expr))
 766                       (eq simplify 'no)
 767                       (let ((simp (math-simplify expr)))
 768                         (or (equal simp expr)
 769                             (progn
 770                               (math-tracing-integral "Trying again after "
 771                                                      "simplification...\n")
 772                               (setq val (math-integral simp 'no t))))))))
 773               (if (eq calc-display-working-message 'lots)
 774                   (message "%s" math-integ-msg)))
 775           (setcar (cdr math-cur-record) (or val
 776                                        (if (or math-enable-subst
 777                                                (not math-any-substs))
 778                                            math-integ-level
 779                                          'cancelled)))))
 780     (setq val math-cur-record)
 781     (while (vectorp (nth 1 val))
 782       (setq val (aref (nth 1 val) 0)))
 783     (setq val (if (memq (nth 1 val) '(parts parts2))
 784                   (progn
 785                     (setcar (cdr val) 'parts2)
 786                     (list 'var 'PARTS val))
 787                 (and (consp (nth 1 val))
 788                      (nth 1 val))))
 789     (math-tracing-integral "Integral of "
 790                            (math-format-value expr 1000)
 791                            "  is  "
 792                            (math-format-value val 1000)
 793                            "\n")
 794     val))
 795
 796 (defun math-integral-contains-parts (expr)
 797   (if (Math-primp expr)
 798       (and (eq (car-safe expr) 'var)
 799            (eq (nth 1 expr) 'PARTS)
 800            (listp (nth 2 expr)))
 801     (while (and (setq expr (cdr expr))
 802                 (not (math-integral-contains-parts (car expr)))))
 803     expr))
 804
 805 (defun math-replace-integral-parts (expr)
 806   (or (Math-primp expr)
 807       (while (setq expr (cdr expr))
 808         (and (consp (car expr))
 809              (if (eq (car (car expr)) 'var)
 810                  (and (eq (nth 1 (car expr)) 'PARTS)
 811                       (consp (nth 2 (car expr)))
 812                       (if (listp (nth 1 (nth 2 (car expr))))
 813                           (progn
 814                             (setcar expr (nth 1 (nth 2 (car expr))))
 815                             (math-replace-integral-parts (cons 'foo expr)))
 816                         (setcar (cdr math-cur-record) 'cancelled)))
 817                (math-replace-integral-parts (car expr)))))))
 818
 819 (defvar math-linear-subst-tried t
 820   "Non-nil means that a linear substitution has been tried.")
 821
 822 ;; The variable math-has-rules is a local variable for math-try-integral,
 823 ;; but is used by math-do-integral, which is called (non-directly) by
 824 ;; math-try-integral.
 825 (defvar math-has-rules)
 826
 827 ;; math-old-integ is a local variable for math-do-integral, but is
 828 ;; used by math-sub-integration.
 829 (defvar math-old-integ)
 830
 831 ;; The variables math-t1, math-t2 and math-t3 are local to
 832 ;; math-do-integral, math-try-solve-for and math-decompose-poly, but
 833 ;; are used by functions they call (directly or indirectly);
 834 ;; math-do-integral calls math-do-integral-methods;
 835 ;; math-try-solve-for calls math-try-solve-prod,
 836 ;; math-solve-find-root-term and math-solve-find-root-in-prod;
 837 ;; math-decompose-poly calls math-solve-poly-funny-powers and
 838 ;; math-solve-crunch-poly.
 839 (defvar math-t1)
 840 (defvar math-t2)
 841 (defvar math-t3)
 842
 843 (defun math-do-integral (expr)
 844   (let ((math-linear-subst-tried nil)
 845         math-t1 math-t2)
 846     (or (cond ((not (math-expr-contains expr math-integ-var))
 847                (math-mul expr math-integ-var))
 848               ((equal expr math-integ-var)
 849                (math-div (math-sqr expr) 2))
 850               ((eq (car expr) '+)
 851                (and (setq math-t1 (math-integral (nth 1 expr)))
 852                     (setq math-t2 (math-integral (nth 2 expr)))
 853                     (math-add math-t1 math-t2)))
 854               ((eq (car expr) '-)
 855                (and (setq math-t1 (math-integral (nth 1 expr)))
 856                     (setq math-t2 (math-integral (nth 2 expr)))
 857                     (math-sub math-t1 math-t2)))
 858               ((eq (car expr) 'neg)
 859                (and (setq math-t1 (math-integral (nth 1 expr)))
 860                     (math-neg math-t1)))
 861               ((eq (car expr) '*)
 862                (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
 863                       (and (setq math-t1 (math-integral (nth 2 expr)))
 864                            (math-mul (nth 1 expr) math-t1)))
 865                      ((not (math-expr-contains (nth 2 expr) math-integ-var))
 866                       (and (setq math-t1 (math-integral (nth 1 expr)))
 867                            (math-mul math-t1 (nth 2 expr))))
 868                      ((memq (car-safe (nth 1 expr)) '(+ -))
 869                       (math-integral (list (car (nth 1 expr))
 870                                            (math-mul (nth 1 (nth 1 expr))
 871                                                      (nth 2 expr))
 872                                            (math-mul (nth 2 (nth 1 expr))
 873                                                      (nth 2 expr)))
 874                                      'yes t))
 875                      ((memq (car-safe (nth 2 expr)) '(+ -))
 876                       (math-integral (list (car (nth 2 expr))
 877                                            (math-mul (nth 1 (nth 2 expr))
 878                                                      (nth 1 expr))
 879                                            (math-mul (nth 2 (nth 2 expr))
 880                                                      (nth 1 expr)))
 881                                      'yes t))))
 882               ((eq (car expr) '/)
 883                (cond ((and (not (math-expr-contains (nth 1 expr)
 884                                                     math-integ-var))
 885                            (not (math-equal-int (nth 1 expr) 1)))
 886                       (and (setq math-t1 (math-integral (math-div 1 (nth 2 expr))))
 887                            (math-mul (nth 1 expr) math-t1)))
 888                      ((not (math-expr-contains (nth 2 expr) math-integ-var))
 889                       (and (setq math-t1 (math-integral (nth 1 expr)))
 890                            (math-div math-t1 (nth 2 expr))))
 891                      ((and (eq (car-safe (nth 1 expr)) '*)
 892                            (not (math-expr-contains (nth 1 (nth 1 expr))
 893                                                     math-integ-var)))
 894                       (and (setq math-t1 (math-integral
 895                                      (math-div (nth 2 (nth 1 expr))
 896                                                (nth 2 expr))))
 897                            (math-mul math-t1 (nth 1 (nth 1 expr)))))
 898                      ((and (eq (car-safe (nth 1 expr)) '*)
 899                            (not (math-expr-contains (nth 2 (nth 1 expr))
 900                                                     math-integ-var)))
 901                       (and (setq math-t1 (math-integral
 902                                      (math-div (nth 1 (nth 1 expr))
 903                                                (nth 2 expr))))
 904                            (math-mul math-t1 (nth 2 (nth 1 expr)))))
 905                      ((and (eq (car-safe (nth 2 expr)) '*)
 906                            (not (math-expr-contains (nth 1 (nth 2 expr))
 907                                                     math-integ-var)))
 908                       (and (setq math-t1 (math-integral
 909                                      (math-div (nth 1 expr)
 910                                                (nth 2 (nth 2 expr)))))
 911                            (math-div math-t1 (nth 1 (nth 2 expr)))))
 912                      ((and (eq (car-safe (nth 2 expr)) '*)
 913                            (not (math-expr-contains (nth 2 (nth 2 expr))
 914                                                     math-integ-var)))
 915                       (and (setq math-t1 (math-integral
 916                                      (math-div (nth 1 expr)
 917                                                (nth 1 (nth 2 expr)))))
 918                            (math-div math-t1 (nth 2 (nth 2 expr)))))
 919                      ((eq (car-safe (nth 2 expr)) 'calcFunc-exp)
 920                       (math-integral
 921                        (math-mul (nth 1 expr)
 922                                  (list 'calcFunc-exp
 923                                        (math-neg (nth 1 (nth 2 expr)))))))))
 924               ((eq (car expr) '^)
 925                (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
 926                       (or (and (setq math-t1 (math-is-polynomial (nth 2 expr)
 927                                                             math-integ-var 1))
 928                                (math-div expr
 929                                          (math-mul (nth 1 math-t1)
 930                                                    (math-normalize
 931                                                     (list 'calcFunc-ln
 932                                                           (nth 1 expr))))))
 933                           (math-integral
 934                            (list 'calcFunc-exp
 935                                  (math-mul (nth 2 expr)
 936                                            (math-normalize
 937                                             (list 'calcFunc-ln
 938                                                   (nth 1 expr)))))
 939                            'yes t)))
 940                      ((not (math-expr-contains (nth 2 expr) math-integ-var))
 941                       (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0))
 942                           (math-integral
 943                            (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr))))
 944                            nil t)
 945                         (or (and (setq math-t1 (math-is-polynomial (nth 1 expr)
 946                                                               math-integ-var
 947                                                               1))
 948                                  (setq math-t2 (math-add (nth 2 expr) 1))
 949                                  (math-div (math-pow (nth 1 expr) math-t2)
 950                                            (math-mul math-t2 (nth 1 math-t1))))
 951                             (and (Math-negp (nth 2 expr))
 952                                  (math-integral
 953                                   (math-div 1
 954                                             (math-pow (nth 1 expr)
 955                                                       (math-neg
 956                                                        (nth 2 expr))))
 957                                   nil t))
 958                             nil))))))
 959
 960         ;; Integral of a polynomial.
 961         (and (setq math-t1 (math-is-polynomial expr math-integ-var 20))
 962              (let ((accum 0)
 963                    (n 1))
 964                (while math-t1
 965                  (if (setq accum (math-add accum
 966                                            (math-div (math-mul (car math-t1)
 967                                                                (math-pow
 968                                                                 math-integ-var
 969                                                                 n))
 970                                                      n))
 971                            math-t1 (cdr math-t1))
 972                      (setq n (1+ n))))
 973                accum))
 974
 975         ;; Try looking it up!
 976         (cond ((= (length expr) 2)
 977                (and (symbolp (car expr))
 978                     (setq math-t1 (get (car expr) 'math-integral))
 979                     (progn
 980                       (while (and math-t1
 981                                   (not (setq math-t2 (funcall (car math-t1)
 982                                                          (nth 1 expr)))))
 983                         (setq math-t1 (cdr math-t1)))
 984                       (and math-t2 (math-normalize math-t2)))))
 985               ((= (length expr) 3)
 986                (and (symbolp (car expr))
 987                     (setq math-t1 (get (car expr) 'math-integral-2))
 988                     (progn
 989                       (while (and math-t1
 990                                   (not (setq math-t2 (funcall (car math-t1)
 991                                                          (nth 1 expr)
 992                                                          (nth 2 expr)))))
 993                         (setq math-t1 (cdr math-t1)))
 994                       (and math-t2 (math-normalize math-t2))))))
 995
 996         ;; Integral of a rational function.
 997         (and (math-ratpoly-p expr math-integ-var)
 998              (setq math-t1 (calcFunc-apart expr math-integ-var))
 999              (not (equal math-t1 expr))
1000              (math-integral math-t1))
1001
1002         ;; Try user-defined integration rules.
1003         (and math-has-rules
1004              (let ((math-old-integ (symbol-function 'calcFunc-integ))
1005                    (input (list 'calcFunc-integtry expr math-integ-var))
1006                    res part)
1007                (unwind-protect
1008                    (progn
1009                      (fset 'calcFunc-integ 'math-sub-integration)
1010                      (setq res (math-rewrite input
1011                                              '(var IntegRules var-IntegRules)
1012                                              1))
1013                      (fset 'calcFunc-integ math-old-integ)
1014                      (and (not (equal res input))
1015                           (if (setq part (math-expr-calls
1016                                           res '(calcFunc-integsubst)))
1017                               (and (memq (length part) '(3 4 5))
1018                                    (let ((parts (mapcar
1019                                                  (function
1020                                                   (lambda (x)
1021                                                     (math-expr-subst
1022                                                      x (nth 2 part)
1023                                                      math-integ-var)))
1024                                                  (cdr part))))
1025                                      (math-integrate-by-substitution
1026                                       expr (car parts) t
1027                                       (or (nth 2 parts)
1028                                           (list 'calcFunc-integfailed
1029                                                 math-integ-var))
1030                                       (nth 3 parts))))
1031                             (if (not (math-expr-calls res
1032                                                       '(calcFunc-integtry
1033                                                         calcFunc-integfailed)))
1034                                 res))))
1035                  (fset 'calcFunc-integ math-old-integ))))
1036
1037         ;; See if the function is a symbolic derivative.
1038         (and (string-match "'" (symbol-name (car expr)))
1039              (let ((name (symbol-name (car expr)))
1040                    (p expr) (n 0) (which nil) (bad nil))
1041                (while (setq n (1+ n) p (cdr p))
1042                  (if (equal (car p) math-integ-var)
1043                      (if which (setq bad t) (setq which n))
1044                    (if (math-expr-contains (car p) math-integ-var)
1045                        (setq bad t))))
1046                (and which (not bad)
1047                     (let ((prime (if (= which 1) "'" (format "'%d" which))))
1048                       (and (string-match (concat prime "\\('['0-9]*\\|$\\)")
1049                                          name)
1050                            (cons (intern
1051                                   (concat
1052                                    (substring name 0 (match-beginning 0))
1053                                    (substring name (+ (match-beginning 0)
1054                                                       (length prime)))))
1055                                  (cdr expr)))))))
1056
1057         ;; Try transformation methods (parts, substitutions).
1058         (and (> math-integ-level 0)
1059              (math-do-integral-methods expr))
1060
1061         ;; Try expanding the function's definition.
1062         (let ((res (math-expand-formula expr)))
1063           (and res
1064                (math-integral res))))))
1065
1066 (defun math-sub-integration (expr &rest rest)
1067   (or (if (or (not rest)
1068               (and (< math-integ-level math-integral-limit)
1069                    (eq (car rest) math-integ-var)))
1070           (math-integral expr)
1071         (let ((res (apply math-old-integ expr rest)))
1072           (and (or (= math-integ-level math-integral-limit)
1073                    (not (math-expr-calls res 'calcFunc-integ)))
1074                res)))
1075       (list 'calcFunc-integfailed expr)))
1076
1077 ;; math-so-far is a local variable for math-do-integral-methods, but
1078 ;; is used by math-integ-try-linear-substitutions and
1079 ;; math-integ-try-substitutions.
1080 (defvar math-so-far)
1081
1082 ;; math-integ-expr is a local variable for math-do-integral-methods,
1083 ;; but is used by math-integ-try-linear-substitutions and
1084 ;; math-integ-try-substitutions.
1085 (defvar math-integ-expr)
1086
1087 (defun math-do-integral-methods (math-integ-expr)
1088   (let ((math-so-far math-integ-var-list-list)
1089         rat-in)
1090
1091     ;; Integration by substitution, for various likely sub-expressions.
1092     ;; (In first pass, we look only for sub-exprs that are linear in X.)
1093     (or (math-integ-try-linear-substitutions math-integ-expr)
1094         (math-integ-try-substitutions math-integ-expr)
1095
1096         ;; If function has sines and cosines, try tan(x/2) substitution.
1097         (and (let ((p (setq rat-in (math-expr-rational-in math-integ-expr))))
1098                (while (and p
1099                            (memq (car (car p)) '(calcFunc-sin
1100                                                  calcFunc-cos
1101                                                  calcFunc-tan
1102                                                  calcFunc-sec
1103                                                  calcFunc-csc
1104                                                  calcFunc-cot))
1105                            (equal (nth 1 (car p)) math-integ-var))
1106                  (setq p (cdr p)))
1107                (null p))
1108              (or (and (math-integ-parts-easy math-integ-expr)
1109                       (math-integ-try-parts math-integ-expr t))
1110                  (math-integrate-by-good-substitution
1111                   math-integ-expr (list 'calcFunc-tan (math-div math-integ-var 2)))))
1112
1113         ;; If function has sinh and cosh, try tanh(x/2) substitution.
1114         (and (let ((p rat-in))
1115                (while (and p
1116                            (memq (car (car p)) '(calcFunc-sinh
1117                                                  calcFunc-cosh
1118                                                  calcFunc-tanh
1119                                                  calcFunc-sech
1120                                                  calcFunc-csch
1121                                                  calcFunc-coth
1122                                                  calcFunc-exp))
1123                            (equal (nth 1 (car p)) math-integ-var))
1124                  (setq p (cdr p)))
1125                (null p))
1126              (or (and (math-integ-parts-easy math-integ-expr)
1127                       (math-integ-try-parts math-integ-expr t))
1128                  (math-integrate-by-good-substitution
1129                   math-integ-expr (list 'calcFunc-tanh (math-div math-integ-var 2)))))
1130
1131         ;; If function has square roots, try sin, tan, or sec substitution.
1132         (and (let ((p rat-in))
1133                (setq math-t1 nil)
1134                (while (and p
1135                            (or (equal (car p) math-integ-var)
1136                                (and (eq (car (car p)) 'calcFunc-sqrt)
1137                                     (setq math-t1 (math-is-polynomial
1138                                               (nth 1 (setq math-t2 (car p)))
1139                                               math-integ-var 2)))))
1140                  (setq p (cdr p)))
1141                (and (null p) math-t1))
1142              (if (cdr (cdr math-t1))
1143                  (if (math-guess-if-neg (nth 2 math-t1))
1144                      (let* ((c (math-sqrt (math-neg (nth 2 math-t1))))
1145                             (d (math-div (nth 1 math-t1) (math-mul -2 c)))
1146                             (a (math-sqrt (math-add (car math-t1) (math-sqr d)))))
1147                        (math-integrate-by-good-substitution
1148                         math-integ-expr (list 'calcFunc-arcsin
1149                                    (math-div-thru
1150                                     (math-add (math-mul c math-integ-var) d)
1151                                     a))))
1152                    (let* ((c (math-sqrt (nth 2 math-t1)))
1153                           (d (math-div (nth 1 math-t1) (math-mul 2 c)))
1154                           (aa (math-sub (car math-t1) (math-sqr d))))
1155                      (if (and nil (not (and (eq d 0) (eq c 1))))
1156                          (math-integrate-by-good-substitution
1157                           math-integ-expr (math-add (math-mul c math-integ-var) d))
1158                        (if (math-guess-if-neg aa)
1159                            (math-integrate-by-good-substitution
1160                             math-integ-expr (list 'calcFunc-arccosh
1161                                        (math-div-thru
1162                                         (math-add (math-mul c math-integ-var)
1163                                                   d)
1164                                         (math-sqrt (math-neg aa)))))
1165                          (math-integrate-by-good-substitution
1166                           math-integ-expr (list 'calcFunc-arcsinh
1167                                      (math-div-thru
1168                                       (math-add (math-mul c math-integ-var)
1169                                                 d)
1170                                       (math-sqrt aa))))))))
1171                (math-integrate-by-good-substitution math-integ-expr math-t2)) )
1172
1173         ;; Try integration by parts.
1174         (math-integ-try-parts math-integ-expr)
1175
1176         ;; Give up.
1177         nil)))
1178
1179 (defun math-integ-parts-easy (expr)
1180   (cond ((Math-primp expr) t)
1181         ((memq (car expr) '(+ - *))
1182          (and (math-integ-parts-easy (nth 1 expr))
1183               (math-integ-parts-easy (nth 2 expr))))
1184         ((eq (car expr) '/)
1185          (and (math-integ-parts-easy (nth 1 expr))
1186               (math-atomic-factorp (nth 2 expr))))
1187         ((eq (car expr) '^)
1188          (and (natnump (nth 2 expr))
1189               (math-integ-parts-easy (nth 1 expr))))
1190         ((eq (car expr) 'neg)
1191          (math-integ-parts-easy (nth 1 expr)))
1192         (t t)))
1193
1194 ;; math-prev-parts-v is local to calcFunc-integ (as well as
1195 ;; math-integrate-by-parts), but is used by math-integ-try-parts.
1196 (defvar math-prev-parts-v)
1197
1198 ;; math-good-parts is local to calcFunc-integ (as well as
1199 ;; math-integ-try-parts), but is used by math-integrate-by-parts.
1200 (defvar math-good-parts)
1201
1202
1203 (defun math-integ-try-parts (expr &optional math-good-parts)
1204   ;; Integration by parts:
1205   ;;   integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
1206   ;;     where h(x) = integ(g(x),x).
1207   (or (let ((exp (calcFunc-expand expr)))
1208         (and (not (equal exp expr))
1209              (math-integral exp)))
1210       (and (eq (car expr) '*)
1211            (let ((first-bad (or (math-polynomial-p (nth 1 expr)
1212                                                    math-integ-var)
1213                                 (equal (nth 2 expr) math-prev-parts-v))))
1214              (or (and first-bad   ; so try this one first
1215                       (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))
1216                  (math-integrate-by-parts (nth 2 expr) (nth 1 expr))
1217                  (and (not first-bad)
1218                       (math-integrate-by-parts (nth 1 expr) (nth 2 expr))))))
1219       (and (eq (car expr) '/)
1220            (math-expr-contains (nth 1 expr) math-integ-var)
1221            (let ((recip (math-div 1 (nth 2 expr))))
1222              (or (math-integrate-by-parts (nth 1 expr) recip)
1223                  (math-integrate-by-parts recip (nth 1 expr)))))
1224       (and (eq (car expr) '^)
1225            (math-integrate-by-parts (math-pow (nth 1 expr)
1226                                               (math-sub (nth 2 expr) 1))
1227                                     (nth 1 expr)))))
1228
1229 (defun math-integrate-by-parts (u vprime)
1230   (let ((math-integ-level (if (or math-good-parts
1231                                   (math-polynomial-p u math-integ-var))
1232                               math-integ-level
1233                             (1- math-integ-level)))
1234         (math-doing-parts t)
1235         v temp)
1236     (and (>= math-integ-level 0)
1237          (unwind-protect
1238              (progn
1239                (setcar (cdr math-cur-record) 'parts)
1240                (math-tracing-integral "Integrating by parts, u = "
1241                                       (math-format-value u 1000)
1242                                       ", v' = "
1243                                       (math-format-value vprime 1000)
1244                                       "\n")
1245                (and (setq v (math-integral vprime))
1246                     (setq temp (calcFunc-deriv u math-integ-var nil t))
1247                     (setq temp (let ((math-prev-parts-v v))
1248                                  (math-integral (math-mul v temp) 'yes)))
1249                     (setq temp (math-sub (math-mul u v) temp))
1250                     (if (eq (nth 1 math-cur-record) 'parts)
1251                         (calcFunc-expand temp)
1252                       (setq v (list 'var 'PARTS math-cur-record)
1253                             temp (let (calc-next-why)
1254                                    (math-simplify-extended
1255                                     (math-solve-for (math-sub v temp) 0 v nil)))
1256                             temp (if (and (eq (car-safe temp) '/)
1257                                           (math-zerop (nth 2 temp)))
1258                                      nil temp)))))
1259            (setcar (cdr math-cur-record) 'busy)))))
1260
1261 ;;; This tries two different formulations, hoping the algebraic simplifier
1262 ;;; will be strong enough to handle at least one.
1263 (defun math-integrate-by-substitution (expr u &optional user uinv uinvprime)
1264   (and (> math-integ-level 0)
1265        (let ((math-integ-level (max (- math-integ-level 2) 0)))
1266          (math-integrate-by-good-substitution expr u user uinv uinvprime))))
1267
1268 (defun math-integrate-by-good-substitution (expr u &optional user
1269                                                  uinv uinvprime)
1270   (let ((math-living-dangerously t)
1271         deriv temp)
1272     (and (setq uinv (if uinv
1273                         (math-expr-subst uinv math-integ-var
1274                                          math-integ-var-2)
1275                       (let (calc-next-why)
1276                         (math-solve-for u
1277                                         math-integ-var-2
1278                                         math-integ-var nil))))
1279          (progn
1280            (math-tracing-integral "Integrating by substitution, u = "
1281                                   (math-format-value u 1000)
1282                                   "\n")
1283            (or (and (setq deriv (calcFunc-deriv u
1284                                                 math-integ-var nil
1285                                                 (not user)))
1286                     (setq temp (math-integral (math-expr-subst
1287                                                (math-expr-subst
1288                                                 (math-expr-subst
1289                                                  (math-div expr deriv)
1290                                                  u
1291                                                  math-integ-var-2)
1292                                                 math-integ-var
1293                                                 uinv)
1294                                                math-integ-var-2
1295                                                math-integ-var)
1296                                               'yes)))
1297                (and (setq deriv (or uinvprime
1298                                     (calcFunc-deriv uinv
1299                                                     math-integ-var-2
1300                                                     math-integ-var
1301                                                     (not user))))
1302                     (setq temp (math-integral (math-mul
1303                                                (math-expr-subst
1304                                                 (math-expr-subst
1305                                                  (math-expr-subst
1306                                                   expr
1307                                                   u
1308                                                   math-integ-var-2)
1309                                                  math-integ-var
1310                                                  uinv)
1311                                                 math-integ-var-2
1312                                                 math-integ-var)
1313                                                deriv)
1314                                               'yes)))))
1315          (math-simplify-extended
1316           (math-expr-subst temp math-integ-var u)))))
1317
1318 ;;; Look for substitutions of the form u = a x + b.
1319 (defun math-integ-try-linear-substitutions (sub-expr)
1320   (setq math-linear-subst-tried t)
1321   (and (not (Math-primp sub-expr))
1322        (or (and (not (memq (car sub-expr) '(+ - * / neg)))
1323                 (not (and (eq (car sub-expr) '^)
1324                           (integerp (nth 2 sub-expr))))
1325                 (math-expr-contains sub-expr math-integ-var)
1326                 (let ((res nil))
1327                   (while (and (setq sub-expr (cdr sub-expr))
1328                               (or (not (math-linear-in (car sub-expr)
1329                                                        math-integ-var))
1330                                   (assoc (car sub-expr) math-so-far)
1331                                   (progn
1332                                     (setq math-so-far (cons (list (car sub-expr))
1333                                                        math-so-far))
1334                                     (not (setq res
1335                                                (math-integrate-by-substitution
1336                                                 math-integ-expr (car sub-expr))))))))
1337                   res))
1338            (let ((res nil))
1339              (while (and (setq sub-expr (cdr sub-expr))
1340                          (not (setq res (math-integ-try-linear-substitutions
1341                                          (car sub-expr))))))
1342              res))))
1343
1344 ;;; Recursively try different substitutions based on various sub-expressions.
1345 (defun math-integ-try-substitutions (sub-expr &optional allow-rat)
1346   (and (not (Math-primp sub-expr))
1347        (not (assoc sub-expr math-so-far))
1348        (math-expr-contains sub-expr math-integ-var)
1349        (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg)))
1350                          (not (and (eq (car sub-expr) '^)
1351                                    (integerp (nth 2 sub-expr)))))
1352                     (setq allow-rat t)
1353                   (prog1 allow-rat (setq allow-rat nil)))
1354                 (not (eq sub-expr math-integ-expr))
1355                 (or (math-integrate-by-substitution math-integ-expr sub-expr)
1356                     (and (eq (car sub-expr) '^)
1357                          (integerp (nth 2 sub-expr))
1358                          (< (nth 2 sub-expr) 0)
1359                          (math-integ-try-substitutions
1360                           (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr)))
1361                           t))))
1362            (let ((res nil))
1363              (setq math-so-far (cons (list sub-expr) math-so-far))
1364              (while (and (setq sub-expr (cdr sub-expr))
1365                          (not (setq res (math-integ-try-substitutions
1366                                          (car sub-expr) allow-rat)))))
1367              res))))
1368
1369 ;; The variable math-expr-parts is local to math-expr-rational-in,
1370 ;; but is used by math-expr-rational-in-rec
1371 (defvar math-expr-parts)
1372
1373 (defun math-expr-rational-in (expr)
1374   (let ((math-expr-parts nil))
1375     (math-expr-rational-in-rec expr)
1376     (mapcar 'car math-expr-parts)))
1377
1378 (defun math-expr-rational-in-rec (expr)
1379   (cond ((Math-primp expr)
1380          (and (equal expr math-integ-var)
1381               (not (assoc expr math-expr-parts))
1382               (setq math-expr-parts (cons (list expr) math-expr-parts))))
1383         ((or (memq (car expr) '(+ - * / neg))
1384              (and (eq (car expr) '^) (integerp (nth 2 expr))))
1385          (math-expr-rational-in-rec (nth 1 expr))
1386          (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr))))
1387         ((and (eq (car expr) '^)
1388               (eq (math-quarter-integer (nth 2 expr)) 2))
1389          (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr))))
1390         (t
1391          (and (not (assoc expr math-expr-parts))
1392               (math-expr-contains expr math-integ-var)
1393               (setq math-expr-parts (cons (list expr) math-expr-parts))))))
1394
1395 (defun math-expr-calls (expr funcs &optional arg-contains)
1396   (if (consp expr)
1397       (if (or (memq (car expr) funcs)
1398               (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt)
1399                    (eq (math-quarter-integer (nth 2 expr)) 2)))
1400           (and (or (not arg-contains)
1401                    (math-expr-contains expr arg-contains))
1402                expr)
1403         (and (not (Math-primp expr))
1404              (let ((res nil))
1405                (while (and (setq expr (cdr expr))
1406                            (not (setq res (math-expr-calls
1407                                            (car expr) funcs arg-contains)))))
1408                res)))))
1409
1410 (defun math-fix-const-terms (expr except-vars)
1411   (cond ((not (math-expr-depends expr except-vars)) 0)
1412         ((Math-primp expr) expr)
1413         ((eq (car expr) '+)
1414          (math-add (math-fix-const-terms (nth 1 expr) except-vars)
1415                    (math-fix-const-terms (nth 2 expr) except-vars)))
1416         ((eq (car expr) '-)
1417          (math-sub (math-fix-const-terms (nth 1 expr) except-vars)
1418                    (math-fix-const-terms (nth 2 expr) except-vars)))
1419         (t expr)))
1420
1421 ;; Command for debugging the Calculator's symbolic integrator.
1422 (defun calc-dump-integral-cache (&optional arg)
1423   (interactive "P")
1424   (let ((buf (current-buffer)))
1425     (unwind-protect
1426         (let ((p math-integral-cache)
1427               math-cur-record)
1428           (display-buffer (get-buffer-create "*Integral Cache*"))
1429           (set-buffer (get-buffer "*Integral Cache*"))
1430           (erase-buffer)
1431           (while p
1432             (setq math-cur-record (car p))
1433             (or arg (math-replace-integral-parts math-cur-record))
1434             (insert (math-format-flat-expr (car math-cur-record) 0)
1435                     " --> "
1436                     (if (symbolp (nth 1 math-cur-record))
1437                         (concat "(" (symbol-name (nth 1 math-cur-record)) ")")
1438                       (math-format-flat-expr (nth 1 math-cur-record) 0))
1439                     "\n")
1440             (setq p (cdr p)))
1441           (goto-char (point-min)))
1442       (set-buffer buf))))
1443
1444 ;; The variable math-max-integral-limit is local to calcFunc-integ,
1445 ;; but is used by math-try-integral.
1446 (defvar math-max-integral-limit)
1447
1448 (defun math-try-integral (expr)
1449   (let ((math-integ-level math-integral-limit)
1450         (math-integ-depth 0)
1451         (math-integ-msg "Working...done")
1452         (math-cur-record nil)   ; a technicality
1453         (math-integrating t)
1454         (calc-prefer-frac t)
1455         (calc-symbolic-mode t)
1456         (math-has-rules (calc-has-rules 'var-IntegRules)))
1457     (or (math-integral expr 'yes)
1458         (and math-any-substs
1459              (setq math-enable-subst t)
1460              (math-integral expr 'yes))
1461         (and (> math-max-integral-limit math-integral-limit)
1462              (setq math-integral-limit math-max-integral-limit
1463                    math-integ-level math-integral-limit)
1464              (math-integral expr 'yes)))))
1465
1466 (defvar var-IntegLimit nil)
1467
1468 (defun calcFunc-integ (expr var &optional low high)
1469   (cond
1470    ;; Do these even if the parts turn out not to be integrable.
1471    ((eq (car-safe expr) '+)
1472     (math-add (calcFunc-integ (nth 1 expr) var low high)
1473               (calcFunc-integ (nth 2 expr) var low high)))
1474    ((eq (car-safe expr) '-)
1475     (math-sub (calcFunc-integ (nth 1 expr) var low high)
1476               (calcFunc-integ (nth 2 expr) var low high)))
1477    ((eq (car-safe expr) 'neg)
1478     (math-neg (calcFunc-integ (nth 1 expr) var low high)))
1479    ((and (eq (car-safe expr) '*)
1480          (not (math-expr-contains (nth 1 expr) var)))
1481     (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high)))
1482    ((and (eq (car-safe expr) '*)
1483          (not (math-expr-contains (nth 2 expr) var)))
1484     (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1485    ((and (eq (car-safe expr) '/)
1486          (not (math-expr-contains (nth 1 expr) var))
1487          (not (math-equal-int (nth 1 expr) 1)))
1488     (math-mul (nth 1 expr)
1489               (calcFunc-integ (math-div 1 (nth 2 expr)) var low high)))
1490    ((and (eq (car-safe expr) '/)
1491          (not (math-expr-contains (nth 2 expr) var)))
1492     (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1493    ((and (eq (car-safe expr) '/)
1494          (eq (car-safe (nth 1 expr)) '*)
1495          (not (math-expr-contains (nth 1 (nth 1 expr)) var)))
1496     (math-mul (nth 1 (nth 1 expr))
1497               (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr))
1498                               var low high)))
1499    ((and (eq (car-safe expr) '/)
1500          (eq (car-safe (nth 1 expr)) '*)
1501          (not (math-expr-contains (nth 2 (nth 1 expr)) var)))
1502     (math-mul (nth 2 (nth 1 expr))
1503               (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr))
1504                               var low high)))
1505    ((and (eq (car-safe expr) '/)
1506          (eq (car-safe (nth 2 expr)) '*)
1507          (not (math-expr-contains (nth 1 (nth 2 expr)) var)))
1508     (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr)))
1509                               var low high)
1510               (nth 1 (nth 2 expr))))
1511    ((and (eq (car-safe expr) '/)
1512          (eq (car-safe (nth 2 expr)) '*)
1513          (not (math-expr-contains (nth 2 (nth 2 expr)) var)))
1514     (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr)))
1515                               var low high)
1516               (nth 2 (nth 2 expr))))
1517    ((eq (car-safe expr) 'vec)
1518     (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high)))
1519                        (cdr expr))))
1520    (t
1521     (let ((state (list calc-angle-mode
1522                        ;;calc-symbolic-mode
1523                        ;;calc-prefer-frac
1524                        calc-internal-prec
1525                        (calc-var-value 'var-IntegRules)
1526                        (calc-var-value 'var-IntegSimpRules))))
1527       (or (equal state math-integral-cache-state)
1528           (setq math-integral-cache-state state
1529                 math-integral-cache nil)))
1530     (let* ((math-max-integral-limit (or (and (natnump var-IntegLimit)
1531                                              var-IntegLimit)
1532                                         3))
1533            (math-integral-limit 1)
1534            (sexpr (math-expr-subst expr var math-integ-var))
1535            (trace-buffer (get-buffer "*Trace*"))
1536            (calc-language (if (eq calc-language 'big) nil calc-language))
1537            (math-any-substs t)
1538            (math-enable-subst nil)
1539            (math-prev-parts-v nil)
1540            (math-doing-parts nil)
1541            (math-good-parts nil)
1542            (res
1543             (if trace-buffer
1544                 (let ((calcbuf (current-buffer))
1545                       (calcwin (selected-window)))
1546                   (unwind-protect
1547                       (progn
1548                         (if (get-buffer-window trace-buffer)
1549                             (select-window (get-buffer-window trace-buffer)))
1550                         (set-buffer trace-buffer)
1551                         (goto-char (point-max))
1552                         (or (assq 'scroll-stop (buffer-local-variables))
1553                             (progn
1554                               (make-local-variable 'scroll-step)
1555                               (setq scroll-step 3)))
1556                         (insert "\n\n\n")
1557                         (set-buffer calcbuf)
1558                         (math-try-integral sexpr))
1559                     (select-window calcwin)
1560                       (set-buffer calcbuf)))
1561               (math-try-integral sexpr))))
1562       (if res
1563           (progn
1564             (if (calc-has-rules 'var-IntegAfterRules)
1565                 (setq res (math-rewrite res '(var IntegAfterRules
1566                                                   var-IntegAfterRules))))
1567             (math-simplify
1568              (if (and low high)
1569                  (math-sub (math-expr-subst res math-integ-var high)
1570                            (math-expr-subst res math-integ-var low))
1571                (setq res (math-fix-const-terms res math-integ-vars))
1572                (if low
1573                    (math-expr-subst res math-integ-var low)
1574                  (math-expr-subst res math-integ-var var)))))
1575         (append (list 'calcFunc-integ expr var)
1576                 (and low (list low))
1577                 (and high (list high))))))))
1578
1579
1580 (math-defintegral calcFunc-inv
1581   (math-integral (math-div 1 u)))
1582
1583 (math-defintegral calcFunc-conj
1584   (let ((int (math-integral u)))
1585     (and int
1586          (list 'calcFunc-conj int))))
1587
1588 (math-defintegral calcFunc-deg
1589   (let ((int (math-integral u)))
1590     (and int
1591          (list 'calcFunc-deg int))))
1592
1593 (math-defintegral calcFunc-rad
1594   (let ((int (math-integral u)))
1595     (and int
1596          (list 'calcFunc-rad int))))
1597
1598 (math-defintegral calcFunc-re
1599   (let ((int (math-integral u)))
1600     (and int
1601          (list 'calcFunc-re int))))
1602
1603 (math-defintegral calcFunc-im
1604   (let ((int (math-integral u)))
1605     (and int
1606          (list 'calcFunc-im int))))
1607
1608 (math-defintegral calcFunc-sqrt
1609   (and (equal u math-integ-var)
1610        (math-mul '(frac 2 3)
1611                  (list 'calcFunc-sqrt (math-pow u 3)))))
1612
1613 (math-defintegral calcFunc-exp
1614   (or (and (equal u math-integ-var)
1615            (list 'calcFunc-exp u))
1616       (let ((p (math-is-polynomial u math-integ-var 2)))
1617         (and (nth 2 p)
1618              (let ((sqa (math-sqrt (math-neg (nth 2 p)))))
1619                (math-div
1620                 (math-mul
1621                  (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi))
1622                                      sqa)
1623                            (math-normalize
1624                             (list 'calcFunc-exp
1625                                   (math-div (math-sub (math-mul (car p)
1626                                                                 (nth 2 p))
1627                                                       (math-div
1628                                                        (math-sqr (nth 1 p))
1629                                                        4))
1630                                             (nth 2 p)))))
1631                  (list 'calcFunc-erf
1632                        (math-sub (math-mul sqa math-integ-var)
1633                                  (math-div (nth 1 p) (math-mul 2 sqa)))))
1634                 2))))))
1635
1636 (math-defintegral calcFunc-ln
1637   (or (and (equal u math-integ-var)
1638            (math-sub (math-mul u (list 'calcFunc-ln u)) u))
1639       (and (eq (car u) '*)
1640            (math-integral (math-add (list 'calcFunc-ln (nth 1 u))
1641                                     (list 'calcFunc-ln (nth 2 u)))))
1642       (and (eq (car u) '/)
1643            (math-integral (math-sub (list 'calcFunc-ln (nth 1 u))
1644                                     (list 'calcFunc-ln (nth 2 u)))))
1645       (and (eq (car u) '^)
1646            (math-integral (math-mul (nth 2 u)
1647                                     (list 'calcFunc-ln (nth 1 u)))))))
1648
1649 (math-defintegral calcFunc-log10
1650   (and (equal u math-integ-var)
1651        (math-sub (math-mul u (list 'calcFunc-ln u))
1652                  (math-div u (list 'calcFunc-ln 10)))))
1653
1654 (math-defintegral-2 calcFunc-log
1655   (math-integral (math-div (list 'calcFunc-ln u)
1656                            (list 'calcFunc-ln v))))
1657
1658 (math-defintegral calcFunc-sin
1659   (or (and (equal u math-integ-var)
1660            (math-neg (math-from-radians-2 (list 'calcFunc-cos u))))
1661       (and (nth 2 (math-is-polynomial u math-integ-var 2))
1662            (math-integral (math-to-exponentials (list 'calcFunc-sin u))))))
1663
1664 (math-defintegral calcFunc-cos
1665   (or (and (equal u math-integ-var)
1666            (math-from-radians-2 (list 'calcFunc-sin u)))
1667       (and (nth 2 (math-is-polynomial u math-integ-var 2))
1668            (math-integral (math-to-exponentials (list 'calcFunc-cos u))))))
1669
1670 (math-defintegral calcFunc-tan
1671   (and (equal u math-integ-var)
1672        (math-from-radians-2
1673         (list 'calcFunc-ln (list 'calcFunc-sec u)))))
1674
1675 (math-defintegral calcFunc-sec
1676   (and (equal u math-integ-var)
1677        (math-from-radians-2
1678         (list 'calcFunc-ln
1679               (math-add
1680                (list 'calcFunc-sec u)
1681                (list 'calcFunc-tan u))))))
1682
1683 (math-defintegral calcFunc-csc
1684   (and (equal u math-integ-var)
1685        (math-from-radians-2
1686         (list 'calcFunc-ln
1687               (math-sub
1688                (list 'calcFunc-csc u)
1689                (list 'calcFunc-cot u))))))
1690
1691 (math-defintegral calcFunc-cot
1692   (and (equal u math-integ-var)
1693        (math-from-radians-2
1694         (list 'calcFunc-ln (list 'calcFunc-sin u)))))
1695
1696 (math-defintegral calcFunc-arcsin
1697   (and (equal u math-integ-var)
1698        (math-add (math-mul u (list 'calcFunc-arcsin u))
1699                  (math-from-radians-2
1700                   (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1701
1702 (math-defintegral calcFunc-arccos
1703   (and (equal u math-integ-var)
1704        (math-sub (math-mul u (list 'calcFunc-arccos u))
1705                  (math-from-radians-2
1706                   (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1707
1708 (math-defintegral calcFunc-arctan
1709   (and (equal u math-integ-var)
1710        (math-sub (math-mul u (list 'calcFunc-arctan u))
1711                  (math-from-radians-2
1712                   (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u)))
1713                             2)))))
1714
1715 (math-defintegral calcFunc-sinh
1716   (and (equal u math-integ-var)
1717        (list 'calcFunc-cosh u)))
1718
1719 (math-defintegral calcFunc-cosh
1720   (and (equal u math-integ-var)
1721        (list 'calcFunc-sinh u)))
1722
1723 (math-defintegral calcFunc-tanh
1724   (and (equal u math-integ-var)
1725        (list 'calcFunc-ln (list 'calcFunc-cosh u))))
1726
1727 (math-defintegral calcFunc-sech
1728   (and (equal u math-integ-var)
1729        (list 'calcFunc-arctan (list 'calcFunc-sinh u))))
1730
1731 (math-defintegral calcFunc-csch
1732   (and (equal u math-integ-var)
1733        (list 'calcFunc-ln (list 'calcFunc-tanh (math-div u 2)))))
1734
1735 (math-defintegral calcFunc-coth
1736   (and (equal u math-integ-var)
1737        (list 'calcFunc-ln (list 'calcFunc-sinh u))))
1738
1739 (math-defintegral calcFunc-arcsinh
1740   (and (equal u math-integ-var)
1741        (math-sub (math-mul u (list 'calcFunc-arcsinh u))
1742                  (list 'calcFunc-sqrt (math-add (math-sqr u) 1)))))
1743
1744 (math-defintegral calcFunc-arccosh
1745   (and (equal u math-integ-var)
1746        (math-sub (math-mul u (list 'calcFunc-arccosh u))
1747                  (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))
1748
1749 (math-defintegral calcFunc-arctanh
1750   (and (equal u math-integ-var)
1751        (math-sub (math-mul u (list 'calcFunc-arctan u))
1752                  (math-div (list 'calcFunc-ln
1753                                  (math-add 1 (math-sqr u)))
1754                            2))))
1755
1756 ;;; (Ax + B) / (ax^2 + bx + c)^n forms.
1757 (math-defintegral-2 /
1758   (math-integral-rational-funcs u v))
1759
1760 (defun math-integral-rational-funcs (u v)
1761   (let ((pu (math-is-polynomial u math-integ-var 1))
1762         (vpow 1) pv)
1763     (and pu
1764          (catch 'int-rat
1765            (if (and (eq (car-safe v) '^) (natnump (nth 2 v)))
1766                (setq vpow (nth 2 v)
1767                      v (nth 1 v)))
1768            (and (setq pv (math-is-polynomial v math-integ-var 2))
1769                 (let ((int (math-mul-thru
1770                             (car pu)
1771                             (math-integral-q02 (car pv) (nth 1 pv)
1772                                                (nth 2 pv) v vpow))))
1773                   (if (cdr pu)
1774                       (setq int (math-add int
1775                                           (math-mul-thru
1776                                            (nth 1 pu)
1777                                            (math-integral-q12
1778                                             (car pv) (nth 1 pv)
1779                                             (nth 2 pv) v vpow)))))
1780                   int))))))
1781
1782 (defun math-integral-q12 (a b c v vpow)
1783   (let (q)
1784     (cond ((not c)
1785            (cond ((= vpow 1)
1786                   (math-sub (math-div math-integ-var b)
1787                             (math-mul (math-div a (math-sqr b))
1788                                       (list 'calcFunc-ln v))))
1789                  ((= vpow 2)
1790                   (math-div (math-add (list 'calcFunc-ln v)
1791                                       (math-div a v))
1792                             (math-sqr b)))
1793                  (t
1794                   (let ((nm1 (math-sub vpow 1))
1795                         (nm2 (math-sub vpow 2)))
1796                     (math-div (math-sub
1797                                (math-div a (math-mul nm1 (math-pow v nm1)))
1798                                (math-div 1 (math-mul nm2 (math-pow v nm2))))
1799                               (math-sqr b))))))
1800           ((math-zerop
1801             (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1802            (let ((part (math-div b (math-mul 2 c))))
1803              (math-mul-thru (math-pow c vpow)
1804                             (math-integral-q12 part 1 nil
1805                                                (math-add math-integ-var part)
1806                                                (* vpow 2)))))
1807           ((= vpow 1)
1808            (and (math-ratp q) (math-negp q)
1809                 (let ((calc-symbolic-mode t))
1810                   (math-ratp (math-sqrt (math-neg q))))
1811                 (throw 'int-rat nil))  ; should have used calcFunc-apart first
1812            (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c))
1813                      (math-mul-thru (math-div b (math-mul 2 c))
1814                                     (math-integral-q02 a b c v 1))))
1815           (t
1816            (let ((n (1- vpow)))
1817              (math-sub (math-neg (math-div
1818                                   (math-add (math-mul b math-integ-var)
1819                                             (math-mul 2 a))
1820                                   (math-mul n (math-mul q (math-pow v n)))))
1821                        (math-mul-thru (math-div (math-mul b (1- (* 2 n)))
1822                                                 (math-mul n q))
1823                                       (math-integral-q02 a b c v n))))))))
1824
1825 (defun math-integral-q02 (a b c v vpow)
1826   (let (q rq part)
1827     (cond ((not c)
1828            (cond ((= vpow 1)
1829                   (math-div (list 'calcFunc-ln v) b))
1830                  (t
1831                   (math-div (math-pow v (- 1 vpow))
1832                             (math-mul (- 1 vpow) b)))))
1833           ((math-zerop
1834             (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1835            (let ((part (math-div b (math-mul 2 c))))
1836              (math-mul-thru (math-pow c vpow)
1837                             (math-integral-q02 part 1 nil
1838                                                (math-add math-integ-var part)
1839                                                (* vpow 2)))))
1840           ((progn
1841              (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b))
1842              (> vpow 1))
1843            (let ((n (1- vpow)))
1844              (math-add (math-div part (math-mul n (math-mul q (math-pow v n))))
1845                        (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c)
1846                                                 (math-mul n q))
1847                                       (math-integral-q02 a b c v n)))))
1848           ((math-guess-if-neg q)
1849            (setq rq (list 'calcFunc-sqrt (math-neg q)))
1850            ;;(math-div-thru (list 'calcFunc-ln
1851            ;;                   (math-div (math-sub part rq)
1852            ;;                             (math-add part rq)))
1853            ;;             rq)
1854            (math-div (math-mul -2 (list 'calcFunc-arctanh
1855                                         (math-div part rq)))
1856                      rq))
1857           (t
1858            (setq rq (list 'calcFunc-sqrt q))
1859            (math-div (math-mul 2 (math-to-radians-2
1860                                   (list 'calcFunc-arctan
1861                                         (math-div part rq))))
1862                      rq)))))
1863
1864
1865 (math-defintegral calcFunc-erf
1866   (and (equal u math-integ-var)
1867        (math-add (math-mul u (list 'calcFunc-erf u))
1868                  (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1869                                        (list 'calcFunc-sqrt
1870                                              '(var pi var-pi)))))))
1871
1872 (math-defintegral calcFunc-erfc
1873   (and (equal u math-integ-var)
1874        (math-sub (math-mul u (list 'calcFunc-erfc u))
1875                  (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1876                                        (list 'calcFunc-sqrt
1877                                              '(var pi var-pi)))))))
1878
1879
1880
1881
1882 (defvar math-tabulate-initial nil)
1883 (defvar math-tabulate-function nil)
1884
1885 ;; The variables calc-low and calc-high are local to calcFunc-table,
1886 ;; but are used by math-scan-for-limits.
1887 (defvar calc-low)
1888 (defvar calc-high)
1889
1890 (defun calcFunc-table (expr 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 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 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) var))
1953          (let* ((calc-next-why nil)
1954                 (low-val (math-solve-for (nth 2 x) 1 var nil))
1955                 (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x)))
1956                                           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 ;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4
3672 ;;; calcalg2.el ends here