Import Upstream version 2c0upstream/2c0

1 files changed, 480 insertions, 0 deletions

diff --git a/cring.scm b/cring.scm
new file mode 100644
index 0000000..3f594bc
--- /dev/null
+++ b/cring.scm
@@ -0,0 +1,480 @@
+;;;"cring.scm" Extend Scheme numerics to any commutative ring.
+;Copyright (C) 1997 Aubrey Jaffer
+;
+;Permission to copy this software, to redistribute it, and to use it
+;for any purpose is granted, subject to the following restrictions and
+;understandings.
+;
+;1.  Any copy made of this software must include this copyright notice
+;in full.
+;
+;2.  I have made no warrantee or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3.  In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+(require 'common-list-functions)
+(require 'relational-database)
+(require 'database-utilities)
+(require 'sort)
+
+(define (symbol-alpha? sym)
+  (char-alphabetic? (string-ref (symbol->string sym) 0)))
+(define (expression-< x y)
+  (cond ((and (number? x) (number? y)) (> x y))	;want negatives last
+	((number? x) #t)
+	((number? y) #f)
+	((and (symbol? x) (symbol? y))
+	 (cond ((eqv? (symbol-alpha? x) (symbol-alpha? y))
+		(string<? (symbol->string x) (symbol->string y)))
+	       (else (symbol-alpha? x))))
+	((symbol? x) #t)
+	((symbol? y) #f)
+	((null? x) #t)
+	((null? y) #f)
+	((expression-< (car x) (car y)) #t)
+	((expression-< (car y) (car x)) #f)
+	(else (expression-< (cdr x) (cdr y)))))
+(define (expression-sort seq) (sort! seq expression-<))
+
+(define cring:db (create-database #f 'alist-table))
+(define-tables cring:db
+  `(operation
+    ((op symbol)
+     (sub-op1 symbol)
+     (sub-op2 symbol))
+    ((reduction expression))
+    (;; This is the distributive rule (* over +)
+     (* + identity
+	,(lambda (exp1 exp2)
+	   ;;(print 'distributing '* '+ exp1 exp2 '==>)
+	   (apply + (map (lambda (trm) (* trm exp2)) (cdr exp1)))))
+     (* - identity
+	,(lambda (exp1 exp2)
+	   ;;(print 'distributing '* '- exp1 exp2 '==>)
+	   (apply - (map (lambda (trm) (* trm exp2)) (cdr exp1)))))
+     (/ + identity
+	,(lambda (exp1 exp2)
+	   ;;(print 'distributing '/ '+ exp1 exp2 '==>)
+	   (apply + (map (lambda (trm) (/ trm exp2)) (cdr exp1)))))
+     (/ - identity
+	,(lambda (exp1 exp2)
+	   ;;(print 'distributing '/ '- exp1 exp2 '==>)
+	   (apply - (map (lambda (trm) (/ trm exp2)) (cdr exp1))))))))
+
+(define cring:op-tab ((cring:db 'open-table) 'operation #t))
+(define cring:rule (cring:op-tab 'get 'reduction))
+(define cring:defrule (cring:op-tab 'row:update))
+(define (cring:define-rule . args) (cring:defrule args))
+
+(define number* *)
+(define number+ +)
+(define number- -)
+(define number/ /)
+(define number^ integer-expt)
+(define is-term-op? (lambda (term op) (and (pair? term) (eq? op (car term)))))
+;;(define (sign x) (if (positive? x) 1 (if (negative? x) -1 0)))
+(define number0? zero?)
+(define (zero? x) (and (number? x) (number0? x)))
+
+(define (make-rat n d)
+  (let* ((g (if (negative? d) (number- (gcd n d)) (gcd n d)))
+	 (n/g (quotient n g))
+	 (d/g (quotient d g)))
+    (case d/g
+      ((1) n/g)
+      (else
+       (case n/g
+	 ((0) 0)
+	 ((1) (list '/ d/g))
+	 (else (list '/ n/g d/g)))))))
+
+(define (rat-number? r)
+  (and (list? r)
+       (<= 2 (length r) 3)
+       (eq? '/ (car r))
+       (every number? (cdr r))))
+
+(define (rat-numerator r)
+  (cond ((number? r) r)
+	((rat-number? r)
+	 (case (length r)
+	   ((2) 1)
+	   ((3) (cadr r))))
+	(else (slib:error 'rat-numerator "of non-rat" r))))
+
+(define (rat-denominator r)
+  (cond ((number? r) 1)
+	((rat-number? r)
+	 (case (length r)
+	   ((2) (cadr r))
+	   ((3) (caddr r))))
+	(else (slib:error 'rat-denominator "of non-rat" r))))
+
+;; To convert to CR internal form, NUMBER-op all the `numbers' in the
+;; argument list and remove them from the argument list.  Collect the
+;; remaining arguments into equivalence classes, keeping track of the
+;; number of arguments in each class.  The returned list is thus:
+;; (<numeric> (<expression1> . <exp1>) ...)
+
+;;; Converts * argument list to CR internal form
+(define (cr*-args->fcts args)
+  ;;(print (cons 'cr*-args->fcts args) '==>)
+  (let loop ((args args) (pow 1) (nums 1) (denoms 1) (arg.exps '()))
+    ;;(print (list 'loop args pow nums denoms arg.exps) '==>)
+    (cond ((null? args) (cons (make-rat nums denoms) arg.exps))
+	  ((number? (car args))
+	   (let ((num^pow (number^ (car args) (abs pow))))
+	     (if (negative? pow)
+		 (loop (cdr args) pow nums (number* num^pow denoms) arg.exps)
+		 (loop (cdr args) pow (number* num^pow nums) denoms arg.exps))))
+	  ((rat-number? (car args))
+	   (let ((num^pow (number^ (rat-numerator (car args)) (abs pow)))
+		 (den^pow (number^ (rat-denominator (car args)) (abs pow))))
+	     (if (negative? pow)
+		 (loop (cdr args) pow
+		       (number* den^pow nums)
+		       (number* num^pow denoms) arg.exps)
+		 (loop (cdr args) pow
+		       (number* num^pow nums)
+		       (number* den^pow denoms) arg.exps))))
+	  ;; Associative Rule
+	  ((is-term-op? (car args) '*) (loop (append (cdar args) (cdr args))
+					     pow
+					     nums denoms
+					     arg.exps))
+	  ;; Do singlet -
+	  ((and (is-term-op? (car args) '-) (= 2 (length (car args))))
+	   ;;(print 'got-here (car args))
+	   (set! arg.exps (loop (cdar args) pow (number- nums) denoms arg.exps))
+	   (loop (cdr args) pow
+		 (rat-numerator (car arg.exps))
+		 (rat-denominator (car arg.exps))
+		 (cdr arg.exps)))
+	  ((and (is-term-op? (car args) '/) (= 2 (length (car args))))
+	   ;; Do singlet /
+	   ;;(print 'got-here=cr+ (car args))
+	   (set! arg.exps (loop (cdar args) (number- pow) nums denoms arg.exps))
+	   (loop (cdr args) pow
+		 (rat-numerator (car arg.exps))
+		 (rat-denominator (car arg.exps))
+		 (cdr arg.exps)))
+	  ((is-term-op? (car args) '/)
+	   ;; Do multi-arg /
+	   ;;(print 'doing '/ (cddar args) (number- pow))
+	   (set! arg.exps
+		 (loop (cddar args) (number- pow) nums denoms arg.exps))
+	   ;;(print 'finishing '/ (cons (cadar args) (cdr args)) pow)
+	   (loop (cons (cadar args) (cdr args))
+		 pow
+		 (rat-numerator (car arg.exps))
+		 (rat-denominator (car arg.exps))
+		 (cdr arg.exps)))
+	  ;; Pull out numeric exponents as powers
+	  ((and (is-term-op? (car args) '^)
+		(= 3 (length (car args)))
+		(number? (caddar args)))
+	   (set! arg.exps (loop (list (cadar args))
+				(number* pow (caddar args))
+				nums denoms
+				arg.exps))
+	   (loop (cdr args) pow
+		 (rat-numerator (car arg.exps))
+		 (rat-denominator (car arg.exps)) (cdr arg.exps)))
+	  ;; combine with same terms
+	  ((assoc (car args) arg.exps)
+	   => (lambda (pair) (set-cdr! pair (number+ pow (cdr pair)))
+		      (loop (cdr args) pow nums denoms arg.exps)))
+	  ;; Add new term to arg.exps
+	  (else (loop (cdr args) pow nums denoms
+		      (cons (cons (car args) pow) arg.exps))))))
+
+;;; Converts + argument list to CR internal form
+(define (cr+-args->trms args)
+  (let loop ((args args) (cof 1) (numbers 0) (arg.exps '()))
+    (cond ((null? args) (cons numbers arg.exps))
+	  ((number? (car args))
+	   (loop (cdr args)
+		 cof
+		 (number+ (number* (car args) cof) numbers)
+		 arg.exps))
+	  ;; Associative Rule
+	  ((is-term-op? (car args) '+) (loop (append (cdar args) (cdr args))
+					     cof
+					     numbers
+					     arg.exps))
+	  ;; Idempotent singlet *
+	  ((and (is-term-op? (car args) '*) (= 2 (length (car args))))
+	   (loop (cons (cadar args) (cdr args))
+		 cof
+		 numbers
+		 arg.exps))
+	  ((and (is-term-op? (car args) '-) (= 2 (length (car args))))
+	   ;; Do singlet -
+	   (set! arg.exps (loop (cdar args) (number- cof) numbers arg.exps))
+	   (loop (cdr args) cof (car arg.exps) (cdr arg.exps)))
+	  ;; Pull out numeric factors as coefficients
+	  ((and (is-term-op? (car args) '*) (some number? (cdar args)))
+	   ;;(print 'got-here (car args) '=> (cons '* (remove-if number? (cdar args))))
+	   (set! arg.exps
+		 (loop (list (cons '* (remove-if number? (cdar args))))
+		       (apply number* cof (remove-if-not number? (cdar args)))
+		       numbers
+		       arg.exps))
+	   (loop (cdr args) cof (car arg.exps) (cdr arg.exps)))
+	  ((is-term-op? (car args) '-)
+	   ;; Do multi-arg -
+	   (set! arg.exps (loop (cddar args) (number- cof) numbers arg.exps))
+	   (loop (cons (cadar args) (cdr args))
+		 cof
+		 (car arg.exps)
+		 (cdr arg.exps)))
+	  ;; combine with same terms
+	  ((assoc (car args) arg.exps)
+	   => (lambda (pair) (set-cdr! pair (number+ cof (cdr pair)))
+		      (loop (cdr args) cof numbers arg.exps)))
+	  ;; Add new term to arg.exps
+	  (else (loop (cdr args) cof numbers
+		      (cons (cons (car args) cof) arg.exps))))))
+
+;;; Converts + or * internal form to Scheme expression
+(define (cr-terms->form op ident inv-op higher-op res.cofs)
+  (define (negative-cof? fct.cof)
+    (negative? (cdr fct.cof)))
+  (define (finish exprs)
+    (if (null? exprs) ident
+	(if (null? (cdr exprs))
+	    (car exprs)
+	    (cons op exprs))))
+  (define (do-terms sign fct.cofs)
+    (expression-sort
+     (map (lambda (fct.cof)
+	    (define cof (number* sign (cdr fct.cof)))
+	    (cond ((eqv? 1 cof) (car fct.cof))
+		  ((number? (car fct.cof)) (number* cof (car fct.cof)))
+		  ((is-term-op? (car fct.cof) higher-op)
+		   (if (eq? higher-op '^)
+		       (list '^ (cadar fct.cof) (* cof (caddar fct.cof)))
+		       (cons higher-op (cons cof (cdar fct.cof)))))
+		  ((eqv? -1 cof) (list inv-op (car fct.cof)))
+		  (else (list higher-op (car fct.cof) cof))))
+	  fct.cofs)))
+  (let* ((all.cofs (remove-if (lambda (fct.cof)
+				(or (zero? (cdr fct.cof))
+				    (eqv? ident (car fct.cof))))
+			      res.cofs))
+	 (cofs (map cdr all.cofs))
+	 (some-positive? (some positive? cofs)))
+    ;;(print op 'positive? some-positive? 'negative? (some negative? cofs) all.cofs)
+    (cond ((and some-positive? (some negative? cofs))
+	   (append (list inv-op
+			 (finish (do-terms
+				  1 (remove-if negative-cof? all.cofs))))
+		   (do-terms -1 (remove-if-not negative-cof? all.cofs))))
+	  (some-positive? (finish (do-terms 1 all.cofs)))
+	  ((not (some negative? cofs)) ident)
+	  (else (list inv-op (finish (do-terms -1 all.cofs)))))))
+
+(define (* . args)
+  (cond
+   ((null? args) 1)
+   ;;((null? (cdr args)) (car args))
+   (else
+    (let ((in (cr*-args->fcts args)))
+      (cond
+       ((zero? (car in)) 0)
+       (else
+	(if (null? (cdr in))
+	    (set-cdr! in (list (cons 1 1))))
+	(let* ((num #f)
+	       (ans (cr-terms->form
+		     '* 1 '/ '^
+		     (apply
+		      (lambda (numeric red.cofs res.cofs)
+			(append
+			 (cond ((number? numeric)
+				(set! num numeric)
+				(list (cons (abs numeric) 1)))
+			       (else
+				(set! num (rat-numerator numeric))
+				(list (cons (abs num) 1)
+				      (cons (rat-denominator numeric) -1))))
+			 red.cofs
+			 res.cofs))
+		      (cr1 '* number* '^ '/ (car in) (cdr in))))))
+	  (if (negative? num) (list '- ans) ans))))))))
+
+(define (+ . args)
+  (cond ((null? args) 0)
+	;;((null? (cdr args)) (car args))
+	(else
+	 (let ((in (cr+-args->trms args)))
+	   (if (null? (cdr in))
+	       (car in)
+	       (cr-terms->form
+		'+ 0 '- '*
+		(apply (lambda (numeric red.cofs res.cofs)
+			 (append
+			  (list (if (and (number? numeric)
+					 (negative? numeric))
+				    (cons (abs numeric) -1)
+				    (cons numeric 1)))
+			  red.cofs
+			  res.cofs))
+		       (cr1 '+ number+ '* '- (car in) (cdr in)))))))))
+
+(define (- arg1 . args)
+  (if (null? args)
+      (if (number? arg1) (number- arg1)
+	  (* -1 arg1)			;(list '- arg1)
+	  )
+      (+ arg1 (* -1 (apply + args)))))
+
+;;(print `(/ ,arg1 ,@args) '=> )
+(define (/ arg1 . args)
+  (if (null? args)
+      (^ arg1 -1)
+      (* arg1 (^ (apply * args) -1))))
+
+(define (^ arg1 arg2)
+  (cond ((and (number? arg2) (integer? arg2))
+	 (* (list '^ arg1 arg2)))
+	(else (list '^ arg1 arg2))))
+
+;; TRY-EACH-PAIR-ONCE algorithm.  I think this does the minimum
+;; number of rule lookups given no information about how to sort
+;; terms.
+
+;; Pick equivalence classes one at a time and move them into the
+;; result set of equivalence classes by searching for rules to
+;; multiply an element of the chosen class by itself (if multiple) and
+;; the element of each class already in the result group.  Each
+;; (multiplicative) term resulting from rule application would be put
+;; in the result class, if that class exists; or put in an argument
+;; class if not.
+
+(define (cr1 op number-op hop inv-op numeric in)
+  (define red.pows '())
+  (define res.pows '())
+  (let loop.arg.pow.s ((arg (caar in)) (pow (cdar in)) (arg.pows (cdr in)))
+    (define (arg-loop arg.pows)
+      (if (null? arg.pows)
+	  (list numeric red.pows res.pows)
+	  (loop.arg.pow.s (caar arg.pows) (cdar arg.pows) (cdr arg.pows))))
+    (define (cring:apply-rule->terms exp1 exp2)
+      ;;(display op)
+      (let ((ans (cring:apply-rule op exp1 exp2)))
+	(cond ((not ans) #f)
+	      ((number? ans) (list ans))
+	      (else (list (cons ans 1))))))
+    (define (cring:apply-inv-rule->terms exp1 exp2)
+      ;;(display inv-op)
+      (let ((ans (cring:apply-rule inv-op exp1 exp2)))
+	(cond ((not ans) #f)
+	      ((number? ans) (list ans))
+	      (else (list (cons ans 1))))))
+    (define (merge-res tmp.pows multiplicity)
+      (cond ((null? tmp.pows))
+	    ((number? (car tmp.pows))
+	     (do ((m (number+ -1 (abs multiplicity)) (number+ -1 m))
+		  (n numeric (number-op n (abs (car tmp.pows)))))
+		 ((negative? m) (set! numeric n)))
+	     (merge-res (cdr tmp.pows) multiplicity))
+	    ((or (assoc (car tmp.pows) res.pows)
+		 (assoc (car tmp.pows) arg.pows))
+	     => (lambda (pair)
+		  (set-cdr! pair (number+
+				  pow (number-op multiplicity (cdar tmp.pows))))
+		  (merge-res (cdr tmp.pows) multiplicity)))
+	    ((assoc (car tmp.pows) red.pows)
+	     => (lambda (pair)
+		  (set! arg.pows
+			(cons (cons (caar tmp.pows)
+				    (number+
+				     (cdr pair)
+				     (number* multiplicity (cdar tmp.pows))))
+			      arg.pows))
+		  (set-cdr! pair 0)
+		  (merge-res (cdr tmp.pows) multiplicity)))
+	    (else (set! arg.pows
+			(cons (cons (caar tmp.pows)
+				    (number* multiplicity (cdar tmp.pows)))
+			      arg.pows))
+		  (merge-res (cdr tmp.pows) multiplicity))))
+    (define (try-fct.pow fct.pow)
+      ;;(print 'try-fct.pow fct.pow op 'arg arg 'pow pow)
+      (cond ((or (zero? (cdr fct.pow)) (number? (car fct.pow))) #f)
+	    ((not (and (number? pow) (number? (cdr fct.pow))
+		       (integer? pow)	;(integer? (cdr fct.pow))
+		       ))
+	     #f)
+	    ;;((zero? pow) (slib:error "Don't try exp-0 terms") #f)
+	    ;;((or (number? arg) (number? (car fct.pow)))
+	    ;; (slib:error 'found-number arg fct.pow) #f)
+	    ((and (positive? pow) (positive? (cdr fct.pow))
+		  (or (cring:apply-rule->terms arg (car fct.pow))
+		      (cring:apply-rule->terms (car fct.pow) arg)))
+	     => (lambda (terms)
+		  ;;(print op op terms)
+		  (let ((multiplicity (min pow (cdr fct.pow))))
+		    (set-cdr! fct.pow (number- (cdr fct.pow) multiplicity))
+		    (set! pow (number- pow multiplicity))
+		    (merge-res terms multiplicity))))
+	    ((and (negative? pow) (negative? (cdr fct.pow))
+		  (or (cring:apply-rule->terms arg (car fct.pow))
+		      (cring:apply-rule->terms (car fct.pow) arg)))
+	     => (lambda (terms)
+		  ;;(print inv-op inv-op terms)
+		  (let ((multiplicity (max pow (cdr fct.pow))))
+		    (set-cdr! fct.pow (number+ (cdr fct.pow) multiplicity))
+		    (set! pow (number+ pow multiplicity))
+		    (merge-res terms multiplicity))))
+	    ((and (positive? pow) (negative? (cdr fct.pow))
+		  (cring:apply-inv-rule->terms arg (car fct.pow)))
+	     => (lambda (terms)
+		  ;;(print op inv-op terms)
+		  (let ((multiplicity (min pow (number- (cdr fct.pow)))))
+		    (set-cdr! fct.pow (number+ (cdr fct.pow) multiplicity))
+		    (set! pow (number- pow multiplicity))
+		    (merge-res terms multiplicity))))
+	    ((and (negative? pow) (positive? (cdr fct.pow))
+		  (cring:apply-inv-rule->terms (car fct.pow) arg))
+	     => (lambda (terms)
+		  ;;(print inv-op op terms)
+		  (let ((multiplicity (max (number- pow) (cdr fct.pow))))
+		    (set-cdr! fct.pow (number- (cdr fct.pow) multiplicity))
+		    (set! pow (number+ pow multiplicity))
+		    (merge-res terms multiplicity))))
+	    (else #f)))
+    ;;(print op numeric 'arg arg 'pow pow 'arg.pows arg.pows 'red.pows red.pows 'res.pows res.pows)
+    ;;(trace arg-loop cring:apply-rule->terms merge-res try-fct.pow) (set! *qp-width* 333)
+    (cond ((or (zero? pow) (number? arg)) (arg-loop arg.pows))
+	  ((assoc arg res.pows) => (lambda (pair)
+				     (set-cdr! pair (number+ pow (cdr pair)))
+				     (arg-loop arg.pows)))
+	  ((and (> (abs pow) 1) (cring:apply-rule->terms arg arg))
+	   => (lambda (terms)
+		(merge-res terms (quotient pow 2))
+		(if (odd? pow)
+		    (loop.arg.pow.s arg 1 arg.pows)
+		    (arg-loop arg.pows))))
+	  ((or (some try-fct.pow res.pows) (some try-fct.pow arg.pows))
+	   (loop.arg.pow.s arg pow arg.pows))
+	  (else (set! res.pows (cons (cons arg pow) res.pows))
+		(arg-loop arg.pows)))))
+
+(define (cring:try-rule op sop1 sop2 exp1 exp2)
+  (let ((rule (cring:rule op sop1 sop2)))
+    (and rule (rule exp1 exp2))))
+
+(define (cring:apply-rule op exp1 exp2)
+  (and (pair? exp1)
+       (or (and (pair? exp2)
+		(cring:try-rule op (car exp1) (car exp2) exp1 exp2))
+	   (cring:try-rule op (car exp1) 'identity exp1 exp2))))
+
+;;(begin (trace cr-terms->form) (set! *qp-width* 333))