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