1;; Copyright 2002-2003 by 2;; Stavros Macrakis (macrakis@alum.mit.edu) and 3;; Barton Willis (willisb@unk.edu) 4 5;; Maxima nset is free software; you can redistribute it and/or 6;; modify it under the terms of the GNU General Public License, 7;; http://www.gnu.org/copyleft/gpl.html. 8 9;; Maxima nset has NO WARRANTY, not even the implied warranty of 10;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 11 12;; A Maxima set package 13 14(in-package :maxima) 15 16(macsyma-module nset) 17 18($put '$nset 1.21 '$version) 19;; Let's remove built-in symbols from list for user-defined properties. 20(setq $props (remove '$nset $props)) 21 22;; Display sets as { .. }. 23 24(defprop $set msize-matchfix grind) 25 26(setf (get '$set 'dissym) '((#\{ ) #\} )) 27(setf (get '$set 'dimension) 'dimension-match) 28 29;; Parse {a, b, c} into set(a, b, c). 30 31(setf (get '$set 'op) "{") 32 33(setf (get '|$}| 'nud) 'delim-err) 34(setf (get '|$}| 'led) 'erb-err) 35 36(def-lbp |$}| 5.) 37 38(setf (get '|${| 'nud) 'parse-matchfix) 39(def-match |${| |$}|) 40(def-lbp |${| 200.) 41;No RBP 42(def-mheader |${| ($set)) 43(def-pos |${| $any) 44;No LPOS 45;No RPOS 46 47(def-operator "{" '$any nil '$any nil nil nil nil '(nud . parse-matchfix) 'msize-matchfix 'dimension-match "}") 48;; Let's remove built-in operators from list for user-defined properties. 49(setq $props (remove "{" $props :test #'equal)) 50(setq $props (remove "}" $props :test #'equal)) 51 52;; DEF-OPERATOR makes "{" map to ${, but it needs to map to $SET. 53;; Just clobber whatever DEF-OPERATOR put into *OPR-TABLE*. 54(putopr "{" '$set) 55 56;; Support for TeXing sets. If your mactex doesn't TeX the empty set 57;; correctly, get the latest mactex.lisp. 58 59(defprop $set tex-matchfix tex) 60(defprop $set (("\\left \\{" ) " \\right \\}") texsym) 61 62(defun require-set (x context-string) 63 (if ($setp x) (cdr x) (merror (intl:gettext "~:M: argument must be a set; found: ~:M") context-string x))) 64 65;; If x is a Maxima list, return a Lisp list of its members; otherwise, 66;; signal an error. Unlike require-set, the function require-list does not 67;; coerce the result to a set. 68 69(defun require-list (x context-string) 70 (if ($listp x) (cdr x) 71 (merror (intl:gettext "~:M: argument must be a list; found: ~:M") context-string x))) 72 73;; If x is a Maxima list or set, return a Lisp list of its members; otherwise, 74;; signal an error. Unlike require-set, the function require-list-or-set 75;; does not coerce the result to a set. 76 77(defun require-list-or-set (x context-string) 78 (if (or ($listp x) ($setp x)) (cdr x) 79 (merror (intl:gettext "~:M: argument must be a list or a set; found: ~:M") context-string x))) 80 81;; When a is a list, return a list of the unique elements of a. 82;; Otherwise just return a. 83 84(defmfun $unique (x) 85 (if ($listp x) 86 `((mlist) ,@(sorted-remove-duplicates (sort (copy-list (cdr x)) '$orderlessp))) 87 x)) 88 89;; When a is a list, setify(a) is equivalent to apply(set, a). When a 90;; isn't a list, signal an error. 91 92(defmfun $setify (a) 93 (simplifya `(($set) ,@(require-list a '$setify)) nil)) 94 95;; When a is a list, convert a and all of its elements that are lists 96;; into sets. When a isn't a list, return a. 97 98(defmfun $fullsetify (a) 99 (cond (($listp a) 100 `(($set) ,@(mapcar '$fullsetify (cdr a)))) 101 (t a))) 102 103;; If a is a set, convert the top-level set to a list; when a isn't a 104;; list, return a. 105 106(defmfun $listify (a) 107 (if ($setp a) `((mlist simp) ,@(cdr a)) a)) 108 109;; full_listify(a) converts all sets in a into lists. 110 111(defmfun $full_listify (a) 112 (setq a ($ratdisrep a)) 113 (cond (($mapatom a) a) 114 (($setp a) (simplify (cons (list 'mlist) (mapcar #'$full_listify (cdr a))))) 115 (t (simplify (cons (car a) (mapcar #'$full_listify (cdr a))))))) 116 117(defprop $set simp-set operators) 118 119;; Simplify a set. 120 121(defun simp-set (a yy z) 122 (declare (ignore yy)) 123 (setq a (mapcar #'(lambda (x) (simplifya x z)) (cdr a))) 124 (setq a (sorted-remove-duplicates (sort a '$orderlessp))) 125 `(($set simp) ,@a)) 126 127;; Return true iff a is an empty set or list 128 129(defmfun $emptyp (a) 130 (or (like a `(($set))) (like a `((mlist))) (and ($matrixp a) (every '$emptyp (margs a))))) 131 132;; Return true iff the operator of a is set. 133 134(defmfun $setp (a) 135 (and (consp a) (consp (car a)) (eq (caar a) '$set))) 136 137;; Return the cardinality of a set. This function works even when $simp is false. 138 139(defmfun $cardinality (a) 140 (if $simp (length (require-set a '$cardinality)) 141 (let (($simp t)) ($cardinality (simplify a))))) 142 143;; Return true iff a is a subset of b. If either argument is a list, first 144;; convert it to a set. Signal an error if a or b aren't lists or sets. 145 146(defmfun $subsetp (a b) 147 (setq a (require-set a '$subsetp)) 148 (setq b (require-set b '$subsetp)) 149 (and (<= (length a) (length b)) (set-subsetp a b))) 150 151;; Return true iff sets a and b are equal; If either argument is a list, first 152;; convert convert it to a set. Signal an error if either a or b aren't lists 153;; or sets. 154 155(defmfun $setequalp (a b) 156 (setq a (require-set a '$setequalp)) 157 (setq b (require-set b '$setequalp)) 158 (and (= (length a) (length b)) (every #'like a b))) 159 160 161;; Adjoin x to the list or set a and return a set. 162 163(defmfun $adjoin (x a) 164 (setq a (require-set a '$adjoin)) 165 (multiple-value-bind (f i b) (b-search-expr x a 0 (length a)) 166 (if (not f) (setq a (prefixconc a i (cons x b)))) 167 `(($set simp) ,@a))) 168 169;; If x is a member of the set or list a, delete x from setify(a); otherwise, return 170;; setify(a). For a set a, disjoin(x,a) == delete(x,a) == setdifference(a,set(x)); 171;; however, disjoin should be the fastest way to delete a member from a set. 172 173(defmfun $disjoin (x a) 174 (setq a (require-set a '$disjoin)) 175 (multiple-value-bind (f i b) (b-search-expr x a 0 (length a)) 176 `(($set simp) ,@(if f (prefixconc a i b) a)))) 177 178;; (previxconc l len rest) is equivalent to (nconc (subseq l len) rest) 179 180(defun prefixconc (l len rest) 181 (do ((res nil (cons (car l) res)) 182 (i len (decf i)) 183 (l l (cdr l))) 184 ((= i 0) (nreconc res rest)) 185 (declare (fixnum i)))) 186 187;; union(a1,a2,...an) returns the union of the sets a1,a2,...,an. 188;; If any argument is a list, convert it to a set. Signal an error 189;; if one of the arguments isn't a list or a set. When union receives 190;; no arguments, it returns the empty set. 191 192(defmfun $union (&rest a) 193 (let ((acc nil)) 194 (dolist (ai a `(($set simp) ,@acc)) 195 (setq acc (set-union acc (require-set ai '$union)))))) 196 197;; Remove elements of b from a. Works on lists or sets. 198 199(defmfun $setdifference (a b) 200 `(($set simp) ,@(sset-difference (require-set a '$setdifference) 201 (require-set b '$setdifference)))) 202 203;; intersection(a1,a2,...an) returns the intersection of the sets 204;; a1,a2,...,an. Signal an error if one of the arguments isn't a 205;; list or a set. intersection must receive at least one argument. 206 207(defmfun $intersection (a &rest b) 208 (let ((acc (require-set a '$intersection))) 209 (cond ((consp b) 210 (dolist (bi b) 211 (setq acc (set-intersect acc (require-set bi '$intersection)))))) 212 `(($set simp) ,@acc))) 213 214;; intersect is an alias for intersection. 215 216(defmfun $intersect (a &rest b) 217 (apply '$intersection (cons a b))) 218 219;; Return true iff x as an element of the set or list a. Use like 220;; to test for equality. Signal an error if a isn't a set or list. 221 222(defmfun $elementp (x a) 223 (setq a (require-set a '$elementp)) 224 (b-search-expr x a 0 (length a))) 225 226;; Return true if and only if the lists or sets a and b are disjoint; 227;; signal an error if a or b aren't lists or sets. 228 229#| 230(defmfun $disjointp-binary-search-version (a b) 231 (setq a (require-set a '$disjointp)) 232 (setq b (require-set b '$disjointp)) 233 (if (> (length a) (length b)) (rotatef a b)) 234 (let ((n (length b))) 235 (catch 'disjoint 236 (dolist (ai a) 237 (if (b-search-expr ai b 0 n) (throw 'disjoint nil))) 238 t))) 239|# 240 241(defmfun $disjointp (a b) 242 (setq a (require-set a '$disjointp)) 243 (setq b (require-set b '$disjointp)) 244 (set-disjointp a b)) 245 246;; Return the set of elements of the list or set a for which the predicate 247;; f evaluates to true; signal an error if a isn't a list or a set. Also, 248;; signal an error if the function f doesn't evaluate to true, false, or 249;; unknown. 250 251(defmfun $subset (a f) 252 (setq a (require-set a '$subset)) 253 (let ((acc nil) (b)) 254 (dolist (x a `(($set simp) ,@(nreverse acc))) 255 (setq b (mfuncall f x)) 256 (cond ((eq t b) (push x acc)) 257 ((not (or (eq b nil) (eq b '$unknown))) 258 (merror (intl:gettext "subset: ~:M(~:M) evaluates to a non-boolean.") f x)))))) 259 260;; Return a list of three sets. The first set is the subset of a for which 261;; the predicate f evaluates to true, the second is the subset of a 262;; for which f evaluates to false, and the third is the subset of a 263;; for which f evaluates to unknown. 264 265(defmfun $partition_set (a f) 266 (setq a (require-set a '$partition_set)) 267 (let ((t-acc) (f-acc) (b)) 268 (dolist (x a `((mlist simp) 269 (($set simp) ,@(nreverse f-acc)) 270 (($set simp) ,@(nreverse t-acc)))) 271 (setq b (mfuncall f x)) 272 (cond ((eq t b) (push x t-acc)) 273 ((or (eq b nil) (eq b '$unknown)) (push x f-acc)) 274 (t 275 (merror (intl:gettext "partition_set: ~:M(~:M) evaluates to a non-boolean.") f x)))))) 276 277;; The symmetric difference of sets, that is (A-B) union (B - A), is associative. 278;; Thus the symmetric difference extends unambiguously to n-arguments. 279 280(defmfun $symmdifference (&rest l) 281 (let ((acc nil)) 282 (dolist (lk l (cons '($set simp) acc)) 283 (setq acc (set-symmetric-difference acc (require-set lk '$symmdifference)))))) 284 285;; Return {x | x in exactly one set l1, l2, ...} 286 287(defmfun $in_exactly_one (&rest l) 288 ;; u = union of l1, l2,... 289 ;; r = members that are in two or more l1, l2, ... 290 (let ((u nil) (r nil)) 291 (dolist (lk l) 292 (setq lk (require-set lk '$in_exactly_one)) 293 (setq r (set-union r (set-intersect u lk))) 294 (setq u (set-union u lk))) 295 (cons '($set simp) (sset-difference u r)))) 296 297;; When k is an integer, return the set of all subsets of the set a 298;; that have exactly k elements; when k is nil, return the power set 299;; of a. Signal an error if the first argument isn't a list or a set. 300 301(defmfun $powerset (a &optional k) 302 (setq a (require-set a "powerset")) 303 (cond ((null k) 304 (cons `($set simp) 305 (mapcar #'(lambda (s) 306 (cons '($set simp) s)) (power-set a)))) 307 ((integerp k) 308 (powerset-subset a k (length a))) 309 (t (merror (intl:gettext "The second argument to powerset must be an integer; found ~:M") k)))) 310 311 312(defun power-set (a) 313 (cond ((null a) `(())) 314 (t 315 (let ((x (car a)) (b (power-set (cdr a)))) 316 (append `(()) (mapcar #'(lambda (u) (cons x u)) b) (cdr b)))))) 317 318(defun powerset-subset (a k n) 319 (let ((s) (b) (acc)) 320 (cond ((= k 0) 321 (setq acc (cons `(($set)) acc))) 322 ((<= k n) 323 (dotimes (i k) 324 (setq s (cons i s))) 325 (setq s (nreverse s)) 326 (while (not (null s)) 327 (setq b nil) 328 (dotimes (i k) 329 (setq b (cons (nth (nth i s) a) b))) 330 (setq acc (cons (cons `($set simp) (nreverse b)) acc)) 331 (setq s (ksubset-lex-successor s k n))))) 332 (cons `($set simp) (nreverse acc)))) 333 334;; This code is based on Algorithm 2.6 "Combinatorial Algorithms Generation, 335;; Enumeration, and Search," CRC Press, 1999 by Donald Kreher and Douglas 336;; Stinson. 337 338(defun ksubset-lex-successor (s k n) 339 (let ((u (copy-list s)) 340 (i (- k 1)) (j) (si (- n k))) 341 (while (and (>= i 0) (= (nth i s) (+ si i))) 342 (decf i)) 343 (cond ((< i 0) 344 nil) 345 (t 346 (setq j i) 347 (setq si (+ 1 (- (nth i s) i))) 348 (while (< j k) 349 (setf (nth j u) (+ si j)) 350 (incf j)) 351 u)))) 352 353;; When the list a is redundant, need-to-simp is set to true; this flag 354;; determines if acc needs to be simplified. Initially, p = (0,1,2,..,n); 355;; the 356 357(defmfun $permutations (a) 358 (cond (($listp a) 359 (setq a (sort (copy-list (cdr a)) '$orderlessp))) 360 (t 361 (setq a (require-set a '$permutations)))) 362 363 (let* ((n (length a)) (p (make-array (+ n 1) :element-type 'fixnum)) 364 (r (make-array (+ n 1) :initial-element 0 :element-type 'fixnum)) 365 (b (make-array (+ n 1) :initial-element 0)) 366 (i) (acc) (q) 367 (need-to-simp (not (= (length a) 368 (length (sorted-remove-duplicates (copy-list a))))))) 369 370 (dotimes (i (+ 1 n)) 371 (setf (aref p i) i)) 372 (dotimes (i n) 373 (setf (aref b (+ i 1)) (nth i a))) 374 375 (cond ((not (null a)) 376 (while (not (null p)) 377 (setq i 1) 378 (setq q nil) 379 (while (<= i n) 380 (setq q (cons (aref b (aref p i)) q)) 381 (incf i)) 382 (setq acc (cons (cons '(mlist simp) (nreverse q)) acc)) 383 (setq p (permutation-lex-successor n p r)))) 384 (t 385 (setq acc `(((mlist)))))) 386 (setq acc (nreverse acc)) 387 (if need-to-simp `(($set) ,@acc) 388 `(($set simp) ,@acc)))) 389 390;; This code is based on Algorithm 2.14 "Combinatorial Algorithms Generation, 391;; Enumeration, and Search," CRC Press, 1999 by Donald Kreher and Douglas 392;; Stinson. 393 394;; The array elements p(1) thru p(n) specify the permutation; the array 395;; r gets used for swapping elements of p. Initially p = (0,1,2,..,n). 396 397(defun permutation-lex-successor (n p r) 398 (declare (type (simple-array fixnum *) p r)) 399 (declare (type fixnum n)) 400 (let ((i (- n 1)) (j n) (m) (tm)) 401 (setf (aref p 0) 0) 402 (while (< (aref p (+ i 1)) (aref p i)) 403 (decf i)) 404 (cond ((= i 0) 405 nil) 406 (t 407 (while (< (aref p j) (aref p i)) 408 (decf j)) 409 (setq tm (aref p j)) 410 (setf (aref p j) (aref p i)) 411 (setf (aref p i) tm) 412 (setq j (+ i 1)) 413 (while (<= j n) 414 (setf (aref r j) (aref p j)) 415 (incf j)) 416 (setq j (+ i 1)) 417 (setq m (+ n i 1)) 418 (while (<= j n) 419 (setf (aref p j) (aref r (- m j))) 420 (incf j)) 421 p)))) 422 423(defmfun $random_permutation (a) 424 (if ($listp a) 425 (setq a (copy-list (cdr a))) 426 (setq a (copy-list (require-set a '$random_permutation)))) 427 428 (let ((n (length a))) 429 (dotimes (i n) 430 (let 431 ((j (+ i ($random (- n i)))) 432 (tmp (nth i a))) 433 (setf (nth i a) (nth j a)) 434 (setf (nth j a) tmp)))) 435 436 `((mlist) ,@a)) 437 438 439#| 440;;; Returns 3 values 441;;; FOUND -- is X in L 442;;; POSITION -- where is X in L; if not in L, position it is before 443;;; REST -- everything after X in L 444 445(defun old-b-search-expr (x l lo len) 446 (declare (fixnum lo len)) 447 (if (= len 0) (values nil lo l) 448 (let ((mid) (midl)) 449 (while (> len 1) 450 (if ($orderlessp x (car (setq midl (nthcdr (setq mid (floor len 2)) l)))) 451 (setq len mid) 452 (setq l midl 453 lo (+ lo mid) 454 len (- len mid)))) 455 (cond (($orderlessp x (nth 0 l)) (values nil lo l)) 456 ((like x (nth 0 l)) (values t lo (cdr l))) 457 (t (values nil (1+ lo) (cdr l))))))) 458|# 459 460;;; Returns 3 values 461;;; FOUND -- is X in L 462;;; POSITION -- where is X in L; if not in L, position it is before 463;;; REST -- everything after X in L 464 465(defun b-search-expr (x l lo len) 466 (declare (fixnum lo len)) 467 (if (= len 0) (values nil lo l) 468 (progn 469 ;; uses great directly instead of $orderlessp; only specrepcheck x once 470 (setq x (specrepcheck x)) 471 (let ((mid) (midl) (midel)) 472 (while (> len 1) 473 (cond 474 ;; Previously, it could hit x and continue searching 475 ;; Since great doesn't guarantee inequality, we need the check for 476 ;; alike1 anyway (hidden inside $orderlessp), so we might as well 477 ;; take advantage of it 478 ((alike1 479 x 480 (setq midel 481 (specrepcheck (car (setq midl (nthcdr (setq mid (floor 482 len 2)) l)))))) 483 (setq l midl 484 lo (+ lo mid) 485 len -1)) 486 487 ((great midel x) 488 (setq len mid)) 489 (t (setq l midl 490 lo (+ lo mid) 491 len (- len mid))))) 492 493 (cond ((= len -1) (values t lo (cdr l))) 494 ((alike1 x (specrepcheck (nth 0 l))) (values t lo (cdr l))) 495 ((great (specrepcheck (nth 0 l)) x) (values nil lo l)) 496 (t (values nil (1+ lo) (cdr l)))))))) 497 498;; Flatten is somewhat difficult to define -- essentially it evaluates an 499;; expression as if its main operator had been declared nary; however, there 500;; is a difference. We have 501 502;; (C2) flatten(f(g(f(f(x))))); 503;; (D2) f(g(f(f(x)))) 504;; (C3) declare(f,nary); 505;; (D3) DONE 506;; (C4) ev(d2); 507;; (D4) f(g(f(x))) 508 509;; Unlike declaring the main operator of an expression to be nary, flatten 510;; doesn't recurse into other function arguments. 511 512;; To successfully flatten an expression, the main operator must be 513;; defined for zero or more arguments; if this isn't the case, 514;; Maxima can halt with an error. So be it. 515 516(defmfun $flatten (e) 517 (cond ((or (specrepp e) (mapatom e)) e) 518 (t (mcons-op-args (mop e) (flattenl-op (margs e) (mop e)))))) 519 520(defun flattenl-op (e op) 521 (mapcan #'(lambda (e) 522 (cond ((or (mapatom e) (not (alike1 (mop e) op))) 523 (list e)) 524 (t (flattenl-op (margs e) op)))) 525 e)) 526 527; doesn't work on f[1](f[1](x)). 528;(defmfun $flatten (e) 529; (if (or (specrepp e) (mapatom e)) e 530; (cons `(,(mop e)) (total-nary e)))) 531 532(defun sorted-remove-duplicates (l) 533 (prog1 l 534 (while (cdr l) 535 (while (and (cdr l) (like (car l) (cadr l)) 536 (rplacd l (cddr l)))) 537 (setq l (cdr l))))) 538 539(defun set-intersect (l1 l2) 540 ;; Only works for lists of sorted by $orderlessp. 541 (with-collector collect 542 (do-merge-symm 543 l1 l2 544 #'like 545 #'$orderlessp 546 #'collect 547 nil))) 548 549(defun set-union (l1 l2) 550 ;; Only works for lists of sorted by $orderlessp. 551 (with-collector collect 552 (do-merge-symm 553 l1 l2 554 #'like 555 #'$orderlessp 556 #'collect 557 #'collect))) 558 559(defun sset-difference (l1 l2) 560 ;; Only works for lists of sorted by $orderlessp. 561 (with-collector collect 562 (do-merge-asym 563 l1 l2 564 #'like 565 #'$orderlessp 566 nil 567 #'collect 568 nil))) 569 570(defun set-subsetp (l1 l2) 571 ;; Is l1 a subset of l2 572 (catch 'subset 573 (do-merge-asym 574 l1 l2 575 #'like 576 #'$orderlessp 577 nil 578 #'(lambda (xx) (declare (ignore xx)) (throw 'subset nil)) 579 nil) 580 t)) 581 582(defun set-symmetric-difference (l1 l2) ; i.e. xor 583 (with-collector collect 584 (do-merge-symm 585 l1 l2 586 #'like 587 #'$orderlessp 588 nil 589 #'collect))) 590 591(defun set-disjointp (l1 l2) 592 (catch 'disjoint 593 (do-merge-symm 594 l1 l2 595 #'like 596 #'$orderlessp 597 #'(lambda (xx) (declare (ignore xx)) (throw 'disjoint nil)) 598 nil) 599 t)) 600 601;; When s = $max, return { x in a | f(x) = maximum of f on a} and 602;; when s = $min, return { x in a | f(x) = minimum of f on a}. 603;; Signal an error when s isn't $max or $min. 604 605(defmfun $extremal_subset (a f s) 606 (setq a (require-set a '$extremal_subset)) 607 (cond ((null a) 608 `(($set simp))) 609 (t 610 (cond ((eq s '$min) 611 (setq s -1)) 612 ((eq s '$max) 613 (setq s 1)) 614 (t 615 (merror (intl:gettext "extremal_subset: third argument must be 'max or 'min; found: ~:M") s))) 616 (let* ((max-subset (nth 0 a)) 617 (mx (mul s (mfuncall f max-subset))) 618 (x)) 619 (setq max-subset `(,max-subset)) 620 (setq a (cdr a)) 621 (dolist (ai a) 622 (setq x (mul s (mfuncall f ai))) 623 (cond ((mevalp_tr (mgrp x mx) t nil) 624 (setq mx x 625 max-subset `(,ai))) 626 ((like x mx) 627 (setq max-subset (cons ai max-subset))))) 628 `(($set simp) ,@(nreverse max-subset)))))) 629 630(defun bool-checked-mfuncall (f x y) 631; (let ((bool (is-boole-check (mfuncall f x y)))) 632; (if (not (or (eq bool 't) (eq bool nil))) 633; (merror "~:M(~:M,~:M) doesn't evaluate to a boolean" f x y) 634; bool))) 635 (let (($prederror nil) (b)) 636 (setq b (mevalp (mfuncall f x y))) 637 (if (or (eq b t) (eq b nil)) b 638 (merror (intl:gettext "equiv_classes: ~:M(~:M, ~:M) evaluates to a non-boolean.") f x y)))) 639 640 641;; Return the set of equivalence classes of f on the set l. The 642;; function f must be an boolean-valued function defined on the 643;; cartesian product of l with l; additionally, f should be an 644;; equivalence relation. 645 646;; The lists acc and tail share structure. 647 648(defmfun $equiv_classes (l f) 649 (setq l (require-set l '$equiv_classes)) 650 (do ((l l (cdr l)) 651 (acc) 652 (tail) 653 (x)) 654 ((null l) (simplify (cons '($set) (mapcar #'(lambda (x) (cons '($set) x)) acc)))) 655 (setq x (car l)) 656 (setq tail (member-if #'(lambda (z) (bool-checked-mfuncall f x (car z))) acc)) 657 (cond ((null tail) 658 (setq acc (cons `(,x) acc))) 659 (t 660 (setf (car tail) (cons x (car tail))))))) 661 662;; cartesian_product(a,b1,b2,...,bn), where a, b1, ..., bn are all sets, 663;; returns the set with members of the form [x0,x1, ..., xn], 664;; where x0 in a, x1 in b1, ... , and xn in bn. 665;; With just one argument cartesian_product(a) returns the 666;; set with members [a1],[a2], ... [an], where a1, ..., an are the members of a. 667;; With no arguments, cartesian_product() returns {[]}. 668 669(defmfun $cartesian_product (&rest b) 670 (if (null b) 671 '(($set) ((mlist))) 672 (if (every #'$setp b) 673 (let ((l (apply #'cartesian-product (mapcar #'cdr b)))) 674 (cons '($set) (mapcar #'(lambda (e) (cons '(mlist) e)) l))) 675 ;; MAYBE JUST PRINT THE LIST OF TYPES OR OPERATORS INSTEAD OF B IN ITS ENTIRETY !! 676 (merror (intl:gettext "cartesian_product: all arguments must be sets; found: ~M") (cons '(mlist) b))))) 677 678;; cartesian_product_list(a,b1,b2,...,bn), where a, b1, ..., bn are all lists, 679;; returns the list with elements of the form [x0,x1, ..., xn], 680;; where x0 in a, x1 in b1, ... , and xn in bn. 681;; With just one argument cartesian_product_list(a) returns the 682;; list with elements [a1],[a2], ... [an], where a1, ..., an are the elements of a. 683;; With no arguments, cartesian_product_list() returns [[]]. 684 685(defmfun $cartesian_product_list (&rest b) 686 (if (null b) 687 '((mlist) ((mlist))) 688 (if (every #'$listp b) 689 (let ((l (apply #'cartesian-product (mapcar #'cdr b)))) 690 (cons '(mlist) (mapcar #'(lambda (e) (cons '(mlist) e)) l))) 691 ;; MAYBE JUST PRINT THE LIST OF TYPES OR OPERATORS INSTEAD OF B IN ITS ENTIRETY !! 692 (merror (intl:gettext "cartesian_product_list: all arguments must be lists; found: ~M") (cons '(mlist) b))))) 693 694;; Assume here that B is nonempty; caller has already handled case B = NIL. 695(defun cartesian-product (&rest b) 696 (setq b (reverse b)) 697 (let 698 ((a) 699 (acc (mapcar #'list (car b)))) 700 (setq b (cdr b)) 701 (dolist (bi b) 702 (setq a nil) 703 (dolist (bij bi (setq acc a)) 704 (setq a (append a (mapcar #'(lambda (x) (cons bij x)) acc))))) 705 acc)) 706 707;; When n is defined, return a set of partitions of the set or list a 708;; into n disjoint subsets. When n isn't defined, return the set of 709;; all partitions. 710 711;; Let S be a set. We say a set P is a partition of S provided 712;; (1) p in P implies p is a set, 713;; (2) p1, p2 in P and p1 # p2 implies p1 and p2 are disjoint, 714;; (3) union(x | x in P) = S. 715;; Thus set() is a partition of set(). 716 717(defmfun $set_partitions (a &optional n-sub) 718 (setq a (require-set a '$set_partitions)) 719 (cond ((and (integerp n-sub) (> n-sub -1)) 720 `(($set) ,@(set-partitions a n-sub))) 721 ((null n-sub) 722 (setq n-sub (length a)) 723 (let ((acc (set-partitions a 0)) (k 1)) 724 (while (<= k n-sub) 725 (setq acc (append acc (set-partitions a k))) 726 (incf k)) 727 `(($set) ,@acc))) 728 (t 729 (merror (intl:gettext "set_partitions: second argument must be a positive integer; found: ~:M") n-sub)))) 730 731(defun set-partitions (a n) 732 (cond ((= n 0) 733 (cond ((null a) 734 (list `(($set)))) 735 (t 736 nil))) 737 ((null a) 738 nil) 739 (t 740 (let ((p) (x) (acc) (w) (s) (z)) 741 (setq x (car a)) 742 (setq p (set-partitions (cdr a) n)) 743 (dolist (pj p) 744 (setq w nil) 745 (setq s (cdr pj)) 746 (while (not (null s)) 747 (setq z (pop s)) 748 (setq acc (cons (simplifya `(($set) ,@w ,($adjoin x z) ,@s) t) acc)) 749 (setq w (cons z w)))) 750 751 (setq x `(($set) ,x)) 752 (setq p (set-partitions (cdr a) (- n 1))) 753 (dolist (pj p acc) 754 (setq acc (cons ($adjoin x pj) acc))))))) 755 756;; Generate the integer partitions in dictionary order. When the optional 757;; argument len is defined, only generate the partitions with exactly len 758;; members, including 0. 759 760(defmfun $integer_partitions (n &optional len) 761 (let ((acc)) 762 (cond ((and (integerp n) (>= n 0)) 763 (setq acc (cond ((= n 0) nil) 764 ((integerp len) (fixed-length-partitions n n len)) 765 (t (integer-partitions n)))) 766 (if (not acc) 767 (setq acc `(((mlist simp)))) 768 (setq acc (mapcar #'(lambda (x) (simplify (cons '(mlist) x))) acc))) 769 `(($set simp) ,@acc)) 770 (t 771 (if len `(($integer_partitions simp) ,n ,len) `(($integer_partitions simp) ,n)))))) 772 773(defun integer-partitions (n) 774 (let ((p `(,n)) (acc nil) (d) (k) (j) (r)) 775 (while (> (car (last p)) 1) 776 (setq acc (cons (copy-list (reverse p)) acc)) 777 (setq p (member t p :key #'(lambda (x) (> x 1)))) 778 (setq k (- (nth 0 p) 1)) 779 (setf (nth 0 p) k) 780 (setq d (- n (reduce #'+ p))) 781 (setq j k) 782 (while (and (> k 0) (> d 0)) 783 (multiple-value-setq (d r) (floor d k)) 784 (setq p (append (make-list d :initial-element k) p)) 785 (setq d r) 786 (decf k))) 787 (setq acc (cons (copy-list (reverse p)) acc)) 788 acc)) 789 790(defun fixed-length-partitions (n b len) 791 (let ((p t) (acc) (i)) 792 (cond ((> n (* b len)) nil) 793 ((= len 1) (if (<= n b) (setq acc `((,n))) nil)) 794 (t 795 (setq len (- len 1)) 796 (setq i (- n (min n b))) 797 (setq n (min n b)) 798 (while (not (null p)) 799 (setq p (mapcar #'(lambda (x) (cons n x)) (fixed-length-partitions i (min i n) len))) 800 (setq acc (append p acc)) 801 (decf n) 802 (incf i)))) 803 acc)) 804 805;; When n is a nonnegative integer, return the number of partitions of n. 806;; If the optional parameter lst has the value "list", return a list of 807;; the number of partitions of 1,2,3, ... , n. If n isn't a nonnegative 808;; integer, return a noun form. 809 810(defmfun $num_partitions (n &optional lst) 811 (cond ((equal n 0) 1) 812 ((and (integerp n) (> n -1)) 813 (let ((p (make-array (+ n 1))) 814 (s (make-array (+ n 1))) 815 (sum) (i) (j)) 816 (setf (aref p 0) 1) 817 (setf (aref p 1) 1) 818 819 (setq i 0) 820 (while (<= i n) 821 (setf (aref s i) (mfuncall '$divsum i 1)) 822 (incf i)) 823 824 (setq i 2) 825 (while (<= i n) 826 (setq sum 0) 827 (setq j 1) 828 (while (<= j i) 829 (setq sum (+ sum (* (aref s j) (aref p (- i j))))) 830 (incf j)) 831 (setf (aref p i) (/ sum i)) 832 (incf i)) 833 (cond ((eq lst '$list) 834 (let ((acc)) 835 (incf n) 836 (dotimes (i n) 837 (push (aref p i) acc)) 838 (setq acc (nreverse acc)) 839 (push '(mlist simp) (cdr acc)))) 840 (t 841 (aref p n))))) 842 (t (if lst `(($num_partitions simp) ,n ,lst) 843 `(($num_partitions simp) ,n))))) 844 845(defmfun $num_distinct_partitions (n &optional lst) 846 (cond ((eql n 0) 1) 847 ((and (integerp n) (> n -1)) 848 (let ((p (make-array (+ n 1))) 849 (s (make-array (+ n 1))) 850 (u (make-array (+ n 1))) 851 (sum) (i) (j)) 852 (setf (aref p 0) 1) 853 (setf (aref p 1) 1) 854 855 (setq i 0) 856 (while (<= i n) 857 (setf (aref s i) (mfuncall '$divsum i 1)) 858 (incf i)) 859 (setq i 0) 860 (while (<= i n) 861 (if (oddp i) 862 (setf (aref u i) (aref s i)) 863 (setf (aref u i) (- (aref s i) (* 2 (aref s (/ i 2)))))) 864 (incf i)) 865 (setq i 2) 866 (while (<= i n) 867 (setq sum 0) 868 (setq j 1) 869 (while (<= j i) 870 (setq sum (+ sum (* (aref u j) (aref p (- i j))))) 871 (incf j)) 872 (setf (aref p i) (/ sum i)) 873 (incf i)) 874 875 (cond ((eq lst '$list) 876 (let ((acc)) 877 (incf n) 878 (dotimes (i n) 879 (push (aref p i) acc)) 880 (setq acc (nreverse acc)) 881 (push '(mlist simp) (cdr acc)))) 882 (t 883 (aref p n))))) 884 (t (if lst `(($num_distinct_partitions simp) ,n ,lst) 885 `(($num_distinct_partitions simp) ,n))))) 886 887;; A n-ary Kronecker delta function: kron_delta(n0,n1, ..., nk) simplifies to 1 if 888;; (meqp ni nj) is true for *all* pairs ni, nj in (n0,n1, ..., nk); it simplifies to 0 if 889;; (mnqp ni nj) is true for *some* pair ni, nj in (n0,n1, ..., nk). Further kron_delta() --> 1 890;; and kron_delta(xxx) --> wrong number of arguments error. Thus 891;; 892;; kron_delta(x0,...,xn) * kron_delta(y0,..., ym) = kron_delta(x0, ..., xn, y0, ..., ym) 893;; 894;; is an identity. 895 896(defprop %kron_delta simp-kron-delta operators) 897(setf (get '$kron_delta 'noun) '%kron_delta) 898(setf (get '%kron_delta 'verb) '$kron_delta) 899(setf (get '$kron_delta 'alias) '%kron_delta) 900(setf (get '%kron_delta 'reversealias) '$kron_delta) 901(defmfun $kron_delta (&rest x) (simplifya `((%kron_delta) ,@x) t)) 902(setf (get '%kron_delta 'real-valued) t) ;; conjugate(kron_delta(xxx)) --> kron_delta(xxx) 903(setf (get '%kron_delta 'integer-valued) t) ;; featurep(kron_delta(xxx), integer) --> true 904(mputprop '%kron_delta t '$scalar) ;; same effect as declare(kron_delta, scalar) 905 906(putprop '%kron_delta #'(lambda (yy) (declare (ignore yy)) (setq sign '$pz)) 'sign-function) 907 908(defun simp-kron-delta (l yy z) 909 (declare (ignore yy)) 910 911 (setq l (cdr l)) ;; remove (($kron_delta simp) 912 (if (and l (null (cdr l))) (wna-err '$kron_delta)) ;; wrong number of arguments error for exactly one argument 913 914 ;; Checking both mnqp and meqp is convenient, but unnecessary. This code misses simplifications that 915 ;; involve three or more arguments. Example: kron_delta(a,b,a+b+1,a-b+5) could (but doesn't) simplify 916 ;; to 0 (the solution set (a = b, a = a+b+1, a=a-b+5) is empty. 917 918 (let ((acc nil) (return-zero nil)) 919 (setq return-zero (catch 'done 920 (dolist (lk l) 921 (setq lk (simpcheck lk z)) 922 (cond ((some #'(lambda (s) (eq t (mnqp s lk))) acc) ;; lk # some member of acc, return zero. 923 (throw 'done t)) 924 ((some #'(lambda (s) (eq t (meqp s lk))) acc)) ;; lk = some member of acc, do nothing 925 (t (push lk acc))));; push lk onto acc 926 nil)) ;; set return-zero to nil 927 (cond (return-zero 0) 928 ((or (null acc) (null (cdr acc))) 1) 929 (t ;; reflection: kron_delta(-a,-b,...) == kron_delta(a,b,...). 930 (let ((neg-acc (sort (mapcar #'neg acc) '$orderlessp))) 931 (setq acc (sort acc '$orderlessp)) 932 `((%kron_delta simp) ,@(if (great (cons '(mlist) neg-acc) (cons '(mlist) acc)) neg-acc acc))))))) 933 934(defprop %kron_delta tex-kron-delta tex) 935 936(defun tex-kron-delta (x l r) 937 (append l `("\\delta_{" ,@(tex-list (cdr x) nil (list "} ") ", ")) r)) 938 939;; Stirling numbers of the first kind. 940 941(defprop $stirling1 simp-stirling1 operators) 942;; Stirling numbers of the first kind. 943 944(defprop $stirling1 simp-stirling1 operators) 945 946;; Stirling1 simplifications; for n,k in Z (Z = set of integers) 947;; (1) stirling1(1,k) = kron_delta(1,k), k >= 0, (http://dlmf.nist.gov/26.8.E2) 948;; (2) stirling1(n,n) = 1, n >= 0 (http://dlmf.nist.gov/26.8.E1) 949;; (3) stirling1(n,n-1) = -binomial(n,2), n >= 1, (http://dlmf.nist.gov/26.8.E16) 950;; (4) stirling1(n,0) = kron_delta(n,0), n >=0 (http://dlmf.nist.gov/26.8.E14 and http://dlmf.nist.gov/26.8.E1) 951;; (5) stirling1(n,1) =(-1)^(n-1) (n-1)!, n >= 1 (http://dlmf.nist.gov/26.8.E14) 952;; (6) stirling1(n,k) = 0, n >= 0 and k > n. 953 954(defun simp-stirling1 (l yy z) 955 (declare (ignore yy)) 956 (let* ((fn (car (pop l))) 957 (n (if l (simplifya (pop l) z) (wna-err fn))) 958 (k (if l (simplifya (pop l) z) (wna-err fn))) 959 (n-is-nonnegative-int (nonnegative-integerp n))) 960 (if l (wna-err fn)) 961 (cond ((and (integerp n) (integerp k) (> n -1) (> k -1)) 962 (integer-stirling1 n k)) 963 ((and (eql n 1) ($featurep k '$integer) (eq t (mgrp k -1))) 964 (take (list '%kron_delta) 1 k)) 965 ((and n-is-nonnegative-int (like n k)) 966 1) 967 ((and (nonnegative-integerp (sub n 1)) (like n (add k 1))) 968 (mul -1 (take (list '%binomial) n 2))) 969 ((and n-is-nonnegative-int (eql k 0)) 970 (take (list '%kron_delta) n 0)) 971 ((and n-is-nonnegative-int (eql k 1) (eq t (mgqp n 1))) 972 (mul (power -1 (sub n 1)) (take (list 'mfactorial) (sub n 1)))) 973 ((and n-is-nonnegative-int ($featurep k '$integer) (eq t (mgrp k n))) 974 0) 975 (t (list (list fn 'simp) n k))))) 976 977;; This code is based on Algorithm 3.17 "Combinatorial Algorithms Generation, 978;; Enumeration, and Search," CRC Press, 1999 by Donald Kreher and Douglas 979;; Stinson. There is a typographical error in Algorithm 3.17; replace i - j 980;; with i - 1. See Theorem 3.14. 981 982(defun integer-stirling1 (m n) 983 (cond ((>= m n) 984 (let ((s (make-array `(,(+ m 1) ,(+ m 1)) :initial-element 0)) 985 (i) (j) (k) (im1)) 986 (setf (aref s 0 0) 1) 987 (setq i 1) 988 (while (<= i m) 989 (setq k (min i n)) 990 (setq j 1) 991 (setq im1 (- i 1)) 992 (while (<= j k) 993 (setf (aref s i j) (- (aref s im1 (- j 1)) 994 (* im1 (aref s im1 j)))) 995 (incf j)) 996 (incf i)) 997 (aref s m n))) 998 (t 0))) 999 1000 1001;; Stirling1 simplifications; for n,k in Z (Z = set of integers) 1002;; (1) stirling2(n,0) = 1, n >= 1 (http://dlmf.nist.gov/26.8.E17 and stirling2(0,0) = 1) 1003;; (2) stirling2(n,n) = 1, n >= 0, (http://dlmf.nist.gov/26.8.E4) 1004;; (3) stirling2(n,1) = 1, n >= 1, (http://dlmf.nist.gov/26.8.E17 and stirling2(0,1) = 0) 1005;; (4) stirling2(n,2) = 2^(n-1) -1 , n >= 1, (http://dlmf.nist.gov/26.8.E17) 1006;; (5) stirling2(n,n-1) = binomial(n,2), n>= 1 (http://dlmf.nist.gov/26.8.E16) 1007;; (6) stirling2(n,k) = 0, n >= 0 and k > n. 1008 1009(defun nonnegative-integerp (e) 1010 (and ($featurep e '$integer) 1011 (member ($sign (specrepcheck e)) `($pos $zero $pz) :test #'eq))) 1012 1013(defprop $stirling2 simp-stirling2 operators) 1014 1015(defun simp-stirling2 (l yy z) 1016 (declare (ignore yy)) 1017 (let* ((fn (car (pop l))) 1018 (n (if l (simplifya (pop l) z) (wna-err fn))) 1019 (k (if l (simplifya (pop l) z) (wna-err fn))) 1020 (n-is-nonnegative-int (nonnegative-integerp n)) 1021 (n-is-positive-int (nonnegative-integerp (sub n 1)))) 1022 (if l (wna-err fn)) 1023 (cond ((and (integerp n) (integerp k) (> n -1) (> k -1)) 1024 (integer-stirling2 n k)) 1025 1026 ((and n-is-positive-int (eql k 0)) 1027 0) 1028 1029 ((and n-is-positive-int (eql k 1)) 1030 1) 1031 1032 ((and n-is-positive-int (eql k 2)) 1033 (add (power 2 (sub n 1)) -1)) 1034 1035 ((and n-is-nonnegative-int (like n k)) 1036 1) 1037 1038 ((and n-is-positive-int (like n (add k 1))) 1039 (take (list '%binomial) n 2)) 1040 1041 ((and n-is-nonnegative-int ($featurep k '$integer) (eq t (mgrp k n))) 1042 0) 1043 1044 (t (list (list fn 'simp) n k))))) 1045 1046;; Stirling2(n,m) = sum((-1)^(m - k) binomial(m k) k^n,i,1,m) / m!. 1047;; See A & S 24.1.4, page 824. 1048 1049(defun integer-stirling2 (n m) 1050 (let ((s (if (= n 0) 1 0)) (i 1) (z) (f 1) (b m)) 1051 (while (<= i m) 1052 (setq z (* b (expt i n))) 1053 (setq f (* f i)) 1054 (setq b (/ (* (- m i) b) (+ i 1))) 1055 (if (oddp i) (setq z (- z))) 1056 (setq s (+ s z)) 1057 (incf i)) 1058 (setq s (/ s f)) 1059 (if (oddp m) (- s) s))) 1060 1061;; Return the Bell number of n; specifically, belln(n) is the 1062;; cardinality of the set of partitions of a set with n elements. 1063 1064(defprop $belln simp-belln operators) 1065 1066;; Simplify the Bell function. Other than evaluation for nonnegative 1067;; integer arguments, there isn't much that can be done. I don't know 1068;; a reasonable extension of the Bell function to non-integers or of 1069;; any simplifications -- we do thread belln over lists, sets, matrices, 1070;; and equalities. 1071 1072(defun simp-belln (n y z) 1073 (oneargcheck n) 1074 (setq y (caar n)) 1075 (setq n (simpcheck (cadr n) z)) 1076 (cond ((and (integerp n) (> n -1)) 1077 (integer-belln n)) 1078 ((or ($listp n) ($setp n) ($matrixp n) (mequalp n)) 1079 (thread y (cdr n) (caar n))) 1080 (t `(($belln simp) ,n)))) 1081 1082(defun integer-belln (n) 1083 (let ((s (if (= n 0) 1 0)) (i 1)) 1084 (while (<= i n) 1085 (setq s (+ s (integer-stirling2 n i))) 1086 (incf i)) 1087 s)) 1088 1089;; The multinomial coefficient; explicitly multinomial_coeff(a1,a2, ... an) = 1090;; (a1 + a2 + ... + an)! / (a1! a2! ... an!). The multinomial coefficient 1091;; gives the number of ways of placing a1 + a2 + ... + an distinct objects 1092;; into n boxes with ak elements in the k-th box. 1093 1094;; multinomial_coeff is symmetric; thus when at least one of its arguments 1095;; is symbolic, we sort them. Additionally any zero element of the 1096;; argument list can be removed without changing the value of 1097;; multinomial_coeff; we make this simplification as well. If 1098;; b is nil following (remove 0 b), something has gone wrong. 1099 1100(defmfun $multinomial_coeff (&rest a) 1101 (let ((n 0) (d 1)) 1102 (dolist (ai a) 1103 (setq n (add n ai)) 1104 (setq d (mult d (simplify `((mfactorial) ,ai))))) 1105 (div (simplify `((mfactorial) ,n)) d))) 1106 1107;; Extend a function f : S x S -> S to n arguments using right associativity. 1108;; Thus rreduce(f,[0,1,2]) -> f(0,f(1,2)). The second argument must be a list. 1109 1110(defmfun $rreduce (f s &optional (init 'no-init)) 1111 (rl-reduce f s t init '$rreduce)) 1112 1113;; Extend a function f : S x S -> S to n arguments using left associativity. 1114;; Thus lreduce(f,[0,1,2]) -> f(f(0,1),2). The second argument must be a list. 1115 1116(defmfun $lreduce (f s &optional (init 'no-init)) 1117 (rl-reduce f s nil init '$lreduce)) 1118 1119(defun rl-reduce (f s left init fn) 1120 (setq s (require-list s fn)) 1121 (cond ((not (equal init 'no-init)) 1122 (reduce #'(lambda (x y) (mfuncall f x y)) s :from-end left 1123 :initial-value init)) 1124 ((null s) 1125 (merror (intl:gettext "~a: either a nonempty set or initial value must be given.") fn)) 1126 (t 1127 (reduce #'(lambda (x y) (mfuncall f x y)) s :from-end left)))) 1128 1129;; Define an operator (signature S x S -> S, for some set S) to be nary and 1130;; define a function for its n-argument reduction. There isn't a user-level 1131;; interface to this mechanism. 1132 1133(defmacro def-nary (fn arg f-body id) 1134 `(setf (get ,fn '$nary) (list #'(lambda ,arg ,f-body) ,id))) 1135 1136(defun xappend (s) 1137 #+(or cmu scl) 1138 (cons '(mlist) (apply 'append (mapcar #'(lambda (x) 1139 (require-list x '$append)) s))) 1140 #-(or cmu scl) 1141 (let ((acc)) 1142 (dolist (si (reverse s) (cons '(mlist) acc)) 1143 (setq acc (append (require-list si '$append) acc))))) 1144 1145(def-nary 'mand (s) (mevalp (cons '(mand) s)) t) 1146(def-nary 'mor (s) (mevalp (cons '(mor) s)) nil) 1147(def-nary 'mplus (s) (simplify (cons '(mplus) s)) 0) 1148(def-nary 'mtimes (s) (simplify (cons '(mtimes) s)) 1) 1149(def-nary '$max (s) (if (null s) '$minf (maximin s '$max)) '$minf) 1150(def-nary '$min (s) (if (null s) '$inf (maximin s '$min)) '$inf) 1151(def-nary '$append (s) (xappend s) '((mlist))) 1152(def-nary '$union (s) ($apply '$union (cons '(mlist) s)) '(($set))) 1153 1154;; Extend a function f : S x S -> S to n arguments. When we 1155;; recognize f as a nary function (associative), if possible we call a Maxima 1156;; function that does the work efficiently -- examples are "+", "min", and "max". 1157;; When there isn't a Maxima function we can call (actually when (get op '$nary) 1158;; returns nil) we give up and use rl-reduce with left-associativity. 1159 1160 1161(defmfun $xreduce (f s &optional (init 'no-init)) 1162 (let* ((op-props (get (if (atom f) ($verbify f) nil) '$nary)) 1163 (opfn (if (consp op-props) (car op-props) nil))) 1164 1165 (cond (opfn 1166 (setq s (require-list-or-set s '$xreduce)) 1167 (if (not (equal init 'no-init)) 1168 (setq s (cons init s))) 1169 1170 (if (null s) 1171 (cadr op-props) ; is this clause really needed? 1172 1173 (funcall opfn s))) 1174 1175 (op-props 1176 ($apply f ($listify s))) 1177 1178 (t 1179 (rl-reduce f ($listify s) nil init '$xreduce))))) 1180 1181 1182;; Extend a function f : S x S -> S to n arguments using a minimum depth tree. 1183;; The function f should be nary (associative); otherwise, the result is somewhat 1184;; difficult to describe -- for an odd number of arguments, we favor the left side of the tree. 1185 1186(defmfun $tree_reduce (f a &optional (init 'no-init)) 1187 (setq a (require-list-or-set a '$tree_reduce)) 1188 (if (not (equal init 'no-init)) (push init a)) 1189 (if (null a) 1190 (merror (intl:gettext "tree_reduce: either a nonempty set or initial value must be given."))) 1191 1192 (let ((acc) (x) (doit nil)) 1193 (while (consp a) 1194 (setq x (pop a)) 1195 (while (consp a) 1196 (push (mfuncall f x (pop a)) acc) 1197 (if (setq doit (consp a)) (setq x (pop a)))) 1198 (if doit (push x acc)) 1199 (setq a (nreverse acc)) 1200 (setq acc nil)) 1201 x)) 1202 1203 1204 1205;; An identity function -- may see some use in things like 1206;; every(identity, [true, true, false, ..]). 1207 1208(defmfun $identity (x) x) 1209 1210;; Maxima 'some' and 'every' functions. The first argument should be 1211;; a predicate (a function that evaluates to true, false, or unknown). 1212;; The functions 'some' and 'every' locally bind $prederror to false. 1213;; Thus within 'some' or 'every,' is(a < b) evaluates to unknown instead 1214;; of signaling an error (as it would when $prederror is true). 1215;; 1216;; Three cases: 1217;; 1218;; (1) some(f, set(a1,...,an)) If any f(ai) evaluates to true, 1219;; 'some' returns true. 'Some' may or may not evaluate all the 1220;; f(ai)'s. Since sets are unordered, 'some' is free to evaluate 1221;; f(ai) in any order. To use 'some' on multiple set arguments, 1222;; they should first be converted to an ordered sequence so that 1223;; their relative alignment becomes well-defined. 1224 1225;; (2) some(f,[a11,...,a1n],[a21,...],...) If any f(ai1,ai2,...) 1226;; evaluates to true, 'some' returns true. 'Some' may or may not 1227;; evaluate all the f(ai)'s. Since sequences are ordered, 'some' 1228;; evaluates in the order of increasing 'i'. 1229 1230;; (3) some(f, matrix([a111,...],[a121,...],[a1n1...]), matrix(...)). 1231;; If any f(a1ij, a2ij, ...) evaluates to true, return true. 'Some' 1232;; may or may not evaluate all the predicates. Since there is no 1233;; natural order for the entries of a matrix, 'some' is free to 1234;; evaluate the predicates in any order. 1235 1236;; Notes: 1237;; (a) 'some' and 'every' automatically apply 'maybe'; thus the following 1238;; work correctly 1239;; 1240;; (C1) some("<",[a,b,5],[1,2,8]); 1241;; (D1) TRUE 1242;; (C2) some("=",[2,3],[2,7]); 1243;; (D2) TRUE 1244;; 1245;; (b) Since 'some' is free to choose the order of evaluation, and 1246;; possibly stop as soon as any one instance returns true, the 1247;; predicate f should not normally have side-effects or signal 1248;; errors. Similarly, 'every' may halt after one instance returns false; 1249;; however, the function 'maybe' is wrapped inside 'errset' This allows 1250;; some things to work that would otherwise signal an error: 1251 1252;; (%i1) some("<",[i,1],[3,12]); 1253;; (%o1) true 1254;; (%i2) every("<",[i,1],[3,12]); 1255;; (%o2) false 1256;; (%i3) maybe(%i < 3); 1257;; `sign' called on an imaginary argument: 1258 1259;; 1260;; (c) The functions 'some' and 'every' effectively use the functions 1261;; 'map' and 'matrixmap' to map the predicate over the arguments. The 1262;; option variable 'maperror' modifies the behavior of 'map' and 1263;; 'matrixmap;' similarly, the value of 'maperror' modifies the behavior 1264;; of 'some' and 'every.' 1265;; 1266;; (d) 'every' behaves similarly to 'some' except that 'every' returns 1267;; true iff every f evaluates to true for all its inputs. 1268;; 1269;; (e) If emptyp(e) is true, then some(f,e) --> false and every(f,e) --> true. 1270;; Thus (provided an error doesn't get signaled), we have the identities: 1271;; 1272;; some(f,s1) or some(f,s2) == some(f, union(s1,s2)), 1273;; every(f,s1) and every(f,s2) == every(f, union(s1,s2)). 1274;; Similarly, some(f) --> false and every(f) --> true. 1275 1276(defun checked-and (x) 1277 (setq x (mfuncall '$maybe `((mand) ,@x))) 1278 (cond ((or (eq x t) (eq x nil) (not $prederror)) x) 1279 ((eq x '$unknown) nil) 1280 (t 1281 ;; FOLLOWING MESSAGE IS UNREACHABLE FROM WHAT I CAN TELL 1282 ;; SINCE MAYBE RETURNS T, NIL, OR '$UNKNOWN 1283 (merror "Predicate isn't true/false valued; maybe you want to set 'prederror' to false")))) 1284 1285(defun checked-or (x) 1286 (setq x (mfuncall '$maybe `((mor) ,@x))) 1287 (cond ((or (eq x t) (eq x nil) (not $prederror)) x) 1288 ((eq x '$unknown) nil) 1289 (t 1290 ;; FOLLOWING MESSAGE IS UNREACHABLE FROM WHAT I CAN TELL 1291 ;; SINCE MAYBE RETURNS T, NIL, OR '$UNKNOWN 1292 (merror "Predicate isn't true/false valued; maybe you want to set 'prederror' to false")))) 1293 1294;; Apply the Maxima function f to x. If an error is signaled, return nil; otherwise 1295;; return (list (mfuncall f x)). 1296 1297(defun ignore-errors-mfuncall (f x) 1298 (let ((errcatch t)) 1299 (declare (special errcatch)) 1300 (errset (mfuncall f x)))) 1301 1302(defmfun $every (f &rest x) 1303 (cond ((or (null x) (and (null (cdr x)) ($emptyp (first x)))) t) 1304 1305 ((or ($listp (first x)) (and ($setp (first x)) (null (cdr x)))) 1306 (setq x (margs (simplify (apply #'map1 (cons f x))))) 1307 (checked-and (mapcar #'car (mapcar #'(lambda (s) (ignore-errors-mfuncall '$maybe s)) x)))) 1308 1309 ((every '$matrixp x) 1310 (let ((fmaplvl 2)) 1311 (setq x (margs (simplify (apply #'fmapl1 (cons f x))))) 1312 (checked-and (mapcar #'(lambda (s) ($every '$identity s)) x)))) 1313 1314 (t 1315 ;; NOT CLEAR FROM PRECEDING CODE WHAT IS "INVALID" HERE 1316 (merror (intl:gettext "every: invalid arguments."))))) 1317 1318(defmfun $some (f &rest x) 1319 (cond ((or (null x) (and (null (cdr x)) ($emptyp (first x)))) nil) 1320 1321 ((or ($listp (first x)) (and ($setp (first x)) (null (cdr x)))) 1322 (setq x (margs (simplify (apply #'map1 (cons f x))))) 1323 (checked-or (mapcar #'car (mapcar #'(lambda (s) (ignore-errors-mfuncall '$maybe s)) x)))) 1324 1325 ((every '$matrixp x) 1326 (let ((fmaplvl 2)) 1327 (setq x (margs (simplify (apply #'fmapl1 (cons f x))))) 1328 (checked-or (mapcar #'(lambda (s) ($some '$identity s)) x)))) 1329 1330 (t 1331 ;; NOT CLEAR FROM PRECEDING CODE WHAT IS "INVALID" HERE 1332 (merror (intl:gettext "some: invalid arguments."))))) 1333 1334(defmspec $makeset (l) 1335 (let* ((fn (car (pop l))) 1336 (f (if l (pop l) (wna-err fn))) 1337 (v (if l (pop l) (wna-err fn))) 1338 (s (if l (pop l) (wna-err fn)))) 1339 (if l (wna-err fn)) 1340 (if (or (not ($listp v)) (not (every #'(lambda (x) (or ($symbolp x) ($subvarp x))) (cdr v)))) 1341 (merror (intl:gettext "makeset: second argument must be a list of symbols; found: ~:M") v)) 1342 (setq s (require-list-or-set (meval s) '$makeset)) 1343 (setq f (list (list 'lambda) v f)) 1344 (setq v (margs v)) 1345 (dolist (sk v) (setq f (subst (gensym) sk f :test #'alike1))) 1346 (simplifya (cons '($set) (mapcar #'(lambda (x) (mfuncall '$apply f x)) s)) t))) 1347 1348;; Thread fn over l and apply op to the resulting list. 1349 1350(defun thread (fn l op) 1351 (simplify (cons `(,op) (mapcar #'(lambda (x) (simplify `((,fn) ,x))) l)))) 1352 1353;; Return a set of the divisors of n. If n isn't a positive integer, 1354;; return a noun form. We consider both 1 and n to be divisors of n. 1355;; The divisors of a negative number are the divisors of its absolute 1356;; value; divisors(0) simplifies to itself. We thread divisors over 1357;; lists, sets, matrices, and equalities. 1358 1359(defprop $divisors simp-divisors operators) 1360 1361(defun simp-divisors (n y z) 1362 (oneargcheck n) 1363 (setq y (caar n)) 1364 (setq n (simpcheck (cadr n) z)) 1365 (cond ((or ($listp n) ($setp n) ($matrixp n) (mequalp n)) 1366 (thread y (cdr n) (caar n))) 1367 ((and (integerp n) (not (zerop n))) 1368 (let (($intfaclim nil) 1369 (n (abs n))) 1370 `(($set simp) ,@(sort (mapcar #'car (divisors (cfactorw n))) #'<)))) 1371 (t `(($divisors simp) ,n)))) 1372 1373;; The Moebius function; it threads over lists, sets, matrices, and equalities. 1374 1375(defprop $moebius simp-moebius operators) 1376 1377(defun simp-moebius (n y z) 1378 (oneargcheck n) 1379 (setq y (caar n)) 1380 (setq n (simpcheck (cadr n) z)) 1381 (cond ((posint n) 1382 (if (= n 1) 1383 1 1384 (let (($intfaclim nil) 1385 (pfl (get-factor-list n))) ; pfl is a list of (prime exponent) pairs 1386 (if (every #'(lambda (x) (= 1 (second x))) pfl) ; if n is not 1387 ; squarefree return 0 1388 (if (evenp (length pfl)) 1 -1) ; else (-1)^(number of prime factors) 1389 0)))) 1390 ((or ($listp n) ($setp n) ($matrixp n) (mequalp n)) 1391 (thread y (cdr n) (caar n))) 1392 (t `(($moebius simp) ,n)))) 1393 1394; Find indices of elements which satisfy a predicate. 1395; Thanks to Bill Wood (william.wood3@comcast.net) for his help. 1396; Released under terms of GNU GPL v2 with Bill's approval. 1397 1398(defmfun $sublist_indices (items pred) 1399 (let ((items (require-list items '$sublist_indices))) 1400 (do ((i 0 (1+ i)) 1401 (xs items (cdr xs)) 1402 (acc '() (if (definitely-so (mfuncall pred (car xs))) (cons (1+ i) acc) acc))) 1403 ((endp xs) `((mlist) ,@(nreverse acc)))))) 1404