[bpt/emacs.git] / lisp / calc / calc-stat.el

;;; calc-stat.el --- statistical functions for Calc

;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
;;   2005, 2006, 2007 Free Software Foundation, Inc.

;; Author: David Gillespie <daveg@synaptics.com>
;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>

;; This file is part of GNU Emacs.

;; GNU Emacs is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2, or (at your option)
;; any later version.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
;; along with GNU Emacs; see the file COPYING.  If not, write to the
;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
;; Boston, MA 02110-1301, USA.

;;; Commentary:

;;; Code:

;; This file is autoloaded from calc-ext.el.

(require 'calc-ext)
(require 'calc-macs)

;;; Statistical operations on vectors.

(defun calc-vector-count (arg)
  (interactive "P")
  (calc-slow-wrapper
   (calc-vector-op "coun" 'calcFunc-vcount arg)))

(defun calc-vector-sum (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (calc-vector-op "vprd" 'calcFunc-vprod arg)
     (calc-vector-op "vsum" 'calcFunc-vsum arg))))

(defun calc-vector-product (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-vector-sum arg))

(defun calc-vector-max (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-vector-op "vmin" 'calcFunc-vmin arg)
     (calc-vector-op "vmax" 'calcFunc-vmax arg))))

(defun calc-vector-min (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-vector-max arg))

(defun calc-vector-mean (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-vector-op "harm" 'calcFunc-vhmean arg)
	 (calc-vector-op "medn" 'calcFunc-vmedian arg))
     (if (calc-is-inverse)
	 (calc-vector-op "meae" 'calcFunc-vmeane arg)
       (calc-vector-op "mean" 'calcFunc-vmean arg)))))

(defun calc-vector-mean-error (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-vector-mean arg))

(defun calc-vector-median (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-vector-mean arg))

(defun calc-vector-harmonic-mean (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-vector-mean arg))

(defun calc-vector-geometric-mean (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (calc-binary-op "geom" 'calcFunc-agmean arg)
     (calc-vector-op "geom" 'calcFunc-vgmean arg))))

(defun calc-vector-sdev (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-vector-op "pvar" 'calcFunc-vpvar arg)
	 (calc-vector-op "var" 'calcFunc-vvar arg))
     (if (calc-is-inverse)
	 (calc-vector-op "psdv" 'calcFunc-vpsdev arg)
       (calc-vector-op "sdev" 'calcFunc-vsdev arg)))))

(defun calc-vector-pop-sdev (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-vector-sdev arg))

(defun calc-vector-variance (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-vector-sdev arg))

(defun calc-vector-pop-variance (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-vector-sdev arg))

(defun calc-vector-covariance (arg)
  (interactive "P")
  (calc-slow-wrapper
   (let ((n (if (eq arg 1) 1 2)))
     (if (calc-is-hyperbolic)
	 (calc-enter-result n "corr" (cons 'calcFunc-vcorr
					   (calc-top-list-n n)))
       (if (calc-is-inverse)
	   (calc-enter-result n "pcov" (cons 'calcFunc-vpcov
					     (calc-top-list-n n)))
	 (calc-enter-result n "cov" (cons 'calcFunc-vcov
					  (calc-top-list-n n))))))))

(defun calc-vector-pop-covariance (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-vector-covariance arg))

(defun calc-vector-correlation (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-vector-covariance arg))

(defun calc-vector-op (name func arg)
  (setq calc-aborted-prefix name
	arg (prefix-numeric-value arg))
  (if (< arg 0)
      (error "Negative arguments not allowed"))
  (calc-enter-result arg name (cons func (calc-top-list-n arg))))


;;; Useful statistical functions

;;; Sum, product, etc., of one or more values or vectors.
;;; Each argument must be either a number or a vector.  Vectors
;;; are flattened, but variables inside are assumed to represent
;;; non-vectors.

(defun calcFunc-vsum (&rest vecs)
  (math-reduce-many-vecs 'calcFunc-add 'calcFunc-vsum vecs 0))

(defun calcFunc-vprod (&rest vecs)
  (math-reduce-many-vecs 'calcFunc-mul 'calcFunc-vprod vecs 1))

(defun calcFunc-vmax (&rest vecs)
  (if (eq (car-safe (car vecs)) 'sdev)
      '(var inf var-inf)
    (if (eq (car-safe (car vecs)) 'intv)
	(nth 3 (math-fix-int-intv (car vecs)))
      (math-reduce-many-vecs 'calcFunc-max 'calcFunc-vmax vecs
			     '(neg (var inf var-inf))))))

(defun calcFunc-vmin (&rest vecs)
  (if (eq (car-safe (car vecs)) 'sdev)
      '(neg (var inf var-inf))
    (if (eq (car-safe (car vecs)) 'intv)
	(nth 2 (math-fix-int-intv (car vecs)))
      (math-reduce-many-vecs 'calcFunc-min 'calcFunc-vmin vecs
			     '(var inf var-inf)))))

(defun math-reduce-many-vecs (func whole-func vecs ident)
  (let ((const-part nil)
	(symb-part nil)
	val vec)
    (let ((calc-internal-prec (+ calc-internal-prec 2)))
      (while vecs
	(setq val (car vecs))
	(and (eq (car-safe val) 'var)
	     (eq (car-safe (calc-var-value (nth 2 val))) 'vec)
	     (setq val (symbol-value (nth 2 val))))
	(cond ((Math-vectorp val)
	       (setq vec (append (and const-part (list const-part))
				 (math-flatten-vector val)))
	       (setq const-part (if vec
				    (calcFunc-reducer
				     (math-calcFunc-to-var func)
				     (cons 'vec vec))
				  ident)))
	      ((or (Math-objectp val) (math-infinitep val))
	       (setq const-part (if const-part
				    (funcall func const-part val)
				  val)))
	      (t
	       (setq symb-part (nconc symb-part (list val)))))
	(setq vecs (cdr vecs))))
    (if const-part
	(progn
	  (setq const-part (math-normalize const-part))
	  (if symb-part
	      (funcall func const-part (cons whole-func symb-part))
	    const-part))
      (if symb-part (cons whole-func symb-part) ident))))


;;; Return the number of data elements among the arguments.
(defun calcFunc-vcount (&rest vecs)
  (let ((count 0))
    (while vecs
      (setq count (if (Math-vectorp (car vecs))
		      (+ count (math-count-elements (car vecs)))
		    (if (Math-objectp (car vecs))
			(1+ count)
		      (if (and (eq (car-safe (car vecs)) 'var)
			       (eq (car-safe (calc-var-value
					      (nth 2 (car vecs))))
				   'vec))
			  (+ count (math-count-elements
				    (symbol-value (nth 2 (car vecs)))))
			(math-reject-arg (car vecs) 'numvecp))))
	    vecs (cdr vecs)))
    count))

(defun math-count-elements (vec)
  (let ((count 0))
    (while (setq vec (cdr vec))
      (setq count (if (Math-vectorp (car vec))
		      (+ count (math-count-elements (car vec)))
		    (1+ count))))
    count))


(defun math-flatten-many-vecs (vecs)
  (let ((p vecs)
	(vec (list 'vec)))
    (while p
      (setq vec (nconc vec
		       (if (Math-vectorp (car p))
			   (math-flatten-vector (car p))
			 (if (Math-objectp (car p))
			     (list (car p))
			   (if (and (eq (car-safe (car p)) 'var)
				    (eq (car-safe (calc-var-value
						   (nth 2 (car p)))) 'vec))
			       (math-flatten-vector (symbol-value
						     (nth 2 (car p))))
			     (math-reject-arg (car p) 'numvecp)))))
	    p (cdr p)))
    vec))

(defun calcFunc-vflat (&rest vecs)
  (math-flatten-many-vecs vecs))

(defun math-split-sdev-vec (vec zero-ok)
  (let ((means (list 'vec))
	(wts (list 'vec))
	(exact nil)
	(p vec))
    (while (and (setq p (cdr p))
		(not (and (consp (car p))
			  (eq (car (car p)) 'sdev)))))
    (if (null p)
	(list vec nil)
      (while (setq vec (cdr vec))
	(if (and (consp (setq p (car vec)))
		 (eq (car p) 'sdev))
	    (or exact
		(setq means (cons (nth 1 p) means)
		      wts (cons (nth 2 p) wts)))
	  (if zero-ok
	      (setq means (cons (nth 1 p) means)
		    wts (cons 0 wts))
	    (or exact
		(setq means (list 'vec)
		      wts nil
		      exact t))
	    (setq means (cons p means)))))
      (list (nreverse means)
	    (and wts (nreverse wts))))))


;;; Return the arithmetic mean of the argument numbers or vectors.
;;; (If numbers are error forms, computes the weighted mean.)
(defun calcFunc-vmean (&rest vecs)
  (let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
	 (means (car split))
	 (wts (nth 1 split))
	 (len (1- (length means))))
    (if (= len 0)
	(math-reject-arg nil "*Must be at least 1 argument")
      (if (and (= len 1) (eq (car-safe (nth 1 means)) 'intv))
	  (let ((x (math-fix-int-intv (nth 1 means))))
	    (calcFunc-vmean (nth 2 x) (nth 3 x)))
	(math-with-extra-prec 2
	  (if (and wts (> len 1))
	      (let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
		     (suminvsqrwts (calcFunc-reduce
				    '(var add var-add)
				    (calcFunc-map '(var div var-div)
						  1 sqrwts))))
		(math-div (calcFunc-reduce '(var add var-add)
					   (calcFunc-map '(var div var-div)
							 means sqrwts))
			  suminvsqrwts))
	    (math-div (calcFunc-reduce '(var add var-add) means) len)))))))

(defun math-fix-int-intv (x)
  (if (math-floatp x)
      x
    (list 'intv 3
	  (if (memq (nth 1 x) '(2 3)) (nth 2 x) (math-add (nth 2 x) 1))
	  (if (memq (nth 1 x) '(1 3)) (nth 3 x) (math-sub (nth 3 x) 1)))))

;;; Compute the mean with an error estimate.
(defun calcFunc-vmeane (&rest vecs)
  (let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
	 (means (car split))
	 (wts (nth 1 split))
	 (len (1- (length means))))
    (if (= len 0)
	(math-reject-arg nil "*Must be at least 1 argument")
      (math-with-extra-prec 2
	(if wts
	    (let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
		   (suminvsqrwts (calcFunc-reduce
				  '(var add var-add)
				  (calcFunc-map '(var div var-div)
						1 sqrwts))))
	      (math-make-sdev
	       (math-div (calcFunc-reduce '(var add var-add)
					  (calcFunc-map '(var div var-div)
							means sqrwts))
			 suminvsqrwts)
	       (list 'calcFunc-sqrt (math-div 1 suminvsqrwts))))
	  (let ((mean (math-div (calcFunc-reduce '(var add var-add) means)
				len)))
	    (math-make-sdev
	     mean
	     (list 'calcFunc-sqrt
		   (math-div (calcFunc-reducer
			      '(var add var-add)
			      (calcFunc-map '(var pow var-pow)
					    (calcFunc-map '(var abs var-abs)
							  (calcFunc-map
							   '(var add var-add)
							   means
							   (math-neg mean)))
					    2))
			     (math-mul len (1- len)))))))))))


;;; Compute the median of a list of values.
(defun calcFunc-vmedian (&rest vecs)
  (let* ((flat (copy-sequence (cdr (math-flatten-many-vecs vecs))))
	 (p flat)
	 (len (length flat))
	 (hlen (/ len 2)))
    (if (= len 0)
	(math-reject-arg nil "*Must be at least 1 argument")
      (if (and (= len 1) (memq (car-safe (car flat)) '(sdev intv)))
	  (calcFunc-vmean (car flat))
	(while p
	  (if (eq (car-safe (car p)) 'sdev)
	      (setcar p (nth 1 (car p))))
	  (or (Math-anglep (car p))
	      (math-reject-arg (car p) 'anglep))
	  (setq p (cdr p)))
	(setq flat (sort flat 'math-lessp))
	(if (= (% len 2) 0)
	    (math-div (math-add (nth (1- hlen) flat) (nth hlen flat)) 2)
	  (nth hlen flat))))))


(defun calcFunc-vgmean (&rest vecs)
  (let* ((flat (math-flatten-many-vecs vecs))
	 (len (1- (length flat))))
    (if (= len 0)
	(math-reject-arg nil "*Must be at least 1 argument")
      (math-with-extra-prec 2
	(let ((x (calcFunc-reduce '(var mul math-mul) flat)))
	  (if (= len 2)
	      (math-sqrt x)
	    (math-pow x (list 'frac 1 len))))))))


(defun calcFunc-agmean (a b)
  (cond ((Math-equal a b) a)
	((math-zerop a) a)
	((math-zerop b) b)
	(calc-symbolic-mode (math-inexact-result))
	((not (Math-realp a)) (math-reject-arg a 'realp))
	((not (Math-realp b)) (math-reject-arg b 'realp))
	(t
	 (math-with-extra-prec 2
	   (setq a (math-float (math-abs a))
		 b (math-float (math-abs b)))
	   (let (mean)
	     (while (not (math-nearly-equal-float a b))
	       (setq mean (math-mul-float (math-add-float a b) '(float 5 -1))
		     b (math-sqrt-float (math-mul-float a b))
		     a mean))
	     a)))))


(defun calcFunc-vhmean (&rest vecs)
  (let* ((flat (math-flatten-many-vecs vecs))
	 (len (1- (length flat))))
    (if (= len 0)
	(math-reject-arg nil "*Must be at least 1 argument")
      (math-with-extra-prec 2
	(math-div len
		  (calcFunc-reduce '(var add math-add)
				   (calcFunc-map '(var inv var-inv) flat)))))))


;;; Compute the sample variance or standard deviation of numbers or vectors.
;;; (If the numbers are error forms, only the mean part of them is used.)
(defun calcFunc-vvar (&rest vecs)
  (if (and (= (length vecs) 1)
	   (memq (car-safe (car vecs)) '(sdev intv)))
      (if (eq (car-safe (car vecs)) 'intv)
	  (math-intv-variance (car vecs) nil)
	(math-sqr (nth 2 (car vecs))))
    (math-covariance vecs nil nil 0)))

(defun calcFunc-vsdev (&rest vecs)
  (if (and (= (length vecs) 1)
	   (memq (car-safe (car vecs)) '(sdev intv)))
      (if (eq (car-safe (car vecs)) 'intv)
	  (if (math-floatp (car vecs))
	      (math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
			(math-sqrt-12))
	    (math-sqrt (calcFunc-vvar (car vecs))))
	(nth 2 (car vecs)))
    (math-sqrt (math-covariance vecs nil nil 0))))

;;; Compute the population variance or std deviation of numbers or vectors.
(defun calcFunc-vpvar (&rest vecs)
  (if (and (= (length vecs) 1)
	   (memq (car-safe (car vecs)) '(sdev intv)))
      (if (eq (car-safe (car vecs)) 'intv)
	  (math-intv-variance (car vecs) t)
	(math-sqr (nth 2 (car vecs))))
    (math-covariance vecs nil t 0)))

(defun calcFunc-vpsdev (&rest vecs)
  (if (and (= (length vecs) 1)
	   (memq (car-safe (car vecs)) '(sdev intv)))
      (if (eq (car-safe (car vecs)) 'intv)
	  (if (math-floatp (car vecs))
	      (math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
			(math-sqrt-12))
	    (math-sqrt (calcFunc-vpvar (car vecs))))
	(nth 2 (car vecs)))
    (math-sqrt (math-covariance vecs nil t 0))))

(defun math-intv-variance (x pop)
  (or (math-constp x) (math-reject-arg x 'constp))
  (if (math-floatp x)
      (math-div (math-sqr (math-sub (nth 3 x) (nth 2 x))) 12)
    (let* ((x (math-fix-int-intv x))
	   (len (math-sub (nth 3 x) (nth 2 x)))
	   (hlen (math-quotient len 2)))
      (math-div (if (math-evenp len)
		    (calcFunc-sum '(^ (var X var-X) 2) '(var X var-X)
				  (math-neg hlen) hlen)
		  (calcFunc-sum '(^ (- (var X var-X) (/ 1 2)) 2)
				'(var X var-X)
				(math-neg hlen) (math-add hlen 1)))
		(if pop (math-add len 1) len)))))

;;; Compute the covariance and linear correlation coefficient.
(defun calcFunc-vcov (vec1 &optional vec2)
  (math-covariance (list vec1) (list vec2) nil 1))

(defun calcFunc-vpcov (vec1 &optional vec2)
  (math-covariance (list vec1) (list vec2) t 1))

(defun calcFunc-vcorr (vec1 &optional vec2)
  (math-covariance (list vec1) (list vec2) nil 2))


(defun math-covariance (vec1 vec2 pop mode)
  (or (car vec2) (= mode 0)
      (progn
	(if (and (eq (car-safe (car vec1)) 'var)
		 (eq (car-safe (calc-var-value (nth 2 (car vec1)))) 'vec))
	    (setq vec1 (symbol-value (nth 2 (car vec1))))
	  (setq vec1 (car vec1)))
	(or (math-matrixp vec1) (math-dimension-error))
	(or (= (length (nth 1 vec1)) 3) (math-dimension-error))
	(setq vec2 (list (math-mat-col vec1 2))
	      vec1 (list (math-mat-col vec1 1)))))
  (math-with-extra-prec 2
    (let* ((split1 (math-split-sdev-vec (math-flatten-many-vecs vec1) nil))
	   (means1 (car split1))
	   (wts1 (nth 1 split1))
	   split2 means2 (wts2 nil)
	   (sqrwts nil)
	   suminvsqrwts
	   (len (1- (length means1))))
      (if (< len (if pop 1 2))
	  (math-reject-arg nil (if pop
				   "*Must be at least 1 argument"
				 "*Must be at least 2 arguments")))
      (if (or wts1 wts2)
	  (setq sqrwts (math-add
			(if wts1
			    (calcFunc-map '(var mul var-mul) wts1 wts1)
			  0)
			(if wts2
			    (calcFunc-map '(var mul var-mul) wts2 wts2)
			  0))
		suminvsqrwts (calcFunc-reduce
			      '(var add var-add)
			      (calcFunc-map '(var div var-div) 1 sqrwts))))
      (or (= mode 0)
	  (progn
	    (setq split2 (math-split-sdev-vec (math-flatten-many-vecs vec2)
					      nil)
		  means2 (car split2)
		  wts2 (nth 2 split1))
	    (or (= len (1- (length means2))) (math-dimension-error))))
      (let* ((diff1 (calcFunc-map
		     '(var add var-add)
		     means1
		     (if sqrwts
			 (math-div (calcFunc-reduce
				    '(var add var-add)
				    (calcFunc-map '(var div var-div)
						  means1 sqrwts))
				   (math-neg suminvsqrwts))
		       (math-div (calcFunc-reducer '(var add var-add) means1)
				 (- len)))))
	     (diff2 (if (= mode 0)
			diff1
		      (calcFunc-map
		       '(var add var-add)
		       means2
		       (if sqrwts
			   (math-div (calcFunc-reduce
				      '(var add var-add)
				      (calcFunc-map '(var div var-div)
						    means2 sqrwts))
				     (math-neg suminvsqrwts))
			 (math-div (calcFunc-reducer '(var add var-add) means2)
				   (- len))))))
	     (covar (calcFunc-map '(var mul var-mul) diff1 diff2)))
	(if sqrwts
	    (setq covar (calcFunc-map '(var div var-div) covar sqrwts)))
	(math-div
	 (calcFunc-reducer '(var add var-add) covar)
	 (if (= mode 2)
	     (let ((var1 (calcFunc-map '(var mul var-mul) diff1 diff1))
		   (var2 (calcFunc-map '(var mul var-mul) diff2 diff2)))
	       (if sqrwts
		   (setq var1 (calcFunc-map '(var div var-div) var1 sqrwts)
			 var2 (calcFunc-map '(var div var-div) var2 sqrwts)))
	       (math-sqrt
		(math-mul (calcFunc-reducer '(var add var-add) var1)
			  (calcFunc-reducer '(var add var-add) var2))))
	   (if sqrwts
	       (if pop
		   suminvsqrwts
		 (math-div (math-mul suminvsqrwts (1- len)) len))
	     (if pop len (1- len)))))))))

(provide 'calc-stat)

;;; arch-tag: 423858e9-8513-489c-9f35-710cd9d9c307
;;; calc-stat.el ends here
Commit	Line	Data
3132f345 CW	1	;;; calc-stat.el --- statistical functions for Calc
3132f345 CW	2
58ba2f8f	3	;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
f0fa15c5	4	;; 2005, 2006, 2007 Free Software Foundation, Inc.
3132f345 CW	5
3132f345 CW	6	;; Author: David Gillespie <daveg@synaptics.com>
e8fff8ed	7	;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
136211a9 EZ	8
	9	;; This file is part of GNU Emacs.
	10
7c671b23 GM	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
136211a9	16	;; GNU Emacs is distributed in the hope that it will be useful,
7c671b23 GM	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.
136211a9	25
3132f345	26	;;; Commentary:
136211a9	27
3132f345	28	;;; Code:
136211a9 EZ	29
136211a9 EZ	30	;; This file is autoloaded from calc-ext.el.
136211a9	31
01a2939c	32	(require 'calc-ext)
136211a9 EZ	33	(require 'calc-macs)
136211a9 EZ	34
136211a9 EZ	35	;;; Statistical operations on vectors.
	36
	37	(defun calc-vector-count (arg)
	38	(interactive "P")
	39	(calc-slow-wrapper
bf77c646	40	(calc-vector-op "coun" 'calcFunc-vcount arg)))
136211a9 EZ	41
	42	(defun calc-vector-sum (arg)
	43	(interactive "P")
	44	(calc-slow-wrapper
	45	(if (calc-is-hyperbolic)
	46	(calc-vector-op "vprd" 'calcFunc-vprod arg)
bf77c646	47	(calc-vector-op "vsum" 'calcFunc-vsum arg))))
136211a9 EZ	48
	49	(defun calc-vector-product (arg)
	50	(interactive "P")
	51	(calc-hyperbolic-func)
bf77c646	52	(calc-vector-sum arg))
136211a9 EZ	53
	54	(defun calc-vector-max (arg)
	55	(interactive "P")
	56	(calc-slow-wrapper
	57	(if (calc-is-inverse)
	58	(calc-vector-op "vmin" 'calcFunc-vmin arg)
bf77c646	59	(calc-vector-op "vmax" 'calcFunc-vmax arg))))
136211a9 EZ	60
	61	(defun calc-vector-min (arg)
	62	(interactive "P")
	63	(calc-invert-func)
bf77c646	64	(calc-vector-max arg))
136211a9 EZ	65
	66	(defun calc-vector-mean (arg)
	67	(interactive "P")
	68	(calc-slow-wrapper
	69	(if (calc-is-hyperbolic)
	70	(if (calc-is-inverse)
	71	(calc-vector-op "harm" 'calcFunc-vhmean arg)
	72	(calc-vector-op "medn" 'calcFunc-vmedian arg))
	73	(if (calc-is-inverse)
	74	(calc-vector-op "meae" 'calcFunc-vmeane arg)
bf77c646	75	(calc-vector-op "mean" 'calcFunc-vmean arg)))))
136211a9 EZ	76
	77	(defun calc-vector-mean-error (arg)
	78	(interactive "P")
	79	(calc-invert-func)
bf77c646	80	(calc-vector-mean arg))
136211a9 EZ	81
	82	(defun calc-vector-median (arg)
	83	(interactive "P")
	84	(calc-hyperbolic-func)
bf77c646	85	(calc-vector-mean arg))
136211a9 EZ	86
	87	(defun calc-vector-harmonic-mean (arg)
	88	(interactive "P")
	89	(calc-invert-func)
	90	(calc-hyperbolic-func)
bf77c646	91	(calc-vector-mean arg))
136211a9 EZ	92
	93	(defun calc-vector-geometric-mean (arg)
	94	(interactive "P")
	95	(calc-slow-wrapper
	96	(if (calc-is-hyperbolic)
	97	(calc-binary-op "geom" 'calcFunc-agmean arg)
bf77c646	98	(calc-vector-op "geom" 'calcFunc-vgmean arg))))
136211a9 EZ	99
	100	(defun calc-vector-sdev (arg)
	101	(interactive "P")
	102	(calc-slow-wrapper
	103	(if (calc-is-hyperbolic)
	104	(if (calc-is-inverse)
	105	(calc-vector-op "pvar" 'calcFunc-vpvar arg)
	106	(calc-vector-op "var" 'calcFunc-vvar arg))
	107	(if (calc-is-inverse)
	108	(calc-vector-op "psdv" 'calcFunc-vpsdev arg)
bf77c646	109	(calc-vector-op "sdev" 'calcFunc-vsdev arg)))))
136211a9 EZ	110
	111	(defun calc-vector-pop-sdev (arg)
	112	(interactive "P")
	113	(calc-invert-func)
bf77c646	114	(calc-vector-sdev arg))
136211a9 EZ	115
	116	(defun calc-vector-variance (arg)
	117	(interactive "P")
	118	(calc-hyperbolic-func)
bf77c646	119	(calc-vector-sdev arg))
136211a9 EZ	120
	121	(defun calc-vector-pop-variance (arg)
	122	(interactive "P")
	123	(calc-invert-func)
	124	(calc-hyperbolic-func)
bf77c646	125	(calc-vector-sdev arg))
136211a9 EZ	126
	127	(defun calc-vector-covariance (arg)
	128	(interactive "P")
	129	(calc-slow-wrapper
	130	(let ((n (if (eq arg 1) 1 2)))
	131	(if (calc-is-hyperbolic)
	132	(calc-enter-result n "corr" (cons 'calcFunc-vcorr
	133	(calc-top-list-n n)))
	134	(if (calc-is-inverse)
	135	(calc-enter-result n "pcov" (cons 'calcFunc-vpcov
	136	(calc-top-list-n n)))
	137	(calc-enter-result n "cov" (cons 'calcFunc-vcov
bf77c646	138	(calc-top-list-n n))))))))
136211a9 EZ	139
	140	(defun calc-vector-pop-covariance (arg)
	141	(interactive "P")
	142	(calc-invert-func)
bf77c646	143	(calc-vector-covariance arg))
136211a9 EZ	144
	145	(defun calc-vector-correlation (arg)
	146	(interactive "P")
	147	(calc-hyperbolic-func)
bf77c646	148	(calc-vector-covariance arg))
136211a9 EZ	149
	150	(defun calc-vector-op (name func arg)
	151	(setq calc-aborted-prefix name
	152	arg (prefix-numeric-value arg))
	153	(if (< arg 0)
	154	(error "Negative arguments not allowed"))
bf77c646	155	(calc-enter-result arg name (cons func (calc-top-list-n arg))))
136211a9 EZ	156
	157
	158
	159
	160	;;; Useful statistical functions
	161
	162	;;; Sum, product, etc., of one or more values or vectors.
	163	;;; Each argument must be either a number or a vector. Vectors
	164	;;; are flattened, but variables inside are assumed to represent
	165	;;; non-vectors.
	166
	167	(defun calcFunc-vsum (&rest vecs)
bf77c646	168	(math-reduce-many-vecs 'calcFunc-add 'calcFunc-vsum vecs 0))
136211a9 EZ	169
136211a9 EZ	170	(defun calcFunc-vprod (&rest vecs)
bf77c646	171	(math-reduce-many-vecs 'calcFunc-mul 'calcFunc-vprod vecs 1))
136211a9 EZ	172
	173	(defun calcFunc-vmax (&rest vecs)
	174	(if (eq (car-safe (car vecs)) 'sdev)
	175	'(var inf var-inf)
	176	(if (eq (car-safe (car vecs)) 'intv)
	177	(nth 3 (math-fix-int-intv (car vecs)))
	178	(math-reduce-many-vecs 'calcFunc-max 'calcFunc-vmax vecs
bf77c646	179	'(neg (var inf var-inf))))))
136211a9 EZ	180
	181	(defun calcFunc-vmin (&rest vecs)
	182	(if (eq (car-safe (car vecs)) 'sdev)
	183	'(neg (var inf var-inf))
	184	(if (eq (car-safe (car vecs)) 'intv)
	185	(nth 2 (math-fix-int-intv (car vecs)))
	186	(math-reduce-many-vecs 'calcFunc-min 'calcFunc-vmin vecs
bf77c646	187	'(var inf var-inf)))))
136211a9 EZ	188
	189	(defun math-reduce-many-vecs (func whole-func vecs ident)
	190	(let ((const-part nil)
	191	(symb-part nil)
	192	val vec)
	193	(let ((calc-internal-prec (+ calc-internal-prec 2)))
	194	(while vecs
	195	(setq val (car vecs))
	196	(and (eq (car-safe val) 'var)
	197	(eq (car-safe (calc-var-value (nth 2 val))) 'vec)
	198	(setq val (symbol-value (nth 2 val))))
	199	(cond ((Math-vectorp val)
	200	(setq vec (append (and const-part (list const-part))
	201	(math-flatten-vector val)))
	202	(setq const-part (if vec
	203	(calcFunc-reducer
	204	(math-calcFunc-to-var func)
	205	(cons 'vec vec))
	206	ident)))
	207	((or (Math-objectp val) (math-infinitep val))
	208	(setq const-part (if const-part
	209	(funcall func const-part val)
	210	val)))
	211	(t
	212	(setq symb-part (nconc symb-part (list val)))))
	213	(setq vecs (cdr vecs))))
	214	(if const-part
	215	(progn
	216	(setq const-part (math-normalize const-part))
	217	(if symb-part
	218	(funcall func const-part (cons whole-func symb-part))
	219	const-part))
bf77c646	220	(if symb-part (cons whole-func symb-part) ident))))
136211a9 EZ	221
	222
	223	;;; Return the number of data elements among the arguments.
	224	(defun calcFunc-vcount (&rest vecs)
	225	(let ((count 0))
	226	(while vecs
	227	(setq count (if (Math-vectorp (car vecs))
	228	(+ count (math-count-elements (car vecs)))
	229	(if (Math-objectp (car vecs))
	230	(1+ count)
	231	(if (and (eq (car-safe (car vecs)) 'var)
	232	(eq (car-safe (calc-var-value
	233	(nth 2 (car vecs))))
	234	'vec))
	235	(+ count (math-count-elements
	236	(symbol-value (nth 2 (car vecs)))))
	237	(math-reject-arg (car vecs) 'numvecp))))
	238	vecs (cdr vecs)))
bf77c646	239	count))
136211a9 EZ	240
	241	(defun math-count-elements (vec)
	242	(let ((count 0))
	243	(while (setq vec (cdr vec))
	244	(setq count (if (Math-vectorp (car vec))
	245	(+ count (math-count-elements (car vec)))
	246	(1+ count))))
bf77c646	247	count))
136211a9 EZ	248
	249
	250	(defun math-flatten-many-vecs (vecs)
	251	(let ((p vecs)
	252	(vec (list 'vec)))
	253	(while p
	254	(setq vec (nconc vec
	255	(if (Math-vectorp (car p))
	256	(math-flatten-vector (car p))
	257	(if (Math-objectp (car p))
	258	(list (car p))
	259	(if (and (eq (car-safe (car p)) 'var)
	260	(eq (car-safe (calc-var-value
	261	(nth 2 (car p)))) 'vec))
	262	(math-flatten-vector (symbol-value
	263	(nth 2 (car p))))
	264	(math-reject-arg (car p) 'numvecp)))))
	265	p (cdr p)))
bf77c646	266	vec))
136211a9 EZ	267
136211a9 EZ	268	(defun calcFunc-vflat (&rest vecs)
bf77c646	269	(math-flatten-many-vecs vecs))
136211a9 EZ	270
	271	(defun math-split-sdev-vec (vec zero-ok)
	272	(let ((means (list 'vec))
	273	(wts (list 'vec))
	274	(exact nil)
	275	(p vec))
	276	(while (and (setq p (cdr p))
	277	(not (and (consp (car p))
	278	(eq (car (car p)) 'sdev)))))
	279	(if (null p)
	280	(list vec nil)
	281	(while (setq vec (cdr vec))
	282	(if (and (consp (setq p (car vec)))
	283	(eq (car p) 'sdev))
	284	(or exact
	285	(setq means (cons (nth 1 p) means)
	286	wts (cons (nth 2 p) wts)))
	287	(if zero-ok
	288	(setq means (cons (nth 1 p) means)
	289	wts (cons 0 wts))
	290	(or exact
	291	(setq means (list 'vec)
	292	wts nil
	293	exact t))
	294	(setq means (cons p means)))))
	295	(list (nreverse means)
bf77c646	296	(and wts (nreverse wts))))))
136211a9 EZ	297
	298
	299	;;; Return the arithmetic mean of the argument numbers or vectors.
	300	;;; (If numbers are error forms, computes the weighted mean.)
	301	(defun calcFunc-vmean (&rest vecs)
	302	(let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
	303	(means (car split))
	304	(wts (nth 1 split))
	305	(len (1- (length means))))
	306	(if (= len 0)
	307	(math-reject-arg nil "*Must be at least 1 argument")
	308	(if (and (= len 1) (eq (car-safe (nth 1 means)) 'intv))
	309	(let ((x (math-fix-int-intv (nth 1 means))))
	310	(calcFunc-vmean (nth 2 x) (nth 3 x)))
	311	(math-with-extra-prec 2
	312	(if (and wts (> len 1))
	313	(let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
	314	(suminvsqrwts (calcFunc-reduce
	315	'(var add var-add)
	316	(calcFunc-map '(var div var-div)
	317	1 sqrwts))))
	318	(math-div (calcFunc-reduce '(var add var-add)
	319	(calcFunc-map '(var div var-div)
	320	means sqrwts))
	321	suminvsqrwts))
bf77c646	322	(math-div (calcFunc-reduce '(var add var-add) means) len)))))))
136211a9 EZ	323
	324	(defun math-fix-int-intv (x)
	325	(if (math-floatp x)
	326	x
	327	(list 'intv 3
	328	(if (memq (nth 1 x) '(2 3)) (nth 2 x) (math-add (nth 2 x) 1))
bf77c646	329	(if (memq (nth 1 x) '(1 3)) (nth 3 x) (math-sub (nth 3 x) 1)))))
136211a9 EZ	330
	331	;;; Compute the mean with an error estimate.
	332	(defun calcFunc-vmeane (&rest vecs)
	333	(let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
	334	(means (car split))
	335	(wts (nth 1 split))
	336	(len (1- (length means))))
	337	(if (= len 0)
	338	(math-reject-arg nil "*Must be at least 1 argument")
	339	(math-with-extra-prec 2
	340	(if wts
	341	(let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
	342	(suminvsqrwts (calcFunc-reduce
	343	'(var add var-add)
	344	(calcFunc-map '(var div var-div)
	345	1 sqrwts))))
	346	(math-make-sdev
	347	(math-div (calcFunc-reduce '(var add var-add)
	348	(calcFunc-map '(var div var-div)
	349	means sqrwts))
	350	suminvsqrwts)
	351	(list 'calcFunc-sqrt (math-div 1 suminvsqrwts))))
	352	(let ((mean (math-div (calcFunc-reduce '(var add var-add) means)
	353	len)))
	354	(math-make-sdev
	355	mean
	356	(list 'calcFunc-sqrt
	357	(math-div (calcFunc-reducer
	358	'(var add var-add)
	359	(calcFunc-map '(var pow var-pow)
	360	(calcFunc-map '(var abs var-abs)
	361	(calcFunc-map
	362	'(var add var-add)
	363	means
	364	(math-neg mean)))
	365	2))
bf77c646	366	(math-mul len (1- len)))))))))))
136211a9 EZ	367
	368
	369	;;; Compute the median of a list of values.
	370	(defun calcFunc-vmedian (&rest vecs)
	371	(let* ((flat (copy-sequence (cdr (math-flatten-many-vecs vecs))))
	372	(p flat)
	373	(len (length flat))
	374	(hlen (/ len 2)))
	375	(if (= len 0)
	376	(math-reject-arg nil "*Must be at least 1 argument")
	377	(if (and (= len 1) (memq (car-safe (car flat)) '(sdev intv)))
	378	(calcFunc-vmean (car flat))
	379	(while p
	380	(if (eq (car-safe (car p)) 'sdev)
	381	(setcar p (nth 1 (car p))))
	382	(or (Math-anglep (car p))
	383	(math-reject-arg (car p) 'anglep))
	384	(setq p (cdr p)))
	385	(setq flat (sort flat 'math-lessp))
	386	(if (= (% len 2) 0)
	387	(math-div (math-add (nth (1- hlen) flat) (nth hlen flat)) 2)
bf77c646	388	(nth hlen flat))))))
136211a9 EZ	389
	390
	391	(defun calcFunc-vgmean (&rest vecs)
	392	(let* ((flat (math-flatten-many-vecs vecs))
	393	(len (1- (length flat))))
	394	(if (= len 0)
	395	(math-reject-arg nil "*Must be at least 1 argument")
	396	(math-with-extra-prec 2
	397	(let ((x (calcFunc-reduce '(var mul math-mul) flat)))
	398	(if (= len 2)
	399	(math-sqrt x)
bf77c646	400	(math-pow x (list 'frac 1 len))))))))
136211a9 EZ	401
	402
	403	(defun calcFunc-agmean (a b)
	404	(cond ((Math-equal a b) a)
	405	((math-zerop a) a)
	406	((math-zerop b) b)
	407	(calc-symbolic-mode (math-inexact-result))
	408	((not (Math-realp a)) (math-reject-arg a 'realp))
	409	((not (Math-realp b)) (math-reject-arg b 'realp))
	410	(t
	411	(math-with-extra-prec 2
	412	(setq a (math-float (math-abs a))
	413	b (math-float (math-abs b)))
	414	(let (mean)
	415	(while (not (math-nearly-equal-float a b))
	416	(setq mean (math-mul-float (math-add-float a b) '(float 5 -1))
	417	b (math-sqrt-float (math-mul-float a b))
	418	a mean))
bf77c646	419	a)))))
136211a9 EZ	420
	421
	422	(defun calcFunc-vhmean (&rest vecs)
	423	(let* ((flat (math-flatten-many-vecs vecs))
	424	(len (1- (length flat))))
	425	(if (= len 0)
	426	(math-reject-arg nil "*Must be at least 1 argument")
	427	(math-with-extra-prec 2
	428	(math-div len
	429	(calcFunc-reduce '(var add math-add)
bf77c646	430	(calcFunc-map '(var inv var-inv) flat)))))))
136211a9 EZ	431
	432
	433
	434	;;; Compute the sample variance or standard deviation of numbers or vectors.
	435	;;; (If the numbers are error forms, only the mean part of them is used.)
	436	(defun calcFunc-vvar (&rest vecs)
	437	(if (and (= (length vecs) 1)
	438	(memq (car-safe (car vecs)) '(sdev intv)))
	439	(if (eq (car-safe (car vecs)) 'intv)
	440	(math-intv-variance (car vecs) nil)
	441	(math-sqr (nth 2 (car vecs))))
bf77c646	442	(math-covariance vecs nil nil 0)))
136211a9 EZ	443
	444	(defun calcFunc-vsdev (&rest vecs)
	445	(if (and (= (length vecs) 1)
	446	(memq (car-safe (car vecs)) '(sdev intv)))
	447	(if (eq (car-safe (car vecs)) 'intv)
	448	(if (math-floatp (car vecs))
	449	(math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
	450	(math-sqrt-12))
	451	(math-sqrt (calcFunc-vvar (car vecs))))
	452	(nth 2 (car vecs)))
bf77c646	453	(math-sqrt (math-covariance vecs nil nil 0))))
136211a9 EZ	454
	455	;;; Compute the population variance or std deviation of numbers or vectors.
	456	(defun calcFunc-vpvar (&rest vecs)
	457	(if (and (= (length vecs) 1)
	458	(memq (car-safe (car vecs)) '(sdev intv)))
	459	(if (eq (car-safe (car vecs)) 'intv)
	460	(math-intv-variance (car vecs) t)
	461	(math-sqr (nth 2 (car vecs))))
bf77c646	462	(math-covariance vecs nil t 0)))
136211a9 EZ	463
	464	(defun calcFunc-vpsdev (&rest vecs)
	465	(if (and (= (length vecs) 1)
	466	(memq (car-safe (car vecs)) '(sdev intv)))
	467	(if (eq (car-safe (car vecs)) 'intv)
	468	(if (math-floatp (car vecs))
	469	(math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
	470	(math-sqrt-12))
	471	(math-sqrt (calcFunc-vpvar (car vecs))))
	472	(nth 2 (car vecs)))
bf77c646	473	(math-sqrt (math-covariance vecs nil t 0))))
136211a9 EZ	474
	475	(defun math-intv-variance (x pop)
	476	(or (math-constp x) (math-reject-arg x 'constp))
	477	(if (math-floatp x)
	478	(math-div (math-sqr (math-sub (nth 3 x) (nth 2 x))) 12)
	479	(let* ((x (math-fix-int-intv x))
	480	(len (math-sub (nth 3 x) (nth 2 x)))
	481	(hlen (math-quotient len 2)))
	482	(math-div (if (math-evenp len)
	483	(calcFunc-sum '(^ (var X var-X) 2) '(var X var-X)
	484	(math-neg hlen) hlen)
	485	(calcFunc-sum '(^ (- (var X var-X) (/ 1 2)) 2)
	486	'(var X var-X)
	487	(math-neg hlen) (math-add hlen 1)))
bf77c646	488	(if pop (math-add len 1) len)))))
136211a9 EZ	489
	490	;;; Compute the covariance and linear correlation coefficient.
	491	(defun calcFunc-vcov (vec1 &optional vec2)
bf77c646	492	(math-covariance (list vec1) (list vec2) nil 1))
136211a9 EZ	493
136211a9 EZ	494	(defun calcFunc-vpcov (vec1 &optional vec2)
bf77c646	495	(math-covariance (list vec1) (list vec2) t 1))
136211a9 EZ	496
136211a9 EZ	497	(defun calcFunc-vcorr (vec1 &optional vec2)
bf77c646	498	(math-covariance (list vec1) (list vec2) nil 2))
136211a9 EZ	499
	500
	501	(defun math-covariance (vec1 vec2 pop mode)
	502	(or (car vec2) (= mode 0)
	503	(progn
	504	(if (and (eq (car-safe (car vec1)) 'var)
	505	(eq (car-safe (calc-var-value (nth 2 (car vec1)))) 'vec))
	506	(setq vec1 (symbol-value (nth 2 (car vec1))))
	507	(setq vec1 (car vec1)))
	508	(or (math-matrixp vec1) (math-dimension-error))
	509	(or (= (length (nth 1 vec1)) 3) (math-dimension-error))
	510	(setq vec2 (list (math-mat-col vec1 2))
	511	vec1 (list (math-mat-col vec1 1)))))
	512	(math-with-extra-prec 2
	513	(let* ((split1 (math-split-sdev-vec (math-flatten-many-vecs vec1) nil))
	514	(means1 (car split1))
	515	(wts1 (nth 1 split1))
	516	split2 means2 (wts2 nil)
	517	(sqrwts nil)
	518	suminvsqrwts
	519	(len (1- (length means1))))
	520	(if (< len (if pop 1 2))
	521	(math-reject-arg nil (if pop
	522	"*Must be at least 1 argument"
	523	"*Must be at least 2 arguments")))
	524	(if (or wts1 wts2)
	525	(setq sqrwts (math-add
	526	(if wts1
	527	(calcFunc-map '(var mul var-mul) wts1 wts1)
	528	0)
	529	(if wts2
	530	(calcFunc-map '(var mul var-mul) wts2 wts2)
	531	0))
	532	suminvsqrwts (calcFunc-reduce
	533	'(var add var-add)
	534	(calcFunc-map '(var div var-div) 1 sqrwts))))
	535	(or (= mode 0)
	536	(progn
	537	(setq split2 (math-split-sdev-vec (math-flatten-many-vecs vec2)
	538	nil)
	539	means2 (car split2)
	540	wts2 (nth 2 split1))
	541	(or (= len (1- (length means2))) (math-dimension-error))))
	542	(let* ((diff1 (calcFunc-map
	543	'(var add var-add)
	544	means1
	545	(if sqrwts
	546	(math-div (calcFunc-reduce
	547	'(var add var-add)
	548	(calcFunc-map '(var div var-div)
	549	means1 sqrwts))
	550	(math-neg suminvsqrwts))
	551	(math-div (calcFunc-reducer '(var add var-add) means1)
	552	(- len)))))
	553	(diff2 (if (= mode 0)
	554	diff1
	555	(calcFunc-map
	556	'(var add var-add)
	557	means2
	558	(if sqrwts
	559	(math-div (calcFunc-reduce
	560	'(var add var-add)
	561	(calcFunc-map '(var div var-div)
	562	means2 sqrwts))
563	(math-neg suminvsqrwts))
564	(math-div (calcFunc-reducer '(var add var-add) means2)
565	(- len))))))
566	(covar (calcFunc-map '(var mul var-mul) diff1 diff2)))
567	(if sqrwts
568	(setq covar (calcFunc-map '(var div var-div) covar sqrwts)))
569	(math-div
570	(calcFunc-reducer '(var add var-add) covar)
571	(if (= mode 2)
572	(let ((var1 (calcFunc-map '(var mul var-mul) diff1 diff1))
573	(var2 (calcFunc-map '(var mul var-mul) diff2 diff2)))
574	(if sqrwts
575	(setq var1 (calcFunc-map '(var div var-div) var1 sqrwts)
576	var2 (calcFunc-map '(var div var-div) var2 sqrwts)))
577	(math-sqrt
578	(math-mul (calcFunc-reducer '(var add var-add) var1)
579	(calcFunc-reducer '(var add var-add) var2))))
580	(if sqrwts
581	(if pop
582	suminvsqrwts
583	(math-div (math-mul suminvsqrwts (1- len)) len))
bf77c646	584	(if pop len (1- len)))))))))
136211a9	585
01a2939c JB	586	(provide 'calc-stat)
01a2939c JB	587
ab5796a9	588	;;; arch-tag: 423858e9-8513-489c-9f35-710cd9d9c307
bf77c646	589	;;; calc-stat.el ends here