lisp/calc/calc-vec.el

   1 ;;; calc-vec.el --- vector functions for Calc
   2
   3 ;; Copyright (C) 1990-1993, 2001-2011 Free Software Foundation, Inc.
   4
   5 ;; Author: David Gillespie <daveg@synaptics.com>
   6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
   7
   8 ;; This file is part of GNU Emacs.
   9
  10 ;; GNU Emacs is free software: you can redistribute it and/or modify
  11 ;; it under the terms of the GNU General Public License as published by
  12 ;; the Free Software Foundation, either version 3 of the License, or
  13 ;; (at your option) any later version.
  14
  15 ;; GNU Emacs is distributed in the hope that it will be useful,
  16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 ;; GNU General Public License for more details.
  19
  20 ;; You should have received a copy of the GNU General Public License
  21 ;; along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.
  22
  23 ;;; Commentary:
  24
  25 ;;; Code:
  26
  27 ;; This file is autoloaded from calc-ext.el.
  28
  29 (require 'calc-ext)
  30 (require 'calc-macs)
  31
  32 ;; Declare functions which are defined elsewhere.
  33 (declare-function math-read-expr-level "calc-aent" (exp-prec &optional exp-term))
  34
  35
  36 (defun calc-display-strings (n)
  37   (interactive "P")
  38   (calc-wrapper
  39    (message (if (calc-change-mode 'calc-display-strings n t t)
  40                 "Displaying vectors of integers as quoted strings"
  41               "Displaying vectors of integers normally"))))
  42
  43
  44 (defun calc-pack (n)
  45   (interactive "P")
  46   (calc-wrapper
  47    (let* ((nn (if n 1 2))
  48           (mode (if n (prefix-numeric-value n) (calc-top-n 1)))
  49           (mode (if (and (Math-vectorp mode) (cdr mode)) (cdr mode)
  50                   (if (integerp mode) mode
  51                     (error "Packing mode must be an integer or vector of integers"))))
  52           (num (calc-pack-size mode))
  53           (items (calc-top-list num nn)))
  54      (calc-enter-result (+ nn num -1) "pack" (calc-pack-items mode items)))))
  55
  56 (defun calc-pack-size (mode)
  57   (cond ((consp mode)
  58          (let ((size 1))
  59            (while mode
  60              (or (integerp (car mode)) (error "Vector of integers expected"))
  61              (setq size (* size (calc-pack-size (car mode)))
  62                    mode (cdr mode)))
  63            (if (= size 0)
  64                (error "Zero dimensions not allowed")
  65              size)))
  66         ((>= mode 0) mode)
  67         (t (or (cdr (assq mode '((-3 . 3) (-13 . 1) (-14 . 3) (-15 . 6))))
  68                2))))
  69
  70 (defun calc-pack-items (mode items)
  71   (cond ((consp mode)
  72          (if (cdr mode)
  73              (let* ((size (calc-pack-size (cdr mode)))
  74                     (len (length items))
  75                     (new nil)
  76                     p row)
  77                (while (> len 0)
  78                  (setq p (nthcdr (1- size) items)
  79                        row items
  80                        items (cdr p)
  81                        len (- len size))
  82                  (setcdr p nil)
  83                  (setq new (cons (calc-pack-items (cdr mode) row) new)))
  84                (calc-pack-items (car mode) (nreverse new)))
  85            (calc-pack-items (car mode) items)))
  86         ((>= mode 0)
  87          (cons 'vec items))
  88         ((= mode -3)
  89          (if (and (math-objvecp (car items))
  90                   (math-objvecp (nth 1 items))
  91                   (math-objvecp (nth 2 items)))
  92              (if (and (math-num-integerp (car items))
  93                       (math-num-integerp (nth 1 items)))
  94                  (if (math-realp (nth 2 items))
  95                      (cons 'hms items)
  96                    (error "Seconds must be real"))
  97                (error "Hours and minutes must be integers"))
  98            (math-normalize (list '+
  99                                  (list '+
 100                                        (if (eq calc-angle-mode 'rad)
 101                                            (list '* (car items)
 102                                                  '(hms 1 0 0))
 103                                          (car items))
 104                                        (list '* (nth 1 items) '(hms 0 1 0)))
 105                                  (list '* (nth 2 items) '(hms 0 0 1))))))
 106         ((= mode -13)
 107          (if (math-realp (car items))
 108              (cons 'date items)
 109            (if (eq (car-safe (car items)) 'date)
 110                (car items)
 111              (if (math-objvecp (car items))
 112                  (error "Date value must be real")
 113                (cons 'calcFunc-date items)))))
 114         ((memq mode '(-14 -15))
 115          (let ((p items))
 116            (while (and p (math-objvecp (car p)))
 117              (or (math-integerp (car p))
 118                  (error "Components must be integers"))
 119              (setq p (cdr p)))
 120            (if p
 121                (cons 'calcFunc-date items)
 122              (list 'date (math-dt-to-date items)))))
 123         ((or (eq (car-safe (car items)) 'vec)
 124              (eq (car-safe (nth 1 items)) 'vec))
 125          (let* ((x (car items))
 126                 (vx (eq (car-safe x) 'vec))
 127                 (y (nth 1 items))
 128                 (vy (eq (car-safe y) 'vec))
 129                 (z nil)
 130                 (n (1- (length (if vx x y)))))
 131            (and vx vy
 132                 (/= n (1- (length y)))
 133                 (error "Vectors must be the same length"))
 134            (while (>= (setq n (1- n)) 0)
 135              (setq z (cons (calc-pack-items
 136                             mode
 137                             (list (if vx (car (setq x (cdr x))) x)
 138                                   (if vy (car (setq y (cdr y))) y)))
 139                            z)))
 140            (cons 'vec (nreverse z))))
 141         ((= mode -1)
 142          (if (and (math-realp (car items)) (math-realp (nth 1 items)))
 143              (cons 'cplx items)
 144            (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
 145                (error "Components must be real"))
 146            (math-normalize (list '+ (car items)
 147                                  (list '* (nth 1 items) '(cplx 0 1))))))
 148         ((= mode -2)
 149          (if (and (math-realp (car items)) (math-anglep (nth 1 items)))
 150              (cons 'polar items)
 151            (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
 152                (error "Components must be real"))
 153            (math-normalize (list '* (car items)
 154                                  (if (math-anglep (nth 1 items))
 155                                      (list 'polar 1 (nth 1 items))
 156                                    (list 'calcFunc-exp
 157                                          (list '*
 158                                                (math-to-radians-2
 159                                                 (nth 1 items))
 160                                                (list 'polar
 161                                                      1
 162                                                      (math-quarter-circle
 163                                                       nil)))))))))
 164         ((= mode -4)
 165          (let ((x (car items))
 166                (sigma (nth 1 items)))
 167            (if (or (math-scalarp x) (not (math-objvecp x)))
 168                (if (or (math-anglep sigma) (not (math-objvecp sigma)))
 169                    (math-make-sdev x sigma)
 170                  (error "Error component must be real"))
 171              (error "Mean component must be real or complex"))))
 172         ((= mode -5)
 173          (let ((a (car items))
 174                (m (nth 1 items)))
 175            (if (and (math-anglep a) (math-anglep m))
 176                (if (math-posp m)
 177                    (math-make-mod a m)
 178                  (error "Modulus must be positive"))
 179              (if (and (math-objectp a) (math-objectp m))
 180                  (error "Components must be real"))
 181              (list 'calcFunc-makemod a m))))
 182         ((memq mode '(-6 -7 -8 -9))
 183          (let ((lo (car items))
 184                (hi (nth 1 items)))
 185            (if (and (or (math-anglep lo) (eq (car lo) 'date)
 186                         (not (math-objvecp lo)))
 187                     (or (math-anglep hi) (eq (car hi) 'date)
 188                         (not (math-objvecp hi))))
 189                (math-make-intv (+ mode 9) lo hi)
 190              (error "Components must be real"))))
 191         ((eq mode -10)
 192          (if (math-zerop (nth 1 items))
 193              (error "Denominator must not be zero")
 194            (if (and (math-integerp (car items)) (math-integerp (nth 1 items)))
 195                (math-normalize (cons 'frac items))
 196              (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
 197                  (error "Components must be integers"))
 198              (cons 'calcFunc-fdiv items))))
 199         ((memq mode '(-11 -12))
 200          (if (and (math-realp (car items)) (math-integerp (nth 1 items)))
 201              (calcFunc-scf (math-float (car items)) (nth 1 items))
 202            (if (and (math-objectp (car items)) (math-objectp (nth 1 items)))
 203                (error "Components must be integers"))
 204            (math-normalize
 205             (list 'calcFunc-scf
 206                   (list 'calcFunc-float (car items))
 207                   (nth 1 items)))))
 208         (t
 209          (error "Invalid packing mode: %d" mode))))
 210
 211 (defvar calc-unpack-with-type nil)
 212 (defun calc-unpack (mode)
 213   (interactive "P")
 214   (calc-wrapper
 215    (let ((calc-unpack-with-type t))
 216      (calc-pop-push-record-list 1 "unpk" (calc-unpack-item
 217                                           (and mode
 218                                                (prefix-numeric-value mode))
 219                                           (calc-top))))))
 220
 221 (defun calc-unpack-type (item)
 222   (cond ((eq (car-safe item) 'vec)
 223          (1- (length item)))
 224         ((eq (car-safe item) 'intv)
 225          (- (nth 1 item) 9))
 226         (t
 227          (or (cdr (assq (car-safe item) '( (cplx . -1) (polar . -2)
 228                                            (hms . -3) (sdev . -4) (mod . -5)
 229                                            (frac . -10) (float . -11)
 230                                            (date . -13) )))
 231              (error "Argument must be a composite object")))))
 232
 233 (defun calc-unpack-item (mode item)
 234   (cond ((not mode)
 235          (if (or (and (not (memq (car-safe item) '(frac float cplx polar vec
 236                                                         hms date sdev mod
 237                                                         intv)))
 238                       (math-objvecp item))
 239                  (eq (car-safe item) 'var))
 240              (error "Argument must be a composite object or function call"))
 241          (if (eq (car item) 'intv)
 242              (cdr (cdr item))
 243            (cdr item)))
 244         ((> mode 0)
 245          (let ((dims nil)
 246                type new row)
 247            (setq item (list item))
 248            (while (> mode 0)
 249              (setq type (calc-unpack-type (car item))
 250                    dims (cons type dims)
 251                    new (calc-unpack-item nil (car item)))
 252              (while (setq item (cdr item))
 253                (or (= (calc-unpack-type (car item)) type)
 254                    (error "Inconsistent types or dimensions in vector elements"))
 255                (setq new (append new (calc-unpack-item nil (car item)))))
 256              (setq item new
 257                    mode (1- mode)))
 258            (if (cdr dims) (setq dims (list (cons 'vec (nreverse dims)))))
 259            (cond ((eq calc-unpack-with-type 'pair)
 260                   (list (car dims) (cons 'vec item)))
 261                  (calc-unpack-with-type
 262                   (append item dims))
 263                  (t item))))
 264         ((eq calc-unpack-with-type 'pair)
 265          (let ((calc-unpack-with-type nil))
 266            (list mode (cons 'vec (calc-unpack-item mode item)))))
 267         ((= mode -3)
 268          (if (eq (car-safe item) 'hms)
 269              (cdr item)
 270            (error "Argument must be an HMS form")))
 271         ((= mode -13)
 272          (if (eq (car-safe item) 'date)
 273              (cdr item)
 274            (error "Argument must be a date form")))
 275         ((= mode -14)
 276          (if (eq (car-safe item) 'date)
 277              (math-date-to-dt (math-floor (nth 1 item)))
 278            (error "Argument must be a date form")))
 279         ((= mode -15)
 280          (if (eq (car-safe item) 'date)
 281              (append (math-date-to-dt (nth 1 item))
 282                      (and (not (math-integerp (nth 1 item)))
 283                           (list 0 0 0)))
 284            (error "Argument must be a date form")))
 285         ((eq (car-safe item) 'vec)
 286          (let ((x nil)
 287                (y nil)
 288                res)
 289            (while (setq item (cdr item))
 290              (setq res (calc-unpack-item mode (car item))
 291                    x (cons (car res) x)
 292                    y (cons (nth 1 res) y)))
 293            (list (cons 'vec (nreverse x))
 294                  (cons 'vec (nreverse y)))))
 295         ((= mode -1)
 296          (if (eq (car-safe item) 'cplx)
 297              (cdr item)
 298            (if (eq (car-safe item) 'polar)
 299                (cdr (math-complex item))
 300              (if (Math-realp item)
 301                  (list item 0)
 302                (error "Argument must be a complex number")))))
 303         ((= mode -2)
 304          (if (or (memq (car-safe item) '(cplx polar))
 305                  (Math-realp item))
 306              (cdr (math-polar item))
 307            (error "Argument must be a complex number")))
 308         ((= mode -4)
 309          (if (eq (car-safe item) 'sdev)
 310              (cdr item)
 311            (list item 0)))
 312         ((= mode -5)
 313          (if (eq (car-safe item) 'mod)
 314              (cdr item)
 315            (error "Argument must be a modulo form")))
 316         ((memq mode '(-6 -7 -8 -9))
 317          (if (eq (car-safe item) 'intv)
 318              (cdr (cdr item))
 319            (list item item)))
 320         ((= mode -10)
 321          (if (eq (car-safe item) 'frac)
 322              (cdr item)
 323            (if (Math-integerp item)
 324                (list item 1)
 325              (error "Argument must be a rational number"))))
 326         ((= mode -11)
 327          (if (eq (car-safe item) 'float)
 328              (list (nth 1 item) (math-normalize (nth 2 item)))
 329            (error "Expected a floating-point number")))
 330         ((= mode -12)
 331          (if (eq (car-safe item) 'float)
 332              (list (calcFunc-mant item) (calcFunc-xpon item))
 333            (error "Expected a floating-point number")))
 334         (t
 335          (error "Invalid unpacking mode: %d" mode))))
 336
 337 (defun calc-diag (n)
 338   (interactive "P")
 339   (calc-wrapper
 340    (calc-enter-result 1 "diag" (if n
 341                                    (list 'calcFunc-diag (calc-top-n 1)
 342                                          (prefix-numeric-value n))
 343                                  (list 'calcFunc-diag (calc-top-n 1))))))
 344
 345 (defun calc-ident (n)
 346   (interactive "NDimension of identity matrix = ")
 347   (calc-wrapper
 348    (calc-enter-result 0 "idn" (if (eq n 0)
 349                                   '(calcFunc-idn 1)
 350                                 (list 'calcFunc-idn 1
 351                                       (prefix-numeric-value n))))))
 352
 353 (defun calc-index (n &optional stack)
 354   (interactive "NSize of vector = \nP")
 355   (calc-wrapper
 356    (if (consp stack)
 357        (calc-enter-result 3 "indx" (cons 'calcFunc-index (calc-top-list-n 3)))
 358      (calc-enter-result 0 "indx" (list 'calcFunc-index
 359                                        (prefix-numeric-value n))))))
 360
 361 (defun calc-build-vector (n)
 362   (interactive "NSize of vector = ")
 363   (calc-wrapper
 364    (calc-enter-result 1 "bldv" (list 'calcFunc-cvec
 365                                      (calc-top-n 1)
 366                                      (prefix-numeric-value n)))))
 367
 368 (defun calc-cons (arg)
 369   (interactive "P")
 370   (calc-wrapper
 371    (if (calc-is-hyperbolic)
 372        (calc-binary-op "rcns" 'calcFunc-rcons arg)
 373      (calc-binary-op "cons" 'calcFunc-cons arg))))
 374
 375
 376 (defun calc-head (arg)
 377   (interactive "P")
 378   (calc-wrapper
 379    (if (calc-is-inverse)
 380        (if (calc-is-hyperbolic)
 381            (calc-unary-op "rtai" 'calcFunc-rtail arg)
 382          (calc-unary-op "tail" 'calcFunc-tail arg))
 383      (if (calc-is-hyperbolic)
 384          (calc-unary-op "rhed" 'calcFunc-rhead arg)
 385        (calc-unary-op "head" 'calcFunc-head arg)))))
 386
 387 (defun calc-tail (arg)
 388   (interactive "P")
 389   (calc-invert-func)
 390   (calc-head arg))
 391
 392 (defun calc-vlength (arg)
 393   (interactive "P")
 394   (calc-wrapper
 395    (if (calc-is-hyperbolic)
 396        (calc-unary-op "dims" 'calcFunc-mdims arg)
 397      (calc-unary-op "len" 'calcFunc-vlen arg))))
 398
 399 (defun calc-arrange-vector (n)
 400   (interactive "NNumber of columns = ")
 401   (calc-wrapper
 402    (calc-enter-result 1 "arng" (list 'calcFunc-arrange (calc-top-n 1)
 403                                      (prefix-numeric-value n)))))
 404
 405 (defun calc-vector-find (arg)
 406   (interactive "P")
 407   (calc-wrapper
 408    (let ((func (cons 'calcFunc-find (calc-top-list-n 2))))
 409      (calc-enter-result
 410       2 "find"
 411       (if arg (append func (list (prefix-numeric-value arg))) func)))))
 412
 413 (defun calc-subvector ()
 414   (interactive)
 415   (calc-wrapper
 416    (if (calc-is-inverse)
 417        (calc-enter-result 3 "rsvc" (cons 'calcFunc-rsubvec
 418                                          (calc-top-list-n 3)))
 419      (calc-enter-result 3 "svec" (cons 'calcFunc-subvec (calc-top-list-n 3))))))
 420
 421 (defun calc-reverse-vector (arg)
 422   (interactive "P")
 423   (calc-wrapper
 424    (calc-unary-op "rev" 'calcFunc-rev arg)))
 425
 426 (defun calc-mask-vector (arg)
 427   (interactive "P")
 428   (calc-wrapper
 429    (calc-binary-op "vmsk" 'calcFunc-vmask arg)))
 430
 431 (defun calc-expand-vector (arg)
 432   (interactive "P")
 433   (calc-wrapper
 434    (if (calc-is-hyperbolic)
 435        (calc-enter-result 3 "vexp" (cons 'calcFunc-vexp (calc-top-list-n 3)))
 436      (calc-binary-op "vexp" 'calcFunc-vexp arg))))
 437
 438 (defun calc-sort ()
 439   (interactive)
 440   (calc-slow-wrapper
 441    (if (calc-is-inverse)
 442        (calc-enter-result 1 "rsrt" (list 'calcFunc-rsort (calc-top-n 1)))
 443      (calc-enter-result 1 "sort" (list 'calcFunc-sort (calc-top-n 1))))))
 444
 445 (defun calc-grade ()
 446   (interactive)
 447   (calc-slow-wrapper
 448    (if (calc-is-inverse)
 449        (calc-enter-result 1 "rgrd" (list 'calcFunc-rgrade (calc-top-n 1)))
 450      (calc-enter-result 1 "grad" (list 'calcFunc-grade (calc-top-n 1))))))
 451
 452 (defun calc-histogram (n)
 453   (interactive "P")
 454   (unless (natnump n)
 455     (setq n (math-read-expr (read-string "Centers of bins: "))))
 456   (calc-slow-wrapper
 457    (if calc-hyperbolic-flag
 458        (calc-enter-result 2 "hist" (list 'calcFunc-histogram
 459                                          (calc-top-n 2)
 460                                          (calc-top-n 1)
 461                                          n))
 462      (calc-enter-result 1 "hist" (list 'calcFunc-histogram
 463                                        (calc-top-n 1)
 464                                        n)))))
 465
 466 (defun calc-transpose (arg)
 467   (interactive "P")
 468   (calc-wrapper
 469    (calc-unary-op "trn" 'calcFunc-trn arg)))
 470
 471 (defun calc-conj-transpose (arg)
 472   (interactive "P")
 473   (calc-wrapper
 474    (calc-unary-op "ctrn" 'calcFunc-ctrn arg)))
 475
 476 (defun calc-cross (arg)
 477   (interactive "P")
 478   (calc-wrapper
 479    (calc-binary-op "cros" 'calcFunc-cross arg)))
 480
 481 (defun calc-kron (arg)
 482   (interactive "P")
 483   (calc-wrapper
 484    (calc-binary-op "kron" 'calcFunc-kron arg)))
 485
 486 (defun calc-remove-duplicates (arg)
 487   (interactive "P")
 488   (calc-wrapper
 489    (calc-unary-op "rdup" 'calcFunc-rdup arg)))
 490
 491 (defun calc-set-union (arg)
 492   (interactive "P")
 493   (calc-wrapper
 494    (calc-binary-op "unio" 'calcFunc-vunion arg '(vec) 'calcFunc-rdup)))
 495
 496 (defun calc-set-intersect (arg)
 497   (interactive "P")
 498   (calc-wrapper
 499    (calc-binary-op "intr" 'calcFunc-vint arg '(vec) 'calcFunc-rdup)))
 500
 501 (defun calc-set-difference (arg)
 502   (interactive "P")
 503   (calc-wrapper
 504    (calc-binary-op "diff" 'calcFunc-vdiff arg '(vec) 'calcFunc-rdup)))
 505
 506 (defun calc-set-xor (arg)
 507   (interactive "P")
 508   (calc-wrapper
 509    (calc-binary-op "xor" 'calcFunc-vxor arg '(vec) 'calcFunc-rdup)))
 510
 511 (defun calc-set-complement (arg)
 512   (interactive "P")
 513   (calc-wrapper
 514    (calc-unary-op "cmpl" 'calcFunc-vcompl arg)))
 515
 516 (defun calc-set-floor (arg)
 517   (interactive "P")
 518   (calc-wrapper
 519    (calc-unary-op "vflr" 'calcFunc-vfloor arg)))
 520
 521 (defun calc-set-enumerate (arg)
 522   (interactive "P")
 523   (calc-wrapper
 524    (calc-unary-op "enum" 'calcFunc-venum arg)))
 525
 526 (defun calc-set-span (arg)
 527   (interactive "P")
 528   (calc-wrapper
 529    (calc-unary-op "span" 'calcFunc-vspan arg)))
 530
 531 (defun calc-set-cardinality (arg)
 532   (interactive "P")
 533   (calc-wrapper
 534    (calc-unary-op "card" 'calcFunc-vcard arg)))
 535
 536 (defun calc-unpack-bits (arg)
 537   (interactive "P")
 538   (calc-wrapper
 539    (if (calc-is-inverse)
 540        (calc-unary-op "bpck" 'calcFunc-vpack arg)
 541      (calc-unary-op "bupk" 'calcFunc-vunpack arg))))
 542
 543 (defun calc-pack-bits (arg)
 544   (interactive "P")
 545   (calc-invert-func)
 546   (calc-unpack-bits arg))
 547
 548
 549 (defun calc-rnorm (arg)
 550   (interactive "P")
 551   (calc-wrapper
 552    (calc-unary-op "rnrm" 'calcFunc-rnorm arg)))
 553
 554 (defun calc-cnorm (arg)
 555   (interactive "P")
 556   (calc-wrapper
 557    (calc-unary-op "cnrm" 'calcFunc-cnorm arg)))
 558
 559 (defun calc-mrow (n &optional nn)
 560   (interactive "NRow number: \nP")
 561   (calc-wrapper
 562    (if (consp nn)
 563        (calc-enter-result 2 "mrow" (cons 'calcFunc-mrow (calc-top-list-n 2)))
 564      (setq n (prefix-numeric-value n))
 565      (if (= n 0)
 566          (calc-enter-result 1 "getd" (list 'calcFunc-getdiag (calc-top-n 1)))
 567        (if (< n 0)
 568            (calc-enter-result 1 "rrow" (list 'calcFunc-mrrow
 569                                              (calc-top-n 1) (- n)))
 570          (calc-enter-result 1 "mrow" (list 'calcFunc-mrow
 571                                            (calc-top-n 1) n)))))))
 572
 573 (defun calc-mcol (n &optional nn)
 574   (interactive "NColumn number: \nP")
 575   (calc-wrapper
 576    (if (consp nn)
 577        (calc-enter-result 2 "mcol" (cons 'calcFunc-mcol (calc-top-list-n 2)))
 578      (setq n (prefix-numeric-value n))
 579      (if (= n 0)
 580          (calc-enter-result 1 "getd" (list 'calcFunc-getdiag (calc-top-n 1)))
 581        (if (< n 0)
 582            (calc-enter-result 1 "rcol" (list 'calcFunc-mrcol
 583                                              (calc-top-n 1) (- n)))
 584          (calc-enter-result 1 "mcol" (list 'calcFunc-mcol
 585                                            (calc-top-n 1) n)))))))
 586
 587
 588 ;;;; Vectors.
 589
 590 (defun calcFunc-mdims (m)
 591   (or (math-vectorp m)
 592       (math-reject-arg m 'vectorp))
 593   (cons 'vec (math-mat-dimens m)))
 594
 595
 596 ;;; Apply a function elementwise to vector A.  [V X V; N X N] [Public]
 597 (defun math-map-vec (f a)
 598   (if (math-vectorp a)
 599       (cons 'vec (mapcar f (cdr a)))
 600     (funcall f a)))
 601
 602 (defun math-dimension-error ()
 603   (calc-record-why "*Dimension error")
 604   (signal 'wrong-type-argument nil))
 605
 606
 607 ;;; Build a vector out of a list of objects.  [Public]
 608 (defun calcFunc-vec (&rest objs)
 609   (cons 'vec objs))
 610
 611
 612 ;;; Build a constant vector or matrix.  [Public]
 613 (defun calcFunc-cvec (obj &rest dims)
 614   (math-make-vec-dimen obj dims))
 615
 616 (defun math-make-vec-dimen (obj dims)
 617   (if dims
 618       (if (natnump (car dims))
 619           (if (or (cdr dims)
 620                   (not (math-numberp obj)))
 621               (cons 'vec (copy-sequence
 622                           (make-list (car dims)
 623                                      (math-make-vec-dimen obj (cdr dims)))))
 624             (cons 'vec (make-list (car dims) obj)))
 625         (math-reject-arg (car dims) 'fixnatnump))
 626     obj))
 627
 628 (defun calcFunc-head (vec)
 629   (if (and (Math-vectorp vec)
 630            (cdr vec))
 631       (nth 1 vec)
 632     (calc-record-why 'vectorp vec)
 633     (list 'calcFunc-head vec)))
 634
 635 (defun calcFunc-tail (vec)
 636   (if (and (Math-vectorp vec)
 637            (cdr vec))
 638       (cons 'vec (cdr (cdr vec)))
 639     (calc-record-why 'vectorp vec)
 640     (list 'calcFunc-tail vec)))
 641
 642 (defun calcFunc-cons (head tail)
 643   (if (Math-vectorp tail)
 644       (cons 'vec (cons head (cdr tail)))
 645     (calc-record-why 'vectorp tail)
 646     (list 'calcFunc-cons head tail)))
 647
 648 (defun calcFunc-rhead (vec)
 649   (if (and (Math-vectorp vec)
 650            (cdr vec))
 651       (let ((vec (copy-sequence vec)))
 652         (setcdr (nthcdr (- (length vec) 2) vec) nil)
 653         vec)
 654     (calc-record-why 'vectorp vec)
 655     (list 'calcFunc-rhead vec)))
 656
 657 (defun calcFunc-rtail (vec)
 658   (if (and (Math-vectorp vec)
 659            (cdr vec))
 660       (nth (1- (length vec)) vec)
 661     (calc-record-why 'vectorp vec)
 662     (list 'calcFunc-rtail vec)))
 663
 664 (defun calcFunc-rcons (head tail)
 665   (if (Math-vectorp head)
 666       (append head (list tail))
 667     (calc-record-why 'vectorp head)
 668     (list 'calcFunc-rcons head tail)))
 669
 670
 671
 672 ;;; Apply a function elementwise to vectors A and B.  [O X O O] [Public]
 673 (defun math-map-vec-2 (f a b)
 674   (if (math-vectorp a)
 675       (if (math-vectorp b)
 676           (let ((v nil))
 677             (while (setq a (cdr a))
 678               (or (setq b (cdr b))
 679                   (math-dimension-error))
 680               (setq v (cons (funcall f (car a) (car b)) v)))
 681             (if a (math-dimension-error))
 682             (cons 'vec (nreverse v)))
 683         (let ((v nil))
 684           (while (setq a (cdr a))
 685             (setq v (cons (funcall f (car a) b) v)))
 686           (cons 'vec (nreverse v))))
 687     (if (math-vectorp b)
 688         (let ((v nil))
 689           (while (setq b (cdr b))
 690             (setq v (cons (funcall f a (car b)) v)))
 691           (cons 'vec (nreverse v)))
 692       (funcall f a b))))
 693
 694
 695
 696 ;;; "Reduce" a function over a vector (left-associatively).  [O X V] [Public]
 697 (defun math-reduce-vec (f a)
 698   (if (math-vectorp a)
 699       (if (cdr a)
 700           (let ((accum (car (setq a (cdr a)))))
 701             (while (setq a (cdr a))
 702               (setq accum (funcall f accum (car a))))
 703             accum)
 704         0)
 705     a))
 706
 707 ;;; Reduce a function over the columns of matrix A.  [V X V] [Public]
 708 (defun math-reduce-cols (f a)
 709   (if (math-matrixp a)
 710       (cons 'vec (math-reduce-cols-col-step f (cdr a) 1 (length (nth 1 a))))
 711     a))
 712
 713 (defun math-reduce-cols-col-step (f a col cols)
 714   (and (< col cols)
 715        (cons (math-reduce-cols-row-step f (nth col (car a)) col (cdr a))
 716              (math-reduce-cols-col-step f a (1+ col) cols))))
 717
 718 (defun math-reduce-cols-row-step (f tot col a)
 719   (if a
 720       (math-reduce-cols-row-step f
 721                                  (funcall f tot (nth col (car a)))
 722                                  col
 723                                  (cdr a))
 724     tot))
 725
 726
 727
 728 (defun math-dot-product (a b)
 729   (if (setq a (cdr a) b (cdr b))
 730       (let ((accum (math-mul (car a) (car b))))
 731         (while (setq a (cdr a) b (cdr b))
 732           (setq accum (math-add accum (math-mul (car a) (car b)))))
 733         accum)
 734     0))
 735
 736
 737 ;;; Return the number of elements in vector V.  [Public]
 738 (defun calcFunc-vlen (v)
 739   (if (math-vectorp v)
 740       (1- (length v))
 741     (if (math-objectp v)
 742         0
 743       (list 'calcFunc-vlen v))))
 744
 745 ;;; Get the Nth row of a matrix.
 746 (defun calcFunc-mrow (mat n)   ; [Public]
 747   (if (Math-vectorp n)
 748       (math-map-vec (function (lambda (x) (calcFunc-mrow mat x))) n)
 749     (if (and (eq (car-safe n) 'intv) (math-constp n))
 750         (calcFunc-subvec mat
 751                          (math-add (nth 2 n) (if (memq (nth 1 n) '(2 3)) 0 1))
 752                          (math-add (nth 3 n) (if (memq (nth 1 n) '(1 3)) 1 0)))
 753       (or (and (integerp (setq n (math-check-integer n)))
 754                (> n 0))
 755           (math-reject-arg n 'fixposintp))
 756       (or (Math-vectorp mat)
 757           (math-reject-arg mat 'vectorp))
 758       (or (nth n mat)
 759           (math-reject-arg n "*Index out of range")))))
 760
 761 (defun calcFunc-subscr (mat n &optional m)
 762   (setq mat (calcFunc-mrow mat n))
 763   (if m
 764       (if (math-num-integerp n)
 765           (calcFunc-mrow mat m)
 766         (calcFunc-mcol mat m))
 767     mat))
 768
 769 ;;; Get the Nth column of a matrix.
 770 (defun math-mat-col (mat n)
 771   (cons 'vec (mapcar (function (lambda (x) (elt x n))) (cdr mat))))
 772
 773 (defun calcFunc-mcol (mat n)   ; [Public]
 774   (if (Math-vectorp n)
 775       (calcFunc-trn
 776        (math-map-vec (function (lambda (x) (calcFunc-mcol mat x))) n))
 777     (if (and (eq (car-safe n) 'intv) (math-constp n))
 778         (if (math-matrixp mat)
 779             (math-map-vec (function (lambda (x) (calcFunc-mrow x n))) mat)
 780           (calcFunc-mrow mat n))
 781       (or (and (integerp (setq n (math-check-integer n)))
 782                (> n 0))
 783           (math-reject-arg n 'fixposintp))
 784       (or (Math-vectorp mat)
 785           (math-reject-arg mat 'vectorp))
 786       (or (if (math-matrixp mat)
 787               (and (< n (length (nth 1 mat)))
 788                    (math-mat-col mat n))
 789             (nth n mat))
 790           (math-reject-arg n "*Index out of range")))))
 791
 792 ;;; Remove the Nth row from a matrix.
 793 (defun math-mat-less-row (mat n)
 794   (if (<= n 0)
 795       (cdr mat)
 796     (cons (car mat)
 797           (math-mat-less-row (cdr mat) (1- n)))))
 798
 799 (defun calcFunc-mrrow (mat n)   ; [Public]
 800   (and (integerp (setq n (math-check-integer n)))
 801        (> n 0)
 802        (< n (length mat))
 803        (math-mat-less-row mat n)))
 804
 805 ;;; Remove the Nth column from a matrix.
 806 (defun math-mat-less-col (mat n)
 807   (cons 'vec (mapcar (function (lambda (x) (math-mat-less-row x n)))
 808                      (cdr mat))))
 809
 810 (defun calcFunc-mrcol (mat n)   ; [Public]
 811   (and (integerp (setq n (math-check-integer n)))
 812        (> n 0)
 813        (if (math-matrixp mat)
 814            (and (< n (length (nth 1 mat)))
 815                 (math-mat-less-col mat n))
 816          (math-mat-less-row mat n))))
 817
 818 (defun calcFunc-getdiag (mat)   ; [Public]
 819   (if (math-square-matrixp mat)
 820       (cons 'vec (math-get-diag-step (cdr mat) 1))
 821     (calc-record-why 'square-matrixp mat)
 822     (list 'calcFunc-getdiag mat)))
 823
 824 (defun math-get-diag-step (row n)
 825   (and row
 826        (cons (nth n (car row))
 827              (math-get-diag-step (cdr row) (1+ n)))))
 828
 829 (defun math-transpose (mat)   ; [Public]
 830   (let ((m nil)
 831         (col (length (nth 1 mat))))
 832     (while (> (setq col (1- col)) 0)
 833       (setq m (cons (math-mat-col mat col) m)))
 834     (cons 'vec m)))
 835
 836 (defun calcFunc-trn (mat)
 837   (if (math-vectorp mat)
 838       (if (math-matrixp mat)
 839           (math-transpose mat)
 840         (math-col-matrix mat))
 841     (if (math-numberp mat)
 842         mat
 843       (math-reject-arg mat 'matrixp))))
 844
 845 (defun calcFunc-ctrn (mat)
 846   (calcFunc-conj (calcFunc-trn mat)))
 847
 848 (defun calcFunc-pack (mode els)
 849   (or (Math-vectorp els) (math-reject-arg els 'vectorp))
 850   (if (and (Math-vectorp mode) (cdr mode))
 851       (setq mode (cdr mode))
 852     (or (integerp mode) (math-reject-arg mode 'fixnump)))
 853   (condition-case err
 854       (if (= (calc-pack-size mode) (1- (length els)))
 855           (calc-pack-items mode (cdr els))
 856         (math-reject-arg els "*Wrong number of elements"))
 857     (error (math-reject-arg els (nth 1 err)))))
 858
 859 (defun calcFunc-unpack (mode thing)
 860   (or (integerp mode) (math-reject-arg mode 'fixnump))
 861   (condition-case err
 862       (cons 'vec (calc-unpack-item mode thing))
 863     (error (math-reject-arg thing (nth 1 err)))))
 864
 865 (defun calcFunc-unpackt (mode thing)
 866   (let ((calc-unpack-with-type 'pair))
 867     (calcFunc-unpack mode thing)))
 868
 869 (defun calcFunc-arrange (vec cols)   ; [Public]
 870   (setq cols (math-check-fixnum cols t))
 871   (if (math-vectorp vec)
 872       (let* ((flat (math-flatten-vector vec))
 873              (mat (list 'vec))
 874              next)
 875         (if (<= cols 0)
 876             (nconc mat flat)
 877           (while (>= (length flat) cols)
 878             (setq next (nthcdr cols flat))
 879             (setcdr (nthcdr (1- cols) flat) nil)
 880             (setq mat (nconc mat (list (cons 'vec flat)))
 881                   flat next))
 882           (if flat
 883               (setq mat (nconc mat (list (cons 'vec flat)))))
 884           mat))))
 885
 886 (defun math-flatten-vector (vec)   ; [L V]
 887   (if (math-vectorp vec)
 888       (apply 'append (mapcar 'math-flatten-vector (cdr vec)))
 889     (list vec)))
 890
 891 (defun calcFunc-vconcat (a b)
 892   (math-normalize (list '| a b)))
 893
 894 (defun calcFunc-vconcatrev (a b)
 895   (math-normalize (list '| b a)))
 896
 897 (defun calcFunc-append (v1 v2)
 898   (if (and (math-vectorp v1) (math-vectorp v2))
 899       (append v1 (cdr v2))
 900     (list 'calcFunc-append v1 v2)))
 901
 902 (defun calcFunc-appendrev (v1 v2)
 903   (calcFunc-append v2 v1))
 904
 905
 906 ;;; Copy a matrix.  [Public]
 907 (defun math-copy-matrix (m)
 908   (if (math-vectorp (nth 1 m))
 909       (cons 'vec (mapcar 'copy-sequence (cdr m)))
 910     (copy-sequence m)))
 911
 912 ;;; Convert a scalar or vector into an NxN diagonal matrix.  [Public]
 913 (defun calcFunc-diag (a &optional n)
 914   (and n (not (integerp n))
 915        (setq n (math-check-fixnum n)))
 916   (if (math-vectorp a)
 917       (if (and n (/= (length a) (1+ n)))
 918           (list 'calcFunc-diag a n)
 919         (if (math-matrixp a)
 920             (if (and n (/= (length (elt a 1)) (1+ n)))
 921                 (list 'calcFunc-diag a n)
 922               a)
 923           (cons 'vec (math-diag-step (cdr a) 0 (1- (length a))))))
 924     (if n
 925         (cons 'vec (math-diag-step (make-list n a) 0 n))
 926       (list 'calcFunc-diag a))))
 927
 928 (defun calcFunc-idn (a &optional n)
 929   (if n
 930       (if (math-vectorp a)
 931           (math-reject-arg a 'numberp)
 932         (calcFunc-diag a n))
 933     (if (integerp calc-matrix-mode)
 934         (calcFunc-idn a calc-matrix-mode)
 935       (list 'calcFunc-idn a))))
 936
 937 (defun math-mimic-ident (a m)
 938   (if (math-square-matrixp m)
 939       (calcFunc-idn a (1- (length m)))
 940     (if (math-vectorp m)
 941         (if (math-zerop a)
 942             (cons 'vec (mapcar (function (lambda (x)
 943                                            (if (math-vectorp x)
 944                                                (math-mimic-ident a x)
 945                                              a)))
 946                                (cdr m)))
 947           (math-dimension-error))
 948       (calcFunc-idn a))))
 949
 950 (defun math-diag-step (a n m)
 951   (if (< n m)
 952       (cons (cons 'vec
 953                   (nconc (make-list n 0)
 954                          (cons (car a)
 955                                (make-list (1- (- m n)) 0))))
 956             (math-diag-step (cdr a) (1+ n) m))
 957     nil))
 958
 959 ;;; Create a vector of consecutive integers. [Public]
 960 (defun calcFunc-index (n &optional start incr)
 961   (if (math-messy-integerp n)
 962       (math-float (calcFunc-index (math-trunc n) start incr))
 963     (and (not (integerp n))
 964          (setq n (math-check-fixnum n)))
 965     (let ((vec nil))
 966       (if start
 967           (progn
 968             (if (>= n 0)
 969                 (while (>= (setq n (1- n)) 0)
 970                   (setq vec (cons start vec)
 971                         start (math-add start (or incr 1))))
 972               (while (<= (setq n (1+ n)) 0)
 973                 (setq vec (cons start vec)
 974                       start (math-mul start (or incr 2)))))
 975             (setq vec (nreverse vec)))
 976         (if (>= n 0)
 977             (while (> n 0)
 978               (setq vec (cons n vec)
 979                     n (1- n)))
 980           (let ((i -1))
 981             (while (>= i n)
 982               (setq vec (cons i vec)
 983                     i (1- i))))))
 984       (cons 'vec vec))))
 985
 986 ;;; Find an element in a vector.
 987 (defun calcFunc-find (vec x &optional start)
 988   (setq start (if start (math-check-fixnum start t) 1))
 989   (if (< start 1) (math-reject-arg start 'posp))
 990   (setq vec (nthcdr start vec))
 991   (let ((n start))
 992     (while (and vec (not (Math-equal x (car vec))))
 993       (setq n (1+ n)
 994             vec (cdr vec)))
 995     (if vec n 0)))
 996
 997 ;;; Return a subvector of a vector.
 998 (defun calcFunc-subvec (vec start &optional end)
 999   (setq start (math-check-fixnum start t)
1000         end (math-check-fixnum (or end 0) t))
1001   (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
1002   (let ((len (1- (length vec))))
1003     (if (<= start 0)
1004         (setq start (+ len start 1)))
1005     (if (<= end 0)
1006         (setq end (+ len end 1)))
1007     (if (or (> start len)
1008             (<= end start))
1009         '(vec)
1010       (setq vec (nthcdr start vec))
1011       (if (<= end len)
1012           (let ((chop (nthcdr (- end start 1) (setq vec (copy-sequence vec)))))
1013             (setcdr chop nil)))
1014       (cons 'vec vec))))
1015
1016 ;;; Remove a subvector from a vector.
1017 (defun calcFunc-rsubvec (vec start &optional end)
1018   (setq start (math-check-fixnum start t)
1019         end (math-check-fixnum (or end 0) t))
1020   (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
1021   (let ((len (1- (length vec))))
1022     (if (<= start 0)
1023         (setq start (+ len start 1)))
1024     (if (<= end 0)
1025         (setq end (+ len end 1)))
1026     (if (or (> start len)
1027             (<= end start))
1028         vec
1029       (let ((tail (nthcdr end vec))
1030             (chop (nthcdr (1- start) (setq vec (copy-sequence vec)))))
1031         (setcdr chop nil)
1032         (append vec tail)))))
1033
1034 ;;; Reverse the order of the elements of a vector.
1035 (defun calcFunc-rev (vec)
1036   (if (math-vectorp vec)
1037       (cons 'vec (reverse (cdr vec)))
1038     (math-reject-arg vec 'vectorp)))
1039
1040 ;;; Compress a vector according to a mask vector.
1041 (defun calcFunc-vmask (mask vec)
1042   (if (math-numberp mask)
1043       (if (math-zerop mask)
1044           '(vec)
1045         vec)
1046     (or (math-vectorp mask) (math-reject-arg mask 'vectorp))
1047     (or (math-constp mask) (math-reject-arg mask 'constp))
1048     (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
1049     (or (= (length mask) (length vec)) (math-dimension-error))
1050     (let ((new nil))
1051       (while (setq mask (cdr mask) vec (cdr vec))
1052         (or (math-zerop (car mask))
1053             (setq new (cons (car vec) new))))
1054       (cons 'vec (nreverse new)))))
1055
1056 ;;; Expand a vector according to a mask vector.
1057 (defun calcFunc-vexp (mask vec &optional filler)
1058   (or (math-vectorp mask) (math-reject-arg mask 'vectorp))
1059   (or (math-constp mask) (math-reject-arg mask 'constp))
1060   (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
1061   (let ((new nil)
1062         (fvec (and filler (math-vectorp filler))))
1063     (while (setq mask (cdr mask))
1064       (if (math-zerop (car mask))
1065           (setq new (cons (or (if fvec
1066                                   (car (setq filler (cdr filler)))
1067                                 filler)
1068                               (car mask)) new))
1069         (setq vec (cdr vec)
1070               new (cons (or (car vec) (car mask)) new))))
1071     (cons 'vec (nreverse new))))
1072
1073
1074 ;;; Compute the row and column norms of a vector or matrix.  [Public]
1075 (defun calcFunc-rnorm (a)
1076   (if (and (Math-vectorp a)
1077            (math-constp a))
1078       (if (math-matrixp a)
1079           (math-reduce-vec 'math-max (math-map-vec 'calcFunc-cnorm a))
1080         (math-reduce-vec 'math-max (math-map-vec 'math-abs a)))
1081     (calc-record-why 'vectorp a)
1082     (list 'calcFunc-rnorm a)))
1083
1084 (defun calcFunc-cnorm (a)
1085   (if (and (Math-vectorp a)
1086            (math-constp a))
1087       (if (math-matrixp a)
1088           (math-reduce-vec 'math-max
1089                            (math-reduce-cols 'math-add-abs a))
1090         (math-reduce-vec 'math-add-abs a))
1091     (calc-record-why 'vectorp a)
1092     (list 'calcFunc-cnorm a)))
1093
1094 (defun math-add-abs (a b)
1095   (math-add (math-abs a) (math-abs b)))
1096
1097
1098 ;;; Sort the elements of a vector into increasing order.
1099 (defun calcFunc-sort (vec)   ; [Public]
1100   (if (math-vectorp vec)
1101       (cons 'vec (sort (copy-sequence (cdr vec)) 'math-beforep))
1102     (math-reject-arg vec 'vectorp)))
1103
1104 (defun calcFunc-rsort (vec)   ; [Public]
1105   (if (math-vectorp vec)
1106       (cons 'vec (nreverse (sort (copy-sequence (cdr vec)) 'math-beforep)))
1107     (math-reject-arg vec 'vectorp)))
1108
1109 ;; The variable math-grade-vec is local to calcFunc-grade and
1110 ;; calcFunc-rgrade, but is used by math-grade-beforep, which is called
1111 ;; by calcFunc-grade and calcFunc-rgrade.
1112 (defvar math-grade-vec)
1113
1114 (defun calcFunc-grade (math-grade-vec)
1115   (if (math-vectorp math-grade-vec)
1116       (let* ((len (1- (length math-grade-vec))))
1117         (cons 'vec (sort (cdr (calcFunc-index len)) 'math-grade-beforep)))
1118     (math-reject-arg math-grade-vec 'vectorp)))
1119
1120 (defun calcFunc-rgrade (math-grade-vec)
1121   (if (math-vectorp math-grade-vec)
1122       (let* ((len (1- (length math-grade-vec))))
1123         (cons 'vec (nreverse (sort (cdr (calcFunc-index len))
1124                                    'math-grade-beforep))))
1125     (math-reject-arg math-grade-vec 'vectorp)))
1126
1127 (defun math-grade-beforep (i j)
1128   (math-beforep (nth i math-grade-vec) (nth j math-grade-vec)))
1129
1130
1131 ;;; Compile a histogram of data from a vector.
1132 (defun calcFunc-histogram (vec wts &optional n)
1133   (or n (setq n wts wts 1))
1134   (or (Math-vectorp vec)
1135       (math-reject-arg vec 'vectorp))
1136   (if (Math-vectorp wts)
1137       (or (= (length vec) (length wts))
1138           (math-dimension-error)))
1139   (cond ((natnump n)
1140          (let ((res (make-vector n 0))
1141                (vp vec)
1142                (wvec (Math-vectorp wts))
1143                (wp wts)
1144                bin)
1145            (while (setq vp (cdr vp))
1146              (setq bin (car vp))
1147              (or (natnump bin)
1148                  (setq bin (math-floor bin)))
1149             (and (natnump bin)
1150                  (< bin n)
1151                  (aset res bin
1152                        (math-add (aref res bin)
1153                                  (if wvec (car (setq wp (cdr wp))) wts)))))
1154            (cons 'vec (append res nil))))
1155         ((Math-vectorp n) ;; n is a vector of midpoints
1156          (let* ((bds (math-vector-avg n))
1157                 (res (make-vector (1- (length n)) 0))
1158                 (vp (cdr vec))
1159                 (wvec (Math-vectorp wts))
1160                 (wp wts)
1161                 num)
1162            (while vp
1163              (setq num (car vp))
1164              (let ((tbds (cdr bds))
1165                    (i 0))
1166                (while (and tbds (Math-lessp (car tbds) num))
1167                  (setq i (1+ i))
1168                  (setq tbds (cdr tbds)))
1169                (aset res i
1170                      (math-add (aref res i)
1171                                (if wvec (car (setq wp (cdr wp))) wts))))
1172              (setq vp (cdr vp)))
1173            (cons 'vec (append res nil))))
1174         (t
1175          (math-reject-arg n "*Expecting an integer or vector"))))
1176
1177 ;;; Replace a vector [a b c ...] with a vector of averages
1178 ;;; [(a+b)/2 (b+c)/2 ...]
1179 (defun math-vector-avg (vec)
1180   (let ((vp (sort (copy-sequence (cdr vec)) 'math-beforep))
1181         (res nil))
1182     (while (and vp (cdr vp))
1183       (setq res (cons (math-div (math-add (car vp) (cadr vp)) 2) res)
1184             vp (cdr vp)))
1185     (cons 'vec (reverse res))))
1186
1187
1188 ;;; Set operations.
1189
1190 (defun calcFunc-vunion (a b)
1191   (if (Math-objectp a)
1192       (setq a (list 'vec a))
1193     (or (math-vectorp a) (math-reject-arg a 'vectorp)))
1194   (if (Math-objectp b)
1195       (setq b (list b))
1196     (or (math-vectorp b) (math-reject-arg b 'vectorp))
1197     (setq b (cdr b)))
1198   (calcFunc-rdup (append a b)))
1199
1200 (defun calcFunc-vint (a b)
1201   (if (and (math-simple-set a) (math-simple-set b))
1202       (progn
1203         (setq a (cdr (calcFunc-rdup a)))
1204         (setq b (cdr (calcFunc-rdup b)))
1205         (let ((vec (list 'vec)))
1206           (while (and a b)
1207             (if (math-beforep (car a) (car b))
1208                 (setq a (cdr a))
1209               (if (Math-equal (car a) (car b))
1210                   (setq vec (cons (car a) vec)
1211                         a (cdr a)))
1212               (setq b (cdr b))))
1213           (nreverse vec)))
1214     (calcFunc-vcompl (calcFunc-vunion (calcFunc-vcompl a)
1215                                       (calcFunc-vcompl b)))))
1216
1217 (defun calcFunc-vdiff (a b)
1218   (if (and (math-simple-set a) (math-simple-set b))
1219       (progn
1220         (setq a (cdr (calcFunc-rdup a)))
1221         (setq b (cdr (calcFunc-rdup b)))
1222         (let ((vec (list 'vec)))
1223           (while a
1224             (while (and b (math-beforep (car b) (car a)))
1225               (setq b (cdr b)))
1226             (if (and b (Math-equal (car a) (car b)))
1227                 (setq a (cdr a)
1228                       b (cdr b))
1229               (setq vec (cons (car a) vec)
1230                     a (cdr a))))
1231           (nreverse vec)))
1232     (calcFunc-vcompl (calcFunc-vunion (calcFunc-vcompl a) b))))
1233
1234 (defun calcFunc-vxor (a b)
1235   (if (and (math-simple-set a) (math-simple-set b))
1236       (progn
1237         (setq a (cdr (calcFunc-rdup a)))
1238         (setq b (cdr (calcFunc-rdup b)))
1239         (let ((vec (list 'vec)))
1240           (while (or a b)
1241             (if (and a
1242                      (or (not b)
1243                          (math-beforep (car a) (car b))))
1244                 (setq vec (cons (car a) vec)
1245                       a (cdr a))
1246               (if (and a (Math-equal (car a) (car b)))
1247                   (setq a (cdr a))
1248                 (setq vec (cons (car b) vec)))
1249               (setq b (cdr b))))
1250           (nreverse vec)))
1251     (let ((ca (calcFunc-vcompl a))
1252           (cb (calcFunc-vcompl b)))
1253       (calcFunc-vunion (calcFunc-vcompl (calcFunc-vunion ca b))
1254                        (calcFunc-vcompl (calcFunc-vunion a cb))))))
1255
1256 (defun calcFunc-vcompl (a)
1257   (setq a (math-prepare-set a))
1258   (let ((vec (list 'vec))
1259         (prev '(neg (var inf var-inf)))
1260         (closed 2))
1261     (while (setq a (cdr a))
1262       (or (and (equal (nth 2 (car a)) '(neg (var inf var-inf)))
1263                (memq (nth 1 (car a)) '(2 3)))
1264           (setq vec (cons (list 'intv
1265                                 (+ closed
1266                                    (if (memq (nth 1 (car a)) '(0 1)) 1 0))
1267                                 prev
1268                                 (nth 2 (car a)))
1269                           vec)))
1270       (setq prev (nth 3 (car a))
1271             closed (if (memq (nth 1 (car a)) '(0 2)) 2 0)))
1272     (or (and (equal prev '(var inf var-inf))
1273              (= closed 0))
1274         (setq vec (cons (list 'intv (+ closed 1)
1275                               prev '(var inf var-inf))
1276                         vec)))
1277     (math-clean-set (nreverse vec))))
1278
1279 (defun calcFunc-vspan (a)
1280   (setq a (math-prepare-set a))
1281   (if (cdr a)
1282       (let ((last (nth (1- (length a)) a)))
1283         (math-make-intv (+ (logand (nth 1 (nth 1 a)) 2)
1284                            (logand (nth 1 last) 1))
1285                         (nth 2 (nth 1 a))
1286                         (nth 3 last)))
1287     '(intv 2 0 0)))
1288
1289 (defun calcFunc-vfloor (a &optional always-vec)
1290   (setq a (math-prepare-set a))
1291   (let ((vec (list 'vec)) (p a) (prev nil) b mask)
1292     (while (setq p (cdr p))
1293       (setq mask (nth 1 (car p))
1294             a (nth 2 (car p))
1295             b (nth 3 (car p)))
1296       (and (memq mask '(0 1))
1297            (not (math-infinitep a))
1298            (setq mask (logior mask 2))
1299            (math-num-integerp a)
1300            (setq a (math-add a 1)))
1301       (setq a (math-ceiling a))
1302       (and (memq mask '(0 2))
1303            (not (math-infinitep b))
1304            (setq mask (logior mask 1))
1305            (math-num-integerp b)
1306            (setq b (math-sub b 1)))
1307       (setq b (math-floor b))
1308       (if (and prev (Math-equal (math-sub a 1) (nth 3 prev)))
1309           (setcar (nthcdr 3 prev) b)
1310         (or (Math-lessp b a)
1311             (setq vec (cons (setq prev (list 'intv mask a b)) vec)))))
1312     (setq vec (nreverse vec))
1313     (math-clean-set vec always-vec)))
1314
1315 (defun calcFunc-vcard (a)
1316   (setq a (calcFunc-vfloor a t))
1317   (or (math-constp a) (math-reject-arg a "*Set must be finite"))
1318   (let ((count 0))
1319     (while (setq a (cdr a))
1320       (if (eq (car-safe (car a)) 'intv)
1321           (setq count (math-add count (math-sub (nth 3 (car a))
1322                                                 (nth 2 (car a))))))
1323       (setq count (math-add count 1)))
1324     count))
1325
1326 (defun calcFunc-venum (a)
1327   (setq a (calcFunc-vfloor a t))
1328   (or (math-constp a) (math-reject-arg a "*Set must be finite"))
1329   (let* ((prev a) (this (cdr prev)) this-val next this-last)
1330     (while this
1331       (setq next (cdr this)
1332                         this-val (car this))
1333       (if (eq (car-safe this-val) 'intv)
1334                   (progn
1335                         (setq this (cdr (calcFunc-index (math-add
1336                                                                                          (math-sub (nth 3 this-val)
1337                                                                                                            (nth 2 this-val))
1338                                                                                          1)
1339                                                                                         (nth 2 this-val))))
1340                         (setq this-last (last this))
1341                         (setcdr this-last next)
1342                         (setcdr prev this)
1343                         (setq prev this-last))
1344                 (setq prev this))
1345           (setq this next)))
1346   a)
1347
1348 (defun calcFunc-vpack (a)
1349   (setq a (calcFunc-vfloor a t))
1350   (if (and (cdr a)
1351            (math-negp (if (eq (car-safe (nth 1 a)) 'intv)
1352                           (nth 2 (nth 1 a))
1353                         (nth 1 a))))
1354       (math-reject-arg (nth 1 a) 'posp))
1355   (let ((accum 0))
1356     (while (setq a (cdr a))
1357       (if (eq (car-safe (car a)) 'intv)
1358           (if (equal (nth 3 (car a)) '(var inf var-inf))
1359               (setq accum (math-sub accum
1360                                     (math-power-of-2 (nth 2 (car a)))))
1361             (setq accum (math-add accum
1362                                   (math-sub
1363                                    (math-power-of-2 (1+ (nth 3 (car a))))
1364                                    (math-power-of-2 (nth 2 (car a)))))))
1365         (setq accum (math-add accum (math-power-of-2 (car a))))))
1366     accum))
1367
1368 (defun calcFunc-vunpack (a &optional w)
1369   (or (math-num-integerp a) (math-reject-arg a 'integerp))
1370   (if w (setq a (math-clip a w)))
1371   (if (math-messy-integerp a) (setq a (math-trunc a)))
1372   (let* ((calc-number-radix 2)
1373          (calc-twos-complement-mode nil)
1374          (neg (math-negp a))
1375          (aa (if neg (math-sub -1 a) a))
1376          (str (if (eq aa 0)
1377                   ""
1378                 (if (consp aa)
1379                     (math-format-bignum-binary (cdr aa))
1380                   (math-format-binary aa))))
1381          (zero (if neg ?1 ?0))
1382          (one (if neg ?0 ?1))
1383          (len (length str))
1384          (vec (list 'vec))
1385          (pos (1- len)) pos2)
1386     (while (>= pos 0)
1387       (if (eq (aref str pos) zero)
1388           (setq pos (1- pos))
1389         (setq pos2 pos)
1390         (while (and (>= pos 0) (eq (aref str pos) one))
1391           (setq pos (1- pos)))
1392         (setq vec (cons (if (= pos (1- pos2))
1393                             (- len pos2 1)
1394                           (list 'intv 3 (- len pos2 1) (- len pos 2)))
1395                         vec))))
1396     (if neg
1397         (setq vec (cons (list 'intv 2 len '(var inf var-inf)) vec)))
1398     (math-clean-set (nreverse vec))))
1399
1400 (defun calcFunc-rdup (a)
1401   (if (math-simple-set a)
1402       (progn
1403         (and (Math-objectp a) (setq a (list 'vec a)))
1404         (or (math-vectorp a) (math-reject-arg a 'vectorp))
1405         (setq a (sort (copy-sequence (cdr a)) 'math-beforep))
1406         (let ((p a))
1407           (while (cdr p)
1408             (if (Math-equal (car p) (nth 1 p))
1409                 (setcdr p (cdr (cdr p)))
1410               (setq p (cdr p)))))
1411         (cons 'vec a))
1412     (math-clean-set (math-prepare-set a))))
1413
1414 (defun math-prepare-set (a)
1415   (if (Math-objectp a)
1416       (setq a (list 'vec a))
1417     (or (math-vectorp a) (math-reject-arg a 'vectorp))
1418     (setq a (cons 'vec (sort (copy-sequence (cdr a)) 'math-beforep))))
1419   (let ((p a) res)
1420
1421     ;; Convert all elements to non-empty intervals.
1422     (while (cdr p)
1423       (if (eq (car-safe (nth 1 p)) 'intv)
1424           (if (math-intv-constp (nth 1 p))
1425               (if (and (memq (nth 1 (nth 1 p)) '(0 1 2))
1426                        (Math-equal (nth 2 (nth 1 p)) (nth 3 (nth 1 p))))
1427                   (setcdr p (cdr (cdr p)))
1428                 (setq p (cdr p)))
1429             (math-reject-arg (nth 1 p) 'constp))
1430         (or (Math-anglep (nth 1 p))
1431             (eq (car (nth 1 p)) 'date)
1432             (equal (nth 1 p) '(var inf var-inf))
1433             (equal (nth 1 p) '(neg (var inf var-inf)))
1434             (math-reject-arg (nth 1 p) 'realp))
1435         (setcar (cdr p) (list 'intv 3 (nth 1 p) (nth 1 p)))
1436         (setq p (cdr p))))
1437
1438     ;; Combine redundant intervals.
1439     (setq p a)
1440     (while (cdr (cdr p))
1441       (if (or (memq (setq res (math-compare (nth 3 (nth 1 p))
1442                                             (nth 2 (nth 2 p))))
1443                     '(-1 2))
1444               (and (eq res 0)
1445                    (memq (nth 1 (nth 1 p)) '(0 2))
1446                    (memq (nth 1 (nth 2 p)) '(0 1))))
1447           (setq p (cdr p))
1448         (setq res (math-compare (nth 3 (nth 1 p)) (nth 3 (nth 2 p))))
1449         (setcdr p (cons (list 'intv
1450                               (+ (logand (logior (nth 1 (nth 1 p))
1451                                                  (if (Math-equal
1452                                                       (nth 2 (nth 1 p))
1453                                                       (nth 2 (nth 2 p)))
1454                                                      (nth 1 (nth 2 p))
1455                                                    0))
1456                                          2)
1457                                  (logand (logior (if (memq res '(1 0 2))
1458                                                      (nth 1 (nth 1 p)) 0)
1459                                                  (if (memq res '(-1 0 2))
1460                                                      (nth 1 (nth 2 p)) 0))
1461                                          1))
1462                               (nth 2 (nth 1 p))
1463                               (if (eq res 1)
1464                                   (nth 3 (nth 1 p))
1465                                 (nth 3 (nth 2 p))))
1466                         (cdr (cdr (cdr p))))))))
1467   a)
1468
1469 (defun math-clean-set (a &optional always-vec)
1470   (let ((p a) res)
1471     (while (cdr p)
1472       (if (and (eq (car-safe (nth 1 p)) 'intv)
1473                (Math-equal (nth 2 (nth 1 p)) (nth 3 (nth 1 p))))
1474           (setcar (cdr p) (nth 2 (nth 1 p))))
1475       (setq p (cdr p)))
1476     (if (and (not (cdr (cdr a)))
1477              (eq (car-safe (nth 1 a)) 'intv)
1478              (not always-vec))
1479         (nth 1 a)
1480       a)))
1481
1482 (defun math-simple-set (a)
1483   (or (and (Math-objectp a)
1484            (not (eq (car-safe a) 'intv)))
1485       (and (Math-vectorp a)
1486            (progn
1487              (while (and (setq a (cdr a))
1488                          (not (eq (car-safe (car a)) 'intv))))
1489              (null a)))))
1490
1491
1492
1493
1494 ;;; Compute a right-handed vector cross product.  [O O O] [Public]
1495 (defun calcFunc-cross (a b)
1496   (if (and (eq (car-safe a) 'vec)
1497            (= (length a) 4))
1498       (if (and (eq (car-safe b) 'vec)
1499                (= (length b) 4))
1500           (list 'vec
1501                 (math-sub (math-mul (nth 2 a) (nth 3 b))
1502                           (math-mul (nth 3 a) (nth 2 b)))
1503                 (math-sub (math-mul (nth 3 a) (nth 1 b))
1504                           (math-mul (nth 1 a) (nth 3 b)))
1505                 (math-sub (math-mul (nth 1 a) (nth 2 b))
1506                           (math-mul (nth 2 a) (nth 1 b))))
1507         (math-reject-arg b "*Three-vector expected"))
1508     (math-reject-arg a "*Three-vector expected")))
1509
1510
1511 ;;; Compute a Kronecker product
1512 (defun calcFunc-kron (x y &optional nocheck)
1513   "The Kronecker product of objects X and Y.
1514 The objects X and Y may be scalars, vectors or matrices.
1515 The type of the result depends on the types of the operands;
1516 the product of two scalars is a scalar,
1517 of one scalar and a vector is a vector,
1518 of two vectors is a vector.
1519 of one vector and a matrix is a matrix,
1520 of two matrices is a matrix."
1521   (unless nocheck
1522     (cond ((or (math-matrixp x)
1523                (math-matrixp y))
1524            (unless (math-matrixp x)
1525              (setq x (if (math-vectorp x)
1526                          (list 'vec x)
1527                        (list 'vec (list 'vec x)))))
1528            (unless (math-matrixp y)
1529              (setq y (if (math-vectorp y)
1530                          (list 'vec y)
1531                        (list 'vec (list 'vec y))))))
1532           ((or (math-vectorp x)
1533                (math-vectorp y))
1534            (unless (math-vectorp x)
1535              (setq x (list 'vec x)))
1536            (unless (math-vectorp y)
1537              (setq y (list 'vec y))))))
1538   (if (math-vectorp x)
1539       (let (ret)
1540         (dolist (v (cdr x))
1541           (dolist (w (cdr y))
1542             (setq ret (cons (calcFunc-kron v w t) ret))))
1543         (cons 'vec (nreverse ret)))
1544     (math-mul x y)))
1545
1546
1547 ;; The variable math-rb-close is local to math-read-brackets, but
1548 ;; is used by math-read-vector, which is called (directly and
1549 ;; indirectly) by math-read-brackets.
1550 (defvar math-rb-close)
1551
1552 ;; The next few variables are local to math-read-exprs in calc-aent.el
1553 ;; and math-read-expr in calc-ext.el, but are set in functions they call.
1554 (defvar math-exp-pos)
1555 (defvar math-exp-str)
1556 (defvar math-exp-old-pos)
1557 (defvar math-exp-token)
1558 (defvar math-exp-keep-spaces)
1559 (defvar math-expr-data)
1560
1561 (defun math-read-brackets (space-sep math-rb-close)
1562   (and space-sep (setq space-sep (not (math-check-for-commas))))
1563   (math-read-token)
1564   (while (eq math-exp-token 'space)
1565     (math-read-token))
1566   (if (or (equal math-expr-data math-rb-close)
1567           (eq math-exp-token 'end))
1568       (progn
1569         (math-read-token)
1570         '(vec))
1571     (let ((save-exp-pos math-exp-pos)
1572           (save-exp-old-pos math-exp-old-pos)
1573           (save-exp-token math-exp-token)
1574           (save-exp-data math-expr-data)
1575           (vals (let ((math-exp-keep-spaces space-sep))
1576                   (if (or (equal math-expr-data "\\dots")
1577                           (equal math-expr-data "\\ldots"))
1578                       '(vec (neg (var inf var-inf)))
1579                     (catch 'syntax (math-read-vector))))))
1580       (if (stringp vals)
1581           (if space-sep
1582               (let ((error-exp-pos math-exp-pos)
1583                     (error-exp-old-pos math-exp-old-pos)
1584                     vals2)
1585                 (setq math-exp-pos save-exp-pos
1586                       math-exp-old-pos save-exp-old-pos
1587                       math-exp-token save-exp-token
1588                       math-expr-data save-exp-data)
1589                 (let ((math-exp-keep-spaces nil))
1590                   (setq vals2 (catch 'syntax (math-read-vector))))
1591                 (if (and (not (stringp vals2))
1592                          (or (assoc math-expr-data '(("\\ldots") ("\\dots") (";")))
1593                              (equal math-expr-data math-rb-close)
1594                              (eq math-exp-token 'end)))
1595                     (setq space-sep nil
1596                           vals vals2)
1597                   (setq math-exp-pos error-exp-pos
1598                         math-exp-old-pos error-exp-old-pos)
1599                   (throw 'syntax vals)))
1600             (throw 'syntax vals)))
1601       (if (or (equal math-expr-data "\\dots")
1602               (equal math-expr-data "\\ldots"))
1603           (progn
1604             (math-read-token)
1605             (setq vals (if (> (length vals) 2)
1606                            (cons 'calcFunc-mul (cdr vals)) (nth 1 vals)))
1607             (let ((exp2 (if (or (equal math-expr-data math-rb-close)
1608                                 (equal math-expr-data ")")
1609                                 (eq math-exp-token 'end))
1610                             '(var inf var-inf)
1611                           (math-read-expr-level 0))))
1612               (setq vals
1613                     (list 'intv
1614                           (if (equal math-expr-data ")") 2 3)
1615                           vals
1616                           exp2)))
1617             (if (not (or (equal math-expr-data math-rb-close)
1618                          (equal math-expr-data ")")
1619                          (eq math-exp-token 'end)))
1620                 (throw 'syntax "Expected `]'")))
1621         (if (equal math-expr-data ";")
1622             (let ((math-exp-keep-spaces space-sep))
1623               (setq vals (cons 'vec (math-read-matrix (list vals))))))
1624         (if (not (or (equal math-expr-data math-rb-close)
1625                      (eq math-exp-token 'end)))
1626             (throw 'syntax "Expected `]'")))
1627       (or (eq math-exp-token 'end)
1628           (math-read-token))
1629       vals)))
1630
1631 (defun math-check-for-commas (&optional balancing)
1632   (let ((count 0)
1633         (pos (1- math-exp-pos)))
1634     (while (and (>= count 0)
1635                 (setq pos (string-match
1636                            (if balancing "[],[{}()<>]" "[],[{}()]")
1637                            math-exp-str (1+ pos)))
1638                 (or (/= (aref math-exp-str pos) ?,) (> count 0) balancing))
1639       (cond ((memq (aref math-exp-str pos) '(?\[ ?\{ ?\( ?\<))
1640              (setq count (1+ count)))
1641             ((memq (aref math-exp-str pos) '(?\] ?\} ?\) ?\>))
1642              (setq count (1- count)))))
1643     (if balancing
1644         pos
1645       (and pos (= (aref math-exp-str pos) ?,)))))
1646
1647 (defun math-read-vector ()
1648   (let* ((val (list (math-read-expr-level 0)))
1649          (last val))
1650     (while (progn
1651              (while (eq math-exp-token 'space)
1652                (math-read-token))
1653              (and (not (eq math-exp-token 'end))
1654                   (not (equal math-expr-data ";"))
1655                   (not (equal math-expr-data math-rb-close))
1656                   (not (equal math-expr-data "\\dots"))
1657                   (not (equal math-expr-data "\\ldots"))))
1658       (if (equal math-expr-data ",")
1659           (math-read-token))
1660       (while (eq math-exp-token 'space)
1661         (math-read-token))
1662       (let ((rest (list (math-read-expr-level 0))))
1663         (setcdr last rest)
1664         (setq last rest)))
1665     (cons 'vec val)))
1666
1667 (defun math-read-matrix (mat)
1668   (while (equal math-expr-data ";")
1669     (math-read-token)
1670     (while (eq math-exp-token 'space)
1671       (math-read-token))
1672     (setq mat (nconc mat (list (math-read-vector)))))
1673   mat)
1674
1675 (provide 'calc-vec)
1676
1677 ;;; calc-vec.el ends here