1;;; calc-arith.el --- arithmetic functions for Calc
2
3;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4;;   2005, 2006, 2007 Free Software Foundation, Inc.
5
6;; Author: David Gillespie <daveg@synaptics.com>
7;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8
9;; This file is part of GNU Emacs.
10
11;; GNU Emacs is free software; you can redistribute it and/or modify
12;; it under the terms of the GNU General Public License as published by
13;; the Free Software Foundation; either version 2, or (at your option)
14;; any later version.
15
16;; GNU Emacs is distributed in the hope that it will be useful,
17;; but WITHOUT ANY WARRANTY; without even the implied warranty of
18;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19;; GNU General Public License for more details.
20
21;; You should have received a copy of the GNU General Public License
22;; along with GNU Emacs; see the file COPYING.  If not, write to the
23;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24;; Boston, MA 02110-1301, USA.
25
26;;; Commentary:
27
28;;; Code:
29
30;; This file is autoloaded from calc-ext.el.
31
32(require 'calc-ext)
33(require 'calc-macs)
34
35;;; The following lists are not exhaustive.
36(defvar math-scalar-functions '(calcFunc-det
37				calcFunc-cnorm calcFunc-rnorm
38				calcFunc-vlen calcFunc-vcount
39				calcFunc-vsum calcFunc-vprod
40				calcFunc-vmin calcFunc-vmax))
41
42(defvar math-nonscalar-functions '(vec calcFunc-idn calcFunc-diag
43				       calcFunc-cvec calcFunc-index
44				       calcFunc-trn
45				       | calcFunc-append
46				       calcFunc-cons calcFunc-rcons
47				       calcFunc-tail calcFunc-rhead))
48
49(defvar math-scalar-if-args-functions '(+ - * / neg))
50
51(defvar math-real-functions '(calcFunc-arg
52			      calcFunc-re calcFunc-im
53			      calcFunc-floor calcFunc-ceil
54			      calcFunc-trunc calcFunc-round
55			      calcFunc-rounde calcFunc-roundu
56			      calcFunc-ffloor calcFunc-fceil
57			      calcFunc-ftrunc calcFunc-fround
58			      calcFunc-frounde calcFunc-froundu))
59
60(defvar math-positive-functions '())
61
62(defvar math-nonnegative-functions '(calcFunc-cnorm calcFunc-rnorm
63				     calcFunc-vlen calcFunc-vcount))
64
65(defvar math-real-scalar-functions '(% calcFunc-idiv calcFunc-abs
66				       calcFunc-choose calcFunc-perm
67				       calcFunc-eq calcFunc-neq
68				       calcFunc-lt calcFunc-gt
69				       calcFunc-leq calcFunc-geq
70				       calcFunc-lnot
71				       calcFunc-max calcFunc-min))
72
73(defvar math-real-if-arg-functions '(calcFunc-sin calcFunc-cos
74				     calcFunc-tan calcFunc-sec
75                                     calcFunc-csc calcFunc-cot
76                                     calcFunc-arctan
77				     calcFunc-sinh calcFunc-cosh
78				     calcFunc-tanh calcFunc-sech
79                                     calcFunc-csch calcFunc-coth
80                                     calcFunc-exp
81				     calcFunc-gamma calcFunc-fact))
82
83(defvar math-integer-functions '(calcFunc-idiv
84				 calcFunc-isqrt calcFunc-ilog
85				 calcFunc-vlen calcFunc-vcount))
86
87(defvar math-num-integer-functions '())
88
89(defvar math-rounding-functions '(calcFunc-floor
90				  calcFunc-ceil
91				  calcFunc-round calcFunc-trunc
92				  calcFunc-rounde calcFunc-roundu))
93
94(defvar math-float-rounding-functions '(calcFunc-ffloor
95					calcFunc-fceil
96					calcFunc-fround calcFunc-ftrunc
97					calcFunc-frounde calcFunc-froundu))
98
99(defvar math-integer-if-args-functions '(+ - * % neg calcFunc-abs
100					   calcFunc-min calcFunc-max
101					   calcFunc-choose calcFunc-perm))
102
103
104;;; Arithmetic.
105
106(defun calc-min (arg)
107  (interactive "P")
108  (calc-slow-wrapper
109   (calc-binary-op "min" 'calcFunc-min arg '(var inf var-inf))))
110
111(defun calc-max (arg)
112  (interactive "P")
113  (calc-slow-wrapper
114   (calc-binary-op "max" 'calcFunc-max arg '(neg (var inf var-inf)))))
115
116(defun calc-abs (arg)
117  (interactive "P")
118  (calc-slow-wrapper
119   (calc-unary-op "abs" 'calcFunc-abs arg)))
120
121
122(defun calc-idiv (arg)
123  (interactive "P")
124  (calc-slow-wrapper
125   (calc-binary-op "\\" 'calcFunc-idiv arg 1)))
126
127
128(defun calc-floor (arg)
129  (interactive "P")
130  (calc-slow-wrapper
131   (if (calc-is-inverse)
132       (if (calc-is-hyperbolic)
133	   (calc-unary-op "ceil" 'calcFunc-fceil arg)
134	 (calc-unary-op "ceil" 'calcFunc-ceil arg))
135     (if (calc-is-hyperbolic)
136	 (calc-unary-op "flor" 'calcFunc-ffloor arg)
137       (calc-unary-op "flor" 'calcFunc-floor arg)))))
138
139(defun calc-ceiling (arg)
140  (interactive "P")
141  (calc-invert-func)
142  (calc-floor arg))
143
144(defun calc-round (arg)
145  (interactive "P")
146  (calc-slow-wrapper
147   (if (calc-is-inverse)
148       (if (calc-is-hyperbolic)
149	   (calc-unary-op "trnc" 'calcFunc-ftrunc arg)
150	 (calc-unary-op "trnc" 'calcFunc-trunc arg))
151     (if (calc-is-hyperbolic)
152	 (calc-unary-op "rond" 'calcFunc-fround arg)
153       (calc-unary-op "rond" 'calcFunc-round arg)))))
154
155(defun calc-trunc (arg)
156  (interactive "P")
157  (calc-invert-func)
158  (calc-round arg))
159
160(defun calc-mant-part (arg)
161  (interactive "P")
162  (calc-slow-wrapper
163   (calc-unary-op "mant" 'calcFunc-mant arg)))
164
165(defun calc-xpon-part (arg)
166  (interactive "P")
167  (calc-slow-wrapper
168   (calc-unary-op "xpon" 'calcFunc-xpon arg)))
169
170(defun calc-scale-float (arg)
171  (interactive "P")
172  (calc-slow-wrapper
173   (calc-binary-op "scal" 'calcFunc-scf arg)))
174
175(defun calc-abssqr (arg)
176  (interactive "P")
177  (calc-slow-wrapper
178   (calc-unary-op "absq" 'calcFunc-abssqr arg)))
179
180(defun calc-sign (arg)
181  (interactive "P")
182  (calc-slow-wrapper
183   (calc-unary-op "sign" 'calcFunc-sign arg)))
184
185(defun calc-increment (arg)
186  (interactive "p")
187  (calc-wrapper
188   (calc-enter-result 1 "incr" (list 'calcFunc-incr (calc-top-n 1) arg))))
189
190(defun calc-decrement (arg)
191  (interactive "p")
192  (calc-wrapper
193   (calc-enter-result 1 "decr" (list 'calcFunc-decr (calc-top-n 1) arg))))
194
195
196(defun math-abs-approx (a)
197  (cond ((Math-negp a)
198	 (math-neg a))
199	((Math-anglep a)
200	 a)
201	((eq (car a) 'cplx)
202	 (math-add (math-abs (nth 1 a)) (math-abs (nth 2 a))))
203	((eq (car a) 'polar)
204	 (nth 1 a))
205	((eq (car a) 'sdev)
206	 (math-abs-approx (nth 1 a)))
207	((eq (car a) 'intv)
208	 (math-max (math-abs (nth 2 a)) (math-abs (nth 3 a))))
209	((eq (car a) 'date)
210	 a)
211	((eq (car a) 'vec)
212	 (math-reduce-vec 'math-add-abs-approx a))
213	((eq (car a) 'calcFunc-abs)
214	 (car a))
215	(t a)))
216
217(defun math-add-abs-approx (a b)
218  (math-add (math-abs-approx a) (math-abs-approx b)))
219
220
221;;;; Declarations.
222
223(defvar math-decls-cache-tag nil)
224(defvar math-decls-cache nil)
225(defvar math-decls-all nil)
226
227;;; Math-decls-cache is an a-list where each entry is a list of the form:
228;;;   (VAR TYPES RANGE)
229;;; where VAR is a variable name (with var- prefix) or function name;
230;;;       TYPES is a list of type symbols (any, int, frac, ...)
231;;;	  RANGE is a sorted vector of intervals describing the range.
232
233(defvar math-super-types
234  '((int numint rat real number)
235    (numint real number)
236    (frac rat real number)
237    (rat real number)
238    (float real number)
239    (real number)
240    (number)
241    (scalar)
242    (sqmatrix matrix vector)
243    (matrix vector)
244    (vector)
245    (const)))
246
247(defun math-setup-declarations ()
248  (or (eq math-decls-cache-tag (calc-var-value 'var-Decls))
249      (let ((p (calc-var-value 'var-Decls))
250	    vec type range)
251	(setq math-decls-cache-tag p
252	      math-decls-cache nil)
253	(and (eq (car-safe p) 'vec)
254	     (while (setq p (cdr p))
255	       (and (eq (car-safe (car p)) 'vec)
256		    (setq vec (nth 2 (car p)))
257		    (condition-case err
258			(let ((v (nth 1 (car p))))
259			  (setq type nil range nil)
260			  (or (eq (car-safe vec) 'vec)
261			      (setq vec (list 'vec vec)))
262			  (while (and (setq vec (cdr vec))
263				      (not (Math-objectp (car vec))))
264			    (and (eq (car-safe (car vec)) 'var)
265				 (let ((st (assq (nth 1 (car vec))
266						 math-super-types)))
267				   (cond (st (setq type (append type st)))
268					 ((eq (nth 1 (car vec)) 'pos)
269					  (setq type (append type
270							     '(real number))
271						range
272						'(intv 1 0 (var inf var-inf))))
273					 ((eq (nth 1 (car vec)) 'nonneg)
274					  (setq type (append type
275							     '(real number))
276						range
277						'(intv 3 0
278						       (var inf var-inf))))))))
279			  (if vec
280			      (setq type (append type '(real number))
281				    range (math-prepare-set (cons 'vec vec))))
282			  (setq type (list type range))
283			  (or (eq (car-safe v) 'vec)
284			      (setq v (list 'vec v)))
285			  (while (setq v (cdr v))
286			    (if (or (eq (car-safe (car v)) 'var)
287				    (not (Math-primp (car v))))
288				(setq math-decls-cache
289				      (cons (cons (if (eq (car (car v)) 'var)
290						      (nth 2 (car v))
291						    (car (car v)))
292						  type)
293					    math-decls-cache)))))
294		      (error nil)))))
295	(setq math-decls-all (assq 'var-All math-decls-cache)))))
296
297(defun math-known-scalarp (a &optional assume-scalar)
298  (math-setup-declarations)
299  (if (if calc-matrix-mode
300	  (eq calc-matrix-mode 'scalar)
301	assume-scalar)
302      (not (math-check-known-matrixp a))
303    (math-check-known-scalarp a)))
304
305(defun math-known-matrixp (a)
306  (and (not (Math-scalarp a))
307       (not (math-known-scalarp a t))))
308
309(defun math-known-square-matrixp (a)
310  (and (math-known-matrixp a)
311       (math-check-known-square-matrixp a)))
312
313;;; Try to prove that A is a scalar (i.e., a non-vector).
314(defun math-check-known-scalarp (a)
315  (cond ((Math-objectp a) t)
316	((memq (car a) math-scalar-functions)
317	 t)
318	((memq (car a) math-real-scalar-functions)
319	 t)
320	((memq (car a) math-scalar-if-args-functions)
321	 (while (and (setq a (cdr a))
322		     (math-check-known-scalarp (car a))))
323	 (null a))
324	((eq (car a) '^)
325	 (math-check-known-scalarp (nth 1 a)))
326	((math-const-var a) t)
327	(t
328	 (let ((decl (if (eq (car a) 'var)
329			 (or (assq (nth 2 a) math-decls-cache)
330			     math-decls-all)
331		       (assq (car a) math-decls-cache)))
332               val)
333           (cond
334            ((memq 'scalar (nth 1 decl))
335             t)
336            ((and (eq (car a) 'var)
337                  (symbolp (nth 2 a))
338                  (boundp (nth 2 a))
339                  (setq val (symbol-value (nth 2 a))))
340             (math-check-known-scalarp val))
341            (t
342             nil))))))
343
344;;; Try to prove that A is *not* a scalar.
345(defun math-check-known-matrixp (a)
346  (cond ((Math-objectp a) nil)
347	((memq (car a) math-nonscalar-functions)
348	 t)
349	((memq (car a) math-scalar-if-args-functions)
350	 (while (and (setq a (cdr a))
351		     (not (math-check-known-matrixp (car a)))))
352	 a)
353	((eq (car a) '^)
354	 (math-check-known-matrixp (nth 1 a)))
355	((math-const-var a) nil)
356	(t
357	 (let ((decl (if (eq (car a) 'var)
358			 (or (assq (nth 2 a) math-decls-cache)
359			     math-decls-all)
360		       (assq (car a) math-decls-cache)))
361               val)
362           (cond
363            ((memq 'matrix (nth 1 decl))
364             t)
365            ((and (eq (car a) 'var)
366                  (symbolp (nth 2 a))
367                  (boundp (nth 2 a))
368                  (setq val (symbol-value (nth 2 a))))
369             (math-check-known-matrixp val))
370            (t
371             nil))))))
372
373;;; Given that A is a matrix, try to prove that it is a square matrix.
374(defun math-check-known-square-matrixp (a)
375  (cond ((math-square-matrixp a)
376         t)
377        ((eq (car-safe a) '^)
378         (math-check-known-square-matrixp (nth 1 a)))
379        ((or
380          (eq (car-safe a) '*)
381          (eq (car-safe a) '+)
382          (eq (car-safe a) '-))
383         (and
384          (math-check-known-square-matrixp (nth 1 a))
385          (math-check-known-square-matrixp (nth 2 a))))
386        (t
387         (let ((decl (if (eq (car a) 'var)
388                         (or (assq (nth 2 a) math-decls-cache)
389                             math-decls-all)
390                       (assq (car a) math-decls-cache)))
391               val)
392           (cond
393            ((memq 'sqmatrix (nth 1 decl))
394             t)
395            ((and (eq (car a) 'var)
396                  (boundp (nth 2 a))
397                  (setq val (symbol-value (nth 2 a))))
398             (math-check-known-square-matrixp val))
399            ((and (or
400                   (integerp calc-matrix-mode)
401                   (eq calc-matrix-mode 'sqmatrix))
402                  (eq (car-safe a) 'var))
403             t)
404            ((memq 'matrix (nth 1 decl))
405             nil)
406            (t
407             nil))))))
408
409;;; Try to prove that A is a real (i.e., not complex).
410(defun math-known-realp (a)
411  (< (math-possible-signs a) 8))
412
413;;; Try to prove that A is real and positive.
414(defun math-known-posp (a)
415  (eq (math-possible-signs a) 4))
416
417;;; Try to prove that A is real and negative.
418(defun math-known-negp (a)
419  (eq (math-possible-signs a) 1))
420
421;;; Try to prove that A is real and nonnegative.
422(defun math-known-nonnegp (a)
423  (memq (math-possible-signs a) '(2 4 6)))
424
425;;; Try to prove that A is real and nonpositive.
426(defun math-known-nonposp (a)
427  (memq (math-possible-signs a) '(1 2 3)))
428
429;;; Try to prove that A is nonzero.
430(defun math-known-nonzerop (a)
431  (memq (math-possible-signs a) '(1 4 5 8 9 12 13)))
432
433;;; Return true if A is negative, or looks negative but we don't know.
434(defun math-guess-if-neg (a)
435  (let ((sgn (math-possible-signs a)))
436    (if (memq sgn '(1 3))
437	t
438      (if (memq sgn '(2 4 6))
439	  nil
440	(math-looks-negp a)))))
441
442;;; Find the possible signs of A, assuming A is a number of some kind.
443;;; Returns an integer with bits:  1  may be negative,
444;;;				   2  may be zero,
445;;;				   4  may be positive,
446;;;				   8  may be nonreal.
447
448(defun math-possible-signs (a &optional origin)
449  (cond ((Math-objectp a)
450	 (if origin (setq a (math-sub a origin)))
451	 (cond ((Math-posp a) 4)
452	       ((Math-negp a) 1)
453	       ((Math-zerop a) 2)
454	       ((eq (car a) 'intv)
455		(cond
456                 ((math-known-posp (nth 2 a)) 4)
457                 ((math-known-negp (nth 3 a)) 1)
458                 ((Math-zerop (nth 2 a)) 6)
459                 ((Math-zerop (nth 3 a)) 3)
460                 (t 7)))
461	       ((eq (car a) 'sdev)
462		(if (math-known-realp (nth 1 a)) 7 15))
463	       (t 8)))
464	((memq (car a) '(+ -))
465	 (cond ((Math-realp (nth 1 a))
466		(if (eq (car a) '-)
467		    (math-neg-signs
468		     (math-possible-signs (nth 2 a)
469					  (if origin
470					      (math-add origin (nth 1 a))
471					    (nth 1 a))))
472		  (math-possible-signs (nth 2 a)
473				       (if origin
474					   (math-sub origin (nth 1 a))
475					 (math-neg (nth 1 a))))))
476	       ((Math-realp (nth 2 a))
477		(let ((org (if (eq (car a) '-)
478			       (nth 2 a)
479			     (math-neg (nth 2 a)))))
480		  (math-possible-signs (nth 1 a)
481				       (if origin
482					   (math-add origin org)
483					 org))))
484	       (t
485		(let ((s1 (math-possible-signs (nth 1 a) origin))
486		      (s2 (math-possible-signs (nth 2 a))))
487		  (if (eq (car a) '-) (setq s2 (math-neg-signs s2)))
488		  (cond ((eq s1 s2) s1)
489			((eq s1 2) s2)
490			((eq s2 2) s1)
491			((>= s1 8) 15)
492			((>= s2 8) 15)
493			((and (eq s1 4) (eq s2 6)) 4)
494			((and (eq s2 4) (eq s1 6)) 4)
495			((and (eq s1 1) (eq s2 3)) 1)
496			((and (eq s2 1) (eq s1 3)) 1)
497			(t 7))))))
498	((eq (car a) 'neg)
499	 (math-neg-signs (math-possible-signs
500			  (nth 1 a)
501			  (and origin (math-neg origin)))))
502	((and origin (Math-zerop origin) (setq origin nil)
503	      nil))
504	((and (or (eq (car a) '*)
505		  (and (eq (car a) '/) origin))
506	      (Math-realp (nth 1 a)))
507	 (let ((s (if (eq (car a) '*)
508		      (if (Math-zerop (nth 1 a))
509			  (math-possible-signs 0 origin)
510			(math-possible-signs (nth 2 a)
511					     (math-div (or origin 0)
512						       (nth 1 a))))
513		    (math-neg-signs
514		     (math-possible-signs (nth 2 a)
515					  (math-div (nth 1 a)
516						    origin))))))
517	   (if (Math-negp (nth 1 a)) (math-neg-signs s) s)))
518	((and (memq (car a) '(* /)) (Math-realp (nth 2 a)))
519	 (let ((s (math-possible-signs (nth 1 a)
520				       (if (eq (car a) '*)
521					   (math-mul (or origin 0) (nth 2 a))
522					 (math-div (or origin 0) (nth 2 a))))))
523	   (if (Math-negp (nth 2 a)) (math-neg-signs s) s)))
524	((eq (car a) 'vec)
525	 (let ((signs 0))
526	   (while (and (setq a (cdr a)) (< signs 15))
527	     (setq signs (logior signs (math-possible-signs
528					(car a) origin))))
529	   signs))
530	(t (let ((sign
531		  (cond
532		   ((memq (car a) '(* /))
533		    (let ((s1 (math-possible-signs (nth 1 a)))
534			  (s2 (math-possible-signs (nth 2 a))))
535		      (cond ((>= s1 8) 15)
536			    ((>= s2 8) 15)
537			    ((and (eq (car a) '/) (memq s2 '(2 3 6 7))) 15)
538			    (t
539			     (logior (if (memq s1 '(4 5 6 7)) s2 0)
540				     (if (memq s1 '(2 3 6 7)) 2 0)
541				     (if (memq s1 '(1 3 5 7))
542					 (math-neg-signs s2) 0))))))
543		   ((eq (car a) '^)
544		    (let ((s1 (math-possible-signs (nth 1 a)))
545			  (s2 (math-possible-signs (nth 2 a))))
546		      (cond ((>= s1 8) 15)
547			    ((>= s2 8) 15)
548			    ((eq s1 4) 4)
549			    ((eq s1 2) (if (eq s2 4) 2 15))
550			    ((eq s2 2) (if (memq s1 '(1 5)) 2 15))
551			    ((Math-integerp (nth 2 a))
552			     (if (math-evenp (nth 2 a))
553				 (if (memq s1 '(3 6 7)) 6 4)
554			       s1))
555			    ((eq s1 6) (if (eq s2 4) 6 15))
556			    (t 7))))
557		   ((eq (car a) '%)
558		    (let ((s2 (math-possible-signs (nth 2 a))))
559		      (cond ((>= s2 8) 7)
560			    ((eq s2 2) 2)
561			    ((memq s2 '(4 6)) 6)
562			    ((memq s2 '(1 3)) 3)
563			    (t 7))))
564		   ((and (memq (car a) '(calcFunc-abs calcFunc-abssqr))
565			 (= (length a) 2))
566		    (let ((s1 (math-possible-signs (nth 1 a))))
567		      (cond ((eq s1 2) 2)
568			    ((memq s1 '(1 4 5)) 4)
569			    (t 6))))
570		   ((and (eq (car a) 'calcFunc-exp) (= (length a) 2))
571		    (let ((s1 (math-possible-signs (nth 1 a))))
572		      (if (>= s1 8)
573			  15
574			(if (or (not origin) (math-negp origin))
575			    4
576			  (setq origin (math-sub (or origin 0) 1))
577			  (if (Math-zerop origin) (setq origin nil))
578			  s1))))
579		   ((or (and (memq (car a) '(calcFunc-ln calcFunc-log10))
580			     (= (length a) 2))
581			(and (eq (car a) 'calcFunc-log)
582			     (= (length a) 3)
583			     (math-known-posp (nth 2 a))))
584		    (if (math-known-nonnegp (nth 1 a))
585			(math-possible-signs (nth 1 a) 1)
586		      15))
587		   ((and (eq (car a) 'calcFunc-sqrt) (= (length a) 2))
588		    (let ((s1 (math-possible-signs (nth 1 a))))
589		      (if (memq s1 '(2 4 6)) s1 15)))
590		   ((memq (car a) math-nonnegative-functions) 6)
591		   ((memq (car a) math-positive-functions) 4)
592		   ((memq (car a) math-real-functions) 7)
593		   ((memq (car a) math-real-scalar-functions) 7)
594		   ((and (memq (car a) math-real-if-arg-functions)
595			 (= (length a) 2))
596		    (if (math-known-realp (nth 1 a)) 7 15)))))
597	     (cond (sign
598		    (if origin
599			(+ (logand sign 8)
600			   (if (Math-posp origin)
601			       (if (memq sign '(1 2 3 8 9 10 11)) 1 7)
602			     (if (memq sign '(2 4 6 8 10 12 14)) 4 7)))
603		      sign))
604		   ((math-const-var a)
605		    (cond ((eq (nth 2 a) 'var-pi)
606			   (if origin
607			       (math-possible-signs (math-pi) origin)
608			     4))
609			  ((eq (nth 2 a) 'var-e)
610			   (if origin
611			       (math-possible-signs (math-e) origin)
612			     4))
613			  ((eq (nth 2 a) 'var-inf) 4)
614			  ((eq (nth 2 a) 'var-uinf) 13)
615			  ((eq (nth 2 a) 'var-i) 8)
616			  (t 15)))
617		   (t
618		    (math-setup-declarations)
619		    (let ((decl (if (eq (car a) 'var)
620				    (or (assq (nth 2 a) math-decls-cache)
621					math-decls-all)
622				  (assq (car a) math-decls-cache))))
623		      (if (and origin
624			       (memq 'int (nth 1 decl))
625			       (not (Math-num-integerp origin)))
626			  5
627			(if (nth 2 decl)
628			    (math-possible-signs (nth 2 decl) origin)
629			  (if (memq 'real (nth 1 decl))
630			      7
631			    15))))))))))
632
633(defun math-neg-signs (s1)
634  (if (>= s1 8)
635      (+ 8 (math-neg-signs (- s1 8)))
636    (+ (if (memq s1 '(1 3 5 7)) 4 0)
637       (if (memq s1 '(2 3 6 7)) 2 0)
638       (if (memq s1 '(4 5 6 7)) 1 0))))
639
640
641;;; Try to prove that A is an integer.
642(defun math-known-integerp (a)
643  (eq (math-possible-types a) 1))
644
645(defun math-known-num-integerp (a)
646  (<= (math-possible-types a t) 3))
647
648(defun math-known-imagp (a)
649  (= (math-possible-types a) 16))
650
651
652;;; Find the possible types of A.
653;;; Returns an integer with bits:  1  may be integer.
654;;;				   2  may be integer-valued float.
655;;;				   4  may be fraction.
656;;;				   8  may be non-integer-valued float.
657;;;				  16  may be imaginary.
658;;;				  32  may be non-real, non-imaginary.
659;;; Real infinities count as integers for the purposes of this function.
660(defun math-possible-types (a &optional num)
661  (cond ((Math-objectp a)
662	 (cond ((Math-integerp a) (if num 3 1))
663	       ((Math-messy-integerp a) (if num 3 2))
664	       ((eq (car a) 'frac) (if num 12 4))
665	       ((eq (car a) 'float) (if num 12 8))
666	       ((eq (car a) 'intv)
667		(if (equal (nth 2 a) (nth 3 a))
668		    (math-possible-types (nth 2 a))
669		  15))
670	       ((eq (car a) 'sdev)
671		(if (math-known-realp (nth 1 a)) 15 63))
672	       ((eq (car a) 'cplx)
673		(if (math-zerop (nth 1 a)) 16 32))
674	       ((eq (car a) 'polar)
675		(if (or (Math-equal (nth 2 a) (math-quarter-circle nil))
676			(Math-equal (nth 2 a)
677				    (math-neg (math-quarter-circle nil))))
678		    16 48))
679	       (t 63)))
680	((eq (car a) '/)
681	 (let* ((t1 (math-possible-types (nth 1 a) num))
682		(t2 (math-possible-types (nth 2 a) num))
683		(t12 (logior t1 t2)))
684	   (if (< t12 16)
685	       (if (> (logand t12 10) 0)
686		   10
687		 (if (or (= t1 4) (= t2 4) calc-prefer-frac)
688		     5
689		   15))
690	     (if (< t12 32)
691		 (if (= t1 16)
692		     (if (= t2 16) 15
693		       (if (< t2 16) 16 31))
694		   (if (= t2 16)
695		       (if (< t1 16) 16 31)
696		     31))
697	       63))))
698	((memq (car a) '(+ - * %))
699	 (let* ((t1 (math-possible-types (nth 1 a) num))
700		(t2 (math-possible-types (nth 2 a) num))
701		(t12 (logior t1 t2)))
702	   (if (eq (car a) '%)
703	       (setq t1 (logand t1 15) t2 (logand t2 15) t12 (logand t12 15)))
704	   (if (< t12 16)
705	       (let ((mask (if (<= t12 3)
706			       1
707			     (if (and (or (and (<= t1 3) (= (logand t2 3) 0))
708					  (and (<= t2 3) (= (logand t1 3) 0)))
709				      (memq (car a) '(+ -)))
710				 4
711			       5))))
712		 (if num
713		     (* mask 3)
714		   (logior (if (and (> (logand t1 5) 0) (> (logand t2 5) 0))
715			       mask 0)
716			   (if (> (logand t12 10) 0)
717			       (* mask 2) 0))))
718	     (if (< t12 32)
719		 (if (eq (car a) '*)
720		     (if (= t1 16)
721			 (if (= t2 16) 15
722			   (if (< t2 16) 16 31))
723		       (if (= t2 16)
724			   (if (< t1 16) 16 31)
725			 31))
726		   (if (= t12 16) 16
727		     (if (or (and (= t1 16) (< t2 16))
728			     (and (= t2 16) (< t1 16))) 32 63)))
729	       63))))
730	((eq (car a) 'neg)
731	 (math-possible-types (nth 1 a)))
732	((eq (car a) '^)
733	 (let* ((t1 (math-possible-types (nth 1 a) num))
734		(t2 (math-possible-types (nth 2 a) num))
735		(t12 (logior t1 t2)))
736	   (if (and (<= t2 3) (math-known-nonnegp (nth 2 a)) (< t1 16))
737	       (let ((mask (logior (if (> (logand t1 3) 0) 1 0)
738				   (logand t1 4)
739				   (if (> (logand t1 12) 0) 5 0))))
740		 (if num
741		     (* mask 3)
742		   (logior (if (and (> (logand t1 5) 0) (> (logand t2 5) 0))
743			       mask 0)
744			   (if (> (logand t12 10) 0)
745			       (* mask 2) 0))))
746	     (if (and (math-known-nonnegp (nth 1 a))
747		      (math-known-posp (nth 2 a)))
748		 15
749	       63))))
750	((eq (car a) 'calcFunc-sqrt)
751	 (let ((t1 (math-possible-signs (nth 1 a))))
752	   (logior (if (> (logand t1 2) 0) 3 0)
753		   (if (> (logand t1 1) 0) 16 0)
754		   (if (> (logand t1 4) 0) 15 0)
755		   (if (> (logand t1 8) 0) 32 0))))
756	((eq (car a) 'vec)
757	 (let ((types 0))
758	   (while (and (setq a (cdr a)) (< types 63))
759	     (setq types (logior types (math-possible-types (car a) t))))
760	   types))
761	((or (memq (car a) math-integer-functions)
762	     (and (memq (car a) math-rounding-functions)
763		  (math-known-nonnegp (or (nth 2 a) 0))))
764	 1)
765	((or (memq (car a) math-num-integer-functions)
766	     (and (memq (car a) math-float-rounding-functions)
767		  (math-known-nonnegp (or (nth 2 a) 0))))
768	 2)
769	((eq (car a) 'calcFunc-frac)
770	 5)
771	((and (eq (car a) 'calcFunc-float) (= (length a) 2))
772	 (let ((t1 (math-possible-types (nth 1 a))))
773	   (logior (if (> (logand t1 3) 0) 2 0)
774		   (if (> (logand t1 12) 0) 8 0)
775		   (logand t1 48))))
776	((and (memq (car a) '(calcFunc-abs calcFunc-abssqr))
777	      (= (length a) 2))
778	 (let ((t1 (math-possible-types (nth 1 a))))
779	   (if (>= t1 16)
780	       15
781	     t1)))
782	((math-const-var a)
783	 (cond ((memq (nth 2 a) '(var-e var-pi var-phi var-gamma)) 8)
784	       ((eq (nth 2 a) 'var-inf) 1)
785	       ((eq (nth 2 a) 'var-i) 16)
786	       (t 63)))
787	(t
788	 (math-setup-declarations)
789	 (let ((decl (if (eq (car a) 'var)
790			 (or (assq (nth 2 a) math-decls-cache)
791			     math-decls-all)
792		       (assq (car a) math-decls-cache))))
793	   (cond ((memq 'int (nth 1 decl))
794		  1)
795		 ((memq 'numint (nth 1 decl))
796		  3)
797		 ((memq 'frac (nth 1 decl))
798		  4)
799		 ((memq 'rat (nth 1 decl))
800		  5)
801		 ((memq 'float (nth 1 decl))
802		  10)
803		 ((nth 2 decl)
804		  (math-possible-types (nth 2 decl)))
805		 ((memq 'real (nth 1 decl))
806		  15)
807		 (t 63))))))
808
809(defun math-known-evenp (a)
810  (cond ((Math-integerp a)
811	 (math-evenp a))
812	((Math-messy-integerp a)
813	 (or (> (nth 2 a) 0)
814	     (math-evenp (math-trunc a))))
815	((eq (car a) '*)
816	 (if (math-known-evenp (nth 1 a))
817	     (math-known-num-integerp (nth 2 a))
818	   (if (math-known-num-integerp (nth 1 a))
819	       (math-known-evenp (nth 2 a)))))
820	((memq (car a) '(+ -))
821	 (or (and (math-known-evenp (nth 1 a))
822		  (math-known-evenp (nth 2 a)))
823	     (and (math-known-oddp (nth 1 a))
824		  (math-known-oddp (nth 2 a)))))
825	((eq (car a) 'neg)
826	 (math-known-evenp (nth 1 a)))))
827
828(defun math-known-oddp (a)
829  (cond ((Math-integerp a)
830	 (math-oddp a))
831	((Math-messy-integerp a)
832	 (and (<= (nth 2 a) 0)
833	      (math-oddp (math-trunc a))))
834	((memq (car a) '(+ -))
835	 (or (and (math-known-evenp (nth 1 a))
836		  (math-known-oddp (nth 2 a)))
837	     (and (math-known-oddp (nth 1 a))
838		  (math-known-evenp (nth 2 a)))))
839	((eq (car a) 'neg)
840	 (math-known-oddp (nth 1 a)))))
841
842
843(defun calcFunc-dreal (expr)
844  (let ((types (math-possible-types expr)))
845    (if (< types 16) 1
846      (if (= (logand types 15) 0) 0
847	(math-reject-arg expr 'realp 'quiet)))))
848
849(defun calcFunc-dimag (expr)
850  (let ((types (math-possible-types expr)))
851    (if (= types 16) 1
852      (if (= (logand types 16) 0) 0
853	(math-reject-arg expr "Expected an imaginary number")))))
854
855(defun calcFunc-dpos (expr)
856  (let ((signs (math-possible-signs expr)))
857    (if (eq signs 4) 1
858      (if (memq signs '(1 2 3)) 0
859	(math-reject-arg expr 'posp 'quiet)))))
860
861(defun calcFunc-dneg (expr)
862  (let ((signs (math-possible-signs expr)))
863    (if (eq signs 1) 1
864      (if (memq signs '(2 4 6)) 0
865	(math-reject-arg expr 'negp 'quiet)))))
866
867(defun calcFunc-dnonneg (expr)
868  (let ((signs (math-possible-signs expr)))
869    (if (memq signs '(2 4 6)) 1
870      (if (eq signs 1) 0
871	(math-reject-arg expr 'posp 'quiet)))))
872
873(defun calcFunc-dnonzero (expr)
874  (let ((signs (math-possible-signs expr)))
875    (if (memq signs '(1 4 5 8 9 12 13)) 1
876      (if (eq signs 2) 0
877	(math-reject-arg expr 'nonzerop 'quiet)))))
878
879(defun calcFunc-dint (expr)
880  (let ((types (math-possible-types expr)))
881    (if (= types 1) 1
882      (if (= (logand types 1) 0) 0
883	(math-reject-arg expr 'integerp 'quiet)))))
884
885(defun calcFunc-dnumint (expr)
886  (let ((types (math-possible-types expr t)))
887    (if (<= types 3) 1
888      (if (= (logand types 3) 0) 0
889	(math-reject-arg expr 'integerp 'quiet)))))
890
891(defun calcFunc-dnatnum (expr)
892  (let ((res (calcFunc-dint expr)))
893    (if (eq res 1)
894	(calcFunc-dnonneg expr)
895      res)))
896
897(defun calcFunc-deven (expr)
898  (if (math-known-evenp expr)
899      1
900    (if (or (math-known-oddp expr)
901	    (= (logand (math-possible-types expr) 3) 0))
902	0
903      (math-reject-arg expr "Can't tell if expression is odd or even"))))
904
905(defun calcFunc-dodd (expr)
906  (if (math-known-oddp expr)
907      1
908    (if (or (math-known-evenp expr)
909	    (= (logand (math-possible-types expr) 3) 0))
910	0
911      (math-reject-arg expr "Can't tell if expression is odd or even"))))
912
913(defun calcFunc-drat (expr)
914  (let ((types (math-possible-types expr)))
915    (if (memq types '(1 4 5)) 1
916      (if (= (logand types 5) 0) 0
917	(math-reject-arg expr "Rational number expected")))))
918
919(defun calcFunc-drange (expr)
920  (math-setup-declarations)
921  (let (range)
922    (if (Math-realp expr)
923	(list 'vec expr)
924      (if (eq (car-safe expr) 'intv)
925	  expr
926	(if (eq (car-safe expr) 'var)
927	    (setq range (nth 2 (or (assq (nth 2 expr) math-decls-cache)
928				   math-decls-all)))
929	  (setq range (nth 2 (assq (car-safe expr) math-decls-cache))))
930	(if range
931	    (math-clean-set (copy-sequence range))
932	  (setq range (math-possible-signs expr))
933	  (if (< range 8)
934	      (aref [(vec)
935		     (intv 2 (neg (var inf var-inf)) 0)
936		     (vec 0)
937		     (intv 3 (neg (var inf var-inf)) 0)
938		     (intv 1 0 (var inf var-inf))
939		     (vec (intv 2 (neg (var inf var-inf)) 0)
940			  (intv 1 0 (var inf var-inf)))
941		     (intv 3 0 (var inf var-inf))
942		     (intv 3 (neg (var inf var-inf)) (var inf var-inf))] range)
943	    (math-reject-arg expr 'realp 'quiet)))))))
944
945(defun calcFunc-dscalar (a)
946  (if (math-known-scalarp a) 1
947    (if (math-known-matrixp a) 0
948      (math-reject-arg a 'objectp 'quiet))))
949
950
951;;;; Arithmetic.
952
953(defsubst calcFunc-neg (a)
954  (math-normalize (list 'neg a)))
955
956(defun math-neg-fancy (a)
957  (cond ((eq (car a) 'polar)
958	 (list 'polar
959	       (nth 1 a)
960	       (if (math-posp (nth 2 a))
961		   (math-sub (nth 2 a) (math-half-circle nil))
962		 (math-add (nth 2 a) (math-half-circle nil)))))
963	((eq (car a) 'mod)
964	 (if (math-zerop (nth 1 a))
965	     a
966	   (list 'mod (math-sub (nth 2 a) (nth 1 a)) (nth 2 a))))
967	((eq (car a) 'sdev)
968	 (list 'sdev (math-neg (nth 1 a)) (nth 2 a)))
969	((eq (car a) 'intv)
970	 (math-make-intv (aref [0 2 1 3] (nth 1 a))
971			 (math-neg (nth 3 a))
972			 (math-neg (nth 2 a))))
973	((and math-simplify-only
974	      (not (equal a math-simplify-only)))
975	 (list 'neg a))
976	((eq (car a) '+)
977	 (math-sub (math-neg (nth 1 a)) (nth 2 a)))
978	((eq (car a) '-)
979	 (math-sub (nth 2 a) (nth 1 a)))
980	((and (memq (car a) '(* /))
981	      (math-okay-neg (nth 1 a)))
982	 (list (car a) (math-neg (nth 1 a)) (nth 2 a)))
983	((and (memq (car a) '(* /))
984	      (math-okay-neg (nth 2 a)))
985	 (list (car a) (nth 1 a) (math-neg (nth 2 a))))
986	((and (memq (car a) '(* /))
987	      (or (math-objectp (nth 1 a))
988		  (and (eq (car (nth 1 a)) '*)
989		       (math-objectp (nth 1 (nth 1 a))))))
990	 (list (car a) (math-neg (nth 1 a)) (nth 2 a)))
991	((and (eq (car a) '/)
992	      (or (math-objectp (nth 2 a))
993		  (and (eq (car (nth 2 a)) '*)
994		       (math-objectp (nth 1 (nth 2 a))))))
995	 (list (car a) (nth 1 a) (math-neg (nth 2 a))))
996	((and (eq (car a) 'var) (memq (nth 2 a) '(var-uinf var-nan)))
997	 a)
998	((eq (car a) 'neg)
999	 (nth 1 a))
1000	(t (list 'neg a))))
1001
1002(defun math-okay-neg (a)
1003  (or (math-looks-negp a)
1004      (eq (car-safe a) '-)))
1005
1006(defun math-neg-float (a)
1007  (list 'float (Math-integer-neg (nth 1 a)) (nth 2 a)))
1008
1009
1010(defun calcFunc-add (&rest rest)
1011  (if rest
1012      (let ((a (car rest)))
1013	(while (setq rest (cdr rest))
1014	  (setq a (list '+ a (car rest))))
1015	(math-normalize a))
1016    0))
1017
1018(defun calcFunc-sub (&rest rest)
1019  (if rest
1020      (let ((a (car rest)))
1021	(while (setq rest (cdr rest))
1022	  (setq a (list '- a (car rest))))
1023	(math-normalize a))
1024    0))
1025
1026(defun math-add-objects-fancy (a b)
1027  (cond ((and (Math-numberp a) (Math-numberp b))
1028	 (let ((aa (math-complex a))
1029	       (bb (math-complex b)))
1030	   (math-normalize
1031	    (let ((res (list 'cplx
1032			     (math-add (nth 1 aa) (nth 1 bb))
1033			     (math-add (nth 2 aa) (nth 2 bb)))))
1034	      (if (math-want-polar a b)
1035		  (math-polar res)
1036		res)))))
1037	((or (Math-vectorp a) (Math-vectorp b))
1038	 (math-map-vec-2 'math-add a b))
1039	((eq (car-safe a) 'sdev)
1040	 (if (eq (car-safe b) 'sdev)
1041	     (math-make-sdev (math-add (nth 1 a) (nth 1 b))
1042			     (math-hypot (nth 2 a) (nth 2 b)))
1043	   (and (or (Math-scalarp b)
1044		    (not (Math-objvecp b)))
1045		(math-make-sdev (math-add (nth 1 a) b) (nth 2 a)))))
1046	((and (eq (car-safe b) 'sdev)
1047	      (or (Math-scalarp a)
1048		  (not (Math-objvecp a))))
1049	 (math-make-sdev (math-add a (nth 1 b)) (nth 2 b)))
1050	((eq (car-safe a) 'intv)
1051	 (if (eq (car-safe b) 'intv)
1052	     (math-make-intv (logior (logand (nth 1 a) (nth 1 b))
1053				     (if (equal (nth 2 a)
1054						'(neg (var inf var-inf)))
1055					 (logand (nth 1 a) 2) 0)
1056				     (if (equal (nth 2 b)
1057						'(neg (var inf var-inf)))
1058					 (logand (nth 1 b) 2) 0)
1059				     (if (equal (nth 3 a) '(var inf var-inf))
1060					 (logand (nth 1 a) 1) 0)
1061				     (if (equal (nth 3 b) '(var inf var-inf))
1062					 (logand (nth 1 b) 1) 0))
1063			     (math-add (nth 2 a) (nth 2 b))
1064			     (math-add (nth 3 a) (nth 3 b)))
1065	   (and (or (Math-anglep b)
1066		    (eq (car b) 'date)
1067		    (not (Math-objvecp b)))
1068		(math-make-intv (nth 1 a)
1069				(math-add (nth 2 a) b)
1070				(math-add (nth 3 a) b)))))
1071	((and (eq (car-safe b) 'intv)
1072	      (or (Math-anglep a)
1073		  (eq (car a) 'date)
1074		  (not (Math-objvecp a))))
1075	 (math-make-intv (nth 1 b)
1076			 (math-add a (nth 2 b))
1077			 (math-add a (nth 3 b))))
1078	((eq (car-safe a) 'date)
1079	 (cond ((eq (car-safe b) 'date)
1080		(math-add (nth 1 a) (nth 1 b)))
1081	       ((eq (car-safe b) 'hms)
1082		(let ((parts (math-date-parts (nth 1 a))))
1083		  (list 'date
1084			(math-add (car parts)   ; this minimizes roundoff
1085				  (math-div (math-add
1086					     (math-add (nth 1 parts)
1087						       (nth 2 parts))
1088					     (math-add
1089					      (math-mul (nth 1 b) 3600)
1090					      (math-add (math-mul (nth 2 b) 60)
1091							(nth 3 b))))
1092					    86400)))))
1093	       ((Math-realp b)
1094		(list 'date (math-add (nth 1 a) b)))
1095	       (t nil)))
1096	((eq (car-safe b) 'date)
1097	 (math-add-objects-fancy b a))
1098	((and (eq (car-safe a) 'mod)
1099	      (eq (car-safe b) 'mod)
1100	      (equal (nth 2 a) (nth 2 b)))
1101	 (math-make-mod (math-add (nth 1 a) (nth 1 b)) (nth 2 a)))
1102	((and (eq (car-safe a) 'mod)
1103	      (Math-anglep b))
1104	 (math-make-mod (math-add (nth 1 a) b) (nth 2 a)))
1105	((and (eq (car-safe b) 'mod)
1106	      (Math-anglep a))
1107	 (math-make-mod (math-add a (nth 1 b)) (nth 2 b)))
1108	((and (or (eq (car-safe a) 'hms) (eq (car-safe b) 'hms))
1109	      (and (Math-anglep a) (Math-anglep b)))
1110	 (or (eq (car-safe a) 'hms) (setq a (math-to-hms a)))
1111	 (or (eq (car-safe b) 'hms) (setq b (math-to-hms b)))
1112	 (math-normalize
1113	  (if (math-negp a)
1114	      (math-neg (math-add (math-neg a) (math-neg b)))
1115	    (if (math-negp b)
1116		(let* ((s (math-add (nth 3 a) (nth 3 b)))
1117		       (m (math-add (nth 2 a) (nth 2 b)))
1118		       (h (math-add (nth 1 a) (nth 1 b))))
1119		  (if (math-negp s)
1120		      (setq s (math-add s 60)
1121			    m (math-add m -1)))
1122		  (if (math-negp m)
1123		      (setq m (math-add m 60)
1124			    h (math-add h -1)))
1125		  (if (math-negp h)
1126		      (math-add b a)
1127		    (list 'hms h m s)))
1128	      (let* ((s (math-add (nth 3 a) (nth 3 b)))
1129		     (m (math-add (nth 2 a) (nth 2 b)))
1130		     (h (math-add (nth 1 a) (nth 1 b))))
1131		(list 'hms h m s))))))
1132	(t (calc-record-why "*Incompatible arguments for +" a b))))
1133
1134(defun math-add-symb-fancy (a b)
1135  (or (and math-simplify-only
1136	   (not (equal a math-simplify-only))
1137	   (list '+ a b))
1138      (and (eq (car-safe b) '+)
1139	   (math-add (math-add a (nth 1 b))
1140		     (nth 2 b)))
1141      (and (eq (car-safe b) '-)
1142	   (math-sub (math-add a (nth 1 b))
1143		     (nth 2 b)))
1144      (and (eq (car-safe b) 'neg)
1145	   (eq (car-safe (nth 1 b)) '+)
1146	   (math-sub (math-sub a (nth 1 (nth 1 b)))
1147		     (nth 2 (nth 1 b))))
1148      (and (or (and (Math-vectorp a) (math-known-scalarp b))
1149	       (and (Math-vectorp b) (math-known-scalarp a)))
1150	   (math-map-vec-2 'math-add a b))
1151      (let ((inf (math-infinitep a)))
1152	(cond
1153	 (inf
1154	  (let ((inf2 (math-infinitep b)))
1155	    (if inf2
1156		(if (or (memq (nth 2 inf) '(var-uinf var-nan))
1157			(memq (nth 2 inf2) '(var-uinf var-nan)))
1158		    '(var nan var-nan)
1159		  (let ((dir (math-infinite-dir a inf))
1160			(dir2 (math-infinite-dir b inf2)))
1161		    (if (and (Math-objectp dir) (Math-objectp dir2))
1162			(if (Math-equal dir dir2)
1163			    a
1164			  '(var nan var-nan)))))
1165	      (if (and (equal a '(var inf var-inf))
1166		       (eq (car-safe b) 'intv)
1167		       (memq (nth 1 b) '(2 3))
1168		       (equal (nth 2 b) '(neg (var inf var-inf))))
1169		  (list 'intv 3 (nth 2 b) a)
1170		(if (and (equal a '(neg (var inf var-inf)))
1171			 (eq (car-safe b) 'intv)
1172			 (memq (nth 1 b) '(1 3))
1173			 (equal (nth 3 b) '(var inf var-inf)))
1174		    (list 'intv 3 a (nth 3 b))
1175		  a)))))
1176	 ((math-infinitep b)
1177	  (if (eq (car-safe a) 'intv)
1178	      (math-add b a)
1179	    b))
1180	 ((eq (car-safe a) '+)
1181	  (let ((temp (math-combine-sum (nth 2 a) b nil nil t)))
1182	    (and temp
1183		 (math-add (nth 1 a) temp))))
1184	 ((eq (car-safe a) '-)
1185	  (let ((temp (math-combine-sum (nth 2 a) b t nil t)))
1186	    (and temp
1187		 (math-add (nth 1 a) temp))))
1188	 ((and (Math-objectp a) (Math-objectp b))
1189	  nil)
1190	 (t
1191	  (math-combine-sum a b nil nil nil))))
1192      (and (Math-looks-negp b)
1193	   (list '- a (math-neg b)))
1194      (and (Math-looks-negp a)
1195	   (list '- b (math-neg a)))
1196      (and (eq (car-safe a) 'calcFunc-idn)
1197	   (= (length a) 2)
1198	   (or (and (eq (car-safe b) 'calcFunc-idn)
1199		    (= (length b) 2)
1200		    (list 'calcFunc-idn (math-add (nth 1 a) (nth 1 b))))
1201	       (and (math-square-matrixp b)
1202		    (math-add (math-mimic-ident (nth 1 a) b) b))
1203	       (and (math-known-scalarp b)
1204		    (math-add (nth 1 a) b))))
1205      (and (eq (car-safe b) 'calcFunc-idn)
1206	   (= (length b) 2)
1207	   (or (and (math-square-matrixp a)
1208		    (math-add a (math-mimic-ident (nth 1 b) a)))
1209	       (and (math-known-scalarp a)
1210		    (math-add a (nth 1 b)))))
1211      (list '+ a b)))
1212
1213
1214(defun calcFunc-mul (&rest rest)
1215  (if rest
1216      (let ((a (car rest)))
1217	(while (setq rest (cdr rest))
1218	  (setq a (list '* a (car rest))))
1219	(math-normalize a))
1220    1))
1221
1222(defun math-mul-objects-fancy (a b)
1223  (cond ((and (Math-numberp a) (Math-numberp b))
1224	 (math-normalize
1225	  (if (math-want-polar a b)
1226	      (let ((a (math-polar a))
1227		    (b (math-polar b)))
1228		(list 'polar
1229		      (math-mul (nth 1 a) (nth 1 b))
1230		      (math-fix-circular (math-add (nth 2 a) (nth 2 b)))))
1231	    (setq a (math-complex a)
1232		  b (math-complex b))
1233	    (list 'cplx
1234		  (math-sub (math-mul (nth 1 a) (nth 1 b))
1235			    (math-mul (nth 2 a) (nth 2 b)))
1236		  (math-add (math-mul (nth 1 a) (nth 2 b))
1237			    (math-mul (nth 2 a) (nth 1 b)))))))
1238	((Math-vectorp a)
1239	 (if (Math-vectorp b)
1240	     (if (math-matrixp a)
1241		 (if (math-matrixp b)
1242		     (if (= (length (nth 1 a)) (length b))
1243			 (math-mul-mats a b)
1244		       (math-dimension-error))
1245		   (if (= (length (nth 1 a)) 2)
1246		       (if (= (length a) (length b))
1247			   (math-mul-mats a (list 'vec b))
1248			 (math-dimension-error))
1249		     (if (= (length (nth 1 a)) (length b))
1250			 (math-mul-mat-vec a b)
1251		       (math-dimension-error))))
1252	       (if (math-matrixp b)
1253		   (if (= (length a) (length b))
1254		       (nth 1 (math-mul-mats (list 'vec a) b))
1255		     (math-dimension-error))
1256		 (if (= (length a) (length b))
1257		     (math-dot-product a b)
1258		   (math-dimension-error))))
1259	   (math-map-vec-2 'math-mul a b)))
1260	((Math-vectorp b)
1261	 (math-map-vec-2 'math-mul a b))
1262	((eq (car-safe a) 'sdev)
1263	 (if (eq (car-safe b) 'sdev)
1264	     (math-make-sdev (math-mul (nth 1 a) (nth 1 b))
1265			     (math-hypot (math-mul (nth 2 a) (nth 1 b))
1266					 (math-mul (nth 2 b) (nth 1 a))))
1267	   (and (or (Math-scalarp b)
1268		    (not (Math-objvecp b)))
1269		(math-make-sdev (math-mul (nth 1 a) b)
1270				(math-mul (nth 2 a) b)))))
1271	((and (eq (car-safe b) 'sdev)
1272	      (or (Math-scalarp a)
1273		  (not (Math-objvecp a))))
1274	 (math-make-sdev (math-mul a (nth 1 b)) (math-mul a (nth 2 b))))
1275	((and (eq (car-safe a) 'intv) (Math-anglep b))
1276	 (if (Math-negp b)
1277	     (math-neg (math-mul a (math-neg b)))
1278	   (math-make-intv (nth 1 a)
1279			   (math-mul (nth 2 a) b)
1280			   (math-mul (nth 3 a) b))))
1281	((and (eq (car-safe b) 'intv) (Math-anglep a))
1282	 (math-mul b a))
1283	((and (eq (car-safe a) 'intv) (math-intv-constp a)
1284	      (eq (car-safe b) 'intv) (math-intv-constp b))
1285	 (let ((lo (math-mul a (nth 2 b)))
1286	       (hi (math-mul a (nth 3 b))))
1287	   (or (eq (car-safe lo) 'intv)
1288	       (setq lo (list 'intv (if (memq (nth 1 b) '(2 3)) 3 0) lo lo)))
1289	   (or (eq (car-safe hi) 'intv)
1290	       (setq hi (list 'intv (if (memq (nth 1 b) '(1 3)) 3 0) hi hi)))
1291	   (math-combine-intervals
1292	    (nth 2 lo) (and (or (memq (nth 1 b) '(2 3))
1293				(math-infinitep (nth 2 lo)))
1294			    (memq (nth 1 lo) '(2 3)))
1295	    (nth 3 lo) (and (or (memq (nth 1 b) '(2 3))
1296				(math-infinitep (nth 3 lo)))
1297			    (memq (nth 1 lo) '(1 3)))
1298	    (nth 2 hi) (and (or (memq (nth 1 b) '(1 3))
1299				(math-infinitep (nth 2 hi)))
1300			    (memq (nth 1 hi) '(2 3)))
1301	    (nth 3 hi) (and (or (memq (nth 1 b) '(1 3))
1302				(math-infinitep (nth 3 hi)))
1303			    (memq (nth 1 hi) '(1 3))))))
1304	((and (eq (car-safe a) 'mod)
1305	      (eq (car-safe b) 'mod)
1306	      (equal (nth 2 a) (nth 2 b)))
1307	 (math-make-mod (math-mul (nth 1 a) (nth 1 b)) (nth 2 a)))
1308	((and (eq (car-safe a) 'mod)
1309	      (Math-anglep b))
1310	 (math-make-mod (math-mul (nth 1 a) b) (nth 2 a)))
1311	((and (eq (car-safe b) 'mod)
1312	      (Math-anglep a))
1313	 (math-make-mod (math-mul a (nth 1 b)) (nth 2 b)))
1314	((and (eq (car-safe a) 'hms) (Math-realp b))
1315	 (math-with-extra-prec 2
1316	   (math-to-hms (math-mul (math-from-hms a 'deg) b) 'deg)))
1317	((and (eq (car-safe b) 'hms) (Math-realp a))
1318	 (math-mul b a))
1319	(t (calc-record-why "*Incompatible arguments for *" a b))))
1320
1321;;; Fast function to multiply floating-point numbers.
1322(defun math-mul-float (a b)   ; [F F F]
1323  (math-make-float (math-mul (nth 1 a) (nth 1 b))
1324		   (+ (nth 2 a) (nth 2 b))))
1325
1326(defun math-sqr-float (a)   ; [F F]
1327  (math-make-float (math-mul (nth 1 a) (nth 1 a))
1328		   (+ (nth 2 a) (nth 2 a))))
1329
1330(defun math-intv-constp (a &optional finite)
1331  (and (or (Math-anglep (nth 2 a))
1332	   (and (equal (nth 2 a) '(neg (var inf var-inf)))
1333		(or (not finite)
1334		    (memq (nth 1 a) '(0 1)))))
1335       (or (Math-anglep (nth 3 a))
1336	   (and (equal (nth 3 a) '(var inf var-inf))
1337		(or (not finite)
1338		    (memq (nth 1 a) '(0 2)))))))
1339
1340(defun math-mul-zero (a b)
1341  (if (math-known-matrixp b)
1342      (if (math-vectorp b)
1343	  (math-map-vec-2 'math-mul a b)
1344	(math-mimic-ident 0 b))
1345    (if (math-infinitep b)
1346	'(var nan var-nan)
1347      (let ((aa nil) (bb nil))
1348	(if (and (eq (car-safe b) 'intv)
1349		 (progn
1350		   (and (equal (nth 2 b) '(neg (var inf var-inf)))
1351			(memq (nth 1 b) '(2 3))
1352			(setq aa (nth 2 b)))
1353		   (and (equal (nth 3 b) '(var inf var-inf))
1354			(memq (nth 1 b) '(1 3))
1355			(setq bb (nth 3 b)))
1356		   (or aa bb)))
1357	    (if (or (math-posp a)
1358		    (and (math-zerop a)
1359			 (or (memq calc-infinite-mode '(-1 1))
1360			     (setq aa '(neg (var inf var-inf))
1361				   bb '(var inf var-inf)))))
1362		(list 'intv 3 (or aa 0) (or bb 0))
1363	      (if (math-negp a)
1364		  (math-neg (list 'intv 3 (or aa 0) (or bb 0)))
1365		'(var nan var-nan)))
1366	  (if (or (math-floatp a) (math-floatp b)) '(float 0 0) 0))))))
1367
1368
1369(defun math-mul-symb-fancy (a b)
1370  (or (and math-simplify-only
1371	   (not (equal a math-simplify-only))
1372	   (list '* a b))
1373      (and (Math-equal-int a 1)
1374	   b)
1375      (and (Math-equal-int a -1)
1376	   (math-neg b))
1377      (and (or (and (Math-vectorp a) (math-known-scalarp b))
1378	       (and (Math-vectorp b) (math-known-scalarp a)))
1379	   (math-map-vec-2 'math-mul a b))
1380      (and (Math-objectp b) (not (Math-objectp a))
1381	   (math-mul b a))
1382      (and (eq (car-safe a) 'neg)
1383	   (math-neg (math-mul (nth 1 a) b)))
1384      (and (eq (car-safe b) 'neg)
1385	   (math-neg (math-mul a (nth 1 b))))
1386      (and (eq (car-safe a) '*)
1387	   (math-mul (nth 1 a)
1388		     (math-mul (nth 2 a) b)))
1389      (and (eq (car-safe a) '^)
1390	   (Math-looks-negp (nth 2 a))
1391	   (not (and (eq (car-safe b) '^) (Math-looks-negp (nth 2 b))))
1392	   (math-known-scalarp b t)
1393	   (math-div b (math-normalize
1394			(list '^ (nth 1 a) (math-neg (nth 2 a))))))
1395      (and (eq (car-safe b) '^)
1396	   (Math-looks-negp (nth 2 b))
1397	   (not (and (eq (car-safe a) '^) (Math-looks-negp (nth 2 a))))
1398           (not (math-known-matrixp (nth 1 b)))
1399	   (math-div a (math-normalize
1400			(list '^ (nth 1 b) (math-neg (nth 2 b))))))
1401      (and (eq (car-safe a) '/)
1402	   (or (math-known-scalarp a t) (math-known-scalarp b t))
1403	   (let ((temp (math-combine-prod (nth 2 a) b t nil t)))
1404	     (if temp
1405		 (math-mul (nth 1 a) temp)
1406	       (math-div (math-mul (nth 1 a) b) (nth 2 a)))))
1407      (and (eq (car-safe b) '/)
1408	   (math-div (math-mul a (nth 1 b)) (nth 2 b)))
1409      (and (eq (car-safe b) '+)
1410	   (Math-numberp a)
1411	   (or (Math-numberp (nth 1 b))
1412	       (Math-numberp (nth 2 b)))
1413	   (math-add (math-mul a (nth 1 b))
1414		     (math-mul a (nth 2 b))))
1415      (and (eq (car-safe b) '-)
1416	   (Math-numberp a)
1417	   (or (Math-numberp (nth 1 b))
1418	       (Math-numberp (nth 2 b)))
1419	   (math-sub (math-mul a (nth 1 b))
1420		     (math-mul a (nth 2 b))))
1421      (and (eq (car-safe b) '*)
1422	   (Math-numberp (nth 1 b))
1423	   (not (Math-numberp a))
1424	   (math-mul (nth 1 b) (math-mul a (nth 2 b))))
1425      (and (eq (car-safe a) 'calcFunc-idn)
1426	   (= (length a) 2)
1427	   (or (and (eq (car-safe b) 'calcFunc-idn)
1428		    (= (length b) 2)
1429		    (list 'calcFunc-idn (math-mul (nth 1 a) (nth 1 b))))
1430	       (and (math-known-scalarp b)
1431		    (list 'calcFunc-idn (math-mul (nth 1 a) b)))
1432	       (and (math-known-matrixp b)
1433		    (math-mul (nth 1 a) b))))
1434      (and (eq (car-safe b) 'calcFunc-idn)
1435	   (= (length b) 2)
1436	   (or (and (math-known-scalarp a)
1437		    (list 'calcFunc-idn (math-mul a (nth 1 b))))
1438	       (and (math-known-matrixp a)
1439		    (math-mul a (nth 1 b)))))
1440      (and (math-identity-matrix-p a t)
1441           (or (and (eq (car-safe b) 'calcFunc-idn)
1442                    (= (length b) 2)
1443                    (list 'calcFunc-idn (math-mul
1444                                         (nth 1 (nth 1 a))
1445                                         (nth 1 b))
1446                          (1- (length a))))
1447               (and (math-known-scalarp b)
1448                    (list 'calcFunc-idn (math-mul
1449                                         (nth 1 (nth 1 a)) b)
1450                          (1- (length a))))
1451               (and (math-known-matrixp b)
1452                    (math-mul (nth 1 (nth 1 a)) b))))
1453      (and (math-identity-matrix-p b t)
1454           (or (and (eq (car-safe a) 'calcFunc-idn)
1455                    (= (length a) 2)
1456                    (list 'calcFunc-idn (math-mul (nth 1 a)
1457                                                  (nth 1 (nth 1 b)))
1458                          (1- (length b))))
1459               (and (math-known-scalarp a)
1460                    (list 'calcFunc-idn (math-mul a (nth 1 (nth 1 b)))
1461                          (1- (length b))))
1462               (and (math-known-matrixp a)
1463                    (math-mul a (nth 1 (nth 1 b))))))
1464      (and (math-looks-negp b)
1465	   (math-mul (math-neg a) (math-neg b)))
1466      (and (eq (car-safe b) '-)
1467	   (math-looks-negp a)
1468	   (math-mul (math-neg a) (math-neg b)))
1469      (cond
1470       ((eq (car-safe b) '*)
1471	(let ((temp (math-combine-prod a (nth 1 b) nil nil t)))
1472	  (and temp
1473	       (math-mul temp (nth 2 b)))))
1474       (t
1475	(math-combine-prod a b nil nil nil)))
1476      (and (equal a '(var nan var-nan))
1477	   a)
1478      (and (equal b '(var nan var-nan))
1479	   b)
1480      (and (equal a '(var uinf var-uinf))
1481	   a)
1482      (and (equal b '(var uinf var-uinf))
1483	   b)
1484      (and (equal b '(var inf var-inf))
1485	   (let ((s1 (math-possible-signs a)))
1486	     (cond ((eq s1 4)
1487		    b)
1488		   ((eq s1 6)
1489		    '(intv 3 0 (var inf var-inf)))
1490		   ((eq s1 1)
1491		    (math-neg b))
1492		   ((eq s1 3)
1493		    '(intv 3 (neg (var inf var-inf)) 0))
1494		   ((and (eq (car a) 'intv) (math-intv-constp a))
1495		    '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))
1496		   ((and (eq (car a) 'cplx)
1497			 (math-zerop (nth 1 a)))
1498		    (list '* (list 'cplx 0 (calcFunc-sign (nth 2 a))) b))
1499		   ((eq (car a) 'polar)
1500		    (list '* (list 'polar 1 (nth 2 a)) b)))))
1501      (and (equal a '(var inf var-inf))
1502	   (math-mul b a))
1503      (list '* a b)))
1504
1505
1506(defun calcFunc-div (a &rest rest)
1507  (while rest
1508    (setq a (list '/ a (car rest))
1509	  rest (cdr rest)))
1510  (math-normalize a))
1511
1512(defun math-div-objects-fancy (a b)
1513  (cond ((and (Math-numberp a) (Math-numberp b))
1514	 (math-normalize
1515	  (cond ((math-want-polar a b)
1516		 (let ((a (math-polar a))
1517		       (b (math-polar b)))
1518		   (list 'polar
1519			 (math-div (nth 1 a) (nth 1 b))
1520			 (math-fix-circular (math-sub (nth 2 a)
1521						      (nth 2 b))))))
1522		((Math-realp b)
1523		 (setq a (math-complex a))
1524		 (list 'cplx (math-div (nth 1 a) b)
1525		       (math-div (nth 2 a) b)))
1526		(t
1527		 (setq a (math-complex a)
1528		       b (math-complex b))
1529		 (math-div
1530		  (list 'cplx
1531			(math-add (math-mul (nth 1 a) (nth 1 b))
1532				  (math-mul (nth 2 a) (nth 2 b)))
1533			(math-sub (math-mul (nth 2 a) (nth 1 b))
1534				  (math-mul (nth 1 a) (nth 2 b))))
1535		  (math-add (math-sqr (nth 1 b))
1536			    (math-sqr (nth 2 b))))))))
1537	((math-matrixp b)
1538	 (if (math-square-matrixp b)
1539	     (let ((n1 (length b)))
1540	       (if (Math-vectorp a)
1541		   (if (math-matrixp a)
1542		       (if (= (length a) n1)
1543			   (math-lud-solve (math-matrix-lud b) a b)
1544			 (if (= (length (nth 1 a)) n1)
1545			     (math-transpose
1546			      (math-lud-solve (math-matrix-lud
1547					       (math-transpose b))
1548					      (math-transpose a) b))
1549			   (math-dimension-error)))
1550		     (if (= (length a) n1)
1551			 (math-mat-col (math-lud-solve (math-matrix-lud b)
1552						       (math-col-matrix a) b)
1553				       1)
1554		       (math-dimension-error)))
1555		 (if (Math-equal-int a 1)
1556		     (calcFunc-inv b)
1557		   (math-mul a (calcFunc-inv b)))))
1558	   (math-reject-arg b 'square-matrixp)))
1559	((and (Math-vectorp a) (Math-objectp b))
1560	 (math-map-vec-2 'math-div a b))
1561	((eq (car-safe a) 'sdev)
1562	 (if (eq (car-safe b) 'sdev)
1563	     (let ((x (math-div (nth 1 a) (nth 1 b))))
1564	       (math-make-sdev x
1565			       (math-div (math-hypot (nth 2 a)
1566						     (math-mul (nth 2 b) x))
1567					 (nth 1 b))))
1568	   (if (or (Math-scalarp b)
1569		   (not (Math-objvecp b)))
1570	       (math-make-sdev (math-div (nth 1 a) b) (math-div (nth 2 a) b))
1571	     (math-reject-arg 'realp b))))
1572	((and (eq (car-safe b) 'sdev)
1573	      (or (Math-scalarp a)
1574		  (not (Math-objvecp a))))
1575	 (let ((x (math-div a (nth 1 b))))
1576	   (math-make-sdev x
1577			   (math-div (math-mul (nth 2 b) x) (nth 1 b)))))
1578	((and (eq (car-safe a) 'intv) (Math-anglep b))
1579	 (if (Math-negp b)
1580	     (math-neg (math-div a (math-neg b)))
1581	   (math-make-intv (nth 1 a)
1582			   (math-div (nth 2 a) b)
1583			   (math-div (nth 3 a) b))))
1584	((and (eq (car-safe b) 'intv) (Math-anglep a))
1585	 (if (or (Math-posp (nth 2 b))
1586		 (and (Math-zerop (nth 2 b)) (or (memq (nth 1 b) '(0 1))
1587						 calc-infinite-mode)))
1588	     (if (Math-negp a)
1589		 (math-neg (math-div (math-neg a) b))
1590	       (let ((calc-infinite-mode 1))
1591		 (math-make-intv (aref [0 2 1 3] (nth 1 b))
1592				 (math-div a (nth 3 b))
1593				 (math-div a (nth 2 b)))))
1594	   (if (or (Math-negp (nth 3 b))
1595		   (and (Math-zerop (nth 3 b)) (or (memq (nth 1 b) '(0 2))
1596						   calc-infinite-mode)))
1597	       (math-neg (math-div a (math-neg b)))
1598	     (if calc-infinite-mode
1599		 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1600	       (math-reject-arg b "*Division by zero")))))
1601	((and (eq (car-safe a) 'intv) (math-intv-constp a)
1602	      (eq (car-safe b) 'intv) (math-intv-constp b))
1603	 (if (or (Math-posp (nth 2 b))
1604		 (and (Math-zerop (nth 2 b)) (or (memq (nth 1 b) '(0 1))
1605						 calc-infinite-mode)))
1606	     (let* ((calc-infinite-mode 1)
1607		    (lo (math-div a (nth 2 b)))
1608		    (hi (math-div a (nth 3 b))))
1609	       (or (eq (car-safe lo) 'intv)
1610		   (setq lo (list 'intv (if (memq (nth 1 b) '(2 3)) 3 0)
1611				  lo lo)))
1612	       (or (eq (car-safe hi) 'intv)
1613		   (setq hi (list 'intv (if (memq (nth 1 b) '(1 3)) 3 0)
1614				  hi hi)))
1615	       (math-combine-intervals
1616		(nth 2 lo) (and (or (memq (nth 1 b) '(2 3))
1617				    (and (math-infinitep (nth 2 lo))
1618					 (not (math-zerop (nth 2 b)))))
1619				(memq (nth 1 lo) '(2 3)))
1620		(nth 3 lo) (and (or (memq (nth 1 b) '(2 3))
1621				    (and (math-infinitep (nth 3 lo))
1622					 (not (math-zerop (nth 2 b)))))
1623				(memq (nth 1 lo) '(1 3)))
1624		(nth 2 hi) (and (or (memq (nth 1 b) '(1 3))
1625				    (and (math-infinitep (nth 2 hi))
1626					 (not (math-zerop (nth 3 b)))))
1627				(memq (nth 1 hi) '(2 3)))
1628		(nth 3 hi) (and (or (memq (nth 1 b) '(1 3))
1629				    (and (math-infinitep (nth 3 hi))
1630					 (not (math-zerop (nth 3 b)))))
1631				(memq (nth 1 hi) '(1 3)))))
1632	   (if (or (Math-negp (nth 3 b))
1633		   (and (Math-zerop (nth 3 b)) (or (memq (nth 1 b) '(0 2))
1634						   calc-infinite-mode)))
1635	       (math-neg (math-div a (math-neg b)))
1636	     (if calc-infinite-mode
1637		 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1638	       (math-reject-arg b "*Division by zero")))))
1639	((and (eq (car-safe a) 'mod)
1640	      (eq (car-safe b) 'mod)
1641	      (equal (nth 2 a) (nth 2 b)))
1642	 (math-make-mod (math-div-mod (nth 1 a) (nth 1 b) (nth 2 a))
1643			(nth 2 a)))
1644	((and (eq (car-safe a) 'mod)
1645	      (Math-anglep b))
1646	 (math-make-mod (math-div-mod (nth 1 a) b (nth 2 a)) (nth 2 a)))
1647	((and (eq (car-safe b) 'mod)
1648	      (Math-anglep a))
1649	 (math-make-mod (math-div-mod a (nth 1 b) (nth 2 b)) (nth 2 b)))
1650	((eq (car-safe a) 'hms)
1651	 (if (eq (car-safe b) 'hms)
1652	     (math-with-extra-prec 1
1653	       (math-div (math-from-hms a 'deg)
1654			 (math-from-hms b 'deg)))
1655	   (math-with-extra-prec 2
1656	     (math-to-hms (math-div (math-from-hms a 'deg) b) 'deg))))
1657	(t (calc-record-why "*Incompatible arguments for /" a b))))
1658
1659(defun math-div-by-zero (a b)
1660  (if (math-infinitep a)
1661      (if (or (equal a '(var nan var-nan))
1662	      (equal b '(var uinf var-uinf))
1663	      (memq calc-infinite-mode '(-1 1)))
1664	  a
1665	'(var uinf var-uinf))
1666    (if calc-infinite-mode
1667	(if (math-zerop a)
1668	    '(var nan var-nan)
1669	  (if (eq calc-infinite-mode 1)
1670	      (math-mul a '(var inf var-inf))
1671	    (if (eq calc-infinite-mode -1)
1672		(math-mul a '(neg (var inf var-inf)))
1673	      (if (eq (car-safe a) 'intv)
1674		  '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1675		'(var uinf var-uinf)))))
1676      (math-reject-arg a "*Division by zero"))))
1677
1678(defun math-div-zero (a b)
1679  (if (math-known-matrixp b)
1680      (if (math-vectorp b)
1681	  (math-map-vec-2 'math-div a b)
1682	(math-mimic-ident 0 b))
1683    (if (equal b '(var nan var-nan))
1684	b
1685      (if (and (eq (car-safe b) 'intv) (math-intv-constp b)
1686	       (not (math-posp b)) (not (math-negp b)))
1687	  (if calc-infinite-mode
1688	      (list 'intv 3
1689		    (if (and (math-zerop (nth 2 b))
1690			     (memq calc-infinite-mode '(1 -1)))
1691			(nth 2 b) '(neg (var inf var-inf)))
1692		    (if (and (math-zerop (nth 3 b))
1693			     (memq calc-infinite-mode '(1 -1)))
1694			(nth 3 b) '(var inf var-inf)))
1695	    (math-reject-arg b "*Division by zero"))
1696	a))))
1697
1698;; For math-div-symb-fancy
1699(defvar math-trig-inverses
1700  '((calcFunc-sin . calcFunc-csc)
1701    (calcFunc-cos . calcFunc-sec)
1702    (calcFunc-tan . calcFunc-cot)
1703    (calcFunc-sec . calcFunc-cos)
1704    (calcFunc-csc . calcFunc-sin)
1705    (calcFunc-cot . calcFunc-tan)
1706    (calcFunc-sinh . calcFunc-csch)
1707    (calcFunc-cosh . calcFunc-sech)
1708    (calcFunc-tanh . calcFunc-coth)
1709    (calcFunc-sech . calcFunc-cosh)
1710    (calcFunc-csch . calcFunc-sinh)
1711    (calcFunc-coth . calcFunc-tanh)))
1712
1713(defvar math-div-trig)
1714(defvar math-div-non-trig)
1715
1716(defun math-div-new-trig (tr)
1717  (if math-div-trig
1718      (setq math-div-trig
1719            (list '* tr math-div-trig))
1720    (setq math-div-trig tr)))
1721
1722(defun math-div-new-non-trig (ntr)
1723  (if math-div-non-trig
1724      (setq math-div-non-trig
1725            (list '* ntr math-div-non-trig))
1726    (setq math-div-non-trig ntr)))
1727
1728(defun math-div-isolate-trig (expr)
1729  (if (eq (car-safe expr) '*)
1730      (progn
1731        (math-div-isolate-trig-term (nth 1 expr))
1732        (math-div-isolate-trig (nth 2 expr)))
1733    (math-div-isolate-trig-term expr)))
1734
1735(defun math-div-isolate-trig-term (term)
1736  (let ((fn (assoc (car-safe term) math-trig-inverses)))
1737    (if fn
1738        (math-div-new-trig
1739         (cons (cdr fn) (cdr term)))
1740      (math-div-new-non-trig term))))
1741
1742(defun math-div-symb-fancy (a b)
1743  (or (and (math-known-matrixp b)
1744           (math-mul a (math-pow b -1)))
1745      (and math-simplify-only
1746	   (not (equal a math-simplify-only))
1747	   (list '/ a b))
1748      (and (Math-equal-int b 1) a)
1749      (and (Math-equal-int b -1) (math-neg a))
1750      (and (Math-vectorp a) (math-known-scalarp b)
1751	   (math-map-vec-2 'math-div a b))
1752      (and (eq (car-safe b) '^)
1753	   (or (Math-looks-negp (nth 2 b)) (Math-equal-int a 1))
1754	   (math-mul a (math-normalize
1755			(list '^ (nth 1 b) (math-neg (nth 2 b))))))
1756      (and (eq (car-safe a) 'neg)
1757	   (math-neg (math-div (nth 1 a) b)))
1758      (and (eq (car-safe b) 'neg)
1759	   (math-neg (math-div a (nth 1 b))))
1760      (and (eq (car-safe a) '/)
1761	   (math-div (nth 1 a) (math-mul (nth 2 a) b)))
1762      (and (eq (car-safe b) '/)
1763	   (or (math-known-scalarp (nth 1 b) t)
1764	       (math-known-scalarp (nth 2 b) t))
1765	   (math-div (math-mul a (nth 2 b)) (nth 1 b)))
1766      (and (eq (car-safe b) 'frac)
1767	   (math-mul (math-make-frac (nth 2 b) (nth 1 b)) a))
1768      (and (eq (car-safe a) '+)
1769	   (or (Math-numberp (nth 1 a))
1770	       (Math-numberp (nth 2 a)))
1771	   (Math-numberp b)
1772	   (math-add (math-div (nth 1 a) b)
1773		     (math-div (nth 2 a) b)))
1774      (and (eq (car-safe a) '-)
1775	   (or (Math-numberp (nth 1 a))
1776	       (Math-numberp (nth 2 a)))
1777	   (Math-numberp b)
1778	   (math-sub (math-div (nth 1 a) b)
1779		     (math-div (nth 2 a) b)))
1780      (and (or (eq (car-safe a) '-)
1781	       (math-looks-negp a))
1782	   (math-looks-negp b)
1783	   (math-div (math-neg a) (math-neg b)))
1784      (and (eq (car-safe b) '-)
1785	   (math-looks-negp a)
1786	   (math-div (math-neg a) (math-neg b)))
1787      (and (eq (car-safe a) 'calcFunc-idn)
1788	   (= (length a) 2)
1789	   (or (and (eq (car-safe b) 'calcFunc-idn)
1790		    (= (length b) 2)
1791		    (list 'calcFunc-idn (math-div (nth 1 a) (nth 1 b))))
1792	       (and (math-known-scalarp b)
1793		    (list 'calcFunc-idn (math-div (nth 1 a) b)))
1794	       (and (math-known-matrixp b)
1795		    (math-div (nth 1 a) b))))
1796      (and (eq (car-safe b) 'calcFunc-idn)
1797	   (= (length b) 2)
1798	   (or (and (math-known-scalarp a)
1799		    (list 'calcFunc-idn (math-div a (nth 1 b))))
1800	       (and (math-known-matrixp a)
1801		    (math-div a (nth 1 b)))))
1802      (and math-simplifying
1803           (let ((math-div-trig nil)
1804                 (math-div-non-trig nil))
1805             (math-div-isolate-trig b)
1806             (if math-div-trig
1807                 (if math-div-non-trig
1808                     (math-div (math-mul a math-div-trig) math-div-non-trig)
1809                   (math-mul a math-div-trig))
1810               nil)))
1811      (if (and calc-matrix-mode
1812	       (or (math-known-matrixp a) (math-known-matrixp b)))
1813	  (math-combine-prod a b nil t nil)
1814	(if (eq (car-safe a) '*)
1815	    (if (eq (car-safe b) '*)
1816		(let ((c (math-combine-prod (nth 1 a) (nth 1 b) nil t t)))
1817		  (and c
1818		       (math-div (math-mul c (nth 2 a)) (nth 2 b))))
1819	      (let ((c (math-combine-prod (nth 1 a) b nil t t)))
1820		(and c
1821		     (math-mul c (nth 2 a)))))
1822	  (if (eq (car-safe b) '*)
1823	      (let ((c (math-combine-prod a (nth 1 b) nil t t)))
1824		(and c
1825		     (math-div c (nth 2 b))))
1826	    (math-combine-prod a b nil t nil))))
1827      (and (math-infinitep a)
1828	   (if (math-infinitep b)
1829	       '(var nan var-nan)
1830	     (if (or (equal a '(var nan var-nan))
1831		     (equal a '(var uinf var-uinf)))
1832		 a
1833	       (if (equal a '(var inf var-inf))
1834		   (if (or (math-posp b)
1835			   (and (eq (car-safe b) 'intv)
1836				(math-zerop (nth 2 b))))
1837		       (if (and (eq (car-safe b) 'intv)
1838				(not (math-intv-constp b t)))
1839			   '(intv 3 0 (var inf var-inf))
1840			 a)
1841		     (if (or (math-negp b)
1842			     (and (eq (car-safe b) 'intv)
1843			      (math-zerop (nth 3 b))))
1844			 (if (and (eq (car-safe b) 'intv)
1845				  (not (math-intv-constp b t)))
1846			     '(intv 3 (neg (var inf var-inf)) 0)
1847			   (math-neg a))
1848		       (if (and (eq (car-safe b) 'intv)
1849				(math-negp (nth 2 b)) (math-posp (nth 3 b)))
1850			   '(intv 3 (neg (var inf var-inf))
1851				  (var inf var-inf)))))))))
1852      (and (math-infinitep b)
1853	   (if (equal b '(var nan var-nan))
1854	       b
1855	     (let ((calc-infinite-mode 1))
1856	       (math-mul-zero b a))))
1857      (list '/ a b)))
1858
1859;;; Division from the left.
1860(defun calcFunc-ldiv (a b)
1861  (if (math-known-scalarp a)
1862      (math-div b a)
1863    (math-mul (math-pow a -1) b)))
1864
1865(defun calcFunc-mod (a b)
1866  (math-normalize (list '% a b)))
1867
1868(defun math-mod-fancy (a b)
1869  (cond ((equal b '(var inf var-inf))
1870	 (if (or (math-posp a) (math-zerop a))
1871	     a
1872	   (if (math-negp a)
1873	       b
1874	     (if (eq (car-safe a) 'intv)
1875		 (if (math-negp (nth 2 a))
1876		     '(intv 3 0 (var inf var-inf))
1877		   a)
1878	       (list '% a b)))))
1879	((and (eq (car-safe a) 'mod) (Math-realp b) (math-posp b))
1880	 (math-make-mod (nth 1 a) b))
1881	((and (eq (car-safe a) 'intv) (math-intv-constp a t) (math-posp b))
1882	 (math-mod-intv a b))
1883	(t
1884	 (if (Math-anglep a)
1885	     (calc-record-why 'anglep b)
1886	   (calc-record-why 'anglep a))
1887	 (list '% a b))))
1888
1889
1890(defun calcFunc-pow (a b)
1891  (math-normalize (list '^ a b)))
1892
1893(defun math-pow-of-zero (a b)
1894  "Raise A to the power of B, where A is a form of zero."
1895  (if (math-floatp b) (setq a (math-float a)))
1896  (cond
1897   ;; 0^0 = 1
1898   ((eq b 0)
1899    1)
1900   ;; 0^0.0, etc., are undetermined
1901   ((Math-zerop b)
1902    (if calc-infinite-mode
1903        '(var nan var-nan)
1904      (math-reject-arg (list '^ a b) "*Indeterminate form")))
1905   ;; 0^positive = 0
1906   ((math-known-posp b)
1907    a)
1908   ;; 0^negative is undefined (let math-div handle it)
1909   ((math-known-negp b)
1910    (math-div 1 a))
1911   ;; 0^infinity is undefined
1912   ((math-infinitep b)
1913    '(var nan var-nan))
1914   ;; Some intervals
1915   ((and (eq (car b) 'intv)
1916         calc-infinite-mode
1917         (math-negp (nth 2 b))
1918         (math-posp (nth 3 b)))
1919    '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))
1920   ;; If none of the above, leave it alone.
1921   (t
1922    (list '^ a b))))
1923
1924(defun math-pow-zero (a b)
1925  (if (eq (car-safe a) 'mod)
1926      (math-make-mod 1 (nth 2 a))
1927    (if (math-known-matrixp a)
1928	(math-mimic-ident 1 a)
1929      (if (math-infinitep a)
1930	  '(var nan var-nan)
1931	(if (and (eq (car a) 'intv) (math-intv-constp a)
1932		 (or (and (not (math-posp a)) (not (math-negp a)))
1933		     (not (math-intv-constp a t))))
1934	    '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1935	  (if (or (math-floatp a) (math-floatp b))
1936	      '(float 1 0) 1))))))
1937
1938(defun math-pow-fancy (a b)
1939  (cond ((and (Math-numberp a) (Math-numberp b))
1940	 (or (if (memq (math-quarter-integer b) '(1 2 3))
1941		 (let ((sqrt (math-sqrt (if (math-floatp b)
1942					    (math-float a) a))))
1943		   (and (Math-numberp sqrt)
1944			(math-pow sqrt (math-mul 2 b))))
1945	       (and (eq (car b) 'frac)
1946		    (integerp (nth 2 b))
1947		    (<= (nth 2 b) 10)
1948		    (let ((root (math-nth-root a (nth 2 b))))
1949		      (and root (math-ipow root (nth 1 b))))))
1950	     (and (or (eq a 10) (equal a '(float 1 1)))
1951		  (math-num-integerp b)
1952		  (calcFunc-scf '(float 1 0) b))
1953	     (and calc-symbolic-mode
1954		  (list '^ a b))
1955	     (math-with-extra-prec 2
1956	       (math-exp-raw
1957		(math-float (math-mul b (math-ln-raw (math-float a))))))))
1958	((or (not (Math-objvecp a))
1959	     (not (Math-objectp b)))
1960	 (let (temp)
1961	   (cond ((and math-simplify-only
1962		       (not (equal a math-simplify-only)))
1963		  (list '^ a b))
1964                 ((and (eq (car-safe a) '*)
1965                       (or
1966                        (and
1967                         (math-known-matrixp (nth 1 a))
1968                         (math-known-matrixp (nth 2 a)))
1969                        (and
1970                         calc-matrix-mode
1971                         (not (eq calc-matrix-mode 'scalar))
1972                         (and (not (math-known-scalarp (nth 1 a)))
1973                              (not (math-known-scalarp (nth 2 a)))))))
1974                  (if (and (= b -1)
1975                           (math-known-square-matrixp (nth 1 a))
1976                           (math-known-square-matrixp (nth 2 a)))
1977                      (math-mul (math-pow-fancy (nth 2 a) -1)
1978                                (math-pow-fancy (nth 1 a) -1))
1979                    (list '^ a b)))
1980		 ((and (eq (car-safe a) '*)
1981		       (or (math-known-num-integerp b)
1982			   (math-known-nonnegp (nth 1 a))
1983			   (math-known-nonnegp (nth 2 a))))
1984		  (math-mul (math-pow (nth 1 a) b)
1985			    (math-pow (nth 2 a) b)))
1986		 ((and (eq (car-safe a) '/)
1987		       (or (math-known-num-integerp b)
1988			   (math-known-nonnegp (nth 2 a))))
1989		  (math-div (math-pow (nth 1 a) b)
1990			    (math-pow (nth 2 a) b)))
1991		 ((and (eq (car-safe a) '/)
1992		       (math-known-nonnegp (nth 1 a))
1993		       (not (math-equal-int (nth 1 a) 1)))
1994		  (math-mul (math-pow (nth 1 a) b)
1995			    (math-pow (math-div 1 (nth 2 a)) b)))
1996		 ((and (eq (car-safe a) '^)
1997		       (or (math-known-num-integerp b)
1998			   (math-known-nonnegp (nth 1 a))))
1999		  (math-pow (nth 1 a) (math-mul (nth 2 a) b)))
2000		 ((and (eq (car-safe a) 'calcFunc-sqrt)
2001		       (or (math-known-num-integerp b)
2002			   (math-known-nonnegp (nth 1 a))))
2003		  (math-pow (nth 1 a) (math-div b 2)))
2004		 ((and (eq (car-safe a) '^)
2005		       (math-known-evenp (nth 2 a))
2006		       (memq (math-quarter-integer b) '(1 2 3))
2007		       (math-known-realp (nth 1 a)))
2008		  (math-abs (math-pow (nth 1 a) (math-mul (nth 2 a) b))))
2009		 ((and (math-looks-negp a)
2010		       (math-known-integerp b)
2011		       (setq temp (or (and (math-known-evenp b)
2012					   (math-pow (math-neg a) b))
2013				      (and (math-known-oddp b)
2014					   (math-neg (math-pow (math-neg a)
2015							       b))))))
2016		  temp)
2017		 ((and (eq (car-safe a) 'calcFunc-abs)
2018		       (math-known-realp (nth 1 a))
2019		       (math-known-evenp b))
2020		  (math-pow (nth 1 a) b))
2021		 ((math-infinitep a)
2022		  (cond ((equal a '(var nan var-nan))
2023			 a)
2024			((eq (car a) 'neg)
2025			 (math-mul (math-pow -1 b) (math-pow (nth 1 a) b)))
2026			((math-posp b)
2027			 a)
2028			((math-negp b)
2029			 (if (math-floatp b) '(float 0 0) 0))
2030			((and (eq (car-safe b) 'intv)
2031			      (math-intv-constp b))
2032			 '(intv 3 0 (var inf var-inf)))
2033			(t
2034			 '(var nan var-nan))))
2035		 ((math-infinitep b)
2036		  (let (scale)
2037		    (cond ((math-negp b)
2038			   (math-pow (math-div 1 a) (math-neg b)))
2039			  ((not (math-posp b))
2040			   '(var nan var-nan))
2041			  ((math-equal-int (setq scale (calcFunc-abssqr a)) 1)
2042			   '(var nan var-nan))
2043			  ((Math-lessp scale 1)
2044			   (if (math-floatp a) '(float 0 0) 0))
2045			  ((Math-lessp 1 a)
2046			   b)
2047			  ((Math-lessp a -1)
2048			   '(var uinf var-uinf))
2049			  ((and (eq (car a) 'intv)
2050				(math-intv-constp a))
2051			   (if (Math-lessp -1 a)
2052			       (if (math-equal-int (nth 3 a) 1)
2053				   '(intv 3 0 1)
2054				 '(intv 3 0 (var inf var-inf)))
2055			     '(intv 3 (neg (var inf var-inf))
2056				    (var inf var-inf))))
2057			  (t (list '^ a b)))))
2058		 ((and (eq (car-safe a) 'calcFunc-idn)
2059		       (= (length a) 2)
2060		       (math-known-num-integerp b))
2061		  (list 'calcFunc-idn (math-pow (nth 1 a) b)))
2062		 (t (if (Math-objectp a)
2063			(calc-record-why 'objectp b)
2064		      (calc-record-why 'objectp a))
2065		    (list '^ a b)))))
2066	((and (eq (car-safe a) 'sdev) (eq (car-safe b) 'sdev))
2067	 (if (and (math-constp a) (math-constp b))
2068	     (math-with-extra-prec 2
2069	       (let* ((ln (math-ln-raw (math-float (nth 1 a))))
2070		      (pow (math-exp-raw
2071			    (math-float (math-mul (nth 1 b) ln)))))
2072		 (math-make-sdev
2073		  pow
2074		  (math-mul
2075		   pow
2076		   (math-hypot (math-mul (nth 2 a)
2077					 (math-div (nth 1 b) (nth 1 a)))
2078			       (math-mul (nth 2 b) ln))))))
2079	   (let ((pow (math-pow (nth 1 a) (nth 1 b))))
2080	     (math-make-sdev
2081	      pow
2082	      (math-mul pow
2083			(math-hypot (math-mul (nth 2 a)
2084					      (math-div (nth 1 b) (nth 1 a)))
2085				    (math-mul (nth 2 b) (calcFunc-ln
2086							 (nth 1 a)))))))))
2087	((and (eq (car-safe a) 'sdev) (Math-numberp b))
2088	 (if (math-constp a)
2089	     (math-with-extra-prec 2
2090	       (let ((pow (math-pow (nth 1 a) (math-sub b 1))))
2091		 (math-make-sdev (math-mul pow (nth 1 a))
2092				 (math-mul pow (math-mul (nth 2 a) b)))))
2093	   (math-make-sdev (math-pow (nth 1 a) b)
2094			   (math-mul (math-pow (nth 1 a) (math-add b -1))
2095				     (math-mul (nth 2 a) b)))))
2096	((and (eq (car-safe b) 'sdev) (Math-numberp a))
2097	 (math-with-extra-prec 2
2098	   (let* ((ln (math-ln-raw (math-float a)))
2099		  (pow (calcFunc-exp (math-mul (nth 1 b) ln))))
2100	     (math-make-sdev pow (math-mul pow (math-mul (nth 2 b) ln))))))
2101	((and (eq (car-safe a) 'intv) (math-intv-constp a)
2102	      (Math-realp b)
2103	      (or (Math-natnump b)
2104		  (Math-posp (nth 2 a))
2105		  (and (math-zerop (nth 2 a))
2106		       (or (Math-posp b)
2107			   (and (Math-integerp b) calc-infinite-mode)))
2108		  (Math-negp (nth 3 a))
2109		  (and (math-zerop (nth 3 a))
2110		       (or (Math-posp b)
2111			   (and (Math-integerp b) calc-infinite-mode)))))
2112	 (if (math-evenp b)
2113	     (setq a (math-abs a)))
2114	 (let ((calc-infinite-mode (if (math-zerop (nth 3 a)) -1 1)))
2115	   (math-sort-intv (nth 1 a)
2116			   (math-pow (nth 2 a) b)
2117			   (math-pow (nth 3 a) b))))
2118	((and (eq (car-safe b) 'intv) (math-intv-constp b)
2119	      (Math-realp a) (Math-posp a))
2120	 (math-sort-intv (nth 1 b)
2121			 (math-pow a (nth 2 b))
2122			 (math-pow a (nth 3 b))))
2123	((and (eq (car-safe a) 'intv) (math-intv-constp a)
2124	      (eq (car-safe b) 'intv) (math-intv-constp b)
2125	      (or (and (not (Math-negp (nth 2 a)))
2126		       (not (Math-negp (nth 2 b))))
2127		  (and (Math-posp (nth 2 a))
2128		       (not (Math-posp (nth 3 b))))))
2129	 (let ((lo (math-pow a (nth 2 b)))
2130	       (hi (math-pow a (nth 3 b))))
2131	   (or (eq (car-safe lo) 'intv)
2132	       (setq lo (list 'intv (if (memq (nth 1 b) '(2 3)) 3 0) lo lo)))
2133	   (or (eq (car-safe hi) 'intv)
2134	       (setq hi (list 'intv (if (memq (nth 1 b) '(1 3)) 3 0) hi hi)))
2135	   (math-combine-intervals
2136	    (nth 2 lo) (and (or (memq (nth 1 b) '(2 3))
2137				(math-infinitep (nth 2 lo)))
2138			    (memq (nth 1 lo) '(2 3)))
2139	    (nth 3 lo) (and (or (memq (nth 1 b) '(2 3))
2140				(math-infinitep (nth 3 lo)))
2141			    (memq (nth 1 lo) '(1 3)))
2142	    (nth 2 hi) (and (or (memq (nth 1 b) '(1 3))
2143				(math-infinitep (nth 2 hi)))
2144			    (memq (nth 1 hi) '(2 3)))
2145	    (nth 3 hi) (and (or (memq (nth 1 b) '(1 3))
2146				(math-infinitep (nth 3 hi)))
2147			    (memq (nth 1 hi) '(1 3))))))
2148	((and (eq (car-safe a) 'mod) (eq (car-safe b) 'mod)
2149	      (equal (nth 2 a) (nth 2 b)))
2150	 (math-make-mod (math-pow-mod (nth 1 a) (nth 1 b) (nth 2 a))
2151			(nth 2 a)))
2152	((and (eq (car-safe a) 'mod) (Math-anglep b))
2153	 (math-make-mod (math-pow-mod (nth 1 a) b (nth 2 a)) (nth 2 a)))
2154	((and (eq (car-safe b) 'mod) (Math-anglep a))
2155	 (math-make-mod (math-pow-mod a (nth 1 b) (nth 2 b)) (nth 2 b)))
2156	((not (Math-numberp a))
2157	 (math-reject-arg a 'numberp))
2158	(t
2159	 (math-reject-arg b 'numberp))))
2160
2161(defun math-quarter-integer (x)
2162  (if (Math-integerp x)
2163      0
2164    (if (math-negp x)
2165	(progn
2166	  (setq x (math-quarter-integer (math-neg x)))
2167	  (and x (- 4 x)))
2168      (if (eq (car x) 'frac)
2169	  (if (eq (nth 2 x) 2)
2170	      2
2171	    (and (eq (nth 2 x) 4)
2172		 (progn
2173		   (setq x (nth 1 x))
2174		   (% (if (consp x) (nth 1 x) x) 4))))
2175	(if (eq (car x) 'float)
2176	    (if (>= (nth 2 x) 0)
2177		0
2178	      (if (= (nth 2 x) -1)
2179		  (progn
2180		    (setq x (nth 1 x))
2181		    (and (= (% (if (consp x) (nth 1 x) x) 10) 5) 2))
2182		(if (= (nth 2 x) -2)
2183		    (progn
2184		      (setq x (nth 1 x)
2185			    x (% (if (consp x) (nth 1 x) x) 100))
2186		      (if (= x 25) 1
2187			(if (= x 75) 3)))))))))))
2188
2189;;; This assumes A < M and M > 0.
2190(defun math-pow-mod (a b m)   ; [R R R R]
2191  (if (and (Math-integerp a) (Math-integerp b) (Math-integerp m))
2192      (if (Math-negp b)
2193	  (math-div-mod 1 (math-pow-mod a (Math-integer-neg b) m) m)
2194	(if (eq m 1)
2195	    0
2196	  (math-pow-mod-step a b m)))
2197    (math-mod (math-pow a b) m)))
2198
2199(defun math-pow-mod-step (a n m)   ; [I I I I]
2200  (math-working "pow" a)
2201  (let ((val (cond
2202	      ((eq n 0) 1)
2203	      ((eq n 1) a)
2204	      (t
2205	       (let ((rest (math-pow-mod-step
2206			    (math-imod (math-mul a a) m)
2207			    (math-div2 n)
2208			    m)))
2209		 (if (math-evenp n)
2210		     rest
2211		   (math-mod (math-mul a rest) m)))))))
2212    (math-working "pow" val)
2213    val))
2214
2215
2216;;; Compute the minimum of two real numbers.  [R R R] [Public]
2217(defun math-min (a b)
2218  (if (and (consp a) (eq (car a) 'intv))
2219      (if (and (consp b) (eq (car b) 'intv))
2220	  (let ((lo (nth 2 a))
2221		(lom (memq (nth 1 a) '(2 3)))
2222		(hi (nth 3 a))
2223		(him (memq (nth 1 a) '(1 3)))
2224		res)
2225	    (if (= (setq res (math-compare (nth 2 b) lo)) -1)
2226		(setq lo (nth 2 b) lom (memq (nth 1 b) '(2 3)))
2227	      (if (= res 0)
2228		  (setq lom (or lom (memq (nth 1 b) '(2 3))))))
2229	    (if (= (setq res (math-compare (nth 3 b) hi)) -1)
2230		(setq hi (nth 3 b) him (memq (nth 1 b) '(1 3)))
2231	      (if (= res 0)
2232		  (setq him (or him (memq (nth 1 b) '(1 3))))))
2233	    (math-make-intv (+ (if lom 2 0) (if him 1 0)) lo hi))
2234	(math-min a (list 'intv 3 b b)))
2235    (if (and (consp b) (eq (car b) 'intv))
2236	(math-min (list 'intv 3 a a) b)
2237      (let ((res (math-compare a b)))
2238	(if (= res 1)
2239	    b
2240	  (if (= res 2)
2241	      '(var nan var-nan)
2242	    a))))))
2243
2244(defun calcFunc-min (&optional a &rest b)
2245  (if (not a)
2246      '(var inf var-inf)
2247    (if (not (or (Math-anglep a) (eq (car a) 'date)
2248		 (and (eq (car a) 'intv) (math-intv-constp a))
2249		 (math-infinitep a)))
2250	(math-reject-arg a 'anglep))
2251    (math-min-list a b)))
2252
2253(defun math-min-list (a b)
2254  (if b
2255      (if (or (Math-anglep (car b)) (eq (car b) 'date)
2256	      (and (eq (car (car b)) 'intv) (math-intv-constp (car b)))
2257	      (math-infinitep (car b)))
2258	  (math-min-list (math-min a (car b)) (cdr b))
2259	(math-reject-arg (car b) 'anglep))
2260    a))
2261
2262;;; Compute the maximum of two real numbers.  [R R R] [Public]
2263(defun math-max (a b)
2264  (if (or (and (consp a) (eq (car a) 'intv))
2265	  (and (consp b) (eq (car b) 'intv)))
2266      (math-neg (math-min (math-neg a) (math-neg b)))
2267    (let ((res (math-compare a b)))
2268      (if (= res -1)
2269	  b
2270	(if (= res 2)
2271	      '(var nan var-nan)
2272	  a)))))
2273
2274(defun calcFunc-max (&optional a &rest b)
2275  (if (not a)
2276      '(neg (var inf var-inf))
2277    (if (not (or (Math-anglep a) (eq (car a) 'date)
2278		 (and (eq (car a) 'intv) (math-intv-constp a))
2279		 (math-infinitep a)))
2280	(math-reject-arg a 'anglep))
2281    (math-max-list a b)))
2282
2283(defun math-max-list (a b)
2284  (if b
2285      (if (or (Math-anglep (car b)) (eq (car b) 'date)
2286	      (and (eq (car (car b)) 'intv) (math-intv-constp (car b)))
2287	      (math-infinitep (car b)))
2288	  (math-max-list (math-max a (car b)) (cdr b))
2289	(math-reject-arg (car b) 'anglep))
2290    a))
2291
2292
2293;;; Compute the absolute value of A.  [O O; r r] [Public]
2294(defun math-abs (a)
2295  (cond ((Math-negp a)
2296	 (math-neg a))
2297	((Math-anglep a)
2298	 a)
2299	((eq (car a) 'cplx)
2300	 (math-hypot (nth 1 a) (nth 2 a)))
2301	((eq (car a) 'polar)
2302	 (nth 1 a))
2303	((eq (car a) 'vec)
2304	 (if (cdr (cdr (cdr a)))
2305	     (math-sqrt (calcFunc-abssqr a))
2306	   (if (cdr (cdr a))
2307	       (math-hypot (nth 1 a) (nth 2 a))
2308	     (if (cdr a)
2309		 (math-abs (nth 1 a))
2310	       a))))
2311	((eq (car a) 'sdev)
2312	 (list 'sdev (math-abs (nth 1 a)) (nth 2 a)))
2313	((and (eq (car a) 'intv) (math-intv-constp a))
2314	 (if (Math-posp a)
2315	     a
2316	   (let* ((nlo (math-neg (nth 2 a)))
2317		  (res (math-compare nlo (nth 3 a))))
2318	     (cond ((= res 1)
2319		    (math-make-intv (if (memq (nth 1 a) '(0 1)) 2 3) 0 nlo))
2320		   ((= res 0)
2321		    (math-make-intv (if (eq (nth 1 a) 0) 2 3) 0 nlo))
2322		   (t
2323		    (math-make-intv (if (memq (nth 1 a) '(0 2)) 2 3)
2324				    0 (nth 3 a)))))))
2325	((math-looks-negp a)
2326	 (list 'calcFunc-abs (math-neg a)))
2327	((let ((signs (math-possible-signs a)))
2328	   (or (and (memq signs '(2 4 6)) a)
2329	       (and (memq signs '(1 3)) (math-neg a)))))
2330	((let ((inf (math-infinitep a)))
2331	   (and inf
2332		(if (equal inf '(var nan var-nan))
2333		    inf
2334		  '(var inf var-inf)))))
2335	(t (calc-record-why 'numvecp a)
2336	   (list 'calcFunc-abs a))))
2337
2338(defalias 'calcFunc-abs 'math-abs)
2339
2340(defun math-float-fancy (a)
2341  (cond ((eq (car a) 'intv)
2342	 (cons (car a) (cons (nth 1 a) (mapcar 'math-float (nthcdr 2 a)))))
2343	((and (memq (car a) '(* /))
2344	      (math-numberp (nth 1 a)))
2345	 (list (car a) (math-float (nth 1 a))
2346	       (list 'calcFunc-float (nth 2 a))))
2347	((and (eq (car a) '/)
2348	      (eq (car (nth 1 a)) '*)
2349	      (math-numberp (nth 1 (nth 1 a))))
2350	 (list '* (math-float (nth 1 (nth 1 a)))
2351	       (list 'calcFunc-float (list '/ (nth 2 (nth 1 a)) (nth 2 a)))))
2352	((math-infinitep a) a)
2353	((eq (car a) 'calcFunc-float) a)
2354	((let ((func (assq (car a) '((calcFunc-floor  . calcFunc-ffloor)
2355				     (calcFunc-ceil   . calcFunc-fceil)
2356				     (calcFunc-trunc  . calcFunc-ftrunc)
2357				     (calcFunc-round  . calcFunc-fround)
2358				     (calcFunc-rounde . calcFunc-frounde)
2359				     (calcFunc-roundu . calcFunc-froundu)))))
2360	   (and func (cons (cdr func) (cdr a)))))
2361	(t (math-reject-arg a 'objectp))))
2362
2363(defalias 'calcFunc-float 'math-float)
2364
2365;; The variable math-trunc-prec is local to math-trunc in calc-misc.el,
2366;; but used by math-trunc-fancy which is called by math-trunc.
2367(defvar math-trunc-prec)
2368
2369(defun math-trunc-fancy (a)
2370  (cond ((eq (car a) 'frac) (math-quotient (nth 1 a) (nth 2 a)))
2371	((eq (car a) 'cplx) (math-trunc (nth 1 a)))
2372	((eq (car a) 'polar) (math-trunc (math-complex a)))
2373	((eq (car a) 'hms) (list 'hms (nth 1 a) 0 0))
2374	((eq (car a) 'date) (list 'date (math-trunc (nth 1 a))))
2375	((eq (car a) 'mod)
2376	 (if (math-messy-integerp (nth 2 a))
2377	     (math-trunc (math-make-mod (nth 1 a) (math-trunc (nth 2 a))))
2378	   (math-make-mod (math-trunc (nth 1 a)) (nth 2 a))))
2379	((eq (car a) 'intv)
2380	 (math-make-intv (+ (if (and (equal (nth 2 a) '(neg (var inf var-inf)))
2381				     (memq (nth 1 a) '(0 1)))
2382				0 2)
2383			    (if (and (equal (nth 3 a) '(var inf var-inf))
2384				     (memq (nth 1 a) '(0 2)))
2385				0 1))
2386			 (if (and (Math-negp (nth 2 a))
2387				  (Math-num-integerp (nth 2 a))
2388				  (memq (nth 1 a) '(0 1)))
2389			     (math-add (math-trunc (nth 2 a)) 1)
2390			   (math-trunc (nth 2 a)))
2391			 (if (and (Math-posp (nth 3 a))
2392				  (Math-num-integerp (nth 3 a))
2393				  (memq (nth 1 a) '(0 2)))
2394			     (math-add (math-trunc (nth 3 a)) -1)
2395			   (math-trunc (nth 3 a)))))
2396	((math-provably-integerp a) a)
2397	((Math-vectorp a)
2398	 (math-map-vec (function (lambda (x) (math-trunc x math-trunc-prec))) a))
2399	((math-infinitep a)
2400	 (if (or (math-posp a) (math-negp a))
2401	     a
2402	   '(var nan var-nan)))
2403	((math-to-integer a))
2404	(t (math-reject-arg a 'numberp))))
2405
2406(defun math-trunc-special (a prec)
2407  (if (Math-messy-integerp prec)
2408      (setq prec (math-trunc prec)))
2409  (or (integerp prec)
2410      (math-reject-arg prec 'fixnump))
2411  (if (and (<= prec 0)
2412	   (math-provably-integerp a))
2413      a
2414    (calcFunc-scf (math-trunc (let ((calc-prefer-frac t))
2415				(calcFunc-scf a prec)))
2416		  (- prec))))
2417
2418(defun math-to-integer (a)
2419  (let ((func (assq (car-safe a) '((calcFunc-ffloor  . calcFunc-floor)
2420				   (calcFunc-fceil   . calcFunc-ceil)
2421				   (calcFunc-ftrunc  . calcFunc-trunc)
2422				   (calcFunc-fround  . calcFunc-round)
2423				   (calcFunc-frounde . calcFunc-rounde)
2424				   (calcFunc-froundu . calcFunc-roundu)))))
2425    (and func (= (length a) 2)
2426	 (cons (cdr func) (cdr a)))))
2427
2428(defun calcFunc-ftrunc (a &optional prec)
2429  (if (and (Math-messy-integerp a)
2430	   (or (not prec) (and (integerp prec)
2431			       (<= prec 0))))
2432      a
2433    (math-float (math-trunc a prec))))
2434
2435;; The variable math-floor-prec is local to math-floor in calc-misc.el,
2436;; but used by math-floor-fancy which is called by math-floor.
2437(defvar math-floor-prec)
2438
2439(defun math-floor-fancy (a)
2440  (cond ((math-provably-integerp a) a)
2441	((eq (car a) 'hms)
2442	 (if (or (math-posp a)
2443		 (and (math-zerop (nth 2 a))
2444		      (math-zerop (nth 3 a))))
2445	     (math-trunc a)
2446	   (math-add (math-trunc a) -1)))
2447	((eq (car a) 'date) (list 'date (math-floor (nth 1 a))))
2448	((eq (car a) 'intv)
2449	 (math-make-intv (+ (if (and (equal (nth 2 a) '(neg (var inf var-inf)))
2450				     (memq (nth 1 a) '(0 1)))
2451				0 2)
2452			    (if (and (equal (nth 3 a) '(var inf var-inf))
2453				     (memq (nth 1 a) '(0 2)))
2454				0 1))
2455			 (math-floor (nth 2 a))
2456			 (if (and (Math-num-integerp (nth 3 a))
2457				  (memq (nth 1 a) '(0 2)))
2458			     (math-add (math-floor (nth 3 a)) -1)
2459			   (math-floor (nth 3 a)))))
2460	((Math-vectorp a)
2461	 (math-map-vec (function (lambda (x) (math-floor x math-floor-prec))) a))
2462	((math-infinitep a)
2463	 (if (or (math-posp a) (math-negp a))
2464	     a
2465	   '(var nan var-nan)))
2466	((math-to-integer a))
2467	(t (math-reject-arg a 'anglep))))
2468
2469(defun math-floor-special (a prec)
2470  (if (Math-messy-integerp prec)
2471      (setq prec (math-trunc prec)))
2472  (or (integerp prec)
2473      (math-reject-arg prec 'fixnump))
2474  (if (and (<= prec 0)
2475	   (math-provably-integerp a))
2476      a
2477    (calcFunc-scf (math-floor (let ((calc-prefer-frac t))
2478				(calcFunc-scf a prec)))
2479		  (- prec))))
2480
2481(defun calcFunc-ffloor (a &optional prec)
2482  (if (and (Math-messy-integerp a)
2483	   (or (not prec) (and (integerp prec)
2484			       (<= prec 0))))
2485      a
2486    (math-float (math-floor a prec))))
2487
2488;;; Coerce A to be an integer (by truncation toward plus infinity).  [I N]
2489(defun math-ceiling (a &optional prec)   ;  [Public]
2490  (cond (prec
2491	 (if (Math-messy-integerp prec)
2492	     (setq prec (math-trunc prec)))
2493	 (or (integerp prec)
2494	     (math-reject-arg prec 'fixnump))
2495	 (if (and (<= prec 0)
2496		  (math-provably-integerp a))
2497	     a
2498	   (calcFunc-scf (math-ceiling (let ((calc-prefer-frac t))
2499					 (calcFunc-scf a prec)))
2500			 (- prec))))
2501	((Math-integerp a) a)
2502	((Math-messy-integerp a) (math-trunc a))
2503	((Math-realp a)
2504	 (if (Math-posp a)
2505	     (math-add (math-trunc a) 1)
2506	   (math-trunc a)))
2507	((math-provably-integerp a) a)
2508	((eq (car a) 'hms)
2509	 (if (or (math-negp a)
2510		 (and (math-zerop (nth 2 a))
2511		      (math-zerop (nth 3 a))))
2512	     (math-trunc a)
2513	   (math-add (math-trunc a) 1)))
2514	((eq (car a) 'date) (list 'date (math-ceiling (nth 1 a))))
2515	((eq (car a) 'intv)
2516	 (math-make-intv (+ (if (and (equal (nth 2 a) '(neg (var inf var-inf)))
2517				     (memq (nth 1 a) '(0 1)))
2518				0 2)
2519			    (if (and (equal (nth 3 a) '(var inf var-inf))
2520				     (memq (nth 1 a) '(0 2)))
2521				0 1))
2522			 (if (and (Math-num-integerp (nth 2 a))
2523				  (memq (nth 1 a) '(0 1)))
2524			     (math-add (math-floor (nth 2 a)) 1)
2525			   (math-ceiling (nth 2 a)))
2526			 (math-ceiling (nth 3 a))))
2527	((Math-vectorp a)
2528	 (math-map-vec (function (lambda (x) (math-ceiling x prec))) a))
2529	((math-infinitep a)
2530	 (if (or (math-posp a) (math-negp a))
2531	     a
2532	   '(var nan var-nan)))
2533	((math-to-integer a))
2534	(t (math-reject-arg a 'anglep))))
2535
2536(defalias 'calcFunc-ceil 'math-ceiling)
2537
2538(defun calcFunc-fceil (a &optional prec)
2539  (if (and (Math-messy-integerp a)
2540	   (or (not prec) (and (integerp prec)
2541			       (<= prec 0))))
2542      a
2543    (math-float (math-ceiling a prec))))
2544
2545(defvar math-rounding-mode nil)
2546
2547;;; Coerce A to be an integer (by rounding to nearest integer).  [I N] [Public]
2548(defun math-round (a &optional prec)
2549  (cond (prec
2550	 (if (Math-messy-integerp prec)
2551	     (setq prec (math-trunc prec)))
2552	 (or (integerp prec)
2553	     (math-reject-arg prec 'fixnump))
2554	 (if (and (<= prec 0)
2555		  (math-provably-integerp a))
2556	     a
2557	   (calcFunc-scf (math-round (let ((calc-prefer-frac t))
2558				       (calcFunc-scf a prec)))
2559			 (- prec))))
2560	((Math-anglep a)
2561	 (if (Math-num-integerp a)
2562	     (math-trunc a)
2563	   (if (and (Math-negp a) (not (eq math-rounding-mode 'up)))
2564	       (math-neg (math-round (math-neg a)))
2565	     (setq a (let ((calc-angle-mode 'deg))   ; in case of HMS forms
2566		       (math-add a (if (Math-ratp a)
2567				       '(frac 1 2)
2568				     '(float 5 -1)))))
2569	     (if (and (Math-num-integerp a) (eq math-rounding-mode 'even))
2570		 (progn
2571		   (setq a (math-floor a))
2572		   (or (math-evenp a)
2573		       (setq a (math-sub a 1)))
2574		   a)
2575	       (math-floor a)))))
2576	((math-provably-integerp a) a)
2577	((eq (car a) 'date) (list 'date (math-round (nth 1 a))))
2578	((eq (car a) 'intv)
2579	 (math-floor (math-add a '(frac 1 2))))
2580	((Math-vectorp a)
2581	 (math-map-vec (function (lambda (x) (math-round x prec))) a))
2582	((math-infinitep a)
2583	 (if (or (math-posp a) (math-negp a))
2584	     a
2585	   '(var nan var-nan)))
2586	((math-to-integer a))
2587	(t (math-reject-arg a 'anglep))))
2588
2589(defalias 'calcFunc-round 'math-round)
2590
2591(defsubst calcFunc-rounde (a &optional prec)
2592  (let ((math-rounding-mode 'even))
2593    (math-round a prec)))
2594
2595(defsubst calcFunc-roundu (a &optional prec)
2596  (let ((math-rounding-mode 'up))
2597    (math-round a prec)))
2598
2599(defun calcFunc-fround (a &optional prec)
2600  (if (and (Math-messy-integerp a)
2601	   (or (not prec) (and (integerp prec)
2602			       (<= prec 0))))
2603      a
2604    (math-float (math-round a prec))))
2605
2606(defsubst calcFunc-frounde (a &optional prec)
2607  (let ((math-rounding-mode 'even))
2608    (calcFunc-fround a prec)))
2609
2610(defsubst calcFunc-froundu (a &optional prec)
2611  (let ((math-rounding-mode 'up))
2612    (calcFunc-fround a prec)))
2613
2614;;; Pull floating-point values apart into mantissa and exponent.
2615(defun calcFunc-mant (x)
2616  (if (Math-realp x)
2617      (if (or (Math-ratp x)
2618	      (eq (nth 1 x) 0))
2619	  x
2620	(list 'float (nth 1 x) (- 1 (math-numdigs (nth 1 x)))))
2621    (calc-record-why 'realp x)
2622    (list 'calcFunc-mant x)))
2623
2624(defun calcFunc-xpon (x)
2625  (if (Math-realp x)
2626      (if (or (Math-ratp x)
2627	      (eq (nth 1 x) 0))
2628	  0
2629	(math-normalize (+ (nth 2 x) (1- (math-numdigs (nth 1 x))))))
2630    (calc-record-why 'realp x)
2631    (list 'calcFunc-xpon x)))
2632
2633(defun calcFunc-scf (x n)
2634  (if (integerp n)
2635      (cond ((eq n 0)
2636	     x)
2637	    ((Math-integerp x)
2638	     (if (> n 0)
2639		 (math-scale-int x n)
2640	       (math-div x (math-scale-int 1 (- n)))))
2641	    ((eq (car x) 'frac)
2642	     (if (> n 0)
2643		 (math-make-frac (math-scale-int (nth 1 x) n) (nth 2 x))
2644	       (math-make-frac (nth 1 x) (math-scale-int (nth 2 x) (- n)))))
2645	    ((eq (car x) 'float)
2646	     (math-make-float (nth 1 x) (+ (nth 2 x) n)))
2647	    ((memq (car x) '(cplx sdev))
2648	     (math-normalize
2649	      (list (car x)
2650		    (calcFunc-scf (nth 1 x) n)
2651		    (calcFunc-scf (nth 2 x) n))))
2652	    ((memq (car x) '(polar mod))
2653	     (math-normalize
2654	      (list (car x)
2655		    (calcFunc-scf (nth 1 x) n)
2656		    (nth 2 x))))
2657	    ((eq (car x) 'intv)
2658	     (math-normalize
2659	      (list (car x)
2660		    (nth 1 x)
2661		    (calcFunc-scf (nth 2 x) n)
2662		    (calcFunc-scf (nth 3 x) n))))
2663	    ((eq (car x) 'vec)
2664	     (math-map-vec (function (lambda (x) (calcFunc-scf x n))) x))
2665	    ((math-infinitep x)
2666	     x)
2667	    (t
2668	     (calc-record-why 'realp x)
2669	     (list 'calcFunc-scf x n)))
2670    (if (math-messy-integerp n)
2671	(if (< (nth 2 n) 10)
2672	    (calcFunc-scf x (math-trunc n))
2673	  (math-overflow n))
2674      (if (math-integerp n)
2675	  (math-overflow n)
2676	(calc-record-why 'integerp n)
2677	(list 'calcFunc-scf x n)))))
2678
2679
2680(defun calcFunc-incr (x &optional step relative-to)
2681  (or step (setq step 1))
2682  (cond ((not (Math-integerp step))
2683	 (math-reject-arg step 'integerp))
2684	((Math-integerp x)
2685	 (math-add x step))
2686	((eq (car x) 'float)
2687	 (if (and (math-zerop x)
2688		  (eq (car-safe relative-to) 'float))
2689	     (math-mul step
2690		       (calcFunc-scf relative-to (- 1 calc-internal-prec)))
2691	   (math-add-float x (math-make-float
2692			      step
2693			      (+ (nth 2 x)
2694				 (- (math-numdigs (nth 1 x))
2695				    calc-internal-prec))))))
2696	((eq (car x) 'date)
2697	 (if (Math-integerp (nth 1 x))
2698	     (math-add x step)
2699	   (math-add x (list 'hms 0 0 step))))
2700	(t
2701	 (math-reject-arg x 'realp))))
2702
2703(defsubst calcFunc-decr (x &optional step relative-to)
2704  (calcFunc-incr x (math-neg (or step 1)) relative-to))
2705
2706(defun calcFunc-percent (x)
2707  (if (math-objectp x)
2708      (let ((calc-prefer-frac nil))
2709	(math-div x 100))
2710    (list 'calcFunc-percent x)))
2711
2712(defun calcFunc-relch (x y)
2713  (if (and (math-objectp x) (math-objectp y))
2714      (math-div (math-sub y x) x)
2715    (list 'calcFunc-relch x y)))
2716
2717;;; Compute the absolute value squared of A.  [F N] [Public]
2718(defun calcFunc-abssqr (a)
2719  (cond ((Math-realp a)
2720	 (math-mul a a))
2721	((eq (car a) 'cplx)
2722	 (math-add (math-sqr (nth 1 a))
2723		   (math-sqr (nth 2 a))))
2724	((eq (car a) 'polar)
2725	 (math-sqr (nth 1 a)))
2726	((and (memq (car a) '(sdev intv)) (math-constp a))
2727	 (math-sqr (math-abs a)))
2728	((eq (car a) 'vec)
2729	 (math-reduce-vec 'math-add (math-map-vec 'calcFunc-abssqr a)))
2730	((math-known-realp a)
2731	 (math-pow a 2))
2732	((let ((inf (math-infinitep a)))
2733	   (and inf
2734		(math-mul (calcFunc-abssqr (math-infinite-dir a inf)) inf))))
2735	(t (calc-record-why 'numvecp a)
2736	   (list 'calcFunc-abssqr a))))
2737
2738(defsubst math-sqr (a)
2739  (math-mul a a))
2740
2741;;;; Number theory.
2742
2743(defun calcFunc-idiv (a b)   ; [I I I] [Public]
2744  (cond ((and (Math-natnump a) (Math-natnump b) (not (eq b 0)))
2745	 (math-quotient a b))
2746	((Math-realp a)
2747	 (if (Math-realp b)
2748	     (let ((calc-prefer-frac t))
2749	       (math-floor (math-div a b)))
2750	   (math-reject-arg b 'realp)))
2751	((eq (car-safe a) 'hms)
2752	 (if (eq (car-safe b) 'hms)
2753	     (let ((calc-prefer-frac t))
2754	       (math-floor (math-div a b)))
2755	   (math-reject-arg b 'hmsp)))
2756	((and (or (eq (car-safe a) 'intv) (Math-realp a))
2757	      (or (eq (car-safe b) 'intv) (Math-realp b)))
2758	 (math-floor (math-div a b)))
2759	((or (math-infinitep a)
2760	     (math-infinitep b))
2761	 (math-div a b))
2762	(t (math-reject-arg a 'anglep))))
2763
2764
2765;;; Combine two terms being added, if possible.
2766(defun math-combine-sum (a b nega negb scalar-okay)
2767  (if (and scalar-okay (Math-objvecp a) (Math-objvecp b))
2768      (math-add-or-sub a b nega negb)
2769    (let ((amult 1) (bmult 1))
2770      (and (consp a)
2771	   (cond ((and (eq (car a) '*)
2772		       (Math-objectp (nth 1 a)))
2773		  (setq amult (nth 1 a)
2774			a (nth 2 a)))
2775		 ((and (eq (car a) '/)
2776		       (Math-objectp (nth 2 a)))
2777		  (setq amult (if (Math-integerp (nth 2 a))
2778				  (list 'frac 1 (nth 2 a))
2779				(math-div 1 (nth 2 a)))
2780			a (nth 1 a)))
2781		 ((eq (car a) 'neg)
2782		  (setq amult -1
2783			a (nth 1 a)))))
2784      (and (consp b)
2785	   (cond ((and (eq (car b) '*)
2786		       (Math-objectp (nth 1 b)))
2787		  (setq bmult (nth 1 b)
2788			b (nth 2 b)))
2789		 ((and (eq (car b) '/)
2790		       (Math-objectp (nth 2 b)))
2791		  (setq bmult (if (Math-integerp (nth 2 b))
2792				  (list 'frac 1 (nth 2 b))
2793				(math-div 1 (nth 2 b)))
2794			b (nth 1 b)))
2795		 ((eq (car b) 'neg)
2796		  (setq bmult -1
2797			b (nth 1 b)))))
2798      (and (if math-simplifying
2799	       (Math-equal a b)
2800	     (equal a b))
2801	   (progn
2802	     (if nega (setq amult (math-neg amult)))
2803	     (if negb (setq bmult (math-neg bmult)))
2804	     (setq amult (math-add amult bmult))
2805	     (math-mul amult a))))))
2806
2807(defun math-add-or-sub (a b aneg bneg)
2808  (if aneg (setq a (math-neg a)))
2809  (if bneg (setq b (math-neg b)))
2810  (if (or (Math-vectorp a) (Math-vectorp b))
2811      (math-normalize (list '+ a b))
2812    (math-add a b)))
2813
2814(defvar math-combine-prod-e '(var e var-e))
2815
2816;;; The following is expanded out four ways for speed.
2817
2818;; math-unit-prefixes is defined in calc-units.el,
2819;; but used here.
2820(defvar math-unit-prefixes)
2821
2822(defun math-combine-prod (a b inva invb scalar-okay)
2823  (cond
2824   ((or (and inva (Math-zerop a))
2825	(and invb (Math-zerop b)))
2826    nil)
2827   ((and scalar-okay (Math-objvecp a) (Math-objvecp b))
2828    (setq a (math-mul-or-div a b inva invb))
2829    (and (Math-objvecp a)
2830	 a))
2831   ((and (eq (car-safe a) '^)
2832	 inva
2833	 (math-looks-negp (nth 2 a)))
2834    (math-mul (math-pow (nth 1 a) (math-neg (nth 2 a))) b))
2835   ((and (eq (car-safe b) '^)
2836	 invb
2837	 (math-looks-negp (nth 2 b)))
2838    (math-mul a (math-pow (nth 1 b) (math-neg (nth 2 b)))))
2839   ((and math-simplifying
2840         (math-combine-prod-trig a b)))
2841   (t (let ((apow 1) (bpow 1))
2842	(and (consp a)
2843	     (cond ((and (eq (car a) '^)
2844			 (or math-simplifying
2845			     (Math-numberp (nth 2 a))))
2846		    (setq apow (nth 2 a)
2847			  a (nth 1 a)))
2848		   ((eq (car a) 'calcFunc-sqrt)
2849		    (setq apow '(frac 1 2)
2850			  a (nth 1 a)))
2851		   ((and (eq (car a) 'calcFunc-exp)
2852			 (or math-simplifying
2853			     (Math-numberp (nth 1 a))))
2854		    (setq apow (nth 1 a)
2855			  a math-combine-prod-e))))
2856	(and (consp a) (eq (car a) 'frac)
2857	     (Math-lessp (nth 1 a) (nth 2 a))
2858	     (setq a (math-div 1 a) apow (math-neg apow)))
2859	(and (consp b)
2860	     (cond ((and (eq (car b) '^)
2861			 (or math-simplifying
2862			     (Math-numberp (nth 2 b))))
2863		    (setq bpow (nth 2 b)
2864			  b (nth 1 b)))
2865		   ((eq (car b) 'calcFunc-sqrt)
2866		    (setq bpow '(frac 1 2)
2867			  b (nth 1 b)))
2868		   ((and (eq (car b) 'calcFunc-exp)
2869			 (or math-simplifying
2870			     (Math-numberp (nth 1 b))))
2871		    (setq bpow (nth 1 b)
2872			  b math-combine-prod-e))))
2873	(and (consp b) (eq (car b) 'frac)
2874	     (Math-lessp (nth 1 b) (nth 2 b))
2875	     (setq b (math-div 1 b) bpow (math-neg bpow)))
2876	(if inva (setq apow (math-neg apow)))
2877	(if invb (setq bpow (math-neg bpow)))
2878	(or (and (if math-simplifying
2879		     (math-commutative-equal a b)
2880		   (equal a b))
2881		 (let ((sumpow (math-add apow bpow)))
2882		   (and (or (not (Math-integerp a))
2883			    (Math-zerop sumpow)
2884			    (eq (eq (car-safe apow) 'frac)
2885				(eq (car-safe bpow) 'frac)))
2886			(progn
2887			  (and (math-looks-negp sumpow)
2888			       (Math-ratp a) (Math-posp a)
2889			       (setq a (math-div 1 a)
2890				     sumpow (math-neg sumpow)))
2891			  (cond ((equal sumpow '(frac 1 2))
2892				 (list 'calcFunc-sqrt a))
2893				((equal sumpow '(frac -1 2))
2894				 (math-div 1 (list 'calcFunc-sqrt a)))
2895				((and (eq a math-combine-prod-e)
2896				      (eq a b))
2897				 (list 'calcFunc-exp sumpow))
2898				(t
2899				 (condition-case err
2900				     (math-pow a sumpow)
2901				   (inexact-result (list '^ a sumpow)))))))))
2902	    (and math-simplifying-units
2903		 math-combining-units
2904		 (let* ((ua (math-check-unit-name a))
2905			ub)
2906		   (and ua
2907			(eq ua (setq ub (math-check-unit-name b)))
2908			(progn
2909			  (setq ua (if (eq (nth 1 a) (car ua))
2910				       1
2911				     (nth 1 (assq (aref (symbol-name (nth 1 a))
2912							0)
2913						  math-unit-prefixes)))
2914				ub (if (eq (nth 1 b) (car ub))
2915				       1
2916				     (nth 1 (assq (aref (symbol-name (nth 1 b))
2917							0)
2918						  math-unit-prefixes))))
2919			  (if (Math-lessp ua ub)
2920			      (let (temp)
2921				(setq temp a a b b temp
2922				      temp ua ua ub ub temp
2923				      temp apow apow bpow bpow temp)))
2924			  (math-mul (math-pow (math-div ua ub) apow)
2925				    (math-pow b (math-add apow bpow)))))))
2926	    (and (equal apow bpow)
2927		 (Math-natnump a) (Math-natnump b)
2928		 (cond ((equal apow '(frac 1 2))
2929			(list 'calcFunc-sqrt (math-mul a b)))
2930		       ((equal apow '(frac -1 2))
2931			(math-div 1 (list 'calcFunc-sqrt (math-mul a b))))
2932		       (t
2933			(setq a (math-mul a b))
2934			(condition-case err
2935			    (math-pow a apow)
2936			  (inexact-result (list '^ a apow)))))))))))
2937
2938(defun math-combine-prod-trig (a b)
2939  (cond
2940   ((and (eq (car-safe a) 'calcFunc-sin)
2941         (eq (car-safe b) 'calcFunc-csc)
2942         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2943    1)
2944   ((and (eq (car-safe a) 'calcFunc-sin)
2945         (eq (car-safe b) 'calcFunc-sec)
2946         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2947    (cons 'calcFunc-tan (cdr a)))
2948   ((and (eq (car-safe a) 'calcFunc-sin)
2949         (eq (car-safe b) 'calcFunc-cot)
2950         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2951    (cons 'calcFunc-cos (cdr a)))
2952   ((and (eq (car-safe a) 'calcFunc-cos)
2953         (eq (car-safe b) 'calcFunc-sec)
2954         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2955    1)
2956   ((and (eq (car-safe a) 'calcFunc-cos)
2957         (eq (car-safe b) 'calcFunc-csc)
2958         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2959    (cons 'calcFunc-cot (cdr a)))
2960   ((and (eq (car-safe a) 'calcFunc-cos)
2961         (eq (car-safe b) 'calcFunc-tan)
2962         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2963    (cons 'calcFunc-sin (cdr a)))
2964   ((and (eq (car-safe a) 'calcFunc-tan)
2965         (eq (car-safe b) 'calcFunc-cot)
2966         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2967    1)
2968   ((and (eq (car-safe a) 'calcFunc-tan)
2969         (eq (car-safe b) 'calcFunc-csc)
2970         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2971    (cons 'calcFunc-sec (cdr a)))
2972   ((and (eq (car-safe a) 'calcFunc-sec)
2973         (eq (car-safe b) 'calcFunc-cot)
2974         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2975    (cons 'calcFunc-csc (cdr a)))
2976   ((and (eq (car-safe a) 'calcFunc-sinh)
2977         (eq (car-safe b) 'calcFunc-csch)
2978         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2979    1)
2980   ((and (eq (car-safe a) 'calcFunc-sinh)
2981         (eq (car-safe b) 'calcFunc-sech)
2982         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2983    (cons 'calcFunc-tanh (cdr a)))
2984   ((and (eq (car-safe a) 'calcFunc-sinh)
2985         (eq (car-safe b) 'calcFunc-coth)
2986         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2987    (cons 'calcFunc-cosh (cdr a)))
2988   ((and (eq (car-safe a) 'calcFunc-cosh)
2989         (eq (car-safe b) 'calcFunc-sech)
2990         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2991    1)
2992   ((and (eq (car-safe a) 'calcFunc-cosh)
2993         (eq (car-safe b) 'calcFunc-csch)
2994         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2995    (cons 'calcFunc-coth (cdr a)))
2996   ((and (eq (car-safe a) 'calcFunc-cosh)
2997         (eq (car-safe b) 'calcFunc-tanh)
2998         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
2999    (cons 'calcFunc-sinh (cdr a)))
3000   ((and (eq (car-safe a) 'calcFunc-tanh)
3001         (eq (car-safe b) 'calcFunc-coth)
3002         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
3003    1)
3004   ((and (eq (car-safe a) 'calcFunc-tanh)
3005         (eq (car-safe b) 'calcFunc-csch)
3006         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
3007    (cons 'calcFunc-sech (cdr a)))
3008   ((and (eq (car-safe a) 'calcFunc-sech)
3009         (eq (car-safe b) 'calcFunc-coth)
3010         (= 0 (math-simplify (math-sub (cdr a) (cdr b)))))
3011    (cons 'calcFunc-csch (cdr a)))
3012   (t
3013    nil)))
3014
3015(defun math-mul-or-div (a b ainv binv)
3016  (if (or (Math-vectorp a) (Math-vectorp b))
3017      (math-normalize
3018       (if ainv
3019	   (if binv
3020	       (list '/ (math-div 1 a) b)
3021	     (list '/ b a))
3022	 (if binv
3023	     (list '/ a b)
3024	   (list '* a b))))
3025    (if ainv
3026	(if binv
3027	    (math-div (math-div 1 a) b)
3028	  (math-div b a))
3029      (if binv
3030	  (math-div a b)
3031	(math-mul a b)))))
3032
3033;; The variable math-com-bterms is local to math-commutative-equal,
3034;; but is used by math-commutative collect, which is called by
3035;; math-commutative-equal.
3036(defvar math-com-bterms)
3037
3038(defun math-commutative-equal (a b)
3039  (if (memq (car-safe a) '(+ -))
3040      (and (memq (car-safe b) '(+ -))
3041	   (let ((math-com-bterms nil) aterms p)
3042	     (math-commutative-collect b nil)
3043	     (setq aterms math-com-bterms math-com-bterms nil)
3044	     (math-commutative-collect a nil)
3045	     (and (= (length aterms) (length math-com-bterms))
3046		  (progn
3047		    (while (and aterms
3048				(progn
3049				  (setq p math-com-bterms)
3050				  (while (and p (not (equal (car aterms)
3051							    (car p))))
3052				    (setq p (cdr p)))
3053				  p))
3054		      (setq math-com-bterms (delq (car p) math-com-bterms)
3055			    aterms (cdr aterms)))
3056		    (not aterms)))))
3057    (equal a b)))
3058
3059(defun math-commutative-collect (b neg)
3060  (if (eq (car-safe b) '+)
3061      (progn
3062	(math-commutative-collect (nth 1 b) neg)
3063	(math-commutative-collect (nth 2 b) neg))
3064    (if (eq (car-safe b) '-)
3065	(progn
3066	  (math-commutative-collect (nth 1 b) neg)
3067	  (math-commutative-collect (nth 2 b) (not neg)))
3068      (setq math-com-bterms (cons (if neg (math-neg b) b) math-com-bterms)))))
3069
3070(provide 'calc-arith)
3071
3072;;; arch-tag: 6c396b5b-14c6-40ed-bb2a-7cc2e8111465
3073;;; calc-arith.el ends here
3074