cgit.sxemacs.org Git - packages/blob - xemacs-packages/calc/calc-vec.el

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