cgit.sxemacs.org Git - packages/blob - xemacs-packages/calc/calc-alg-2.el

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