I have reported work on the functions TMS und TIMESIN to improve the
simplification of the multiplication. I have improved the code further.
The following shows what we have for 17 examples. Only 5 examples
simplifies correctly. 12 examples have problems.
********************** Problem 1 ***************
Input:
2
-------
sqrt(2)
Result:
2
-------
sqrt(2)
This differed from the expected result:
sqrt(2)
********************** Problem 2 ***************
Input:
sqrt(2)
-------
2
Result:
sqrt(2)
-------
2
This differed from the expected result:
1
-------
sqrt(2)
********************** Problem 3 ***************
Input:
sqrt(2)
-------
2
Result:
1
-------
sqrt(2)
... Which was correct.
********************** Problem 4 ***************
Input:
a sqrt(2)
---------
2
Result:
- 1/2
2 a
... Which was correct.
********************** Problem 5 ***************
Input:
1 sqrt(8)
------- - -------
sqrt(8) 8
Result:
1 1
--------- - ----
2 sqrt(2) 3/2
2
This differed from the expected result:
0
********************** Problem 6 ***************
Input:
%pi 3 %pi
integrate(sin(t), t, ---, -----)
4 4
Result:
sqrt(2) 1
------- + -------
2 sqrt(2)
This differed from the expected result:
1 1
------- + -------
sqrt(2) sqrt(2)
********************** Problem 7 ***************
Input:
k
2 a 2
Result:
k
2 a 2
This differed from the expected result:
1 + k
a 2
********************** Problem 8 ***************
Input:
k
a 2 2
Result:
k + 1
a 2
... Which was correct.
********************** Problem 9 ***************
Input:
! sqrt(3) %i 1 !
! ---------- - - !
! %i 1 1/3 - sqrt(3) %i 1 2 2 !
!(--------- - -) (------------ - -) + --------------------!
! 6 sqrt(3) 6 2 2 %i 1 1/3!
! 3 (--------- - -) !
! 6 sqrt(3) 6 !
Result:
- 1/2 5 %pi 5 %pi
3 sin(-----) sin(-----)
18 18 2
sqrt((----------------- - ----------)
2 2 sqrt(3)
5 %pi - 1/2 5
%pi
cos(-----) 3
cos(-----)
5 %pi 18
18 2
+ (sin(-----) - ---------- -
-----------------) )
18 2 sqrt(3) 2
This differed from the expected result:
5 %pi
cos(-----)
5 %pi 18
sin(-----) - ----------
18 sqrt(3)
********************** Problem 10 ***************
Input:
k k k k
a 2 b c 2 b 3 3 a 2 3 a 2 3 c
Result:
3 2 2 2 k + 1 2 k + 2
2 a b c 2 3
This differed from the expected result:
3 2 2 2 + 2 k 2 + 2 k
a b c 2 3
********************** Problem 11 ***************
Input:
1 k k
(-) 2 3
2
---------
3
Result:
k - 1 k
2 3
---------
3
This differed from the expected result:
k - 1 k - 1
2 3
********************** Problem 12 ***************
Input:
k k
(2 3 1) 1
-----------
2
-----------
3
Result:
k - 1 k - 1
2 3
... Which was correct.
********************** Problem 13 ***************
Input:
1 k k
(-) 2 3
2
---------
- 3
Result:
k - 1 k
2 3
- ---------
3
This differed from the expected result:
k - 1 k - 1
- 2 3
********************** Problem 14 ***************
Input:
1 k
(-) a b 2 1
2 k
(------------) (- 1) 3
- 3
-----------------------
2 3
Result:
k - 1 k
a b 2 3
-------------
18
This differed from the expected result:
k - 2 k - 2
a b 2 3
********************** Problem 15 ***************
Input:
k k
12 2 3
Result:
k + 1 k + 2
3 2
... Which was correct.
********************** Problem 16 ***************
Input:
k k
2 3 12
Result:
k + 2 k
(3 2 ) 3
This differed from the expected result:
2 + k 1 + k
2 3
********************** Problem 17 ***************
Input:
k k 2
2 (2 3 12) 2
Result:
2 k + 2 2 (k + 2) + 1
2 3 2
This differed from the expected result:
2 + 2 (2 + k) 2 (1 + k)
2 3
5/17 tests passed.
The following 12 problems failed: (1 2 5 6 7 9 10 11 13 14 16 17)
(%o1) [rtest_simptimes.mac, 1, 2, 5, 6, 7, 9, 10, 11, 13, 14, 16, 17]
The following shows the results with the improved code. All example
including the examples from bug the bug reports related to this topic
will work:
********************** Problem 1 ***************
Input:
2
-------
sqrt(2)
Result:
sqrt(2)
... Which was correct.
********************** Problem 2 ***************
Input:
sqrt(2)
-------
2
Result:
1
-------
sqrt(2)
... Which was correct.
********************** Problem 3 ***************
Input:
sqrt(2)
-------
2
Result:
1
-------
sqrt(2)
... Which was correct.
********************** Problem 4 ***************
Input:
a sqrt(2)
---------
2
Result:
a
-------
sqrt(2)
... Which was correct.
********************** Problem 5 ***************
Input:
1 sqrt(8)
------- - -------
sqrt(8) 8
Result:
0
... Which was correct.
********************** Problem 6 ***************
Input:
%pi 3 %pi
integrate(sin(t), t, ---, -----)
4 4
Result:
2
-------
sqrt(2)
... Which was correct.
********************** Problem 7 ***************
Input:
k
2 a 2
Result:
k + 1
a 2
... Which was correct.
********************** Problem 8 ***************
Input:
k
a 2 2
Result:
k + 1
a 2
... Which was correct.
********************** Problem 9 ***************
Input:
! sqrt(3) %i 1 !
! ---------- - - !
! %i 1 1/3 - sqrt(3) %i 1 2 2 !
!(--------- - -) (------------ - -) + --------------------!
! 6 sqrt(3) 6 2 2 %i 1 1/3!
! 3 (--------- - -) !
! 6 sqrt(3) 6 !
Result:
5 %pi
cos(-----)
5 %pi 18
sin(-----) - ----------
18 sqrt(3)
... Which was correct.
********************** Problem 10 ***************
Input:
k k k k
a 2 b c 2 b 3 3 a 2 3 a 2 3 c
Result:
3 2 2 2 k + 2 2 k + 2
a b c 2 3
... Which was correct.
********************** Problem 11 ***************
Input:
1 k k
(-) 2 3
2
---------
3
Result:
k - 1 k - 1
2 3
... Which was correct.
********************** Problem 12 ***************
Input:
k k
(2 3 1) 1
-----------
2
-----------
3
Result:
k - 1 k - 1
2 3
... Which was correct.
********************** Problem 13 ***************
Input:
1 k k
(-) 2 3
2
---------
- 3
Result:
k - 1 k - 1
- 2 3
... Which was correct.
********************** Problem 14 ***************
Input:
1 k
(-) a b 2 1
2 k
(------------) (- 1) 3
- 3
-----------------------
2 3
Result:
k - 2 k - 2
a b 2 3
... Which was correct.
********************** Problem 15 ***************
Input:
k k
12 2 3
Result:
k + 1 k + 2
3 2
... Which was correct.
********************** Problem 16 ***************
Input:
k k
2 3 12
Result:
k + 1 k + 2
3 2
... Which was correct.
********************** Problem 17 ***************
Input:
k k 2
2 (2 3 12) 2
Result:
2 (k + 1) 2 (k + 2) + 2
3 2
... Which was correct.
17/17 tests passed.
(%o3) []
The examples show that the simplification will work commutative, e.g.
12*2^k*3^k gives the same result as 2^k*3^k*12. Therefore a lot of known
problems will vanish.
Most impressive is the problem 9. The complicated expression will
simplify to an very simple expression without any other help.
There is one change:
Expressions like 3*sqrt(3) simplifies always to 3^(3/2). This behavior
is necessary to get consistent simplifications in much more complicated
expressions too.
There is one problem (not a new problem):
Expression like 2*sqrt(2)+2*sqrt(2) simplifies to 2*2^(3/2) and not
2^(5/2). The reason is that the routine PLUSIN does not simplify
products accordingly. This can be improved too.
The testsuite and the share_testsuite give only some errors. In almost
all cases the expressions simplify differently.
I have appended the modified code of TMS and TIMESIN to this posting.
Dieter Kaiser
-------------- next part --------------
;;;; Created on 2009-06-14 15:53:01
(defvar *debug-timesin* nil)
(defun tms (factor power product &aux tem)
(let ((rulesw nil)
(z nil))
(when (mplusp product) (setq product (list '(mtimes simp) product)))
(cond ((zerop1 factor)
(cond ((mnegp power)
(if errorsw
(throw 'errorsw t)
(merror "Division by 0")))
(t factor)))
((and (null product)
(or (and (mtimesp factor) (equal power 1))
(and (setq product (list '(mtimes) 1)) nil)))
(setq tem (append '((mtimes)) (if (mnump (cadr factor)) nil '(1))
(cdr factor) nil))
(if (= (length tem) 1)
(setq tem (copy-list tem))
tem))
; We do not handle a numeric factor in tms but in timesin
; ((mnump factor)
; ;; We need to do something better here. Look through
; ;; product to see if there are any terms of the form
; ;; factor^k, and adjust the exponent.
;
; ;;(format t "tms mnump factor = ~A~%" factor)
; ;;(format t "tms product = ~A~%" product)
; (let ((expo nil))
; (do ((p (cdr product) (cdr p)))
; ((or (null p) expo))
; ;;(format t "p = ~A~%" p)
; (when (and (mexptp (car p))
; (integerp (second (car p)))
; ;;(integerp factor)
; (setf expo (exponent-of factor (second (car p)))))
; (let* ((q (div factor (power (second (car p)) expo)))
; (temp (mul q (list '(mexpt)
; (second (car p))
; (add expo (third (car p)))))))
; ;;(format t "temp = ~A~%" temp)
; ;;(format t "p = ~A~%" p)
; ;;(format t "cdr p = ~A~%" (cdr p))
; (setf temp (append (list temp) (cdr p)))
; ;;(format t "new temp = ~A~%" temp)
; ;;(rplaca p temp)
; (rplacd p (cdr temp))
; (rplaca p (car temp))
; ;;(format t "mod p = ~A~%" p)
; )))
; (unless expo
; (rplaca (cdr product) (timesk (cadr product) (expta factor power)))))
; (if (and (mtimesp product)
; (mtimesp (cadr product)))
; (rplacd product (append (cdadr product) (cddr product))))
; product)
((mtimesp factor)
; We do not hanlde a numeric factor in tms but in timesin
; (when (mnump (cadr factor))
; (setq factor (cdr factor))
; (rplaca (cdr product)
; (timesk (cadr product) (expta (car factor) power))))
;; Multiply each factor into procduct
(do ((factor-list (cdr factor) (cdr factor-list)))
((or (null factor-list) (zerop1 product)) product)
(setq z (timesin (car factor-list) (cdr product) power))
(when rulesw
(setq rulesw nil)
(setq product (tms-format-product z)))))
(t
;; Multiply factor^power into product
(setq z (timesin factor (cdr product) power))
(if rulesw
(tms-format-product z)
product)))))
(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 initialize the work of timesin
(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)
(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)))
start
;; Go through the list of terms in fm and look what is to do
(cond ((null (cdr fm))
;; The list of terms is empty. The loop is finshed.
(when *debug-timesin* (format t "~& GOTO LESS~%"))
(go less))
((and (mnump temp)
(not (or (integerp temp)
(ratnump temp))))
;; Stop the loop for a float or bigfloat number.
(go less))
((mexptp (cadr fm))
(when *debug-timesin*
(format t "~& in START (cond mexptp) with ~A~%" (cadr fm)))
(cond ((alike1 (car x) (cadadr fm))
(when *debug-timesin*
(format t "~& in START (cond alike1)~%"))
(cond ((zerop1 (setq w (plsk (caddr (cadr fm)) w)))
(go del))
((and (mnump w)
(or (mnump (car x))
(eq (car x) '$%i)))
(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))
(go const))
((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))))
(sgn (signum numerator)))
(setf expo (exponent-of (abs numerator) (second (cadr fm))))
(when *debug-timesin*
(format t "~& : numerator = ~A~%" numerator)
(format t "~& : denomerator = ~A~%" denom)
(format t "~& : exponent = ~A~%" expo)
(format t "~& : sign = ~A~%" sgn))
(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.
(setq x (mul sgn
(car y)
(div (div (mul sgn numerator)
(power (second (cadr fm))
expo))
denom)))
(setf y (rplaca y 1))
; (setf y
; (rplaca y
; (mul sgn
; (car fm)
; (div (div (mul sgn numerator)
; (power (second (cadr fm))
; expo))
; denom))))
;; Add in the a^(m+k) term.
(rplacd fm (cddr fm))
(rplacd fm (cons temp (cdr fm)))
(setq temp x
x (list x 1)
w 1
fm y)
(go start)
; (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.
(setq x (mul (car y)
(div numerator
(div denom
(power (second (cadr fm))
(- expo))))))
(setf y (rplaca y 1))
; (setf y (rplaca y
; (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)))
(setq temp x
x (list x 1)
w 1
fm y)
(go start)
; (return (cdr fm))
)
;; Next term in list of products.
(setq fm (cdr fm))
(go start)
))
((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 a number goto start and look in the next term.
(when *debug-timesin*
(format t "~& in TIMESIN start: mnump. GOTO START~%"))
(setq fm (cdr fm))
(go start))
((maxima-constantp (car x))
;; (progn
;; (format t "start: maxima-constantp~%")
;; (format t " temp = ~A~%" temp))
(when (great temp (cadr fm))
;;(format t " go gr~%")
(go gr)))
((great (car x) (cadr fm))
;;(format t "greater, go gr~%")
(go gr)))
less
(cond ((mnump temp)
;; Multiply a number into the list of products.
(when *debug-timesin*
(format t "~& TIMESIN in less: (mnump temp) with ~A~%" temp))
(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 x^w * rational or integer.
(when *debug-timesin*
(format t "~& TIMESIN in less: (mexptp temp) with ~A~%" temp))
(let* ((numerator (if (integerp (car y))
(car y)
(second (car y))))
(denom (if (integerp (car y))
1
(third (car y))))
(sgn (signum numerator)))
(setf expo (exponent-of (abs numerator) (car x)))
(when *debug-timesin*
(format t "~& : numerator = ~A~%" numerator)
(format t "~& : denomerator = ~A~%" denom)
(format t "~& : exponent = ~A~%" expo)
(format t "~& : sign = ~A~%" sgn))
(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.
(setq x (div (div numerator
(power (car x) expo))
denom))
(setf y (rplaca y 1))
; (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))
(setq temp x
x (list x 1)
w 1
fm y)
(go start))
(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.
(setq x (div numerator
(div denom
(power (car x)
(- expo)))))
(setf y (rplaca y 1))
; (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))
(setq temp x
x (list x 1)
w 1
fm y)
(go start))
;; 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))
(merror "timesin: We should never reach this point")
;; 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)))
(merror "timesin: We should never reach this point")
;; 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)))
(go del))
((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)))
del
(return (rplacd fm (cddr fm)))
%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))))))))