path: root/cring.scm

;;;"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))