module/slib/cring.scm

   1 ;;;"cring.scm" Extend Scheme numerics to any commutative ring.
   2 ;Copyright (C) 1997, 1998 Aubrey Jaffer
   3 ;
   4 ;Permission to copy this software, to redistribute it, and to use it
   5 ;for any purpose is granted, subject to the following restrictions and
   6 ;understandings.
   7 ;
   8 ;1.  Any copy made of this software must include this copyright notice
   9 ;in full.
  10 ;
  11 ;2.  I have made no warrantee or representation that the operation of
  12 ;this software will be error-free, and I am under no obligation to
  13 ;provide any services, by way of maintenance, update, or otherwise.
  14 ;
  15 ;3.  In conjunction with products arising from the use of this
  16 ;material, there shall be no use of my name in any advertising,
  17 ;promotional, or sales literature without prior written consent in
  18 ;each case.
  19
  20 (require 'common-list-functions)
  21 (require 'relational-database)
  22 (require 'database-utilities)
  23 (require 'sort)
  24
  25 (define cring:db (create-database #f 'alist-table))
  26 (define (make-ruleset . rules)
  27   (define name #f)
  28   (cond ((and (not (null? rules)) (symbol? (car rules)))
  29          (set! name (car rules))
  30          (set! rules (cdr rules)))
  31         (else (set! name (gentemp))))
  32   (define-tables cring:db
  33     (list name
  34           '((op symbol)
  35             (sub-op1 symbol)
  36             (sub-op2 symbol))
  37           '((reduction expression))
  38           rules))
  39   (let ((table ((cring:db 'open-table) name #t)))
  40     (and table
  41          (list (table 'get 'reduction)
  42                (table 'row:update)
  43                table))))
  44 (define *ruleset* (make-ruleset 'default))
  45 (define (cring:define-rule . args)
  46   (if *ruleset*
  47       ((cadr *ruleset*) args)
  48       (slib:warn "No ruleset in *ruleset*")))
  49
  50 (define (combined-rulesets . rulesets)
  51   (define name #f)
  52   (cond ((symbol? (car rulesets))
  53          (set! name (car rulesets))
  54          (set! rulesets (cdr rulesets)))
  55         (else (set! name (gentemp))))
  56   (apply make-ruleset name
  57          (apply append
  58                 (map (lambda (ruleset) (((caddr ruleset) 'row:retrieve*)))
  59                      rulesets))))
  60
  61 ;;; Distribute * over + (and -)
  62 (define distribute*
  63   (make-ruleset
  64    'distribute*
  65    `(* + identity
  66        ,(lambda (exp1 exp2)
  67           ;;(print 'distributing '* '+ exp1 exp2 '==>)
  68           (apply + (map (lambda (trm) (* trm exp2)) (cdr exp1)))))
  69    `(* - identity
  70        ,(lambda (exp1 exp2)
  71           ;;(print 'distributing '* '- exp1 exp2 '==>)
  72           (apply - (map (lambda (trm) (* trm exp2)) (cdr exp1)))))))
  73
  74 ;;; Distribute / over + (and -)
  75 (define distribute/
  76   (make-ruleset
  77    'distribute/
  78    `(/ + identity
  79        ,(lambda (exp1 exp2)
  80           ;;(print 'distributing '/ '+ exp1 exp2 '==>)
  81           (apply + (map (lambda (trm) (/ trm exp2)) (cdr exp1)))))
  82    `(/ - identity
  83        ,(lambda (exp1 exp2)
  84           ;;(print 'distributing '/ '- exp1 exp2 '==>)
  85           (apply - (map (lambda (trm) (/ trm exp2)) (cdr exp1)))))))
  86
  87 (define (symbol-alpha? sym)
  88   (char-alphabetic? (string-ref (symbol->string sym) 0)))
  89 (define (expression-< x y)
  90   (cond ((and (number? x) (number? y)) (> x y)) ;want negatives last
  91         ((number? x) #t)
  92         ((number? y) #f)
  93         ((and (symbol? x) (symbol? y))
  94          (cond ((eqv? (symbol-alpha? x) (symbol-alpha? y))
  95                 (string<? (symbol->string x) (symbol->string y)))
  96                (else (symbol-alpha? x))))
  97         ((symbol? x) #t)
  98         ((symbol? y) #f)
  99         ((null? x) #t)
 100         ((null? y) #f)
 101         ((expression-< (car x) (car y)) #t)
 102         ((expression-< (car y) (car x)) #f)
 103         (else (expression-< (cdr x) (cdr y)))))
 104 (define (expression-sort seq) (sort! seq expression-<))
 105
 106 (define number* *)
 107 (define number+ +)
 108 (define number- -)
 109 (define number/ /)
 110 (define number^ integer-expt)
 111 (define is-term-op? (lambda (term op) (and (pair? term) (eq? op (car term)))))
 112 ;;(define (sign x) (if (positive? x) 1 (if (negative? x) -1 0)))
 113 (define number0? zero?)
 114 (define (zero? x) (and (number? x) (number0? x)))
 115
 116 ;; To convert to CR internal form, NUMBER-op all the `numbers' in the
 117 ;; argument list and remove them from the argument list.  Collect the
 118 ;; remaining arguments into equivalence classes, keeping track of the
 119 ;; number of arguments in each class.  The returned list is thus:
 120 ;; (<numeric> (<expression1> . <exp1>) ...)
 121
 122 ;;; Converts * argument list to CR internal form
 123 (define (cr*-args->fcts args)
 124   ;;(print (cons 'cr*-args->fcts args) '==>)
 125   (let loop ((args args) (pow 1) (nums 1) (arg.exps '()))
 126     ;;(print (list 'loop args pow nums denoms arg.exps) '==>)
 127     (cond ((null? args) (cons nums arg.exps))
 128           ((number? (car args))
 129            (let ((num^pow (number^ (car args) (abs pow))))
 130              (if (negative? pow)
 131                  (loop (cdr args) pow (number/ (number* num^pow nums))
 132                        arg.exps)
 133                  (loop (cdr args) pow (number* num^pow nums) arg.exps))))
 134           ;; Associative Rule
 135           ((is-term-op? (car args) '*) (loop (append (cdar args) (cdr args))
 136                                              pow nums arg.exps))
 137           ;; Do singlet -
 138           ((and (is-term-op? (car args) '-) (= 2 (length (car args))))
 139            ;;(print 'got-here (car args))
 140            (set! arg.exps (loop (cdar args) pow (number- nums) arg.exps))
 141            (loop (cdr args) pow
 142                  (car arg.exps)
 143                  (cdr arg.exps)))
 144           ((and (is-term-op? (car args) '/) (= 2 (length (car args))))
 145            ;; Do singlet /
 146            ;;(print 'got-here=cr+ (car args))
 147            (set! arg.exps (loop (cdar args) (number- pow) nums arg.exps))
 148            (loop (cdr args) pow
 149                  (car arg.exps)
 150                  (cdr arg.exps)))
 151           ((is-term-op? (car args) '/)
 152            ;; Do multi-arg /
 153            ;;(print 'doing '/ (cddar args) (number- pow))
 154            (set! arg.exps
 155                  (loop (cddar args) (number- pow) nums arg.exps))
 156            ;;(print 'finishing '/ (cons (cadar args) (cdr args)) pow)
 157            (loop (cons (cadar args) (cdr args))
 158                  pow
 159                  (car arg.exps)
 160                  (cdr arg.exps)))
 161           ;; Pull out numeric exponents as powers
 162           ((and (is-term-op? (car args) '^)
 163                 (= 3 (length (car args)))
 164                 (number? (caddar args)))
 165            (set! arg.exps (loop (list (cadar args))
 166                                 (number* pow (caddar args))
 167                                 nums
 168                                 arg.exps))
 169            (loop (cdr args) pow (car arg.exps) (cdr arg.exps)))
 170           ;; combine with same terms
 171           ((assoc (car args) arg.exps)
 172            => (lambda (pair) (set-cdr! pair (number+ pow (cdr pair)))
 173                       (loop (cdr args) pow nums arg.exps)))
 174           ;; Add new term to arg.exps
 175           (else (loop (cdr args) pow nums
 176                       (cons (cons (car args) pow) arg.exps))))))
 177
 178 ;;; Converts + argument list to CR internal form
 179 (define (cr+-args->trms args)
 180   (let loop ((args args) (cof 1) (numbers 0) (arg.exps '()))
 181     (cond ((null? args) (cons numbers arg.exps))
 182           ((number? (car args))
 183            (loop (cdr args)
 184                  cof
 185                  (number+ (number* (car args) cof) numbers)
 186                  arg.exps))
 187           ;; Associative Rule
 188           ((is-term-op? (car args) '+) (loop (append (cdar args) (cdr args))
 189                                              cof
 190                                              numbers
 191                                              arg.exps))
 192           ;; Idempotent singlet *
 193           ((and (is-term-op? (car args) '*) (= 2 (length (car args))))
 194            (loop (cons (cadar args) (cdr args))
 195                  cof
 196                  numbers
 197                  arg.exps))
 198           ((and (is-term-op? (car args) '-) (= 2 (length (car args))))
 199            ;; Do singlet -
 200            (set! arg.exps (loop (cdar args) (number- cof) numbers arg.exps))
 201            (loop (cdr args) cof (car arg.exps) (cdr arg.exps)))
 202           ;; Pull out numeric factors as coefficients
 203           ((and (is-term-op? (car args) '*) (some number? (cdar args)))
 204            ;;(print 'got-here (car args) '=> (cons '* (remove-if number? (cdar args))))
 205            (set! arg.exps
 206                  (loop (list (cons '* (remove-if number? (cdar args))))
 207                        (apply number* cof (remove-if-not number? (cdar args)))
 208                        numbers
 209                        arg.exps))
 210            (loop (cdr args) cof (car arg.exps) (cdr arg.exps)))
 211           ((is-term-op? (car args) '-)
 212            ;; Do multi-arg -
 213            (set! arg.exps (loop (cddar args) (number- cof) numbers arg.exps))
 214            (loop (cons (cadar args) (cdr args))
 215                  cof
 216                  (car arg.exps)
 217                  (cdr arg.exps)))
 218           ;; combine with same terms
 219           ((assoc (car args) arg.exps)
 220            => (lambda (pair) (set-cdr! pair (number+ cof (cdr pair)))
 221                       (loop (cdr args) cof numbers arg.exps)))
 222           ;; Add new term to arg.exps
 223           (else (loop (cdr args) cof numbers
 224                       (cons (cons (car args) cof) arg.exps))))))
 225
 226 ;;; Converts + or * internal form to Scheme expression
 227 (define (cr-terms->form op ident inv-op higher-op res.cofs)
 228   (define (negative-cof? fct.cof)
 229     (negative? (cdr fct.cof)))
 230   (define (finish exprs)
 231     (if (null? exprs) ident
 232         (if (null? (cdr exprs))
 233             (car exprs)
 234             (cons op exprs))))
 235   (define (do-terms sign fct.cofs)
 236     (expression-sort
 237      (map (lambda (fct.cof)
 238             (define cof (number* sign (cdr fct.cof)))
 239             (cond ((eqv? 1 cof) (car fct.cof))
 240                   ((number? (car fct.cof)) (number* cof (car fct.cof)))
 241                   ((is-term-op? (car fct.cof) higher-op)
 242                    (if (eq? higher-op '^)
 243                        (list '^ (cadar fct.cof) (* cof (caddar fct.cof)))
 244                        (cons higher-op (cons cof (cdar fct.cof)))))
 245                   ((eqv? -1 cof) (list inv-op (car fct.cof)))
 246                   (else (list higher-op (car fct.cof) cof))))
 247           fct.cofs)))
 248   (let* ((all.cofs (remove-if (lambda (fct.cof)
 249                                 (or (zero? (cdr fct.cof))
 250                                     (eqv? ident (car fct.cof))))
 251                               res.cofs))
 252          (cofs (map cdr all.cofs))
 253          (some-positive? (some positive? cofs)))
 254     ;;(print op 'positive? some-positive? 'negative? (some negative? cofs) all.cofs)
 255     (cond ((and some-positive? (some negative? cofs))
 256            (append (list inv-op
 257                          (finish (do-terms
 258                                   1 (remove-if negative-cof? all.cofs))))
 259                    (do-terms -1 (remove-if-not negative-cof? all.cofs))))
 260           (some-positive? (finish (do-terms 1 all.cofs)))
 261           ((not (some negative? cofs)) ident)
 262           (else (list inv-op (finish (do-terms -1 all.cofs)))))))
 263
 264 (define (* . args)
 265   (cond
 266    ((null? args) 1)
 267    ;;This next line is commented out so ^ will collapse numerical expressions.
 268    ;;((null? (cdr args)) (car args))
 269    (else
 270     (let ((in (cr*-args->fcts args)))
 271       (cond
 272        ((zero? (car in)) 0)
 273        (else
 274         (if (null? (cdr in))
 275             (set-cdr! in (list (cons 1 1))))
 276         (let* ((num #f)
 277                (ans (cr-terms->form
 278                      '* 1 '/ '^
 279                      (apply
 280                       (lambda (numeric red.cofs res.cofs)
 281                         (set! num numeric)
 282                         (append
 283                          ;;(list (cons (abs numeric) 1))
 284                          red.cofs
 285                          res.cofs))
 286                       (cr1 '* number* '^ '/ (car in) (cdr in))))))
 287           (cond ((number0? (+ -1 num)) ans)
 288                 ((number? ans) (number* num ans))
 289                 ((number0? (+ 1 num))
 290                  (if (and (list? ans) (= 2 (length ans)) (eq? '- (car ans)))
 291                      (cadr ans)
 292                      (list '- ans)))
 293                 ((not (pair? ans)) (list '* num ans))
 294                 (else
 295                  (case (car ans)
 296                    ((*) (append (list '* num) (cdr ans)))
 297                    ((+) (apply + (map (lambda (mon) (* num mon)) (cdr ans))))
 298                    ((-) (apply - (map (lambda (mon) (* num mon)) (cdr ans))))
 299                    (else (list '* num ans))))))))))))
 300
 301 (define (+ . args)
 302   (cond ((null? args) 0)
 303         ;;((null? (cdr args)) (car args))
 304         (else
 305          (let ((in (cr+-args->trms args)))
 306            (if (null? (cdr in))
 307                (car in)
 308                (cr-terms->form
 309                 '+ 0 '- '*
 310                 (apply (lambda (numeric red.cofs res.cofs)
 311                          (append
 312                           (list (if (and (number? numeric)
 313                                          (negative? numeric))
 314                                     (cons (abs numeric) -1)
 315                                     (cons numeric 1)))
 316                           red.cofs
 317                           res.cofs))
 318                        (cr1 '+ number+ '* '- (car in) (cdr in)))))))))
 319
 320 (define (- arg1 . args)
 321   (if (null? args)
 322       (if (number? arg1) (number- arg1)
 323           (* -1 arg1)                   ;(list '- arg1)
 324           )
 325       (+ arg1 (* -1 (apply + args)))))
 326
 327 ;;(print `(/ ,arg1 ,@args) '=> )
 328 (define (/ arg1 . args)
 329   (if (null? args)
 330       (^ arg1 -1)
 331       (* arg1 (^ (apply * args) -1))))
 332
 333 (define (^ arg1 arg2)
 334   (cond ((and (number? arg2) (integer? arg2))
 335          (* (list '^ arg1 arg2)))
 336         (else (list '^ arg1 arg2))))
 337
 338 ;; TRY-EACH-PAIR-ONCE algorithm.  I think this does the minimum
 339 ;; number of rule lookups given no information about how to sort
 340 ;; terms.
 341
 342 ;; Pick equivalence classes one at a time and move them into the
 343 ;; result set of equivalence classes by searching for rules to
 344 ;; multiply an element of the chosen class by itself (if multiple) and
 345 ;; the element of each class already in the result group.  Each
 346 ;; (multiplicative) term resulting from rule application would be put
 347 ;; in the result class, if that class exists; or put in an argument
 348 ;; class if not.
 349
 350 (define (cr1 op number-op hop inv-op numeric in)
 351   (define red.pows '())
 352   (define res.pows '())
 353   (define (cring:apply-rule->terms exp1 exp2) ;(display op)
 354     (let ((ans (cring:apply-rule op exp1 exp2)))
 355       (cond ((not ans) #f)
 356             ((number? ans) (list ans))
 357             (else (list (cons ans 1))))))
 358   (define (cring:apply-inv-rule->terms exp1 exp2) ;(display inv-op)
 359     (let ((ans (cring:apply-rule inv-op exp1 exp2)))
 360       (cond ((not ans) #f)
 361             ((number? ans) (list ans))
 362             (else (list (cons ans 1))))))
 363   (let loop.arg.pow.s ((arg (caar in)) (pow (cdar in)) (arg.pows (cdr in)))
 364     (define (arg-loop arg.pows)
 365       (cond ((not (null? arg.pows))
 366              (loop.arg.pow.s (caar arg.pows) (cdar arg.pows) (cdr arg.pows)))
 367             (else (list numeric red.pows res.pows)))) ; Actually return!
 368     (define (merge-res tmp.pows multiplicity)
 369       (cond ((null? tmp.pows))
 370             ((number? (car tmp.pows))
 371              (do ((m (number+ -1 (abs multiplicity)) (number+ -1 m))
 372                   (n numeric (number-op n (abs (car tmp.pows)))))
 373                  ((negative? m) (set! numeric n)))
 374              (merge-res (cdr tmp.pows) multiplicity))
 375             ((or (assoc (car tmp.pows) res.pows)
 376                  (assoc (car tmp.pows) arg.pows))
 377              => (lambda (pair)
 378                   (set-cdr! pair (number+
 379                                   pow (number-op multiplicity (cdar tmp.pows))))
 380                   (merge-res (cdr tmp.pows) multiplicity)))
 381             ((assoc (car tmp.pows) red.pows)
 382              => (lambda (pair)
 383                   (set! arg.pows
 384                         (cons (cons (caar tmp.pows)
 385                                     (number+
 386                                      (cdr pair)
 387                                      (number* multiplicity (cdar tmp.pows))))
 388                               arg.pows))
 389                   (set-cdr! pair 0)
 390                   (merge-res (cdr tmp.pows) multiplicity)))
 391             (else (set! arg.pows
 392                         (cons (cons (caar tmp.pows)
 393                                     (number* multiplicity (cdar tmp.pows)))
 394                               arg.pows))
 395                   (merge-res (cdr tmp.pows) multiplicity))))
 396     (define (try-fct.pow fct.pow)
 397       ;;(print 'try-fct.pow fct.pow op 'arg arg 'pow pow)
 398       (cond ((or (zero? (cdr fct.pow)) (number? (car fct.pow))) #f)
 399             ((not (and (number? pow) (number? (cdr fct.pow))
 400                        (integer? pow)   ;(integer? (cdr fct.pow))
 401                        ))
 402              #f)
 403             ;;((zero? pow) (slib:error "Don't try exp-0 terms") #f)
 404             ;;((or (number? arg) (number? (car fct.pow)))
 405             ;; (slib:error 'found-number arg fct.pow) #f)
 406             ((and (positive? pow) (positive? (cdr fct.pow))
 407                   (or (cring:apply-rule->terms arg (car fct.pow))
 408                       (cring:apply-rule->terms (car fct.pow) arg)))
 409              => (lambda (terms)
 410                   ;;(print op op terms)
 411                   (let ((multiplicity (min pow (cdr fct.pow))))
 412                     (set-cdr! fct.pow (number- (cdr fct.pow) multiplicity))
 413                     (set! pow (number- pow multiplicity))
 414                     (merge-res terms multiplicity))))
 415             ((and (negative? pow) (negative? (cdr fct.pow))
 416                   (or (cring:apply-rule->terms arg (car fct.pow))
 417                       (cring:apply-rule->terms (car fct.pow) arg)))
 418              => (lambda (terms)
 419                   ;;(print inv-op inv-op terms)
 420                   (let ((multiplicity (max pow (cdr fct.pow))))
 421                     (set-cdr! fct.pow (number+ (cdr fct.pow) multiplicity))
 422                     (set! pow (number+ pow multiplicity))
 423                     (merge-res terms multiplicity))))
 424             ((and (positive? pow) (negative? (cdr fct.pow))
 425                   (cring:apply-inv-rule->terms arg (car fct.pow)))
 426              => (lambda (terms)
 427                   ;;(print op inv-op terms)
 428                   (let ((multiplicity (min pow (number- (cdr fct.pow)))))
 429                     (set-cdr! fct.pow (number+ (cdr fct.pow) multiplicity))
 430                     (set! pow (number- pow multiplicity))
 431                     (merge-res terms multiplicity))))
 432             ((and (negative? pow) (positive? (cdr fct.pow))
 433                   (cring:apply-inv-rule->terms (car fct.pow) arg))
 434              => (lambda (terms)
 435                   ;;(print inv-op op terms)
 436                   (let ((multiplicity (max (number- pow) (cdr fct.pow))))
 437                     (set-cdr! fct.pow (number- (cdr fct.pow) multiplicity))
 438                     (set! pow (number+ pow multiplicity))
 439                     (merge-res terms multiplicity))))
 440             (else #f)))
 441     ;;(print op numeric 'arg arg 'pow pow 'arg.pows arg.pows 'red.pows red.pows 'res.pows res.pows)
 442     ;;(trace arg-loop cring:apply-rule->terms merge-res try-fct.pow) (set! *qp-width* 333)
 443     (cond ((or (zero? pow) (eqv? 1 arg)) ;(number? arg) arg seems to always be 1
 444            (arg-loop arg.pows))
 445           ((assoc arg res.pows) => (lambda (pair)
 446                                      (set-cdr! pair (number+ pow (cdr pair)))
 447                                      (arg-loop arg.pows)))
 448           ((and (> (abs pow) 1) (cring:apply-rule->terms arg arg))
 449            => (lambda (terms)
 450                 (merge-res terms (quotient pow 2))
 451                 (if (odd? pow)
 452                     (loop.arg.pow.s arg 1 arg.pows)
 453                     (arg-loop arg.pows))))
 454           ((or (some try-fct.pow res.pows) (some try-fct.pow arg.pows))
 455            (loop.arg.pow.s arg pow arg.pows))
 456           (else (set! res.pows (cons (cons arg pow) res.pows))
 457                 (arg-loop arg.pows)))))
 458
 459 (define (cring:try-rule op sop1 sop2 exp1 exp2)
 460   (and *ruleset*
 461        (let ((rule ((car *ruleset*) op sop1 sop2)))
 462          (and rule (rule exp1 exp2)))))
 463
 464 (define (cring:apply-rule op exp1 exp2)
 465   (and (pair? exp1)
 466        (or (and (pair? exp2)
 467                 (cring:try-rule op (car exp1) (car exp2) exp1 exp2))
 468            (cring:try-rule op (car exp1) 'identity exp1 exp2))))
 469
 470 ;;(begin (trace cr-terms->form) (set! *qp-width* 333))