We have a bug report ID: 3247367 "expand returns unsimplified". Two
examples are:
(%i1) (1-sqrt(5))^3-4*(1-sqrt(5))^2+8, expand;
(%o1) 5^(3/2)-5^(3/2)
(%i2) 1/sqrt(2)+1/sqrt(2)+1/sqrt(2);
(%o2) sqrt(2)+1/sqrt(2)
The last result should be 3/sqrt(2). The problem is that we have
simplifications for powers of integers like 2/sqrt(2) -> sqrt(2). This
type of simplifications causes problems in the routine plusin.
The general expression we have to simplify more complete is
v*a^(x+n) + w*a^(x+m) -> (v*a^n + w*a^m) * a^x
The symbols v, w, n, m, a must be integers. x is any expression.
Because of this the factor (v*a^n+w*a^m) is an integer too.
I have done an implementation of this type of simplification. The
algorithm is a bit tricky, when x+n and x+m are rational numbers. The
following shows the complete routine plusin. The additional code can be
further optimized, but at first I have tried to get a working algorithm.
; --------------------------------------------------------------------
(defun plusin (x fm)
(prog (x1 flag check w xnew n m a x2 v)
(setq w 1)
(setq v 1)
(cond ((mtimesp x)
(setq check x)
(if (mnump (cadr x)) (setq w (cadr x) x (cddr x))
(setq x (cdr x))))
(t (setq x (ncons x))))
(setq x1 (if (null (cdr x)) (car x) (cons '(mtimes) x))
xnew (list* '(mtimes) w x))
start
(cond ((null (cdr fm)))
((and (alike1 x1 (cadr fm)) (null (cdr x)))
(go equ))
;; Implement the simplification of
;; (v*a^(x+n)+w*a^(x+m) -> (v*a^n+w*a^m)*a^x
((and (or (mexptp (setq x2 (cadr fm)))
(and (mtimesp x2)
(null (cadddr x2))
(integerp (setq v (cadr x2)))
(mexptp (setq x2 (caddr x2)))))
(integerp (setq a (cadr x2)))
(mexptp x1)
(equal a (cadr x1))
(integerp (sub (caddr x2) (caddr x1))))
(setq n (if (and (mplusp (caddr x2))
(mnump (cadr (caddr x2))))
(cadr (caddr x2))
(if (mnump (caddr x2))
(caddr x2)
0)))
(setq m (if (and (mplusp (caddr x1))
(mnump (cadr (caddr x1))))
(cadr (caddr x1))
(if (mnump (caddr x1))
(caddr x1)
0)))
(cond ((integerp n)
(setq x1 (mul (addk (mul v (power a n))
(mul w (power a m)))
(muln (cons (power a (- m)) x) t)))
(go equt2))
(t
(multiple-value-bind (n1 d1)
(truncate (num1 n) (denom1 n))
(multiple-value-bind (n2 d2)
(truncate (num1 m) (denom1 m))
(cond ((equal d1 d2)
(setq x1 (mul (add (mul v (power a n1))
(mul w (power a n2)))
(power a (div d1 (denom1
n)))))
(go equt2))
((equal d2 -1)
(setq n1 (add n1 (div (sub d1 d2) (denom1
n))))
(setq d1 d2)
(setq x1 (mul (add (mul v (power a n1))
(mul w (power a n2)))
(power a (div d1 (denom1
n)))))
(go equt2))
((equal d1 -1)
(setq n2 (add n2 (div (sub d2 d1) (denom1
n))))
(setq d2 d1)
(setq x1 (mul (add (mul v (power a n1))
(mul w (power a n2)))
(power a (div d1 (denom1
n)))))
(go equt2))
(t (merror "Internal error in simplus.")))
)))))
((mtimesp (cadr fm))
(cond ((alike1 x1 (cadr fm))
(go equt))
((and (mnump (cadadr fm)) (alike x (cddadr fm)))
(setq flag t) ; found common factor
(go equt))
((great xnew (cadr fm)) (go gr))))
((great x1 (cadr fm)) (go gr)))
(setq xnew (eqtest (testt xnew) (or check '((foo)))))
(return (cdr (rplacd fm (cons xnew (cdr fm)))))
gr
(setq fm (cdr fm))
(go start)
equ
(when *debug-simplus*
(format t "~&in PLUSIN - equ~%"))
(rplaca (cdr fm)
(if (equal w -1)
(list* '(mtimes simp) 0 x)
;; Call muln to get a simplified product.
(if (mtimesp (setq x1 (muln (cons (addk 1 w) x) t)))
(testtneg x1)
x1)))
del
(cond ((not (mtimesp (cadr fm)))
(go check))
((onep (cadadr fm))
;; Do this simplification for an integer 1, not for 1.0 and
1.0b0
(rplacd (cadr fm) (cddadr fm))
(return (cdr fm)))
((not (zerop1 (cadadr fm)))
(return (cdr fm)))
;; Handle the multiplication with a zero.
((and (or (not $listarith) (not $doallmxops))
(mxorlistp (caddr (cadr fm))))
(return (rplacd fm
(cons (constmx 0 (caddr (cadr fm))) (cddr
fm))))))
;; (cadadr fm) is zero. If the first term of fm is a number,
;; add it to preserve the type.
(when (mnump (car fm))
(rplaca fm (addk (car fm) (cadadr fm))))
(return (rplacd fm (cddr fm)))
equt
(when *debug-simplus*
(format t "~&in PLUSIN - equt~%"))
;; Call muln to get a simplified product.
(setq x1 (muln (cons (addk w (if flag (cadadr fm) 1)) x) t))
equt2
(rplaca (cdr fm)
(if (zerop1 x1)
(list* '(mtimes) x1 x)
(if (mtimesp x1) (testtneg x1) x1)))
(if (not (mtimesp (cadr fm))) (go check))
(when (and (onep (cadadr fm)) flag (null (cdddr (cadr fm))))
;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
(rplaca (cdr fm) (caddr (cadr fm))) (go check))
(go del)
check
(if (mplusp (cadr fm)) (setq *plusflag* t)) ; A nested mplus
expression
(return (cdr fm))))
; --------------------------------------------------------------------
With this implementation I get the following results. The examples give
the expected results.
(%i1) (1-sqrt(5))^3-4*(1-sqrt(5))^2+8, expand;
(%o1) 0
(%i2) 1/sqrt(2)+1/sqrt(2)+1/sqrt(2);
(%o2) 3/sqrt(2)
(%i3) 2^(9/5)+2^(4/5);
(%o3) 3*2^(4/5)
(%i4) 3*sqrt(2)+2*sqrt(2);
(%o4) 5*sqrt(2)
(%i5) 2*sqrt(2)+3*sqrt(2);
(%o5) 5*sqrt(2)
(%i6) (1-sqrt(5))^3, expand;
(%o6) 16-8*sqrt(5)
(%i7) p : z^3-2^(3/2)*%i*z^2-4*z^2+2^(5/2)*%i*z+2*z;
(%o7) z^3-2^(3/2)*%i*z^2-4*z^2+2^(5/2)*%i*z+2*z
(%i8) divide(p, (z-2-sqrt(2)*%i),z);
(%o8) [z^2+(-sqrt(2)*%i-2)*z,0]
An example with a symbolic power:
(%i9) 2^a + 3*2^(a+1);
(%o9) 7*2^a
The testsuite and share_testsuite have no problems. I am looking more
systematic for further examples and possible errors of the code.
Dieter Kaiser