root/trunk/lisp/calc/calc-vec.el

;;; calc-vec.el --- vector functions for Calc

;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
;;   2005, 2006, 2007, 2008 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 3, 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)

(defun calc-display-strings (n)
  (interactive "P")
  (calc-wrapper
   (message (if (calc-change-mode 'calc-display-strings n t t)
                "Displaying vectors of integers as quoted strings"
              "Displaying vectors of integers normally"))))


(defun calc-pack (n)
  (interactive "P")
  (calc-wrapper
   (let* ((nn (if n 1 2))
          (mode (if n (prefix-numeric-value n) (calc-top-n 1)))
          (mode (if (and (Math-vectorp mode) (cdr mode)) (cdr mode)
                  (if (integerp mode) mode
                    (error "Packing mode must be an integer or vector of integers"))))
          (num (calc-pack-size mode))
          (items (calc-top-list num nn)))
     (calc-enter-result (+ nn num -1) "pack" (calc-pack-items mode items)))))

(defun calc-pack-size (mode)
  (cond ((consp mode)
         (let ((size 1))
           (while mode
             (or (integerp (car mode)) (error "Vector of integers expected"))
             (setq size (* size (calc-pack-size (car mode)))
                   mode (cdr mode)))
           (if (= size 0)
               (error "Zero dimensions not allowed")
             size)))
        ((>= mode 0) mode)
        (t (or (cdr (assq mode '((-3 . 3) (-13 . 1) (-14 . 3) (-15 . 6))))
               2))))

(defun calc-pack-items (mode items)
  (cond ((consp mode)
         (if (cdr mode)
             (let* ((size (calc-pack-size (cdr mode)))
                    (len (length items))
                    (new nil)
                    p row)
               (while (> len 0)
                 (setq p (nthcdr (1- size) items)
                       row items
                       items (cdr p)
                       len (- len size))
                 (setcdr p nil)
                 (setq new (cons (calc-pack-items (cdr mode) row) new)))
               (calc-pack-items (car mode) (nreverse new)))
           (calc-pack-items (car mode) items)))
        ((>= mode 0)
         (cons 'vec items))
        ((= mode -3)
         (if (and (math-objvecp (car items))
                  (math-objvecp (nth 1 items))
                  (math-objvecp (nth 2 items)))
             (if (and (math-num-integerp (car items))
                      (math-num-integerp (nth 1 items)))
                 (if (math-realp (nth 2 items))
                     (cons 'hms items)
                   (error "Seconds must be real"))
               (error "Hours and minutes must be integers"))
           (math-normalize (list '+
                                 (list '+
                                       (if (eq calc-angle-mode 'rad)
                                           (list '* (car items)
                                                 '(hms 1 0 0))
                                         (car items))
                                       (list '* (nth 1 items) '(hms 0 1 0)))
                                 (list '* (nth 2 items) '(hms 0 0 1))))))
        ((= mode -13)
         (if (math-realp (car items))
             (cons 'date items)
           (if (eq (car-safe (car items)) 'date)
               (car items)
             (if (math-objvecp (car items))
                 (error "Date value must be real")
               (cons 'calcFunc-date items)))))
        ((memq mode '(-14 -15))
         (let ((p items))
           (while (and p (math-objvecp (car p)))
             (or (math-integerp (car p))
                 (error "Components must be integers"))
             (setq p (cdr p)))
           (if p
               (cons 'calcFunc-date items)
             (list 'date (math-dt-to-date items)))))
        ((or (eq (car-safe (car items)) 'vec)
             (eq (car-safe (nth 1 items)) 'vec))
         (let* ((x (car items))
                (vx (eq (car-safe x) 'vec))
                (y (nth 1 items))
                (vy (eq (car-safe y) 'vec))
                (z nil)
                (n (1- (length (if vx x y)))))
           (and vx vy
                (/= n (1- (length y)))
                (error "Vectors must be the same length"))
           (while (>= (setq n (1- n)) 0)
             (setq z (cons (calc-pack-items
                            mode
                            (list (if vx (car (setq x (cdr x))) x)
                                  (if vy (car (setq y (cdr y))) y)))
                           z)))
           (cons 'vec (nreverse z))))
        ((= mode -1)
         (if (and (math-realp (car items)) (math-realp (nth 1 items)))
             (cons 'cplx items)
           (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
               (error "Components must be real"))
           (math-normalize (list '+ (car items)
                                 (list '* (nth 1 items) '(cplx 0 1))))))
        ((= mode -2)
         (if (and (math-realp (car items)) (math-anglep (nth 1 items)))
             (cons 'polar items)
           (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
               (error "Components must be real"))
           (math-normalize (list '* (car items)
                                 (if (math-anglep (nth 1 items))
                                     (list 'polar 1 (nth 1 items))
                                   (list 'calcFunc-exp
                                         (list '*
                                               (math-to-radians-2
                                                (nth 1 items))
                                               (list 'polar
                                                     1
                                                     (math-quarter-circle
                                                      nil)))))))))
        ((= mode -4)
         (let ((x (car items))
               (sigma (nth 1 items)))
           (if (or (math-scalarp x) (not (math-objvecp x)))
               (if (or (math-anglep sigma) (not (math-objvecp sigma)))
                   (math-make-sdev x sigma)
                 (error "Error component must be real"))
             (error "Mean component must be real or complex"))))
        ((= mode -5)
         (let ((a (car items))
               (m (nth 1 items)))
           (if (and (math-anglep a) (math-anglep m))
               (if (math-posp m)
                   (math-make-mod a m)
                 (error "Modulus must be positive"))
             (if (and (math-objectp a) (math-objectp m))
                 (error "Components must be real"))
             (list 'calcFunc-makemod a m))))
        ((memq mode '(-6 -7 -8 -9))
         (let ((lo (car items))
               (hi (nth 1 items)))
           (if (and (or (math-anglep lo) (eq (car lo) 'date)
                        (not (math-objvecp lo)))
                    (or (math-anglep hi) (eq (car hi) 'date)
                        (not (math-objvecp hi))))
               (math-make-intv (+ mode 9) lo hi)
             (error "Components must be real"))))
        ((eq mode -10)
         (if (math-zerop (nth 1 items))
             (error "Denominator must not be zero")
           (if (and (math-integerp (car items)) (math-integerp (nth 1 items)))
               (math-normalize (cons 'frac items))
             (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
                 (error "Components must be integers"))
             (cons 'calcFunc-fdiv items))))
        ((memq mode '(-11 -12))
         (if (and (math-realp (car items)) (math-integerp (nth 1 items)))
             (calcFunc-scf (math-float (car items)) (nth 1 items))
           (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
               (error "Components must be integers"))
           (math-normalize
            (list 'calcFunc-scf
                  (list 'calcFunc-float (car items))
                  (nth 1 items)))))
        (t
         (error "Invalid packing mode: %d" mode))))

(defvar calc-unpack-with-type nil)
(defun calc-unpack (mode)
  (interactive "P")
  (calc-wrapper
   (let ((calc-unpack-with-type t))
     (calc-pop-push-record-list 1 "unpk" (calc-unpack-item
                                          (and mode
                                               (prefix-numeric-value mode))
                                          (calc-top))))))

(defun calc-unpack-type (item)
  (cond ((eq (car-safe item) 'vec)
         (1- (length item)))
        ((eq (car-safe item) 'intv)
         (- (nth 1 item) 9))
        (t
         (or (cdr (assq (car-safe item) '( (cplx . -1) (polar . -2)
                                           (hms . -3) (sdev . -4) (mod . -5)
                                           (frac . -10) (float . -11)
                                           (date . -13) )))
             (error "Argument must be a composite object")))))

(defun calc-unpack-item (mode item)
  (cond ((not mode)
         (if (or (and (not (memq (car-safe item) '(frac float cplx polar vec
                                                        hms date sdev mod
                                                        intv)))
                      (math-objvecp item))
                 (eq (car-safe item) 'var))
             (error "Argument must be a composite object or function call"))
         (if (eq (car item) 'intv)
             (cdr (cdr item))
           (cdr item)))
        ((> mode 0)
         (let ((dims nil)
               type new row)
           (setq item (list item))
           (while (> mode 0)
             (setq type (calc-unpack-type (car item))
                   dims (cons type dims)
                   new (calc-unpack-item nil (car item)))
             (while (setq item (cdr item))
               (or (= (calc-unpack-type (car item)) type)
                   (error "Inconsistent types or dimensions in vector elements"))
               (setq new (append new (calc-unpack-item nil (car item)))))
             (setq item new
                   mode (1- mode)))
           (if (cdr dims) (setq dims (list (cons 'vec (nreverse dims)))))
           (cond ((eq calc-unpack-with-type 'pair)
                  (list (car dims) (cons 'vec item)))
                 (calc-unpack-with-type
                  (append item dims))
                 (t item))))
        ((eq calc-unpack-with-type 'pair)
         (let ((calc-unpack-with-type nil))
           (list mode (cons 'vec (calc-unpack-item mode item)))))
        ((= mode -3)
         (if (eq (car-safe item) 'hms)
             (cdr item)
           (error "Argument must be an HMS form")))
        ((= mode -13)
         (if (eq (car-safe item) 'date)
             (cdr item)
           (error "Argument must be a date form")))
        ((= mode -14)
         (if (eq (car-safe item) 'date)
             (math-date-to-dt (math-floor (nth 1 item)))
           (error "Argument must be a date form")))
        ((= mode -15)
         (if (eq (car-safe item) 'date)
             (append (math-date-to-dt (nth 1 item))
                     (and (not (math-integerp (nth 1 item)))
                          (list 0 0 0)))
           (error "Argument must be a date form")))
        ((eq (car-safe item) 'vec)
         (let ((x nil)
               (y nil)
               res)
           (while (setq item (cdr item))
             (setq res (calc-unpack-item mode (car item))
                   x (cons (car res) x)
                   y (cons (nth 1 res) y)))
           (list (cons 'vec (nreverse x))
                 (cons 'vec (nreverse y)))))
        ((= mode -1)
         (if (eq (car-safe item) 'cplx)
             (cdr item)
           (if (eq (car-safe item) 'polar)
               (cdr (math-complex item))
             (if (Math-realp item)
                 (list item 0)
               (error "Argument must be a complex number")))))
        ((= mode -2)
         (if (or (memq (car-safe item) '(cplx polar))
                 (Math-realp item))
             (cdr (math-polar item))
           (error "Argument must be a complex number")))
        ((= mode -4)
         (if (eq (car-safe item) 'sdev)
             (cdr item)
           (list item 0)))
        ((= mode -5)
         (if (eq (car-safe item) 'mod)
             (cdr item)
           (error "Argument must be a modulo form")))
        ((memq mode '(-6 -7 -8 -9))
         (if (eq (car-safe item) 'intv)
             (cdr (cdr item))
           (list item item)))
        ((= mode -10)
         (if (eq (car-safe item) 'frac)
             (cdr item)
           (if (Math-integerp item)
               (list item 1)
             (error "Argument must be a rational number"))))
        ((= mode -11)
         (if (eq (car-safe item) 'float)
             (list (nth 1 item) (math-normalize (nth 2 item)))
           (error "Expected a floating-point number")))
        ((= mode -12)
         (if (eq (car-safe item) 'float)
             (list (calcFunc-mant item) (calcFunc-xpon item))
           (error "Expected a floating-point number")))
        (t
         (error "Invalid unpacking mode: %d" mode))))

(defun calc-diag (n)
  (interactive "P")
  (calc-wrapper
   (calc-enter-result 1 "diag" (if n
                                   (list 'calcFunc-diag (calc-top-n 1)
                                         (prefix-numeric-value n))
                                 (list 'calcFunc-diag (calc-top-n 1))))))

(defun calc-ident (n)
  (interactive "NDimension of identity matrix = ")
  (calc-wrapper
   (calc-enter-result 0 "idn" (if (eq n 0)
                                  '(calcFunc-idn 1)
                                (list 'calcFunc-idn 1
                                      (prefix-numeric-value n))))))

(defun calc-index (n &optional stack)
  (interactive "NSize of vector = \nP")
  (calc-wrapper
   (if (consp stack)
       (calc-enter-result 3 "indx" (cons 'calcFunc-index (calc-top-list-n 3)))
     (calc-enter-result 0 "indx" (list 'calcFunc-index
                                       (prefix-numeric-value n))))))

(defun calc-build-vector (n)
  (interactive "NSize of vector = ")
  (calc-wrapper
   (calc-enter-result 1 "bldv" (list 'calcFunc-cvec
                                     (calc-top-n 1)
                                     (prefix-numeric-value n)))))

(defun calc-cons (arg)
  (interactive "P")
  (calc-wrapper
   (if (calc-is-hyperbolic)
       (calc-binary-op "rcns" 'calcFunc-rcons arg)
     (calc-binary-op "cons" 'calcFunc-cons arg))))


(defun calc-head (arg)
  (interactive "P")
  (calc-wrapper
   (if (calc-is-inverse)
       (if (calc-is-hyperbolic)
           (calc-unary-op "rtai" 'calcFunc-rtail arg)
         (calc-unary-op "tail" 'calcFunc-tail arg))
     (if (calc-is-hyperbolic)
         (calc-unary-op "rhed" 'calcFunc-rhead arg)
       (calc-unary-op "head" 'calcFunc-head arg)))))

(defun calc-tail (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-head arg))

(defun calc-vlength (arg)
  (interactive "P")
  (calc-wrapper
   (if (calc-is-hyperbolic)
       (calc-unary-op "dims" 'calcFunc-mdims arg)
     (calc-unary-op "len" 'calcFunc-vlen arg))))

(defun calc-arrange-vector (n)
  (interactive "NNumber of columns = ")
  (calc-wrapper
   (calc-enter-result 1 "arng" (list 'calcFunc-arrange (calc-top-n 1)
                                     (prefix-numeric-value n)))))

(defun calc-vector-find (arg)
  (interactive "P")
  (calc-wrapper
   (let ((func (cons 'calcFunc-find (calc-top-list-n 2))))
     (calc-enter-result
      2 "find"
      (if arg (append func (list (prefix-numeric-value arg))) func)))))

(defun calc-subvector ()
  (interactive)
  (calc-wrapper
   (if (calc-is-inverse)
       (calc-enter-result 3 "rsvc" (cons 'calcFunc-rsubvec
                                         (calc-top-list-n 3)))
     (calc-enter-result 3 "svec" (cons 'calcFunc-subvec (calc-top-list-n 3))))))

(defun calc-reverse-vector (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "rev" 'calcFunc-rev arg)))

(defun calc-mask-vector (arg)
  (interactive "P")
  (calc-wrapper
   (calc-binary-op "vmsk" 'calcFunc-vmask arg)))

(defun calc-expand-vector (arg)
  (interactive "P")
  (calc-wrapper
   (if (calc-is-hyperbolic)
       (calc-enter-result 3 "vexp" (cons 'calcFunc-vexp (calc-top-list-n 3)))
     (calc-binary-op "vexp" 'calcFunc-vexp arg))))

(defun calc-sort ()
  (interactive)
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-enter-result 1 "rsrt" (list 'calcFunc-rsort (calc-top-n 1)))
     (calc-enter-result 1 "sort" (list 'calcFunc-sort (calc-top-n 1))))))

(defun calc-grade ()
  (interactive)
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-enter-result 1 "rgrd" (list 'calcFunc-rgrade (calc-top-n 1)))
     (calc-enter-result 1 "grad" (list 'calcFunc-grade (calc-top-n 1))))))

(defun calc-histogram (n)
  (interactive "NNumber of bins: ")
  (calc-slow-wrapper
   (if calc-hyperbolic-flag
       (calc-enter-result 2 "hist" (list 'calcFunc-histogram
                                         (calc-top-n 2)
                                         (calc-top-n 1)
                                         (prefix-numeric-value n)))
     (calc-enter-result 1 "hist" (list 'calcFunc-histogram
                                       (calc-top-n 1)
                                       (prefix-numeric-value n))))))

(defun calc-transpose (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "trn" 'calcFunc-trn arg)))

(defun calc-conj-transpose (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "ctrn" 'calcFunc-ctrn arg)))

(defun calc-cross (arg)
  (interactive "P")
  (calc-wrapper
   (calc-binary-op "cros" 'calcFunc-cross arg)))

(defun calc-remove-duplicates (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "rdup" 'calcFunc-rdup arg)))

(defun calc-set-union (arg)
  (interactive "P")
  (calc-wrapper
   (calc-binary-op "unio" 'calcFunc-vunion arg '(vec) 'calcFunc-rdup)))

(defun calc-set-intersect (arg)
  (interactive "P")
  (calc-wrapper
   (calc-binary-op "intr" 'calcFunc-vint arg '(vec) 'calcFunc-rdup)))

(defun calc-set-difference (arg)
  (interactive "P")
  (calc-wrapper
   (calc-binary-op "diff" 'calcFunc-vdiff arg '(vec) 'calcFunc-rdup)))

(defun calc-set-xor (arg)
  (interactive "P")
  (calc-wrapper
   (calc-binary-op "xor" 'calcFunc-vxor arg '(vec) 'calcFunc-rdup)))

(defun calc-set-complement (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "cmpl" 'calcFunc-vcompl arg)))

(defun calc-set-floor (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "vflr" 'calcFunc-vfloor arg)))

(defun calc-set-enumerate (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "enum" 'calcFunc-venum arg)))

(defun calc-set-span (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "span" 'calcFunc-vspan arg)))

(defun calc-set-cardinality (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "card" 'calcFunc-vcard arg)))

(defun calc-unpack-bits (arg)
  (interactive "P")
  (calc-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "bpck" 'calcFunc-vpack arg)
     (calc-unary-op "bupk" 'calcFunc-vunpack arg))))

(defun calc-pack-bits (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-unpack-bits arg))


(defun calc-rnorm (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "rnrm" 'calcFunc-rnorm arg)))

(defun calc-cnorm (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "cnrm" 'calcFunc-cnorm arg)))

(defun calc-mrow (n &optional nn)
  (interactive "NRow number: \nP")
  (calc-wrapper
   (if (consp nn)
       (calc-enter-result 2 "mrow" (cons 'calcFunc-mrow (calc-top-list-n 2)))
     (setq n (prefix-numeric-value n))
     (if (= n 0)
         (calc-enter-result 1 "getd" (list 'calcFunc-getdiag (calc-top-n 1)))
       (if (< n 0)
           (calc-enter-result 1 "rrow" (list 'calcFunc-mrrow
                                             (calc-top-n 1) (- n)))
         (calc-enter-result 1 "mrow" (list 'calcFunc-mrow
                                           (calc-top-n 1) n)))))))

(defun calc-mcol (n &optional nn)
  (interactive "NColumn number: \nP")
  (calc-wrapper
   (if (consp nn)
       (calc-enter-result 2 "mcol" (cons 'calcFunc-mcol (calc-top-list-n 2)))
     (setq n (prefix-numeric-value n))
     (if (= n 0)
         (calc-enter-result 1 "getd" (list 'calcFunc-getdiag (calc-top-n 1)))
       (if (< n 0)
           (calc-enter-result 1 "rcol" (list 'calcFunc-mrcol
                                             (calc-top-n 1) (- n)))
         (calc-enter-result 1 "mcol" (list 'calcFunc-mcol
                                           (calc-top-n 1) n)))))))


;;;; Vectors.

(defun calcFunc-mdims (m)
  (or (math-vectorp m)
      (math-reject-arg m 'vectorp))
  (cons 'vec (math-mat-dimens m)))


;;; Apply a function elementwise to vector A.  [V X V; N X N] [Public]
(defun math-map-vec (f a)
  (if (math-vectorp a)
      (cons 'vec (mapcar f (cdr a)))
    (funcall f a)))

(defun math-dimension-error ()
  (calc-record-why "*Dimension error")
  (signal 'wrong-type-argument nil))


;;; Build a vector out of a list of objects.  [Public]
(defun calcFunc-vec (&rest objs)
  (cons 'vec objs))


;;; Build a constant vector or matrix.  [Public]
(defun calcFunc-cvec (obj &rest dims)
  (math-make-vec-dimen obj dims))

(defun math-make-vec-dimen (obj dims)
  (if dims
      (if (natnump (car dims))
          (if (or (cdr dims)
                  (not (math-numberp obj)))
              (cons 'vec (copy-sequence
                          (make-list (car dims)
                                     (math-make-vec-dimen obj (cdr dims)))))
            (cons 'vec (make-list (car dims) obj)))
        (math-reject-arg (car dims) 'fixnatnump))
    obj))

(defun calcFunc-head (vec)
  (if (and (Math-vectorp vec)
           (cdr vec))
      (nth 1 vec)
    (calc-record-why 'vectorp vec)
    (list 'calcFunc-head vec)))

(defun calcFunc-tail (vec)
  (if (and (Math-vectorp vec)
           (cdr vec))
      (cons 'vec (cdr (cdr vec)))
    (calc-record-why 'vectorp vec)
    (list 'calcFunc-tail vec)))

(defun calcFunc-cons (head tail)
  (if (Math-vectorp tail)
      (cons 'vec (cons head (cdr tail)))
    (calc-record-why 'vectorp tail)
    (list 'calcFunc-cons head tail)))

(defun calcFunc-rhead (vec)
  (if (and (Math-vectorp vec)
           (cdr vec))
      (let ((vec (copy-sequence vec)))
        (setcdr (nthcdr (- (length vec) 2) vec) nil)
        vec)
    (calc-record-why 'vectorp vec)
    (list 'calcFunc-rhead vec)))

(defun calcFunc-rtail (vec)
  (if (and (Math-vectorp vec)
           (cdr vec))
      (nth (1- (length vec)) vec)
    (calc-record-why 'vectorp vec)
    (list 'calcFunc-rtail vec)))

(defun calcFunc-rcons (head tail)
  (if (Math-vectorp head)
      (append head (list tail))
    (calc-record-why 'vectorp head)
    (list 'calcFunc-rcons head tail)))



;;; Apply a function elementwise to vectors A and B.  [O X O O] [Public]
(defun math-map-vec-2 (f a b)
  (if (math-vectorp a)
      (if (math-vectorp b)
          (let ((v nil))
            (while (setq a (cdr a))
              (or (setq b (cdr b))
                  (math-dimension-error))
              (setq v (cons (funcall f (car a) (car b)) v)))
            (if a (math-dimension-error))
            (cons 'vec (nreverse v)))
        (let ((v nil))
          (while (setq a (cdr a))
            (setq v (cons (funcall f (car a) b) v)))
          (cons 'vec (nreverse v))))
    (if (math-vectorp b)
        (let ((v nil))
          (while (setq b (cdr b))
            (setq v (cons (funcall f a (car b)) v)))
          (cons 'vec (nreverse v)))
      (funcall f a b))))



;;; "Reduce" a function over a vector (left-associatively).  [O X V] [Public]
(defun math-reduce-vec (f a)
  (if (math-vectorp a)
      (if (cdr a)
          (let ((accum (car (setq a (cdr a)))))
            (while (setq a (cdr a))
              (setq accum (funcall f accum (car a))))
            accum)
        0)
    a))

;;; Reduce a function over the columns of matrix A.  [V X V] [Public]
(defun math-reduce-cols (f a)
  (if (math-matrixp a)
      (cons 'vec (math-reduce-cols-col-step f (cdr a) 1 (length (nth 1 a))))
    a))

(defun math-reduce-cols-col-step (f a col cols)
  (and (< col cols)
       (cons (math-reduce-cols-row-step f (nth col (car a)) col (cdr a))
             (math-reduce-cols-col-step f a (1+ col) cols))))

(defun math-reduce-cols-row-step (f tot col a)
  (if a
      (math-reduce-cols-row-step f
                                 (funcall f tot (nth col (car a)))
                                 col
                                 (cdr a))
    tot))



(defun math-dot-product (a b)
  (if (setq a (cdr a) b (cdr b))
      (let ((accum (math-mul (car a) (car b))))
        (while (setq a (cdr a) b (cdr b))
          (setq accum (math-add accum (math-mul (car a) (car b)))))
        accum)
    0))


;;; Return the number of elements in vector V.  [Public]
(defun calcFunc-vlen (v)
  (if (math-vectorp v)
      (1- (length v))
    (if (math-objectp v)
        0
      (list 'calcFunc-vlen v))))

;;; Get the Nth row of a matrix.
(defun calcFunc-mrow (mat n)   ; [Public]
  (if (Math-vectorp n)
      (math-map-vec (function (lambda (x) (calcFunc-mrow mat x))) n)
    (if (and (eq (car-safe n) 'intv) (math-constp n))
        (calcFunc-subvec mat
                         (math-add (nth 2 n) (if (memq (nth 1 n) '(2 3)) 0 1))
                         (math-add (nth 3 n) (if (memq (nth 1 n) '(1 3)) 1 0)))
      (or (and (integerp (setq n (math-check-integer n)))
               (> n 0))
          (math-reject-arg n 'fixposintp))
      (or (Math-vectorp mat)
          (math-reject-arg mat 'vectorp))
      (or (nth n mat)
          (math-reject-arg n "*Index out of range")))))

(defun calcFunc-subscr (mat n &optional m)
  (setq mat (calcFunc-mrow mat n))
  (if m
      (if (math-num-integerp n)
          (calcFunc-mrow mat m)
        (calcFunc-mcol mat m))
    mat))

;;; Get the Nth column of a matrix.
(defun math-mat-col (mat n)
  (cons 'vec (mapcar (function (lambda (x) (elt x n))) (cdr mat))))

(defun calcFunc-mcol (mat n)   ; [Public]
  (if (Math-vectorp n)
      (calcFunc-trn
       (math-map-vec (function (lambda (x) (calcFunc-mcol mat x))) n))
    (if (and (eq (car-safe n) 'intv) (math-constp n))
        (if (math-matrixp mat)
            (math-map-vec (function (lambda (x) (calcFunc-mrow x n))) mat)
          (calcFunc-mrow mat n))
      (or (and (integerp (setq n (math-check-integer n)))
               (> n 0))
          (math-reject-arg n 'fixposintp))
      (or (Math-vectorp mat)
          (math-reject-arg mat 'vectorp))
      (or (if (math-matrixp mat)
              (and (< n (length (nth 1 mat)))
                   (math-mat-col mat n))
            (nth n mat))
          (math-reject-arg n "*Index out of range")))))

;;; Remove the Nth row from a matrix.
(defun math-mat-less-row (mat n)
  (if (<= n 0)
      (cdr mat)
    (cons (car mat)
          (math-mat-less-row (cdr mat) (1- n)))))

(defun calcFunc-mrrow (mat n)   ; [Public]
  (and (integerp (setq n (math-check-integer n)))
       (> n 0)
       (< n (length mat))
       (math-mat-less-row mat n)))

;;; Remove the Nth column from a matrix.
(defun math-mat-less-col (mat n)
  (cons 'vec (mapcar (function (lambda (x) (math-mat-less-row x n)))
                     (cdr mat))))

(defun calcFunc-mrcol (mat n)   ; [Public]
  (and (integerp (setq n (math-check-integer n)))
       (> n 0)
       (if (math-matrixp mat)
           (and (< n (length (nth 1 mat)))
                (math-mat-less-col mat n))
         (math-mat-less-row mat n))))

(defun calcFunc-getdiag (mat)   ; [Public]
  (if (math-square-matrixp mat)
      (cons 'vec (math-get-diag-step (cdr mat) 1))
    (calc-record-why 'square-matrixp mat)
    (list 'calcFunc-getdiag mat)))

(defun math-get-diag-step (row n)
  (and row
       (cons (nth n (car row))
             (math-get-diag-step (cdr row) (1+ n)))))

(defun math-transpose (mat)   ; [Public]
  (let ((m nil)
        (col (length (nth 1 mat))))
    (while (> (setq col (1- col)) 0)
      (setq m (cons (math-mat-col mat col) m)))
    (cons 'vec m)))

(defun calcFunc-trn (mat)
  (if (math-vectorp mat)
      (if (math-matrixp mat)
          (math-transpose mat)
        (math-col-matrix mat))
    (if (math-numberp mat)
        mat
      (math-reject-arg mat 'matrixp))))

(defun calcFunc-ctrn (mat)
  (calcFunc-conj (calcFunc-trn mat)))

(defun calcFunc-pack (mode els)
  (or (Math-vectorp els) (math-reject-arg els 'vectorp))
  (if (and (Math-vectorp mode) (cdr mode))
      (setq mode (cdr mode))
    (or (integerp mode) (math-reject-arg mode 'fixnump)))
  (condition-case err
      (if (= (calc-pack-size mode) (1- (length els)))
          (calc-pack-items mode (cdr els))
        (math-reject-arg els "*Wrong number of elements"))
    (error (math-reject-arg els (nth 1 err)))))

(defun calcFunc-unpack (mode thing)
  (or (integerp mode) (math-reject-arg mode 'fixnump))
  (condition-case err
      (cons 'vec (calc-unpack-item mode thing))
    (error (math-reject-arg thing (nth 1 err)))))

(defun calcFunc-unp