lisp/calc/calc-vec.el

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