lisp/calc/calc-vec.el

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