lisp/calc/calc-vec.el

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