summaryrefslogtreecommitdiffstats
path: root/cring.scm
diff options
context:
space:
mode:
Diffstat (limited to 'cring.scm')
-rw-r--r--cring.scm480
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))