\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra sups.spad} \author{Clifton J. Williamson} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{domain ISUPS InnerSparseUnivariatePowerSeries} <>= )abbrev domain ISUPS InnerSparseUnivariatePowerSeries ++ Author: Clifton J. Williamson ++ Date Created: 28 October 1994 ++ Date Last Updated: 9 March 1995 ++ Basic Operations: ++ Related Domains: SparseUnivariateTaylorSeries, SparseUnivariateLaurentSeries ++ SparseUnivariatePuiseuxSeries ++ Also See: ++ AMS Classifications: ++ Keywords: sparse, series ++ Examples: ++ References: ++ Description: InnerSparseUnivariatePowerSeries is an internal domain ++ used for creating sparse Taylor and Laurent series. InnerSparseUnivariatePowerSeries(Coef): Exports == Implementation where Coef : Ring B ==> Boolean COM ==> OrderedCompletion Integer I ==> Integer L ==> List NNI ==> NonNegativeInteger OUT ==> OutputForm PI ==> PositiveInteger REF ==> Reference OrderedCompletion Integer RN ==> Fraction Integer Term ==> Record(k:Integer,c:Coef) SG ==> String ST ==> Stream Term Exports ==> UnivariatePowerSeriesCategory(Coef,Integer) with makeSeries: (REF,ST) -> % ++ makeSeries(refer,str) creates a power series from the reference ++ \spad{refer} and the stream \spad{str}. getRef: % -> REF ++ getRef(f) returns a reference containing the order to which the ++ terms of f have been computed. getStream: % -> ST ++ getStream(f) returns the stream of terms representing the series f. series: ST -> % ++ series(st) creates a series from a stream of non-zero terms, ++ where a term is an exponent-coefficient pair. The terms in the ++ stream should be ordered by increasing order of exponents. monomial?: % -> B ++ monomial?(f) tests if f is a single monomial. multiplyCoefficients: (I -> Coef,%) -> % ++ multiplyCoefficients(fn,f) returns the series ++ \spad{sum(fn(n) * an * x^n,n = n0..)}, ++ where f is the series \spad{sum(an * x^n,n = n0..)}. iExquo: (%,%,B) -> Union(%,"failed") ++ iExquo(f,g,taylor?) is the quotient of the power series f and g. ++ If \spad{taylor?} is \spad{true}, then we must have ++ \spad{order(f) >= order(g)}. taylorQuoByVar: % -> % ++ taylorQuoByVar(a0 + a1 x + a2 x**2 + ...) ++ returns \spad{a1 + a2 x + a3 x**2 + ...} iCompose: (%,%) -> % ++ iCompose(f,g) returns \spad{f(g(x))}. This is an internal function ++ which should only be called for Taylor series \spad{f(x)} and ++ \spad{g(x)} such that the constant coefficient of \spad{g(x)} is zero. seriesToOutputForm: (ST,REF,Symbol,Coef,RN) -> OutputForm ++ seriesToOutputForm(st,refer,var,cen,r) prints the series ++ \spad{f((var - cen)^r)}. if Coef has Algebra Fraction Integer then integrate: % -> % ++ integrate(f(x)) returns an anti-derivative of the power series ++ \spad{f(x)} with constant coefficient 0. ++ Warning: function does not check for a term of degree -1. cPower: (%,Coef) -> % ++ cPower(f,r) computes \spad{f^r}, where f has constant coefficient 1. ++ For use when the coefficient ring is commutative. cRationalPower: (%,RN) -> % ++ cRationalPower(f,r) computes \spad{f^r}. ++ For use when the coefficient ring is commutative. cExp: % -> % ++ cExp(f) computes the exponential of the power series f. ++ For use when the coefficient ring is commutative. cLog: % -> % ++ cLog(f) computes the logarithm of the power series f. ++ For use when the coefficient ring is commutative. cSin: % -> % ++ cSin(f) computes the sine of the power series f. ++ For use when the coefficient ring is commutative. cCos: % -> % ++ cCos(f) computes the cosine of the power series f. ++ For use when the coefficient ring is commutative. cTan: % -> % ++ cTan(f) computes the tangent of the power series f. ++ For use when the coefficient ring is commutative. cCot: % -> % ++ cCot(f) computes the cotangent of the power series f. ++ For use when the coefficient ring is commutative. cSec: % -> % ++ cSec(f) computes the secant of the power series f. ++ For use when the coefficient ring is commutative. cCsc: % -> % ++ cCsc(f) computes the cosecant of the power series f. ++ For use when the coefficient ring is commutative. cAsin: % -> % ++ cAsin(f) computes the arcsine of the power series f. ++ For use when the coefficient ring is commutative. cAcos: % -> % ++ cAcos(f) computes the arccosine of the power series f. ++ For use when the coefficient ring is commutative. cAtan: % -> % ++ cAtan(f) computes the arctangent of the power series f. ++ For use when the coefficient ring is commutative. cAcot: % -> % ++ cAcot(f) computes the arccotangent of the power series f. ++ For use when the coefficient ring is commutative. cAsec: % -> % ++ cAsec(f) computes the arcsecant of the power series f. ++ For use when the coefficient ring is commutative. cAcsc: % -> % ++ cAcsc(f) computes the arccosecant of the power series f. ++ For use when the coefficient ring is commutative. cSinh: % -> % ++ cSinh(f) computes the hyperbolic sine of the power series f. ++ For use when the coefficient ring is commutative. cCosh: % -> % ++ cCosh(f) computes the hyperbolic cosine of the power series f. ++ For use when the coefficient ring is commutative. cTanh: % -> % ++ cTanh(f) computes the hyperbolic tangent of the power series f. ++ For use when the coefficient ring is commutative. cCoth: % -> % ++ cCoth(f) computes the hyperbolic cotangent of the power series f. ++ For use when the coefficient ring is commutative. cSech: % -> % ++ cSech(f) computes the hyperbolic secant of the power series f. ++ For use when the coefficient ring is commutative. cCsch: % -> % ++ cCsch(f) computes the hyperbolic cosecant of the power series f. ++ For use when the coefficient ring is commutative. cAsinh: % -> % ++ cAsinh(f) computes the inverse hyperbolic sine of the power ++ series f. For use when the coefficient ring is commutative. cAcosh: % -> % ++ cAcosh(f) computes the inverse hyperbolic cosine of the power ++ series f. For use when the coefficient ring is commutative. cAtanh: % -> % ++ cAtanh(f) computes the inverse hyperbolic tangent of the power ++ series f. For use when the coefficient ring is commutative. cAcoth: % -> % ++ cAcoth(f) computes the inverse hyperbolic cotangent of the power ++ series f. For use when the coefficient ring is commutative. cAsech: % -> % ++ cAsech(f) computes the inverse hyperbolic secant of the power ++ series f. For use when the coefficient ring is commutative. cAcsch: % -> % ++ cAcsch(f) computes the inverse hyperbolic cosecant of the power ++ series f. For use when the coefficient ring is commutative. Implementation ==> add import REF Rep := Record(%ord: REF,%str: Stream Term) -- when the value of 'ord' is n, this indicates that all non-zero -- terms of order up to and including n have been computed; -- when 'ord' is plusInfinity, all terms have been computed; -- lazy evaluation of 'str' has the side-effect of modifying the value -- of 'ord' --% Local functions makeTerm: (Integer,Coef) -> Term getCoef: Term -> Coef getExpon: Term -> Integer iSeries: (ST,REF) -> ST iExtend: (ST,COM,REF) -> ST iTruncate0: (ST,REF,REF,COM,I,I) -> ST iTruncate: (%,COM,I) -> % iCoefficient: (ST,Integer) -> Coef iOrder: (ST,COM,REF) -> I iMap1: ((Coef,I) -> Coef,I -> I,B,ST,REF,REF,Integer) -> ST iMap2: ((Coef,I) -> Coef,I -> I,B,%) -> % iPlus1: ((Coef,Coef) -> Coef,ST,REF,ST,REF,REF,I) -> ST iPlus2: ((Coef,Coef) -> Coef,%,%) -> % productByTerm: (Coef,I,ST,REF,REF,I) -> ST productLazyEval: (ST,REF,ST,REF,COM) -> Void iTimes: (ST,REF,ST,REF,REF,I) -> ST iDivide: (ST,REF,ST,REF,Coef,I,REF,I) -> ST divide: (%,I,%,I,Coef) -> % compose0: (ST,REF,ST,REF,I,%,%,I,REF,I) -> ST factorials?: () -> Boolean termOutput: (RN,Coef,OUT) -> OUT showAll?: () -> Boolean --% macros makeTerm(exp,coef) == [exp,coef] getCoef term == term.c getExpon term == term.k makeSeries(refer,x) == [refer,x] getRef ups == ups.%ord getStream ups == ups.%str --% creation and destruction of series monomial(coef,expon) == nix : ST := empty() st := zero? coef => nix concat(makeTerm(expon,coef),nix) makeSeries(ref plusInfinity(),st) monomial? ups == (not empty? getStream ups) and (empty? rst getStream ups) coerce(n:I) == n :: Coef :: % coerce(r:Coef) == monomial(r,0) iSeries(x,refer) == empty? x => (setref(refer,plusInfinity()); empty()) setref(refer,(getExpon frst x) :: COM) concat(frst x,iSeries(rst x,refer)) series(x:ST) == empty? x => 0 n := getExpon frst x; refer := ref(n :: COM) makeSeries(refer,iSeries(x,refer)) --% values characteristic == characteristic$Coef 0 == monomial(0,0) 1 == monomial(1,0) iExtend(st,n,refer) == (deref refer) < n => explicitlyEmpty? st => (setref(refer,plusInfinity()); st) explicitEntries? st => iExtend(rst st,n,refer) iExtend(lazyEvaluate st,n,refer) st extend(x,n) == (iExtend(getStream x,n :: COM,getRef x); x) complete x == (iExtend(getStream x,plusInfinity(),getRef x); x) iTruncate0(x,xRefer,refer,minExp,maxExp,n) == delay explicitlyEmpty? x => (setref(refer,plusInfinity()); empty()) nn := n :: COM while (deref xRefer) < nn repeat lazyEvaluate x explicitEntries? x => (nx := getExpon(xTerm := frst x)) > maxExp => (setref(refer,plusInfinity()); empty()) setref(refer,nx :: COM) (nx :: COM) >= minExp => concat(makeTerm(nx,getCoef xTerm),_ iTruncate0(rst x,xRefer,refer,minExp,maxExp,nx + 1)) iTruncate0(rst x,xRefer,refer,minExp,maxExp,nx + 1) -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref xRefer)@I setref(refer,degr :: COM) iTruncate0(x,xRefer,refer,minExp,maxExp,degr + 1) iTruncate(ups,minExp,maxExp) == x := getStream ups; xRefer := getRef ups explicitlyEmpty? x => 0 explicitEntries? x => deg := getExpon frst x refer := ref((deg - 1) :: COM) makeSeries(refer,iTruncate0(x,xRefer,refer,minExp,maxExp,deg)) -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref xRefer)@I refer := ref(degr :: COM) makeSeries(refer,iTruncate0(x,xRefer,refer,minExp,maxExp,degr + 1)) truncate(ups,n) == iTruncate(ups,minusInfinity(),n) truncate(ups,n1,n2) == if n1 > n2 then (n1,n2) := (n2,n1) iTruncate(ups,n1 :: COM,n2) iCoefficient(st,n) == explicitEntries? st => term := frst st (expon := getExpon term) > n => 0 expon = n => getCoef term iCoefficient(rst st,n) 0 coefficient(x,n) == (extend(x,n); iCoefficient(getStream x,n)) elt(x:%,n:Integer) == coefficient(x,n) iOrder(st,n,refer) == explicitlyEmpty? st => error "order: series has infinite order" explicitEntries? st => ((r := getExpon frst st) :: COM) >= n => retract(n)@Integer r -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref refer)@I (degr :: COM) >= n => retract(n)@Integer iOrder(lazyEvaluate st,n,refer) order x == iOrder(getStream x,plusInfinity(),getRef x) order(x,n) == iOrder(getStream x,n :: COM,getRef x) terms x == getStream x --% predicates zero? ups == x := getStream ups; ref := getRef ups whatInfinity(n := deref ref) = 1 => explicitlyEmpty? x count : NNI := _$streamCount$Lisp for i in 1..count repeat explicitlyEmpty? x => return true explicitEntries? x => return false lazyEvaluate x false ups1 = ups2 == zero?(ups1 - ups2) --% arithmetic iMap1(cFcn,eFcn,check?,x,xRefer,refer,n) == delay -- when this function is called, all terms in 'x' of order < n have been -- computed and we compute the eFcn(n)th order coefficient of the result explicitlyEmpty? x => (setref(refer,plusInfinity()); empty()) -- if terms in 'x' up to order n have not been computed, -- apply lazy evaluation nn := n :: COM while (deref xRefer) < nn repeat lazyEvaluate x -- 'x' may now be empty: retest explicitlyEmpty? x => (setref(refer,plusInfinity()); empty()) -- must have nx >= n explicitEntries? x => xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm newCoef := cFcn(xCoef,nx); m := eFcn nx setref(refer,m :: COM) not check? => concat(makeTerm(m,newCoef),_ iMap1(cFcn,eFcn,check?,rst x,xRefer,refer,nx + 1)) zero? newCoef => iMap1(cFcn,eFcn,check?,rst x,xRefer,refer,nx + 1) concat(makeTerm(m,newCoef),_ iMap1(cFcn,eFcn,check?,rst x,xRefer,refer,nx + 1)) -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref xRefer)@I setref(refer,eFcn(degr) :: COM) iMap1(cFcn,eFcn,check?,x,xRefer,refer,degr + 1) iMap2(cFcn,eFcn,check?,ups) == -- 'eFcn' must be a strictly increasing function, -- i.e. i < j => eFcn(i) < eFcn(j) xRefer := getRef ups; x := getStream ups explicitlyEmpty? x => 0 explicitEntries? x => deg := getExpon frst x refer := ref(eFcn(deg - 1) :: COM) makeSeries(refer,iMap1(cFcn,eFcn,check?,x,xRefer,refer,deg)) -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref xRefer)@I refer := ref(eFcn(degr) :: COM) makeSeries(refer,iMap1(cFcn,eFcn,check?,x,xRefer,refer,degr + 1)) map(fcn,x) == iMap2(fcn(#1),#1,true,x) differentiate x == iMap2(#2 * #1,#1 - 1,true,x) multiplyCoefficients(f,x) == iMap2(f(#2) * #1,#1,true,x) multiplyExponents(x,n) == iMap2(#1,n * #1,false,x) iPlus1(op,x,xRefer,y,yRefer,refer,n) == delay -- when this function is called, all terms in 'x' and 'y' of order < n -- have been computed and we are computing the nth order coefficient of -- the result; note the 'op' is either '+' or '-' explicitlyEmpty? x => iMap1(op(0,#1),#1,false,y,yRefer,refer,n) explicitlyEmpty? y => iMap1(op(#1,0),#1,false,x,xRefer,refer,n) -- if terms up to order n have not been computed, -- apply lazy evaluation nn := n :: COM while (deref xRefer) < nn repeat lazyEvaluate x while (deref yRefer) < nn repeat lazyEvaluate y -- 'x' or 'y' may now be empty: retest explicitlyEmpty? x => iMap1(op(0,#1),#1,false,y,yRefer,refer,n) explicitlyEmpty? y => iMap1(op(#1,0),#1,false,x,xRefer,refer,n) -- must have nx >= n, ny >= n -- both x and y have explicit terms explicitEntries?(x) and explicitEntries?(y) => xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm yCoef := getCoef(yTerm := frst y); ny := getExpon yTerm nx = ny => setref(refer,nx :: COM) zero? (coef := op(xCoef,yCoef)) => iPlus1(op,rst x,xRefer,rst y,yRefer,refer,nx + 1) concat(makeTerm(nx,coef),_ iPlus1(op,rst x,xRefer,rst y,yRefer,refer,nx + 1)) nx < ny => setref(refer,nx :: COM) concat(makeTerm(nx,op(xCoef,0)),_ iPlus1(op,rst x,xRefer,y,yRefer,refer,nx + 1)) setref(refer,ny :: COM) concat(makeTerm(ny,op(0,yCoef)),_ iPlus1(op,x,xRefer,rst y,yRefer,refer,ny + 1)) -- y has no term of degree n explicitEntries? x => xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm -- can't have deref(yRefer) = infty unless all terms have been computed (degr := retract(deref yRefer)@I) < nx => setref(refer,deref yRefer) iPlus1(op,x,xRefer,y,yRefer,refer,degr + 1) setref(refer,nx :: COM) concat(makeTerm(nx,op(xCoef,0)),_ iPlus1(op,rst x,xRefer,y,yRefer,refer,nx + 1)) -- x has no term of degree n explicitEntries? y => yCoef := getCoef(yTerm := frst y); ny := getExpon yTerm -- can't have deref(xRefer) = infty unless all terms have been computed (degr := retract(deref xRefer)@I) < ny => setref(refer,deref xRefer) iPlus1(op,x,xRefer,y,yRefer,refer,degr + 1) setref(refer,ny :: COM) concat(makeTerm(ny,op(0,yCoef)),_ iPlus1(op,x,xRefer,rst y,yRefer,refer,ny + 1)) -- neither x nor y has a term of degree n setref(refer,xyRef := min(deref xRefer,deref yRefer)) -- can't have xyRef = infty unless all terms have been computed iPlus1(op,x,xRefer,y,yRefer,refer,retract(xyRef)@I + 1) iPlus2(op,ups1,ups2) == xRefer := getRef ups1; x := getStream ups1 xDeg := explicitlyEmpty? x => return map(op(0$Coef,#1),ups2) explicitEntries? x => (getExpon frst x) - 1 -- can't have deref(xRefer) = infty unless all terms have been computed retract(deref xRefer)@I yRefer := getRef ups2; y := getStream ups2 yDeg := explicitlyEmpty? y => return map(op(#1,0$Coef),ups1) explicitEntries? y => (getExpon frst y) - 1 -- can't have deref(yRefer) = infty unless all terms have been computed retract(deref yRefer)@I deg := min(xDeg,yDeg); refer := ref(deg :: COM) makeSeries(refer,iPlus1(op,x,xRefer,y,yRefer,refer,deg + 1)) x + y == iPlus2(#1 + #2,x,y) x - y == iPlus2(#1 - #2,x,y) - y == iMap2(_-#1,#1,false,y) -- gives correct defaults for I, NNI and PI n:I * x:% == (zero? n => 0; map(n * #1,x)) n:NNI * x:% == (zero? n => 0; map(n * #1,x)) n:PI * x:% == (zero? n => 0; map(n * #1,x)) productByTerm(coef,expon,x,xRefer,refer,n) == iMap1(coef * #1,#1 + expon,true,x,xRefer,refer,n) productLazyEval(x,xRefer,y,yRefer,nn) == explicitlyEmpty?(x) or explicitlyEmpty?(y) => void() explicitEntries? x => explicitEntries? y => void() xDeg := (getExpon frst x) :: COM while (xDeg + deref(yRefer)) < nn repeat lazyEvaluate y explicitEntries? y => yDeg := (getExpon frst y) :: COM while (yDeg + deref(xRefer)) < nn repeat lazyEvaluate x lazyEvaluate x -- if x = y, then y may now have explicit entries if lazy? y then lazyEvaluate y productLazyEval(x,xRefer,y,yRefer,nn) iTimes(x,xRefer,y,yRefer,refer,n) == delay -- when this function is called, we are computing the nth order -- coefficient of the product productLazyEval(x,xRefer,y,yRefer,n :: COM) explicitlyEmpty?(x) or explicitlyEmpty?(y) => (setref(refer,plusInfinity()); empty()) -- must have nx + ny >= n explicitEntries?(x) and explicitEntries?(y) => xCoef := getCoef(xTerm := frst x); xExpon := getExpon xTerm yCoef := getCoef(yTerm := frst y); yExpon := getExpon yTerm expon := xExpon + yExpon setref(refer,expon :: COM) scRefer := ref(expon :: COM) scMult := productByTerm(xCoef,xExpon,rst y,yRefer,scRefer,yExpon + 1) prRefer := ref(expon :: COM) pr := iTimes(rst x,xRefer,y,yRefer,prRefer,expon + 1) sm := iPlus1(#1 + #2,scMult,scRefer,pr,prRefer,refer,expon + 1) zero?(coef := xCoef * yCoef) => sm concat(makeTerm(expon,coef),sm) explicitEntries? x => xExpon := getExpon frst x -- can't have deref(yRefer) = infty unless all terms have been computed degr := retract(deref yRefer)@I setref(refer,(xExpon + degr) :: COM) iTimes(x,xRefer,y,yRefer,refer,xExpon + degr + 1) explicitEntries? y => yExpon := getExpon frst y -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref xRefer)@I setref(refer,(yExpon + degr) :: COM) iTimes(x,xRefer,y,yRefer,refer,yExpon + degr + 1) -- can't have deref(xRefer) = infty unless all terms have been computed xDegr := retract(deref xRefer)@I yDegr := retract(deref yRefer)@I setref(refer,(xDegr + yDegr) :: COM) iTimes(x,xRefer,y,yRefer,refer,xDegr + yDegr + 1) ups1:% * ups2:% == xRefer := getRef ups1; x := getStream ups1 xDeg := explicitlyEmpty? x => return 0 explicitEntries? x => (getExpon frst x) - 1 -- can't have deref(xRefer) = infty unless all terms have been computed retract(deref xRefer)@I yRefer := getRef ups2; y := getStream ups2 yDeg := explicitlyEmpty? y => return 0 explicitEntries? y => (getExpon frst y) - 1 -- can't have deref(yRefer) = infty unless all terms have been computed retract(deref yRefer)@I deg := xDeg + yDeg + 1; refer := ref(deg :: COM) makeSeries(refer,iTimes(x,xRefer,y,yRefer,refer,deg + 1)) iDivide(x,xRefer,y,yRefer,rym,m,refer,n) == delay -- when this function is called, we are computing the nth order -- coefficient of the result explicitlyEmpty? x => (setref(refer,plusInfinity()); empty()) -- if terms up to order n - m have not been computed, -- apply lazy evaluation nm := (n + m) :: COM while (deref xRefer) < nm repeat lazyEvaluate x -- 'x' may now be empty: retest explicitlyEmpty? x => (setref(refer,plusInfinity()); empty()) -- must have nx >= n + m explicitEntries? x => newCoef := getCoef(xTerm := frst x) * rym; nx := getExpon xTerm prodRefer := ref(nx :: COM) prod := productByTerm(-newCoef,nx - m,rst y,yRefer,prodRefer,1) sumRefer := ref(nx :: COM) sum := iPlus1(#1 + #2,rst x,xRefer,prod,prodRefer,sumRefer,nx + 1) setref(refer,(nx - m) :: COM); term := makeTerm(nx - m,newCoef) concat(term,iDivide(sum,sumRefer,y,yRefer,rym,m,refer,nx - m + 1)) -- can't have deref(xRefer) = infty unless all terms have been computed degr := retract(deref xRefer)@I setref(refer,(degr - m) :: COM) iDivide(x,xRefer,y,yRefer,rym,m,refer,degr - m + 1) divide(ups1,deg1,ups2,deg2,r) == xRefer := getRef ups1; x := getStream ups1 yRefer := getRef ups2; y := getStream ups2 refer := ref((deg1 - deg2) :: COM) makeSeries(refer,iDivide(x,xRefer,y,yRefer,r,deg2,refer,deg1 - deg2 + 1)) iExquo(ups1,ups2,taylor?) == xRefer := getRef ups1; x := getStream ups1 yRefer := getRef ups2; y := getStream ups2 n : I := 0 -- try to find first non-zero term in y -- give up after 1000 lazy evaluations while not explicitEntries? y repeat explicitlyEmpty? y => return "failed" lazyEvaluate y (n := n + 1) > 1000 => return "failed" yCoef := getCoef(yTerm := frst y); ny := getExpon yTerm (ry := recip yCoef) case "failed" => "failed" nn := ny :: COM if taylor? then while (deref(xRefer) < nn) repeat explicitlyEmpty? x => return 0 explicitEntries? x => return "failed" lazyEvaluate x -- check if ups2 is a monomial empty? rst y => iMap2(#1 * (ry :: Coef),#1 - ny,false,ups1) explicitlyEmpty? x => 0 nx := explicitEntries? x => ((deg := getExpon frst x) < ny) and taylor? => return "failed" deg - 1 -- can't have deref(xRefer) = infty unless all terms have been computed retract(deref xRefer)@I divide(ups1,nx,ups2,ny,ry :: Coef) taylorQuoByVar ups == iMap2(#1,#1 - 1,false,ups - monomial(coefficient(ups,0),0)) compose0(x,xRefer,y,yRefer,yOrd,y1,yn0,n0,refer,n) == delay -- when this function is called, we are computing the nth order -- coefficient of the composite explicitlyEmpty? x => (setref(refer,plusInfinity()); empty()) -- if terms in 'x' up to order n have not been computed, -- apply lazy evaluation nn := n :: COM; yyOrd := yOrd :: COM while (yyOrd * deref(xRefer)) < nn repeat lazyEvaluate x explicitEntries? x => xCoef := getCoef(xTerm := frst x); n1 := getExpon xTerm zero? n1 => setref(refer,n1 :: COM) concat(makeTerm(n1,xCoef),_ compose0(rst x,xRefer,y,yRefer,yOrd,y1,yn0,n0,refer,n1 + 1)) yn1 := yn0 * y1 ** ((n1 - n0) :: NNI) z := getStream yn1; zRefer := getRef yn1 degr := yOrd * n1; prodRefer := ref((degr - 1) :: COM) prod := iMap1(xCoef * #1,#1,true,z,zRefer,prodRefer,degr) coRefer := ref((degr + yOrd - 1) :: COM) co := compose0(rst x,xRefer,y,yRefer,yOrd,y1,yn1,n1,coRefer,degr + yOrd) setref(refer,(degr - 1) :: COM) iPlus1(#1 + #2,prod,prodRefer,co,coRefer,refer,degr) -- can't have deref(xRefer) = infty unless all terms have been computed degr := yOrd * (retract(deref xRefer)@I + 1) setref(refer,(degr - 1) :: COM) compose0(x,xRefer,y,yRefer,yOrd,y1,yn0,n0,refer,degr) iCompose(ups1,ups2) == x := getStream ups1; xRefer := getRef ups1 y := getStream ups2; yRefer := getRef ups2 -- try to compute the order of 'ups2' n : I := _$streamCount$Lisp for i in 1..n while not explicitEntries? y repeat explicitlyEmpty? y => return coefficient(ups1,0) :: % lazyEvaluate y explicitlyEmpty? y => coefficient(ups1,0) :: % yOrd : I := explicitEntries? y => getExpon frst y retract(deref yRefer)@I compRefer := ref((-1) :: COM) makeSeries(compRefer,_ compose0(x,xRefer,y,yRefer,yOrd,ups2,1,0,compRefer,0)) if Coef has Algebra Fraction Integer then integrate x == iMap2(1/(#2 + 1) * #1,#1 + 1,true,x) --% Fixed point computations Ys ==> Y$ParadoxicalCombinatorsForStreams(Term) integ0: (ST,REF,REF,I) -> ST integ0(x,intRef,ansRef,n) == delay nLess1 := (n - 1) :: COM while (deref intRef) < nLess1 repeat lazyEvaluate x explicitlyEmpty? x => (setref(ansRef,plusInfinity()); empty()) explicitEntries? x => xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm setref(ansRef,(n1 := (nx + 1)) :: COM) concat(makeTerm(n1,inv(n1 :: RN) * xCoef),_ integ0(rst x,intRef,ansRef,n1)) -- can't have deref(intRef) = infty unless all terms have been computed degr := retract(deref intRef)@I; setref(ansRef,(degr + 1) :: COM) integ0(x,intRef,ansRef,degr + 2) integ1: (ST,REF,REF) -> ST integ1(x,intRef,ansRef) == integ0(x,intRef,ansRef,1) lazyInteg: (Coef,() -> ST,REF,REF) -> ST lazyInteg(a,xf,intRef,ansRef) == ansStr : ST := integ1(delay xf,intRef,ansRef) concat(makeTerm(0,a),ansStr) cPower(f,r) == -- computes f^r. f should have constant coefficient 1. fp := differentiate f fInv := iExquo(1,f,false) :: %; y := r * fp * fInv yRef := getRef y; yStr := getStream y intRef := ref((-1) :: COM); ansRef := ref(0 :: COM) ansStr := Ys lazyInteg(1,iTimes(#1,ansRef,yStr,yRef,intRef,0),_ intRef,ansRef) makeSeries(ansRef,ansStr) iExp: (%,Coef) -> % iExp(f,cc) == -- computes exp(f). cc = exp coefficient(f,0) fp := differentiate f fpRef := getRef fp; fpStr := getStream fp intRef := ref((-1) :: COM); ansRef := ref(0 :: COM) ansStr := Ys lazyInteg(cc,iTimes(#1,ansRef,fpStr,fpRef,intRef,0),_ intRef,ansRef) makeSeries(ansRef,ansStr) sincos0: (Coef,Coef,L ST,REF,REF,ST,REF,ST,REF) -> L ST sincos0(sinc,cosc,list,sinRef,cosRef,fpStr,fpRef,fpStr2,fpRef2) == sinStr := first list; cosStr := second list prodRef1 := ref((-1) :: COM); prodRef2 := ref((-1) :: COM) prodStr1 := iTimes(cosStr,cosRef,fpStr,fpRef,prodRef1,0) prodStr2 := iTimes(sinStr,sinRef,fpStr2,fpRef2,prodRef2,0) [lazyInteg(sinc,prodStr1,prodRef1,sinRef),_ lazyInteg(cosc,prodStr2,prodRef2,cosRef)] iSincos: (%,Coef,Coef,I) -> Record(%sin: %, %cos: %) iSincos(f,sinc,cosc,sign) == fp := differentiate f fpRef := getRef fp; fpStr := getStream fp fp2 := (one? sign => fp; -fp) fpRef2 := getRef fp2; fpStr2 := getStream fp2 sinRef := ref(0 :: COM); cosRef := ref(0 :: COM) sincos := Ys(sincos0(sinc,cosc,#1,sinRef,cosRef,fpStr,fpRef,fpStr2,fpRef2),2) sinStr := (zero? sinc => rst first sincos; first sincos) cosStr := (zero? cosc => rst second sincos; second sincos) [makeSeries(sinRef,sinStr),makeSeries(cosRef,cosStr)] tan0: (Coef,ST,REF,ST,REF,I) -> ST tan0(cc,ansStr,ansRef,fpStr,fpRef,sign) == sqRef := ref((-1) :: COM) sqStr := iTimes(ansStr,ansRef,ansStr,ansRef,sqRef,0) one : % := 1; oneStr := getStream one; oneRef := getRef one yRef := ref((-1) :: COM) yStr : ST := one? sign => iPlus1(#1 + #2,oneStr,oneRef,sqStr,sqRef,yRef,0) iPlus1(#1 - #2,oneStr,oneRef,sqStr,sqRef,yRef,0) intRef := ref((-1) :: COM) lazyInteg(cc,iTimes(yStr,yRef,fpStr,fpRef,intRef,0),intRef,ansRef) iTan: (%,%,Coef,I) -> % iTan(f,fp,cc,sign) == -- computes the tangent (and related functions) of f. fpRef := getRef fp; fpStr := getStream fp ansRef := ref(0 :: COM) ansStr := Ys tan0(cc,#1,ansRef,fpStr,fpRef,sign) zero? cc => makeSeries(ansRef,rst ansStr) makeSeries(ansRef,ansStr) --% Error Reporting TRCONST : SG := "series expansion involves transcendental constants" NPOWERS : SG := "series expansion has terms of negative degree" FPOWERS : SG := "series expansion has terms of fractional degree" MAYFPOW : SG := "series expansion may have terms of fractional degree" LOGS : SG := "series expansion has logarithmic term" NPOWLOG : SG := "series expansion has terms of negative degree or logarithmic term" NOTINV : SG := "leading coefficient not invertible" --% Rational powers and transcendental functions orderOrFailed : % -> Union(I,"failed") orderOrFailed uts == -- returns the order of x or "failed" -- if -1 is returned, the series is identically zero x := getStream uts for n in 0..1000 repeat explicitlyEmpty? x => return -1 explicitEntries? x => return getExpon frst x lazyEvaluate x "failed" RATPOWERS : Boolean := Coef has "**": (Coef,RN) -> Coef TRANSFCN : Boolean := Coef has TranscendentalFunctionCategory cRationalPower(uts,r) == (ord0 := orderOrFailed uts) case "failed" => error "**: series with many leading zero coefficients" order := ord0 :: I (n := order exquo denom(r)) case "failed" => error "**: rational power does not exist" cc := coefficient(uts,order) (ccInv := recip cc) case "failed" => error concat("**: ",NOTINV) ccPow := one? cc => cc one? denom r => not negative?(num := numer r) => cc ** (num :: NNI) (ccInv :: Coef) ** ((-num) :: NNI) RATPOWERS => cc ** r error "** rational power of coefficient undefined" uts1 := (ccInv :: Coef) * uts uts2 := uts1 * monomial(1,-order) monomial(ccPow,(n :: I) * numer(r)) * cPower(uts2,r :: Coef) cExp uts == zero?(cc := coefficient(uts,0)) => iExp(uts,1) TRANSFCN => iExp(uts,exp cc) error concat("exp: ",TRCONST) cLog uts == zero?(cc := coefficient(uts,0)) => error "log: constant coefficient should not be 0" one? cc => integrate(differentiate(uts) * (iExquo(1,uts,true) :: %)) TRANSFCN => y := iExquo(1,uts,true) :: % (log(cc) :: %) + integrate(y * differentiate(uts)) error concat("log: ",TRCONST) sincos: % -> Record(%sin: %, %cos: %) sincos uts == zero?(cc := coefficient(uts,0)) => iSincos(uts,0,1,-1) TRANSFCN => iSincos(uts,sin cc,cos cc,-1) error concat("sincos: ",TRCONST) cSin uts == sincos(uts).%sin cCos uts == sincos(uts).%cos cTan uts == zero?(cc := coefficient(uts,0)) => iTan(uts,differentiate uts,0,1) TRANSFCN => iTan(uts,differentiate uts,tan cc,1) error concat("tan: ",TRCONST) cCot uts == zero? uts => error "cot: cot(0) is undefined" zero?(cc := coefficient(uts,0)) => error error concat("cot: ",NPOWERS) TRANSFCN => iTan(uts,-differentiate uts,cot cc,1) error concat("cot: ",TRCONST) cSec uts == zero?(cc := coefficient(uts,0)) => iExquo(1,cCos uts,true) :: % TRANSFCN => cosUts := cCos uts zero? coefficient(cosUts,0) => error concat("sec: ",NPOWERS) iExquo(1,cosUts,true) :: % error concat("sec: ",TRCONST) cCsc uts == zero? uts => error "csc: csc(0) is undefined" TRANSFCN => sinUts := cSin uts zero? coefficient(sinUts,0) => error concat("csc: ",NPOWERS) iExquo(1,sinUts,true) :: % error concat("csc: ",TRCONST) cAsin uts == zero?(cc := coefficient(uts,0)) => integrate(cRationalPower(1 - uts*uts,-1/2) * differentiate(uts)) TRANSFCN => x := 1 - uts * uts cc = 1 or cc = -1 => -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("asin: ",MAYFPOW) (order := ord :: I) = -1 => return asin(cc) :: % odd? order => error concat("asin: ",FPOWERS) c0 := asin(cc) :: % c0 + integrate(cRationalPower(x,-1/2) * differentiate(uts)) c0 := asin(cc) :: % c0 + integrate(cRationalPower(x,-1/2) * differentiate(uts)) error concat("asin: ",TRCONST) cAcos uts == zero? uts => TRANSFCN => acos(0)$Coef :: % error concat("acos: ",TRCONST) TRANSFCN => x := 1 - uts * uts cc := coefficient(uts,0) cc = 1 or cc = -1 => -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("acos: ",MAYFPOW) (order := ord :: I) = -1 => return acos(cc) :: % odd? order => error concat("acos: ",FPOWERS) c0 := acos(cc) :: % c0 + integrate(-cRationalPower(x,-1/2) * differentiate(uts)) c0 := acos(cc) :: % c0 + integrate(-cRationalPower(x,-1/2) * differentiate(uts)) error concat("acos: ",TRCONST) cAtan uts == zero?(cc := coefficient(uts,0)) => y := iExquo(1,(1 :: %) + uts*uts,true) :: % integrate(y * (differentiate uts)) TRANSFCN => (y := iExquo(1,(1 :: %) + uts*uts,true)) case "failed" => error concat("atan: ",LOGS) (atan(cc) :: %) + integrate((y :: %) * (differentiate uts)) error concat("atan: ",TRCONST) cAcot uts == TRANSFCN => (y := iExquo(1,(1 :: %) + uts*uts,true)) case "failed" => error concat("acot: ",LOGS) cc := coefficient(uts,0) (acot(cc) :: %) + integrate(-(y :: %) * (differentiate uts)) error concat("acot: ",TRCONST) cAsec uts == zero?(cc := coefficient(uts,0)) => error "asec: constant coefficient should not be 0" TRANSFCN => x := uts * uts - 1 y := cc = 1 or cc = -1 => -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("asec: ",MAYFPOW) (order := ord :: I) = -1 => return asec(cc) :: % odd? order => error concat("asec: ",FPOWERS) cRationalPower(x,-1/2) * differentiate(uts) cRationalPower(x,-1/2) * differentiate(uts) (z := iExquo(y,uts,true)) case "failed" => error concat("asec: ",NOTINV) (asec(cc) :: %) + integrate(z :: %) error concat("asec: ",TRCONST) cAcsc uts == zero?(cc := coefficient(uts,0)) => error "acsc: constant coefficient should not be 0" TRANSFCN => x := uts * uts - 1 y := cc = 1 or cc = -1 => -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("acsc: ",MAYFPOW) (order := ord :: I) = -1 => return acsc(cc) :: % odd? order => error concat("acsc: ",FPOWERS) -cRationalPower(x,-1/2) * differentiate(uts) -cRationalPower(x,-1/2) * differentiate(uts) (z := iExquo(y,uts,true)) case "failed" => error concat("asec: ",NOTINV) (acsc(cc) :: %) + integrate(z :: %) error concat("acsc: ",TRCONST) sinhcosh: % -> Record(%sinh: %, %cosh: %) sinhcosh uts == zero?(cc := coefficient(uts,0)) => tmp := iSincos(uts,0,1,1) [tmp.%sin,tmp.%cos] TRANSFCN => tmp := iSincos(uts,sinh cc,cosh cc,1) [tmp.%sin,tmp.%cos] error concat("sinhcosh: ",TRCONST) cSinh uts == sinhcosh(uts).%sinh cCosh uts == sinhcosh(uts).%cosh cTanh uts == zero?(cc := coefficient(uts,0)) => iTan(uts,differentiate uts,0,-1) TRANSFCN => iTan(uts,differentiate uts,tanh cc,-1) error concat("tanh: ",TRCONST) cCoth uts == tanhUts := cTanh uts zero? tanhUts => error "coth: coth(0) is undefined" zero? coefficient(tanhUts,0) => error concat("coth: ",NPOWERS) iExquo(1,tanhUts,true) :: % cSech uts == coshUts := cCosh uts zero? coefficient(coshUts,0) => error concat("sech: ",NPOWERS) iExquo(1,coshUts,true) :: % cCsch uts == sinhUts := cSinh uts zero? coefficient(sinhUts,0) => error concat("csch: ",NPOWERS) iExquo(1,sinhUts,true) :: % cAsinh uts == x := 1 + uts * uts zero?(cc := coefficient(uts,0)) => cLog(uts + cRationalPower(x,1/2)) TRANSFCN => (ord := orderOrFailed x) case "failed" => error concat("asinh: ",MAYFPOW) (order := ord :: I) = -1 => return asinh(cc) :: % odd? order => error concat("asinh: ",FPOWERS) -- the argument to 'log' must have a non-zero constant term cLog(uts + cRationalPower(x,1/2)) error concat("asinh: ",TRCONST) cAcosh uts == zero? uts => TRANSFCN => acosh(0)$Coef :: % error concat("acosh: ",TRCONST) TRANSFCN => cc := coefficient(uts,0); x := uts*uts - 1 cc = 1 or cc = -1 => -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("acosh: ",MAYFPOW) (order := ord :: I) = -1 => return acosh(cc) :: % odd? order => error concat("acosh: ",FPOWERS) -- the argument to 'log' must have a non-zero constant term cLog(uts + cRationalPower(x,1/2)) cLog(uts + cRationalPower(x,1/2)) error concat("acosh: ",TRCONST) cAtanh uts == half := inv(2 :: RN) :: Coef zero?(cc := coefficient(uts,0)) => half * (cLog(1 + uts) - cLog(1 - uts)) TRANSFCN => cc = 1 or cc = -1 => error concat("atanh: ",LOGS) half * (cLog(1 + uts) - cLog(1 - uts)) error concat("atanh: ",TRCONST) cAcoth uts == zero? uts => TRANSFCN => acoth(0)$Coef :: % error concat("acoth: ",TRCONST) TRANSFCN => cc := coefficient(uts,0); half := inv(2 :: RN) :: Coef cc = 1 or cc = -1 => error concat("acoth: ",LOGS) half * (cLog(uts + 1) - cLog(uts - 1)) error concat("acoth: ",TRCONST) cAsech uts == zero? uts => error "asech: asech(0) is undefined" TRANSFCN => zero?(cc := coefficient(uts,0)) => error concat("asech: ",NPOWLOG) x := 1 - uts * uts cc = 1 or cc = -1 => -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("asech: ",MAYFPOW) (order := ord :: I) = -1 => return asech(cc) :: % odd? order => error concat("asech: ",FPOWERS) (utsInv := iExquo(1,uts,true)) case "failed" => error concat("asech: ",NOTINV) cLog((1 + cRationalPower(x,1/2)) * (utsInv :: %)) (utsInv := iExquo(1,uts,true)) case "failed" => error concat("asech: ",NOTINV) cLog((1 + cRationalPower(x,1/2)) * (utsInv :: %)) error concat("asech: ",TRCONST) cAcsch uts == zero? uts => error "acsch: acsch(0) is undefined" TRANSFCN => zero?(cc := coefficient(uts,0)) => error concat("acsch: ",NPOWLOG) x := uts * uts + 1 -- compute order of 'x' (ord := orderOrFailed x) case "failed" => error concat("acsc: ",MAYFPOW) (order := ord :: I) = -1 => return acsch(cc) :: % odd? order => error concat("acsch: ",FPOWERS) (utsInv := iExquo(1,uts,true)) case "failed" => error concat("acsch: ",NOTINV) cLog((1 + cRationalPower(x,1/2)) * (utsInv :: %)) error concat("acsch: ",TRCONST) --% Output forms -- check a global Lisp variable factorials?() == false termOutput(k,c,vv) == -- creates a term c * vv ** k k = 0 => c :: OUT mon := (k = 1 => vv; vv ** (k :: OUT)) -- if factorials?() and k > 1 then -- c := factorial(k)$IntegerCombinatoricFunctions * c -- mon := mon / hconcat(k :: OUT,"!" :: OUT) c = 1 => mon c = -1 => -mon (c :: OUT) * mon -- check a global Lisp variable showAll?() == true seriesToOutputForm(st,refer,var,cen,r) == vv := zero? cen => var :: OUT paren(var :: OUT - cen :: OUT) l : L OUT := empty() while explicitEntries? st repeat term := frst st l := concat(termOutput(getExpon(term) * r,getCoef term,vv),l) st := rst st l := explicitlyEmpty? st => l (deg := retractIfCan(deref refer)@Union(I,"failed")) case I => concat(prefix("O" :: OUT,[vv ** ((((deg :: I) + 1) * r) :: OUT)]),l) l empty? l => (0$Coef) :: OUT reduce("+",reverse! l) @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --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. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}