\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra manip.spad} \author{Manuel Bronstein, Robert Sutor} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package FACTFUNC FactoredFunctions} <>= )abbrev package FACTFUNC FactoredFunctions ++ Author: Manuel Bronstein ++ Date Created: 2 Feb 1988 ++ Date Last Updated: 25 Jun 1990 ++ Description: computes various functions on factored arguments. -- not visible to the user FactoredFunctions(M:IntegralDomain): Exports == Implementation where N ==> NonNegativeInteger Exports ==> with nthRoot: (Factored M,N) -> Record(exponent:N,coef:M,radicand:List M) ++ nthRoot(f, n) returns \spad{(p, r, [r1,...,rm])} such that ++ the nth-root of f is equal to \spad{r * pth-root(r1 * ... * rm)}, ++ where r1,...,rm are distinct factors of f, ++ each of which has an exponent smaller than p in f. log : Factored M -> List Record(coef:N, logand:M) ++ log(f) returns \spad{[(a1,b1),...,(am,bm)]} such that ++ the logarithm of f is equal to \spad{a1*log(b1) + ... + am*log(bm)}. Implementation ==> add nthRoot(ff, n) == coeff:M := 1 radi:List(M) := (one? unit ff => empty(); [unit ff]) lf := factors ff d:N := empty? radi => gcd(concat(n, [t.exponent::N for t in lf]))::N 1 n := n quo d for term in lf repeat qr := divide(term.exponent::N quo d, n) coeff := coeff * term.factor ** qr.quotient not zero?(qr.remainder) => radi := concat!(radi, term.factor ** qr.remainder) [n, coeff, radi] log ff == ans := unit ff concat([1, unit ff], [[term.exponent::N, term.factor] for term in factors ff]) @ \section{package POLYROOT PolynomialRoots} <>= )abbrev package POLYROOT PolynomialRoots ++ Author: Manuel Bronstein ++ Date Created: 15 July 1988 ++ Date Last Updated: 10 November 1993 ++ Description: computes n-th roots of quotients of ++ multivariate polynomials -- not visible to the user PolynomialRoots(E, V, R, P, F):Exports == Implementation where E: OrderedAbelianMonoidSup V: OrderedSet R: IntegralDomain P: PolynomialCategory(R, E, V) F: Field with coerce: P -> % numer : $ -> P ++ numer(x) \undocumented denom : $ -> P ++ denom(x) \undocumented N ==> NonNegativeInteger Z ==> Integer Q ==> Fraction Z REC ==> Record(exponent:N, coef:F, radicand:F) Exports ==> with rroot: (R, N) -> REC ++ rroot(f, n) returns \spad{[m,c,r]} such ++ that \spad{f**(1/n) = c * r**(1/m)}. qroot : (Q, N) -> REC ++ qroot(f, n) returns \spad{[m,c,r]} such ++ that \spad{f**(1/n) = c * r**(1/m)}. if R has GcdDomain then froot: (F, N) -> REC ++ froot(f, n) returns \spad{[m,c,r]} such ++ that \spad{f**(1/n) = c * r**(1/m)}. nthr: (P, N) -> Record(exponent:N,coef:P,radicand:List P) ++ nthr(p,n) should be local but conditional Implementation ==> add import FactoredFunctions Z import FactoredFunctions P rsplit: List P -> Record(coef:R, poly:P) zroot : (Z, N) -> Record(exponent:N, coef:Z, radicand:Z) zroot(x, n) == zero? x or one? x => [1, x, 1] s := nthRoot(squareFree x, n) [s.exponent, s.coef, */s.radicand] if R has imaginary: () -> R then czroot: (Z, N) -> REC czroot(x, n) == rec := zroot(x, n) rec.exponent = 2 and negative? rec.radicand => [rec.exponent, rec.coef * imaginary()::P::F, (-rec.radicand)::F] [rec.exponent, rec.coef::F, rec.radicand::F] qroot(x, n) == sn := czroot(numer x, n) sd := czroot(denom x, n) m := lcm(sn.exponent, sd.exponent)::N [m, sn.coef / sd.coef, (sn.radicand ** (m quo sn.exponent)) / (sd.radicand ** (m quo sd.exponent))] else qroot(x, n) == sn := zroot(numer x, n) sd := zroot(denom x, n) m := lcm(sn.exponent, sd.exponent)::N [m, sn.coef::F / sd.coef::F, (sn.radicand ** (m quo sn.exponent))::F / (sd.radicand ** (m quo sd.exponent))::F] if R has RetractableTo Fraction Z then rroot(x, n) == (r := retractIfCan(x)@Union(Fraction Z,"failed")) case "failed" => [n, 1, x::P::F] qroot(r::Q, n) else if R has RetractableTo Z then rroot(x, n) == (r := retractIfCan(x)@Union(Z,"failed")) case "failed" => [n, 1, x::P::F] qroot(r::Z::Q, n) else rroot(x, n) == [n, 1, x::P::F] rsplit l == r := 1$R p := 1$P for q in l repeat if (u := retractIfCan(q)@Union(R, "failed")) case "failed" then p := p * q else r := r * u::R [r, p] if R has GcdDomain then if R has RetractableTo Z then nthr(x, n) == (r := retractIfCan(x)@Union(Z,"failed")) case "failed" => nthRoot(squareFree x, n) rec := zroot(r::Z, n) [rec.exponent, rec.coef::P, [rec.radicand::P]] else nthr(x, n) == nthRoot(squareFree x, n) froot(x, n) == zero? x or one? x => [1, x, 1] sn := nthr(numer x, n) sd := nthr(denom x, n) pn := rsplit(sn.radicand) pd := rsplit(sd.radicand) rn := rroot(pn.coef, sn.exponent) rd := rroot(pd.coef, sd.exponent) m := lcm([rn.exponent, rd.exponent, sn.exponent, sd.exponent])::N [m, (sn.coef::F / sd.coef::F) * (rn.coef / rd.coef), ((rn.radicand ** (m quo rn.exponent)) / (rd.radicand ** (m quo rd.exponent))) * (pn.poly ** (m quo sn.exponent))::F / (pd.poly ** (m quo sd.exponent))::F] @ \section{package ALGMANIP AlgebraicManipulations} <>= )abbrev package ALGMANIP AlgebraicManipulations ++ Author: Manuel Bronstein ++ Date Created: 28 Mar 1988 ++ Date Last Updated: 5 August 1993 ++ Description: ++ AlgebraicManipulations provides functions to simplify and expand ++ expressions involving algebraic operators. ++ Keywords: algebraic, manipulation. AlgebraicManipulations(R, F): Exports == Implementation where R : IntegralDomain F : Join(Field, ExpressionSpace) with numer : $ -> SparseMultivariatePolynomial(R, Kernel $) ++ numer(x) \undocumented denom : $ -> SparseMultivariatePolynomial(R, Kernel $) ++ denom(x) \undocumented coerce : SparseMultivariatePolynomial(R, Kernel $) -> $ ++ coerce(x) \undocumented N ==> NonNegativeInteger Z ==> Integer OP ==> BasicOperator SY ==> Symbol K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) RF ==> Fraction P REC ==> Record(ker:List K, exponent: List Z) Exports ==> with rootSplit: F -> F ++ rootSplit(f) transforms every radical of the form ++ \spad{(a/b)**(1/n)} appearing in f into \spad{a**(1/n) / b**(1/n)}. ++ This transformation is not in general valid for all ++ complex numbers \spad{a} and b. ratDenom : F -> F ++ ratDenom(f) rationalizes the denominators appearing in f ++ by moving all the algebraic quantities into the numerators. ratDenom : (F, F) -> F ++ ratDenom(f, a) removes \spad{a} from the denominators in f ++ if \spad{a} is an algebraic kernel. ratDenom : (F, List F) -> F ++ ratDenom(f, [a1,...,an]) removes the ai's which are ++ algebraic kernels from the denominators in f. ratDenom : (F, List K) -> F ++ ratDenom(f, [a1,...,an]) removes the ai's which are ++ algebraic from the denominators in f. ratPoly : F -> SparseUnivariatePolynomial F ++ ratPoly(f) returns a polynomial p such that p has no ++ algebraic coefficients, and \spad{p(f) = 0}. if R has Join(GcdDomain, RetractableTo Integer) and F has FunctionSpace(R) then rootPower : F -> F ++ rootPower(f) transforms every radical power of the form ++ \spad{(a**(1/n))**m} into a simpler form if \spad{m} and ++ \spad{n} have a common factor. rootProduct: F -> F ++ rootProduct(f) combines every product of the form ++ \spad{(a**(1/n))**m * (a**(1/s))**t} into a single power ++ of a root of \spad{a}, and transforms every radical power ++ of the form \spad{(a**(1/n))**m} into a simpler form. rootSimp : F -> F ++ rootSimp(f) transforms every radical of the form ++ \spad{(a * b**(q*n+r))**(1/n)} appearing in f into ++ \spad{b**q * (a * b**r)**(1/n)}. ++ This transformation is not in general valid for all ++ complex numbers b. rootKerSimp: (OP, F, N) -> F ++ rootKerSimp(op,f,n) should be local but conditional. Implementation ==> add macro ALGOP == '%alg macro NTHR == 'nthRoot import PolynomialCategoryQuotientFunctions(IndexedExponents K,K,R,P,F) innerRF : (F, List K) -> F rootExpand : K -> F algkernels : List K -> List K rootkernels: List K -> List K dummy := kernel(new()$SY)$K ratDenom x == innerRF(x, algkernels tower x) ratDenom(x:F, l:List K):F == innerRF(x, algkernels l) ratDenom(x:F, y:F) == ratDenom(x, [y]) ratDenom(x:F, l:List F) == ratDenom(x, [retract(y)@K for y in l]$List(K)) algkernels l == select!(has?(operator #1, ALGOP), l) rootkernels l == select!(is?(operator #1, NTHR::SY), l) ratPoly x == numer univariate(denom(ratDenom inv(dummy::P::F - x))::F, dummy) rootSplit x == lk := rootkernels tower x eval(x, lk, [rootExpand k for k in lk]) rootExpand k == x := first argument k n := second argument k op := operator k op(numer(x)::F, n) / op(denom(x)::F, n) -- all the kernels in ll must be algebraic innerRF(x, ll) == empty?(l := sort!(#1 > #2, kernels x)$List(K)) or empty? setIntersection(ll, tower x) => x lk := empty()$List(K) k : K while not member?(k := first l, ll) repeat lk := concat(k, lk) empty?(l := rest l) => return eval(x, lk, [map(innerRF(#1, ll), kk) for kk in lk]) q := univariate(eval(x, lk, [map(innerRF(#1, ll), kk) for kk in lk]), k, minPoly k) map(innerRF(#1, ll), q) (map(innerRF(#1, ll), k)) if R has Join(GcdDomain, RetractableTo Integer) and F has FunctionSpace(R) then import PolynomialRoots(IndexedExponents K, K, R, P, F) sroot : K -> F inroot : (OP, F, N) -> F radeval: (P, K) -> F breakup: List K -> List REC if R has RadicalCategory then rootKerSimp(op, x, n) == (r := retractIfCan(x)@Union(R, "failed")) case R => nthRoot(r::R, n)::F inroot(op, x, n) else rootKerSimp(op, x, n) == inroot(op, x, n) -- l is a list of nth-roots, returns a list of records of the form -- [a**(1/n1),a**(1/n2),...], [n1,n2,...]] -- such that the whole list covers l exactly breakup l == empty? l => empty() k := first l a := first(arg := argument(k := first l)) n := retract(second arg)@Z expo := empty()$List(Z) others := same := empty()$List(K) for kk in rest l repeat if (a = first(arg := argument kk)) then same := concat(kk, same) expo := concat(retract(second arg)@Z, expo) else others := concat(kk, others) ll := breakup others concat([concat(k, same), concat(n, expo)], ll) rootProduct x == import P,F,Z import List K import List Z for rec in breakup rootkernels tower x repeat k0 := first(l := rec.ker) nx := numer x; dx := denom x if empty? rest l then x := radeval(nx, k0) / radeval(dx, k0) else n := lcm(rec.exponent) k := kernel(operator k0, [first argument k0, n::F], height k0)$K lv := [monomial(1, k, (n quo m)::N) for m in rec.exponent]$List(P) x := radeval(eval(nx, l, lv), k) / radeval(eval(dx, l, lv), k) x rootPower x == for k in rootkernels tower x repeat x := radeval(numer x, k) / radeval(denom x, k) x -- replaces (a**(1/n))**m in p by a power of a simpler radical of a if -- n and m have a common factor radeval(p, k) == a := first(arg := argument k) n := (retract(second arg)@Integer)::NonNegativeInteger ans:F := 0 q := univariate(p, k) while positive?(d := degree q) repeat term := one?(g := gcd(d, n)) => monomial(1, k, d) monomial(1, kernel(operator k, [a,(n quo g)::F], height k), d quo g) ans := ans + leadingCoefficient(q)::F * term::F q := reductum q leadingCoefficient(q)::F + ans inroot(op, x, n) == one? x => x (x ~= -1) and (one?(num := numer x) or (num = -1)) => inv inroot(op, (num * denom x)::F, n) (u := isExpt(x, op)) case "failed" => kernel(op, [x, n::F]) pr := u::Record(var:K, exponent:Integer) q := pr.exponent /$Fraction(Z) (n * retract(second argument(pr.var))@Z) qr := divide(numer q, denom q) x := first argument(pr.var) x ** qr.quotient * rootKerSimp(op,x,denom(q)::N) ** qr.remainder sroot k == pr := froot(first(arg := argument k),(retract(second arg)@Z)::N) pr.coef * rootKerSimp(operator k, pr.radicand, pr.exponent) rootSimp x == lk := rootkernels tower x eval(x, lk, [sroot k for k in lk]) @ \section{package SIMPAN SimplifyAlgebraicNumberConvertPackage} <>= )abbrev package SIMPAN SimplifyAlgebraicNumberConvertPackage ++ Package to allow simplify to be called on AlgebraicNumbers ++ by converting to EXPR(INT) SimplifyAlgebraicNumberConvertPackage(): with simplify: AlgebraicNumber -> Expression(Integer) ++ simplify(an) applies simplifications to an == add simplify(a:AlgebraicNumber) == simplify(a::Expression(Integer))$TranscendentalManipulations(Integer, Expression Integer) @ \section{package TRMANIP TranscendentalManipulations} <>= )abbrev package TRMANIP TranscendentalManipulations ++ Transformations on transcendental objects ++ Author: Bob Sutor, Manuel Bronstein ++ Date Created: Way back ++ Date Last Updated: 22 January 1996, added simplifyLog MCD. ++ Description: ++ TranscendentalManipulations provides functions to simplify and ++ expand expressions involving transcendental operators. ++ Keywords: transcendental, manipulation. TranscendentalManipulations(R, F): Exports == Implementation where R : GcdDomain F : Join(FunctionSpace R, TranscendentalFunctionCategory) Z ==> Integer K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial P POWER ==> '%power POW ==> Record(val: F,exponent: Z) PRODUCT ==> Record(coef : Z, var : K) FPR ==> Fraction Polynomial R Exports ==> with expand : F -> F ++ expand(f) performs the following expansions on f:\begin{items} ++ \item 1. logs of products are expanded into sums of logs, ++ \item 2. trigonometric and hyperbolic trigonometric functions ++ of sums are expanded into sums of products of trigonometric ++ and hyperbolic trigonometric functions. ++ \item 3. formal powers of the form \spad{(a/b)**c} are expanded into ++ \spad{a**c * b**(-c)}. ++ \end{items} simplify : F -> F ++ simplify(f) performs the following simplifications on f:\begin{items} ++ \item 1. rewrites trigs and hyperbolic trigs in terms ++ of \spad{sin} ,\spad{cos}, \spad{sinh}, \spad{cosh}. ++ \item 2. rewrites \spad{sin**2} and \spad{sinh**2} in terms ++ of \spad{cos} and \spad{cosh}, ++ \item 3. rewrites \spad{exp(a)*exp(b)} as \spad{exp(a+b)}. ++ \item 4. rewrites \spad{(a**(1/n))**m * (a**(1/s))**t} as a single ++ power of a single radical of \spad{a}. ++ \end{items} htrigs : F -> F ++ htrigs(f) converts all the exponentials in f into ++ hyperbolic sines and cosines. simplifyExp: F -> F ++ simplifyExp(f) converts every product \spad{exp(a)*exp(b)} ++ appearing in f into \spad{exp(a+b)}. simplifyLog : F -> F ++ simplifyLog(f) converts every \spad{log(a) - log(b)} appearing in f ++ into \spad{log(a/b)}, every \spad{log(a) + log(b)} into \spad{log(a*b)} ++ and every \spad{n*log(a)} into \spad{log(a^n)}. expandPower: F -> F ++ expandPower(f) converts every power \spad{(a/b)**c} appearing ++ in f into \spad{a**c * b**(-c)}. expandLog : F -> F ++ expandLog(f) converts every \spad{log(a/b)} appearing in f into ++ \spad{log(a) - log(b)}, and every \spad{log(a*b)} into ++ \spad{log(a) + log(b)}.. cos2sec : F -> F ++ cos2sec(f) converts every \spad{cos(u)} appearing in f into ++ \spad{1/sec(u)}. cosh2sech : F -> F ++ cosh2sech(f) converts every \spad{cosh(u)} appearing in f into ++ \spad{1/sech(u)}. cot2trig : F -> F ++ cot2trig(f) converts every \spad{cot(u)} appearing in f into ++ \spad{cos(u)/sin(u)}. coth2trigh : F -> F ++ coth2trigh(f) converts every \spad{coth(u)} appearing in f into ++ \spad{cosh(u)/sinh(u)}. csc2sin : F -> F ++ csc2sin(f) converts every \spad{csc(u)} appearing in f into ++ \spad{1/sin(u)}. csch2sinh : F -> F ++ csch2sinh(f) converts every \spad{csch(u)} appearing in f into ++ \spad{1/sinh(u)}. sec2cos : F -> F ++ sec2cos(f) converts every \spad{sec(u)} appearing in f into ++ \spad{1/cos(u)}. sech2cosh : F -> F ++ sech2cosh(f) converts every \spad{sech(u)} appearing in f into ++ \spad{1/cosh(u)}. sin2csc : F -> F ++ sin2csc(f) converts every \spad{sin(u)} appearing in f into ++ \spad{1/csc(u)}. sinh2csch : F -> F ++ sinh2csch(f) converts every \spad{sinh(u)} appearing in f into ++ \spad{1/csch(u)}. tan2trig : F -> F ++ tan2trig(f) converts every \spad{tan(u)} appearing in f into ++ \spad{sin(u)/cos(u)}. tanh2trigh : F -> F ++ tanh2trigh(f) converts every \spad{tanh(u)} appearing in f into ++ \spad{sinh(u)/cosh(u)}. tan2cot : F -> F ++ tan2cot(f) converts every \spad{tan(u)} appearing in f into ++ \spad{1/cot(u)}. tanh2coth : F -> F ++ tanh2coth(f) converts every \spad{tanh(u)} appearing in f into ++ \spad{1/coth(u)}. cot2tan : F -> F ++ cot2tan(f) converts every \spad{cot(u)} appearing in f into ++ \spad{1/tan(u)}. coth2tanh : F -> F ++ coth2tanh(f) converts every \spad{coth(u)} appearing in f into ++ \spad{1/tanh(u)}. removeCosSq: F -> F ++ removeCosSq(f) converts every \spad{cos(u)**2} appearing in f into ++ \spad{1 - sin(x)**2}, and also reduces higher ++ powers of \spad{cos(u)} with that formula. removeSinSq: F -> F ++ removeSinSq(f) converts every \spad{sin(u)**2} appearing in f into ++ \spad{1 - cos(x)**2}, and also reduces higher powers of ++ \spad{sin(u)} with that formula. removeCoshSq:F -> F ++ removeCoshSq(f) converts every \spad{cosh(u)**2} appearing in f into ++ \spad{1 - sinh(x)**2}, and also reduces higher powers of ++ \spad{cosh(u)} with that formula. removeSinhSq:F -> F ++ removeSinhSq(f) converts every \spad{sinh(u)**2} appearing in f into ++ \spad{1 - cosh(x)**2}, and also reduces higher powers ++ of \spad{sinh(u)} with that formula. if R has PatternMatchable(R) and R has ConvertibleTo(Pattern(R)) and F has ConvertibleTo(Pattern(R)) and F has PatternMatchable R then expandTrigProducts : F -> F ++ expandTrigProducts(e) replaces \axiom{sin(x)*sin(y)} by ++ \spad{(cos(x-y)-cos(x+y))/2}, \axiom{cos(x)*cos(y)} by ++ \spad{(cos(x-y)+cos(x+y))/2}, and \axiom{sin(x)*cos(y)} by ++ \spad{(sin(x-y)+sin(x+y))/2}. Note that this operation uses ++ the pattern matcher and so is relatively expensive. To avoid ++ getting into an infinite loop the transformations are applied ++ at most ten times. Implementation ==> add import FactoredFunctions(P) import PolynomialCategoryLifting(IndexedExponents K, K, R, P, F) import PolynomialCategoryQuotientFunctions(IndexedExponents K,K,R,P,F) smpexp : P -> F termexp : P -> F exlog : P -> F smplog : P -> F smpexpand : P -> F smp2htrigs: P -> F kerexpand : K -> F expandpow : K -> F logexpand : K -> F sup2htrigs: (UP, F) -> F supexp : (UP, F, F, Z) -> F ueval : (F, String, F -> F) -> F ueval2 : (F, String, F -> F) -> F powersimp : (P, List K) -> F t2t : F -> F c2t : F -> F c2s : F -> F s2c : F -> F s2c2 : F -> F th2th : F -> F ch2th : F -> F ch2sh : F -> F sh2ch : F -> F sh2ch2 : F -> F simplify0 : F -> F simplifyLog1 : F -> F logArgs : List F -> F import F import List F if R has PatternMatchable R and R has ConvertibleTo Pattern R and F has ConvertibleTo(Pattern(R)) and F has PatternMatchable R then XX : F := coerce new()$Symbol YY : F := coerce new()$Symbol sinCosRule : RewriteRule(R,R,F) := rule(cos(XX)*sin(YY),(sin(XX+YY)-sin(XX-YY))/2::F) sinSinRule : RewriteRule(R,R,F) := rule(sin(XX)*sin(YY),(cos(XX-YY)-cos(XX+YY))/2::F) cosCosRule : RewriteRule(R,R,F) := rule(cos(XX)*cos(YY),(cos(XX-YY)+cos(XX+YY))/2::F) expandTrigProducts(e:F):F == applyRules([sinCosRule,sinSinRule,cosCosRule],e,10)$ApplyRules(R,R,F) logArgs(l:List F):F == -- This function will take a list of Expressions (implicitly a sum) and -- add them up, combining log terms. It also replaces n*log(x) by -- log(x^n). import K sum : F := 0 arg : F := 1 for term in l repeat is?(term,'log) => arg := arg * simplifyLog(first(argument(first(kernels(term))))) -- Now look for multiples, including negative ones. prod : Union(PRODUCT, "failed") := isMult(term) (prod case PRODUCT) and is?(prod.var,'log) => arg := arg * simplifyLog ((first argument(prod.var))**(prod.coef)) sum := sum+term sum+log(arg) simplifyLog(e:F):F == simplifyLog1(numerator e)/simplifyLog1(denominator e) simplifyLog1(e:F):F == freeOf?(e,'log) => e -- Check for n*log(u) prod : Union(PRODUCT, "failed") := isMult(e) (prod case PRODUCT) and is?(prod.var,'log) => log simplifyLog ((first argument(prod.var))**(prod.coef)) termList : Union(List(F),"failed") := isTimes(e) -- I'm using two variables, termList and terms, to work round a -- bug in the old compiler. not (termList case "failed") => -- We want to simplify each log term in the product and then multiply -- them together. However, if there is a constant or arithmetic -- expression (i.e. somwthing which looks like a Polynomial) we would -- like to combine it with a log term. terms :List F := [simplifyLog(term) for term in termList::List(F)] exprs :List F := [] nterms :List F := [] for term in terms repeat if retractIfCan(term)@Union(FPR,"failed") case FPR then exprs := cons(term, exprs) else nterms := cons(term, nterms) terms := nterms if not empty? exprs then foundLog := false i : NonNegativeInteger := 0 while (not(foundLog) and (i < #terms)) repeat i := i+1 if is?(terms.i,'log) then args : List F := argument(retract(terms.i)@K) setelt(terms,i, log simplifyLog1(first(args)**(*/exprs))) foundLog := true -- The next line deals with a situation which shouldn't occur, -- since we have checked whether we are freeOf log already. if not foundLog then terms := append(exprs,terms) */terms terms : Union(List(F),"failed") := isPlus(e) not (terms case "failed") => logArgs(terms) expt : Union(POW, "failed") := isPower(e) (expt case POW) and not one? expt.exponent => simplifyLog(expt.val)**(expt.exponent) kers : List K := kernels e not(one?(#kers)) => e -- Have a constant kernel(operator first kers,[simplifyLog(u) for u in argument first kers]) if R has RetractableTo Integer then simplify x == rootProduct(simplify0 x)$AlgebraicManipulations(R,F) else simplify x == simplify0 x expandpow k == a := expandPower first(arg := argument k) b := expandPower second arg ne:F := (one? numer a => 1; numer(a)::F ** b) de:F := (one? denom a => 1; denom(a)::F ** (-b)) ne * de termexp p == exponent:F := 0 coef := (leadingCoefficient p)::P lpow := select(is?(#1, POWER)$K, lk := variables p)$List(K) for k in lk repeat d := degree(p, k) if is?(k, 'exp) then exponent := exponent + d * first argument k else if not is?(k, POWER) then -- Expand arguments to functions as well ... MCD 23/1/97 --coef := coef * monomial(1, k, d) coef := coef * monomial(1, kernel(operator k,[simplifyExp u for u in argument k], height k), d) coef::F * exp exponent * powersimp(p, lpow) expandPower f == l := select(is?(#1, POWER)$K, kernels f)$List(K) eval(f, l, [expandpow k for k in l]) -- l is a list of pure powers appearing as kernels in p powersimp(p, l) == empty? l => 1 k := first l -- k = a**b a := first(arg := argument k) exponent := degree(p, k) * second arg empty?(lk := select(a = first argument #1, rest l)) => (a ** exponent) * powersimp(p, rest l) for k0 in lk repeat exponent := exponent + degree(p, k0) * second argument k0 (a ** exponent) * powersimp(p, setDifference(rest l, lk)) t2t x == sin(x) / cos(x) c2t x == cos(x) / sin(x) c2s x == inv sin x s2c x == inv cos x s2c2 x == 1 - cos(x)**2 th2th x == sinh(x) / cosh(x) ch2th x == cosh(x) / sinh(x) ch2sh x == inv sinh x sh2ch x == inv cosh x sh2ch2 x == cosh(x)**2 - 1 ueval(x, s,f) == eval(x, s::Symbol, f) ueval2(x,s,f) == eval(x, s::Symbol, 2, f) cos2sec x == ueval(x, "cos", inv sec #1) sin2csc x == ueval(x, "sin", inv csc #1) csc2sin x == ueval(x, "csc", c2s) sec2cos x == ueval(x, "sec", s2c) tan2cot x == ueval(x, "tan", inv cot #1) cot2tan x == ueval(x, "cot", inv tan #1) tan2trig x == ueval(x, "tan", t2t) cot2trig x == ueval(x, "cot", c2t) cosh2sech x == ueval(x, "cosh", inv sech #1) sinh2csch x == ueval(x, "sinh", inv csch #1) csch2sinh x == ueval(x, "csch", ch2sh) sech2cosh x == ueval(x, "sech", sh2ch) tanh2coth x == ueval(x, "tanh", inv coth #1) coth2tanh x == ueval(x, "coth", inv tanh #1) tanh2trigh x == ueval(x, "tanh", th2th) coth2trigh x == ueval(x, "coth", ch2th) removeCosSq x == ueval2(x, "cos", 1 - (sin #1)**2) removeSinSq x == ueval2(x, "sin", s2c2) removeCoshSq x== ueval2(x, "cosh", 1 + (sinh #1)**2) removeSinhSq x== ueval2(x, "sinh", sh2ch2) expandLog x == smplog(numer x) / smplog(denom x) simplifyExp x == (smpexp numer x) / (smpexp denom x) expand x == (smpexpand numer x) / (smpexpand denom x) smpexpand p == map(kerexpand, #1::F, p) smplog p == map(logexpand, #1::F, p) smp2htrigs p == map(htrigs(#1::F), #1::F, p) htrigs f == (m := mainKernel f) case "failed" => f op := operator(k := m::K) arg := [htrigs x for x in argument k]$List(F) num := univariate(numer f, k) den := univariate(denom f, k) is?(op,'exp) => g1 := cosh(a := first arg) + sinh(a) g2 := cosh(a) - sinh(a) supexp(num,g1,g2,b:= (degree num)::Z quo 2)/supexp(den,g1,g2,b) sup2htrigs(num, g1:= op arg) / sup2htrigs(den, g1) supexp(p, f1, f2, bse) == ans:F := 0 while p ~= 0 repeat g := htrigs(leadingCoefficient(p)::F) if ((d := degree(p)::Z - bse) >= 0) then ans := ans + g * f1 ** d else ans := ans + g * f2 ** (-d) p := reductum p ans sup2htrigs(p, f) == (map(smp2htrigs, p)$SparseUnivariatePolynomialFunctions2(P, F)) f exlog p == +/[r.coef * log(r.logand::F) for r in log squareFree p] logexpand k == nullary?(op := operator k) => k::F is?(op,'log) => exlog(numer(x := expandLog first argument k)) - exlog denom x op [expandLog x for x in argument k]$List(F) kerexpand k == nullary?(op := operator k) => k::F is?(op, POWER) => expandpow k arg := first argument k is?(op,'sec) => inv expand cos arg is?(op,'csc) => inv expand sin arg is?(op,'log) => exlog(numer(x := expand arg)) - exlog denom x num := numer arg den := denom arg (b := (reductum num) / den) ~= 0 => a := (leadingMonomial num) / den is?(op,'exp) => exp(expand a) * expand(exp b) is?(op,'sin) => sin(expand a) * expand(cos b) + cos(expand a) * expand(sin b) is?(op,'cos) => cos(expand a) * expand(cos b) - sin(expand a) * expand(sin b) is?(op,'tan) => ta := tan expand a tb := expand tan b (ta + tb) / (1 - ta * tb) is?(op,'cot) => cta := cot expand a ctb := expand cot b (cta * ctb - 1) / (ctb + cta) op [expand x for x in argument k]$List(F) op [expand x for x in argument k]$List(F) smpexp p == ans:F := 0 while p ~= 0 repeat ans := ans + termexp leadingMonomial p p := reductum p ans -- this now works in 3 passes over the expression: -- pass1 rewrites trigs and htrigs in terms of sin,cos,sinh,cosh -- pass2 rewrites sin**2 and sinh**2 in terms of cos and cosh. -- pass3 groups exponentials together simplify0 x == simplifyExp eval(eval(x, ['tan,'cot,'sec,'csc, 'tanh,'coth,'sech,'csch], [t2t,c2t,s2c,c2s,th2th,ch2th,sh2ch,ch2sh]), ['sin, 'sinh], [2, 2], [s2c2, sh2ch2]) @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. --Copyright (C) 2007-2009, Gabriel Dos Reis. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> -- SPAD files for the functional world should be compiled in the -- following order: -- -- op kl function funcpkgs MANIP combfunc -- algfunc elemntry constant funceval fe <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}