;;; -*- Mode:Common-Lisp; Package:mma; Base:10 -*-
(declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
(eval-when (compile) (load "mma")) ;; need the symbols like Plus
(in-package :mma)
;; RatExpand[(a5+b5)*(c5+d5)]
;; old version changed version
;; (a5 + b5) c5 + (a5 + b5) d5 a5 c5 + b5 c5 + a5 d5 + b5 d5
;; extracted from simp1.lisp
(defun intol* (a b)
(cond ((eql a 0) 0)
((eq a 1) b)
((eq b 1) a)
;; additional line to promote distributed form
((and *expand* (IsPlus a))
(ucons 'Plus (mapcar #'(lambda(r)(intol* r b)) (cdr a))))
(t (ucons 'Times
(uappend (intol*chk a) (intol*chk b))))))
;; extracted from rat1.lisp
(defun fpe-insert (poly exp u &optional (monom nil) &aux ctest)
(labels
((fi1(ul)
;; this routine recurses down the list of factors
(cond ((null ul)
;; Since this factor is nowhere to be seen in
;; this fpe, insert the new factor and its power
;; at the end of the list.
(list (cons poly exp)))
;;; added for stein@uni-paderborn.de 1/5/94
((coefonep (caar ul))(cons (car ul)(fi1 (cdr ul))))
;; fully expand. multiply into first factor that
;; better be of multiplicity 1. Presumably is the
;; only factor, but we can manage if it is not.
(*expand*
(if (= (cdar ul) exp) ;probably 1=exp
(cons (cons (p* poly (caar ul)) exp) (cdr ul))
(error "illegal use of expand")))
;;; end
((eq 'e (setq ctest (pcompare poly (caar ul))))
;; We found this factor. Increment the exponent.
(cons (cons (caar ul) (+ exp (cdar ul)))(cdr ul)))
;; Otherwise we keep looking for the factor.
((eq 'l ctest)
;; this poly is less in ordering than (caar ul),
;; so place it right here.
(cons (cons poly exp) ul))
;; This keep searching
(t (cons (car ul)(fi1 (cdr ul)))))))
(cond ((coefp poly)
;; inserting an integer factor: modify the first element
(cond ((coefonep poly) u)
((coefzerop poly) (make-fpe 0 1))
(t (cons (cons (coef* (caar u)
(coef^ poly exp))
1)
(cdr u)))))
((= exp 0) u) ;;multiplying u by z^0 gives u
((and monom (monomialp poly)(or (> (degree poly) 1)
(not (coefonep (lc poly)))))
;; the poly is something like (3*x^2) ^4
;; insert 81 (i.e. 3^4) into 1*(x)^8.
;; should work recursively for (3*y^2*x^2)^4
(fpe-insert (lc poly) exp (fpe-insert
(vector (svref poly 0) 0 1)
(* (degree poly) exp) u )
monom))
;; we can, more generally, not insist on monomial,
;; but only 0 x^1 term. That is,
;; (a*x^5+b*x^2) --> x^2*(a*x^3+b). or
;; (x^5+x^2)^3 --> x^6*(x^3+1)^3. etc.
;; note that 0 constant term is not useful because
;; x is encoded as (vector n 0 1) and has 0 const term
((and monom
(not *expand*)
(coefzerop(constc poly))
(> (length (the simple-vector poly)) 3)
(coefzerop (svref poly 2)))
(do ((i 2 (1+ i)))
((null (coefzerop (svref poly i)));there's a non-0
(fpe-insert
(vector (svref poly 0) 0 1)
(* exp (1- i)) ;; the degree of the factor
(fpe-insert (polyshift poly (1- i)) exp u)))
(declare (fixnum i))
;; no do-body
))
(t (fi1 u)))))