I had a look at the simplifier to find a solution to problems of the
following type:
BUG ID: 721575 2/sqrt(2) doesn\'t simplify
BUG ID 2029041 a*sqrt(2)/2 unsimplified
BUG ID 1923119 1/sqrt(8)-sqrt(8)/8
BUG ID 1927178 integrate(sin(t),t,%pi/4,3*%pi/4)
BUG ID: 1480562 2*a*2^k isn't simplified to a*2^(k+1)
Following problems are related to these bugs:
(1/2)*sqrt(2) is different from sqrt(2)/2 and
2/sqrt(2) is different from 1/sqrt(2)*(1/2)
1.
I have found the code which will give correct results for all examples
from above and related problems.
2.
The testsuite has only 19 failures with the improved code. Most failures
are due to more simple results with the improved code or a different
order of the terms (more like the expected canonical order). In some
cases the expressions are not fully simplified. These failures might be
simplification errors which have vanished by accident.
Only one problem fails completely. We get again a infinite loop for the
integral (rtestint.mac): integrate(1/(1-x^5), x, 0, inf).
3.
In the share testsuite I have got 5 failures. 4 failures are due to a
different order of terms or a different simplification result. One
failure I have not understand:
zeromatrixp(m1 . transpose(m1) - m . transpose(m)) does no longer work
as expected (This might be an error in MTIMES).
4.
To get the improved code, I have cut out the handling of numbers
completely from the routine TMS and moved it to the routine TIMESIN.
5.
There is one change in the behavior of Maxima. Expressions like e.g.
3*sqrt(3) or sqrt(3)*3 simplifies for all numbers.
So we get always 3*sqrt(3) --> 3^(3/2).
The biggest problem is, that the code we have, distribute the handling
of numbers over the two routines TMS and TIMESIN. The simplification of
numbers in TMS makes it impossible to do a more consistent
simplification in the routine TIMESIN. When we cut out the code in TMS,
we can improve TIMESIN consistently to handle all desired cases. Remark:
It seems to me that TIMESIN never was designed to handle numbers, but if
we add it to TMS a lot of work has to be doubled.
So what do you think? Should we risk an implementation of a new function
TIMESIN. Any suggestions how to go on?
This is the improved routine TIMESIN which handles numbers too and does
the desired simplifications to overcome the bugs (I should do even more
comments in it).
(defvar *debug-timesin* nil)
(defun timesin (x y w) ; Multiply X^W into Y
(prog (fm temp z check u expo)
(if (mexptp x) (setq check x))
top
;; Prepare the factor X^W and store the factor in temp.
(cond ((equal w 1)
(setq temp x))
(t
(setq temp (cons '(mexpt)
(if check
(list (cadr x) (mult (caddr x) w))
(list x w))))
(if (and (not timesinp) (not (eq x '$%i)))
(let ((timesinp t))
(setq temp (simplifya temp t))))))
(setq x (if (mexptp temp)
(cdr temp)
(list temp 1)))
(setq w (cadr x)
fm y)
start
(when *debug-timesin*
(format t "~& in TIMESIN START:")
(format t "~& x = ~A~%" x)
(format t "~& w = ~A~%" w)
(format t "~& fm = ~A~%" fm)
(format t "~& temp = ~A~%" temp)
(format t "~& (car fm) = ~A~%" (car fm))
(format t "~& (car y) = ~A~%" (car y)))
(cond ((null (cdr fm))
(when *debug-timesin* (format t "~& GOTO LESS.~%"))
(go less))
((mexptp (cadr fm))
;; A mexpt expression in the list of products
(when *debug-timesin*
(format t " in TIMESIN: cond mexptp for ~A~%" (cadr fm)))
(cond ((alike1 (car x) (cadadr fm))
;; Multiply x^r*x^w
(cond ((zerop1 (setq w (plsk (caddr (cadr fm)) w)))
;; r+w=0. The factor vanishs. Remove the
factor.
(return (rplacd fm (cddr fm))))
((and (mnump w)
(or (mnump (car x))
(eq (car x) '$%i)))
;; Base of the factor is a Maxima number or %i
;; and the power is Maxima number.
(when *debug-timesin*
(format t "~& in TIMES: mexpt found a
number~%"))
(rplacd fm (cddr fm))
(cond ((mnump (setq x (if (mnump (car x))
(exptrl (car x) w)
(power (car x) w))))
(return (rplaca y (timesk (car y) x))))
((mtimesp x)
(go times))
(t
(setq temp x
x (if (mexptp x) (cdr x) (list x
1)))
(setq w (cadr x)
fm y)
(go start))))
((maxima-constantp (car x))
;; A constant.
(when *debug-timesin*
(format t
"~& in TIMESIN: maxima-constantp. Go
const~%"))
;; A constant. Remove the factor from the
list.
;; Set x and check. w is already set. Start
again.
(rplacd fm (cddr fm))
(setq x (car x) check nil)
(go top))
((onep1 w)
(cond ((mtimesp (car x))
;; A base which is a mtimes expression.
;; Remove the factor from the lists of
products.
(rplacd fm (cddr fm))
;; Multiply the factors of the base
with
;; the list of all remaining products.
(setq rulesw t)
(return (muln (nconc y (cdar x)) t)))
(t (return (rplaca (cdr fm) (car x))))))
(t
(go spcheck))))
;; At this place we have to add code for a rational
number
;; as a factor to the list of products.
((and (onep1 w)
(or (ratnump (car x))
(and (integerp (car x))
(not (onep (car x))))))
;; Multiplying a^k * rational.
(when *debug-timesin*
(format t "~& in times: found ratnump~%"))
(let ((numerator (if (integerp (car x))
(car x)
(second (car x))))
(denom (if (integerp (car x))
1
(third (car x)))))
(setf expo (exponent-of numerator (second (cadr
fm))))
(when *debug-timesin*
(format t "~& : numerator = ~A~%" numerator)
(format t "~& : denomerator = ~A~%" denom)
(format t "~& : exponent = ~A~%" expo))
(when expo
;; We have a^m*a^k.
(setq temp (power (second (cadr fm))
(add (third (cadr fm)) expo)))
(when *debug-timesin*
(format t "~& : temp = ~A~%" temp))
;; Set fm to have 1/denom term.
(setf fm
(rplaca fm
(mul (car fm)
(div (div numerator
(power (second (cadr
fm))
expo))
denom))))
;; Add in the a^(m+k) term.
(rplacd fm (cddr fm))
(rplacd fm (cons temp (cdr fm)))
(return (cdr fm)))
(setf expo (exponent-of (inv denom) (second (cadr
fm))))
(when *debug-timesin*
(format t "~& : numerator = ~A~%" numerator)
(format t "~& : denomerator = ~A~%" denom)
(format t "~& : exponent = ~A~%" expo))
(when expo
;; We have a^(-m)*a^k.
(setq temp (power (second (cadr fm))
(add (third (cadr fm)) expo)))
(when *debug-timesin*
(format t "~& : temp = ~A~%" temp))
;; Set fm to have the numerator term.
(setf fm (rplaca fm
(mul (car fm)
(div numerator
(div denom
(power (second
(cadr fm))
(-
expo)))))))
;; Add in the a^(k-m) term.
(rplacd fm (cddr fm))
(rplacd fm (cons temp (cdr fm)))
(return (cdr fm)))
;; The rational doesn't contain any (simple) powers
of
;; the exponential term. We're done. (This is
;; basically the T case below.)
nil))
((or (maxima-constantp (car x))
(maxima-constantp (cadadr fm)))
(if (great temp (cadr fm))
(go gr)))
((great (car x) (cadadr fm))
(go gr)))
(go less))
((alike1 (car x) (cadr fm))
;;(format t "start: alike1 go equ~%")
(go equ))
((mnump temp)
(when *debug-timesin*
(format t "~& in TIMESIN start: mnump. GOTO START~%"))
(setq fm (cdr fm))
(go start))
((maxima-constantp (car x))
(when *debug-timesin*
(format t "~& in TIMESIN start: maxima-constantp~%"))
(when (great temp (cadr fm))
(when *debug-timesin* (format t "~& GOTO GREAT~&"))
(go gr)))
((great (car x) (cadr fm))
(when *debug-timesin*
(format t "~& in TIMESIN start: GOTO GREAT~&"))
(go gr)))
less
(cond ((mnump temp)
;; Multiply a number into the list of products.
(return (rplaca y (timesk (car y) temp))))
((and (eq (car x) '$%i)
(fixnump w))
(go %i))
((and (eq (car x) '$%e)
$numer
(integerp w))
(return (rplaca y (timesk (car y) (exp w)))))
((and (onep1 w)
(not (constant (car x))))
(go less1))
;; At this point we will insert a mexpt expression,
;; but first we look at the car of the list of products and
;; modify the expression if we found a rational number.
((and (mexptp temp)
(not (onep1 (car y)))
(or (integerp (car y))
(ratnump (car y))))
;; Multiplying a^k * rational.
(when *debug-timesin*
(format t "~& in times: found ratnump in (car y)~%"))
(let ((numerator (if (integerp (car y))
(car y)
(second (car y))))
(denom (if (integerp (car y))
1
(third (car y)))))
(setf expo (exponent-of numerator (car x)))
(when *debug-timesin*
(format t "~& : numerator = ~A~%" numerator)
(format t "~& : denomerator = ~A~%" denom)
(format t "~& : exponent = ~A~%" expo))
(when expo
;; We have a^m*a^k.
(setq temp (power (car x)
(add (cadr x) expo)))
(when *debug-timesin*
(format t "~& : temp = ~A~%" temp))
;; Set fm to have 1/denom term.
(setf y
(rplaca y
(div (div numerator
(power (car x) expo))
denom)))
;; Add in the a^(m+k) term.
(rplacd fm (cons temp (cdr fm)))
(return (cdr fm)))
(setf expo (exponent-of (inv denom) (car x)))
(when *debug-timesin*
(format t "~& : numerator = ~A~%" numerator)
(format t "~& : denomerator = ~A~%" denom)
(format t "~& : exponent = ~A~%" expo))
(when expo
;; We have a^(-m)*a^k.
(setq temp (power (car x)
(add (cadr x) expo)))
(when *debug-timesin*
(format t "~& : temp = ~A~%" temp))
;; Set fm to have the numerator term.
(setf y (rplaca y
(div numerator
(div denom
(power (car x)
(- expo))))))
;; Add in the a^(k-m) term.
(rplacd fm (cons temp (cdr fm)))
(return (cdr fm)))
;; The rational doesn't contain any (simple) powers of
;; the exponential term. We're done. (This is
;; basically the T case below.)
(return (cdr (rplacd fm (cons temp (cdr fm)))))
))
((and (maxima-constantp (car x))
(do ((l (cdr fm) (cdr l)))
((null (cdr l)))
(when (and (mexptp (cadr l))
(alike1 (car x) (cadadr l)))
(setq fm l)
(return t))))
(go start))
((or (and (mnump (car x))
(mnump w))
(and (eq (car x) '$%e)
$%emode
(setq u (%especial w))))
(setq x (cond (u)
((alike (cdr check) x)
check)
(t
(exptrl (car x) w))))
(cond ((mnump x)
(return (rplaca y (timesk (car y) x))))
((mtimesp x)
(go times))
((mexptp x)
(return (cdr (rplacd fm (cons x (cdr fm))))))
(t
(setq temp x
x (list x 1)
w 1
fm y)
(go start))))
((onep1 w)
(go less1))
((ratnump (car fm))
(when *debug-timesin*
(merror "We never should reach this part of TIMESIN."))
;; Multiplying a^k * rational.
;;(format t "timesin a^k * rat~%")
(let ((numerator (second (car fm)))
(denom (third (car fm))))
;; (format t "numerator = ~A~%" numerator)
(setf expo (exponent-of numerator (car x)))
(when expo
;; We have a^m*a^k.
(setq temp (list '(mexpt) (car x) (add w expo)))
;; Set fm to have 1/denom term.
(setf fm (rplaca fm (div (div numerator (power (car x)
expo))
denom)))
;; Add in the a^(m+k) term.
(rplacd fm (cons temp (cdr fm)))
(return (cdr fm)))
(setf expo (exponent-of (inv denom) (car x)))
(when expo
;; We have a^(-m)*a^k.
(setq temp (list '(mexpt) (car x) (add w expo)))
;; Set fm to have the numerator term.
(setf fm (rplaca fm (div numerator
(div denom
(power (car x) (-
expo))))))
;; Add in the a^(k-m) term.
(rplacd fm (cons temp (cdr fm)))
(return (cdr fm)))
;; The rational doesn't contain any (simple) powers of
;; the exponential term. We're done. (This is
;; basically the T case below.)
(setq temp (list '(mexpt) (car x) w))
(setq temp (eqtest temp (or check '((foo)))))
(return (cdr (rplacd fm (cons temp (cdr fm)))))
))
((setf expo (exponent-of (car fm) (car x)))
;; Got something like a*a^k, where a is a number.
;;(format t "go a*a^k~%")
(setq temp (list '(mexpt) (car x) (add w expo)))
(setf fm (rplaca fm (div (car fm) (power (car x) expo))))
(return (cdr (rplacd fm (cons temp (cdr fm))))))
(t
;;(format t "default less cond~%")
(setq temp (list '(mexpt) (car x) w))
(setq temp (eqtest temp (or check '((foo)))))
;;(format t "temp = ~A~%" temp)
;;(format t "fm = ~A~%" fm)
(return (cdr (rplacd fm (cons temp (cdr fm)))))))
less1
(return (cdr (rplacd fm (cons (car x) (cdr fm)))))
gr
(setq fm (cdr fm))
(go start)
equ
(cond ((and (eq (car x) '$%i) (equal w 1))
(rplacd fm (cddr fm))
(return (rplaca y (timesk -1 (car y)))))
((zerop1 (setq w (plsk 1 w)))
(return (rplacd fm (cddr fm))))
((and (mnump (car x)) (mnump w))
(return (rplaca (cdr fm) (exptrl (car x) w))))
((maxima-constantp (car x))
(go const)))
spcheck
(setq z (list '(mexpt) (car x) w))
(cond ((alike1 (setq x (simplifya z t)) z)
(return (rplaca (cdr fm) x)))
(t
(rplacd fm (cddr fm))
(setq rulesw t)
(return (muln (cons x y) t))))
const
(rplacd fm (cddr fm))
(setq x (car x) check nil)
(go top)
times
(setq z (tms x 1 (setq temp (cons '(mtimes) y))))
(return (cond ((eq z temp)
(cdr z))
(t
(setq rulesw t) z)))
%i
(if (minusp (setq w (rem w 4)))
(incf w 4))
(return (cond ((zerop w)
fm)
((= w 2)
(rplaca y (timesk -1 (car y))))
((= w 3)
(rplaca y (timesk -1 (car y)))
(rplacd fm (cons '$%i (cdr fm))))
(t
(rplacd fm (cons '$%i (cdr fm))))))))
Dieter Kaiser