share/contrib/altsimp/altsimp.lisp

   1 ;; Author: Barton Willis with help from Richard Fateman
   2
   3 #|
   4 To simplify a sum with n terms, the standard simplus function calls
   5 great O(n^2) times. By using sorting more effectively, this code
   6 reduces the calls to O(n log_2(n)).
   7
   8 Also, this code tries to be "infinity correct" for addition. By this I
   9 mean inf + inf --> inf, inf + number --> inf, and inf + minf --> und.
  10 Since -1 * inf does not simplify to minf, this code doesn't simplify
  11 inf - inf to und; consequently, this portion of the code is largely
  12 untested. There are other problems too.  For one, this code does
  13 f(inf) - f(inf) --> 0 (comment from Stavros Macrakis). I don't know
  14 how far we can go with such things without making a mess. You could
  15 argue that Maxima should do
  16
  17    f(x) - f(x) --> if finitep(f(x)) then 0 else und
  18
  19 instead of f(x) - f(x) --> 0.
  20
  21 There is a great deal more that could be done. We could tag each
  22 infinity with a unique identifier (idea from RJF). That way we could
  23 have x : inf, x - x --> 0 and x : inf, y : inf, x - y --> und, for
  24 example.
  25
  26 In short, this code is highly experimental; you should not use it for
  27 anything that is important. I place it in /share/contrib because
  28 others encouraged me too; also I thought it might be useful to others
  29 who would like to experiment with similar code. Since a great deal of
  30 work has gone into the current simplus code, I'm not sure that a total
  31 rewrite is the best route.
  32
  33 Maybe the special case dispatch part of this code makes the task of
  34 extending simplus (to intervals, for example) easier. The basic design
  35 of this code is due to Barton Willis.
  36 |#
  37
  38 #|
  39
  40 As of late May 2021, altsimp (with use_extended_real_arithmetic set to false), runs
  41 the testsuite + the share testsuite with nine failures; the failures are
  42
  43  rtest11.mac problems (4 8)
  44  rtest14.mac problem: (62)
  45  rtest_expintegral.mac problems: (173 174)
  46  rtest_dgesv.mac problem: (5)
  47  rtest_to_poly_solve.mac problem: (241)
  48
  49 Additionally, rtest3 #146 fails (errcatch(integrate(exp(2^(3/2)*x^2-sqrt(2)*x^2),x))), but
  50 actually, altsimp fixes this bug. This is due to the fact that with altsimp, exp(2^(3/2)*x^2-sqrt(2)*x^2)
  51 simplifies to %e^(sqrt(2)*x^2).
  52
  53 Reasons for these failures include:
  54
  55 * Possible differences in ordering of addition of floating point numbers.
  56
  57 Using CCL 1.12.1 (64 bit) and altsimp, timings are:
  58
  59 took 752,166,000 microseconds (752.166000 seconds) to run.
  60       12,556,790 microseconds ( 12.556790 seconds, 1.67%) of which was spent in GC.
  61 During that period, and with 4 available CPU cores,
  62      687,062,500 microseconds (687.062500 seconds) were spent in user mode
  63       38,109,375 microseconds ( 38.109375 seconds) were spent in system mode
  64 87,194,237,408 bytes of memory allocated.
  65
  66 And timings for standard Maxima are:
  67
  68 took 689,231,000 microseconds (689.231000 seconds) to run.
  69       11,857,789 microseconds ( 11.857789 seconds, 1.72%) of which was spent in GC.
  70 During that period, and with 4 available CPU cores,
  71      644,234,375 microseconds (644.234400 seconds) were spent in user mode
  72       23,062,500 microseconds ( 23.062500 seconds) were spent in system mode
  73  77,100,115,920 bytes of memory allocated.
  74
  75  Compared to standard Maxima, altsimp runs the testsuites more slowly and uses more
  76  memory. Using sbcl 2.0, the results are similar, but faster and more memory.
  77
  78 Using CCL 1.12.1 (64 bit) and use_extended_real_arithmetic set to true, the
  79 testsuites give 45 failures:
  80
  81 Error(s) found:
  82   rtest11.mac problems:  (4 8)
  83   rtest14.mac problem: (62)
  84   rtest16.mac problems:  (457 459)
  85   rtest3.mac problem: (146)
  86   rtest_expintegral.mac problems: (173 174)
  87   rtest_powerseries.mac problems: (53 63)
  88   solve_rec/rtest_simplify_sum.mac problems: (3 4 5 6 7 10 11 12 13 14 16 17 19
  89   20 21 22 23 28 29 53 67 68 69 70 71 72)
  90   rtest_dgesv.mac problem: (5)
  91   rtest_linalg.mac problems: (92 94)
  92   rtest_abs_integrate.mac problems: (75 76 77 132 140)
  93   rtest_to_poly_solve.mac problem: (241)
  94
  95 Most of these additional failures are due to dispatching sign on und.
  96
  97
  98
  99 Speculation on how to speed up simplification of sums:
 100
 101 The altsimp algorithm uses sorting to speed up simplification of expressions with a
 102 large number of summands, but at least for running the testsuites, simplus is mostly
 103 called with just two summands. If one of the two summands is a sum, say XXX, altsimp
 104 re-examines the terms of XXX for common terms. That is, of course, wasteful. If the
 105 summands were marked in a way that indicated that pairs of terms were not common, that
 106 would possibly speed the code.
 107
 108 |#
 109
 110 (in-package :maxima)
 111 (declaim (optimize (speed 3) (safety 0)))
 112
 113 ;; When the option variable $use_extended_real_arithmetic is true, simplus adds the extended real
 114 ;; numbers (minf,inf,infinity, und, ind) in a mathematically logical way. But since simptimes and
 115 ;; simpexpt do not have an option of doing correct extended real number arithmetic, this feature of
 116 ;; simplus is limited.
 117
 118 (defmvar $use_extended_real_arithmetic nil)
 119
 120 (define-modify-macro mincf (&optional (i 1)) addk)
 121
 122 ;; Return true if z = 0, 0.0, or a bigfloat zero.
 123 (defun mzerop (z)
 124   (and (mnump z)
 125        (or (and (numberp z)(= z 0))
 126                  (and (bigfloatp z)(= (cadr z) 0))))) ;bigfloat zeros may be diff precisions
 127
 128 ;; Return true if x has the form integer^(rational)
 129 (defun surd-p (x)
 130   (and (mexptp x) (integerp (cadr x)) ($ratnump (caddr x))))
 131
 132 (defun generalized-surd-p (x)
 133   (and (mexptp x) (integerp (cadr x)) (mplusp (caddr x))
 134          (or ($ratnump (second (caddr x))) (integerp (second (caddr x))))))
 135
 136 (defun surd-convert (x)
 137    (cond ((or (integerp x) (mnump x))
 138            (cons 1 x))
 139          ((surd-p x)
 140           ;; x = a^(p/q) = a^n * a^(p/q - n), where n = floor(p/q).
 141             (let ((a (cadr x)) (p (caddr x)) (q) (n))
 142               (setq q ($denom p))
 143               (setq p ($num p))
 144               (setq n (floor p q))
 145               (cons (take '(mexpt) a (sub (caddr x) n)) (take '(mexpt) a n))))
 146           ((generalized-surd-p x) ; x = ((mexpt) a ((mplus) n a1 a2 ...)))
 147            (let ((a (second x)) (n (second (third x))) (p) (q))
 148               (setq p ($num n))
 149               (setq q ($denom n))
 150               (setq n (floor p q))
 151               (cons (take '(mexpt) a (sub (third x) n)) (take '(mexpt) a n))))
 152
 153           (t (cons x 1))))
 154
 155 ;; Convert a term (a non sum expression) to the form (e . coefficient), where
 156 ;; coefficient is a rational or floating point number. Factors of the form
 157 ;; integer^(rational) are converted by surd-convert.
 158 (defun convert-to-coeff-form (x)
 159   (let ((c 1) (xx) (qq 1))
 160     (cond ((mnump x) (cons 1 x))
 161           ((or (surd-p x) (generalized-surd-p x))
 162             (surd-convert x))
 163
 164                 ((mtimesp x)
 165                   (pop x)  ;remove (car x) which is (mtimes ..)
 166             (while (or (mnump (car x)) (surd-p (car x)) (generalized-surd-p (car x)))
 167                (setq xx (surd-convert (car x)))
 168                (setq c (mul c (cdr xx)))
 169                (setq qq (mul qq (car xx)))
 170                (pop x))
 171
 172             (when (not (eql qq 1))
 173                 (push qq x))
 174
 175             (when (null x)
 176                (push 1 x))
 177
 178                   (if (null (cdr x)); if only one more item, that's it.
 179                               (cons (car x) c)
 180                               (cons (cons (get 'mtimes 'msimpind) x) c)))
 181           (t (cons x 1)))))
 182
 183 ;; Previously there was a specialized function for multiplying a number times an expression. The
 184 ;; motivation was, I think, speed. But the specialized function was responsible for 22 testsuite
 185 ;; failures (May 2021) and it didn't contribute to running the testsuite any faster. So let us
 186 ;; replace the specialized function with a call to mul.
 187
 188 (defun number-times-expr (cf e)
 189   (mul cf e))
 190
 191 ;; Add an expression x to a list of equalities l.
 192 (defun add-expr-mequal (x l)
 193   (setq l (mapcar 'cdr l))
 194   (push (list x x) l)
 195   (setq l (list (reduce #'add (mapcar 'first l)) (reduce #'add (mapcar 'second l))))
 196   (simplifya (cons '(mequal) l) t))
 197
 198 (defun add-expr-mrat (x l)
 199   (ratf (cons '(mplus) (cons (ratf x) l))))
 200
 201 (defun add-expr-taylor (x l)
 202   ($taylor (cons '(mplus) (cons x l))))
 203
 204 (defun add-expr-mlist (x l)
 205   (setq l (if (cdr l) (reduce 'addmx l) (car l)))
 206   (simplifya (cons (list 'mlist) (mapcar #'(lambda (s) (add x s)) (cdr l))) t))
 207
 208 ;; Simple demo showing how to define addition for a new object.
 209 ;; We could append simplification rules for intervals:
 210 ;;  (a) interval(a,a) --> a,
 211 ;;  (b) if p > q then interval(p,q) --> standardized empty interval?
 212
 213 (defun add-expr-interval (x l)
 214   (setq l (mapcar #'(lambda (s) `((mlist) ,@(cdr s))) l))
 215   (setq l (if (cdr l) (reduce #'addmx l) (car l)))
 216   (simplifya (cons (list '$interval) (mapcar #'(lambda (s) (add x s)) (cdr l))) t))
 217
 218 ;; Add an expression x to a list of matrices l. The Maxima function mxplusc
 219 ;; calls pls. Here is a version that doesn't call pls. I'm not sure I've captured all
 220 ;; features of mxplusc.
 221 (defun add-expr-matrix (x l)
 222   (setq l (if (cdr l) (reduce #'addmx l) (car l)))
 223   (cond ((and $listarith ($matrixp l))
 224           (let ((errset nil)
 225                             (errcatch t)
 226                         ($errormsg nil)
 227                             (*mdebug* nil))
 228                 (setq x (errset (simplifya (cons '($matrix) (cdr (add x ($args l)))) t)))
 229                 (if x (car x) (merror "Attempt to add noncomformable matrices"))))
 230         (t (list (get 'mplus 'msimpind) x l))))
 231
 232 ;; Return a + b, where a, b in {minf, inf, ind, und, infinity}. I should
 233 ;; extend this to allow zeroa and zerob (but I'm not sure zeroa and zerob
 234 ;; are supposed to be allowed outside the limit code). Internally, Maxima
 235 ;; uses prin-inf (see defint.lisp) to represent $inf. We could include
 236 ;; prin-inf as an extended real.
 237
 238 ;;We use a hashtable to represent the addition table--this should be easy
 239 ;; to extend and to modify.
 240 (defvar *extended-real-add-table* (make-hash-table :test #'equal :size 16))
 241
 242 (mapcar #'(lambda (a) (setf (gethash (list (first a) (second a)) *extended-real-add-table*) (third a)))
 243    (list (list '$minf '$minf '$minf)
 244          (list '$minf '$inf '$und)
 245          (list '$minf '$infinity '$und)
 246          (list '$minf '$und '$und)
 247          (list '$minf '$ind '$minf)
 248
 249          (list '$inf '$inf '$inf)
 250          (list '$inf '$infinity '$und)
 251          (list '$inf '$und '$und)
 252          (list '$inf '$ind '$inf)
 253
 254          (list '$infinity '$infinity '$und)
 255          (list '$infinity '$und '$und)
 256          (list '$infinity '$ind '$infinity)
 257
 258          (list '$ind '$ind '$ind)
 259          (list '$ind '$und '$und)
 260
 261          (list '$und '$und '$und)))
 262
 263 (defun add-extended-real(a b)
 264   (gethash (list a b) *extended-real-add-table* (gethash (list b a) *extended-real-add-table* '$und)))
 265
 266 ;; Add an expression x to a list of infinities. We do explicit number + extended real --> extended real,
 267 ;; but for a general expression XXX we do XXX + extended real --> nounform.
 268 (defun add-expr-infinities (x l)
 269   (setq l (if l (reduce #'add-extended-real l) (car l)))
 270   (if (mnump x) l (list (get 'mplus 'msimpind) x l)))
 271
 272 ;; The functions pls & plusin are parts of the standard simplus code. Let's issue
 273 ;; errors when these functions are called. Unfortunately, code in share\affine calls
 274 ;; pls directly, so the affine code will not run using altsimp.
 275 (defun pls (x out)
 276     (mtell "in pls; x = ~M,  out = ~M ~%" x out)
 277     (merror "Error: called pls ~%"))
 278
 279 (defun plusin (x fm)
 280    (mtell "in plusin; x = ~M,  fm = ~M ~%" x fm)
 281    (merror "Error: called plusin ~%"))
 282
 283 ;; I assumed that if a list of distinct members is sorted using great,
 284 ;; then it's still sorted after multiplying each list member by a nonzero
 285 ;; maxima number. I'm not sure this is true.
 286
 287 ;; If l has n summands, simplus calls great O(n log_2(n)) times. All
 288 ;; other spendy functions are called O(n) times. The standard simplus
 289 ;; function calls great O(n^2) times, I think.
 290
 291 ;; The binary64 value of %e.
 292 (defvar %e-float64 (exp 1.0d0))
 293
 294 (defun simplus (l w z)
 295   (declare (ignore w))
 296
 297   (let ((acc nil) (cf) (x) (num-sum 0) (do-over nil) (mequal-terms nil) (mrat-terms nil)
 298         (inf-terms nil) (matrix-terms nil) (mlist-terms nil) (taylor-terms nil) (interval-terms nil)
 299   (op) (atom-hash (make-hash-table :test #'eq :size 8)))
 300
 301   (setq l (margs l))
 302   ;; simplify and flatten
 303     (dolist (li l)
 304         (setq li (simplifya li z))
 305      ;; When numer is true, simplifya converts %pi & %phi to their binary64 value,
 306      ;; but only converts %e when both numer & %enumer are true. Here we convert
 307      ;; %e to its binary64 value.
 308      (when (and $numer (atom li) (eq li '$%e)) ; $%e --> 2.718...
 309         (setq li  %e-float64))
 310         (if (mplusp li) (setq acc (append acc (cdr li))) (push li acc)))
 311
 312     (setq l acc)
 313     (setq acc nil)
 314     (dolist (li l)
 315       (cond ((mnump li) (mincf num-sum li))
 316                   ;; factor out infrequent cases.
 317                 ((and (consp li) (consp (car li)) (member (caar li) '(mequal mrat $matrix mlist $interval)))
 318                         (setq op (caar li))
 319                         (cond ((eq op 'mequal)
 320                                       (push li mequal-terms))
 321                                     (($taylorp li)
 322                                     (push li taylor-terms))
 323                                     ((eq op 'mrat)
 324                                        (push li mrat-terms))
 325                                     ((eq op '$matrix)
 326                                        (push li matrix-terms))
 327                                     ((eq op '$interval)
 328                                        (push li interval-terms))
 329                                     ((eq op 'mlist)
 330                                       (if $listarith (push li mlist-terms) (push (convert-to-coeff-form li) acc)))))
 331
 332                     ;; Put non-infinite atoms into a hashtable; push infinite atoms into inf-terms.
 333                     ((atom li)
 334                         (if (and $use_extended_real_arithmetic (member li '($minf $inf $infinity $und $ind)))
 335                                      (push li inf-terms)
 336                             (progn
 337                                       (setq cf (gethash li atom-hash))
 338                                       (setf (gethash li atom-hash) (if cf (1+ cf) 1)))))
 339
 340                 (t (push (convert-to-coeff-form li) acc))))
 341
 342      ;; push atoms in the hashtable into the accumulator acc; sort acc.
 343     (maphash #'(lambda (cf a) (push (cons cf a) acc)) atom-hash)
 344     (setq l (sort acc 'great :key 'car))
 345
 346     ;; common term crunch: when the new coefficient is -1 or 1 (for example, 5*a - 4*a),
 347     ;; set the "do-over" flag to true. In this case, the sum needs to be re-simplified.
 348     ;; Without the do over flag, a + 5*a - 4*a --> a + a. Last I checked, the testsuite
 349     ;; does not test the do-over scheme.
 350
 351     (setq acc nil)
 352     (while l
 353       (setq x (pop l))
 354       (setq cf (cdr x))
 355       (setq x (car x))
 356       (while (and l (like x (caar l)))
 357         (mincf cf (cdr (pop l))))
 358         (if (and (or (eql cf 1) (eql cf -1)) (mplusp x)) (setq do-over t))
 359         (setq x (number-times-expr cf x))
 360         (cond ((mnump x) (mincf num-sum x))
 361                     ((not (mzerop x)) (push x acc))))
 362
 363     ;; Do x + 0 --> x, x + 0.0 --> x, and x + 0.0b0 --> x.
 364     (if (not (mzerop num-sum)) (push num-sum acc))
 365
 366     (setq acc
 367           (cond (do-over (simplifya `((mplus) ,@acc) nil))
 368                 ((null acc) num-sum)
 369                       ((null (cdr acc)) (car acc))
 370                       (t (cons '(mplus simp) acc))))
 371
 372     ;; special case dispatch
 373     (when mequal-terms
 374              (setq acc (add-expr-mequal acc mequal-terms)))
 375     (when taylor-terms
 376         (setq acc (add-expr-taylor acc taylor-terms)))
 377     (when mrat-terms
 378               (setq acc (add-expr-mrat acc mrat-terms)))
 379     (when mlist-terms
 380               (setq acc (add-expr-mlist acc mlist-terms)))
 381     (when interval-terms
 382         (setq acc (add-expr-interval acc interval-terms)))
 383     (when matrix-terms
 384               (setq acc (add-expr-matrix acc matrix-terms)))
 385     (when inf-terms
 386         (setq acc (add-expr-infinities acc inf-terms)))
 387     acc))