lisp/calc/calc-vec.el

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