\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra puiseux.spad} \author{Clifton J. Williamson, Scott C. Morrison} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{category UPXSCCA UnivariatePuiseuxSeriesConstructorCategory} <>= )abbrev category UPXSCCA UnivariatePuiseuxSeriesConstructorCategory ++ Author: Clifton J. Williamson ++ Date Created: 6 February 1990 ++ Date Last Updated: 22 March 1990 ++ Basic Operations: ++ Related Domains: ++ Also See: ++ AMS Classifications: ++ Keywords: series, Puiseux, Laurent ++ Examples: ++ References: ++ Description: ++ This is a category of univariate Puiseux series constructed ++ from univariate Laurent series. A Puiseux series is represented ++ by a pair \spad{[r,f(x)]}, where r is a positive rational number and ++ \spad{f(x)} is a Laurent series. This pair represents the Puiseux ++ series \spad{f(x^r)}. UnivariatePuiseuxSeriesConstructorCategory(Coef,ULS):_ Category == Definition where Coef : Ring ULS : UnivariateLaurentSeriesCategory Coef I ==> Integer RN ==> Fraction Integer Definition ==> Join(UnivariatePuiseuxSeriesCategory(Coef),_ RetractableTo ULS,CoercibleFrom ULS) with puiseux: (RN,ULS) -> % ++ \spad{puiseux(r,f(x))} returns \spad{f(x^r)}. rationalPower: % -> RN ++ \spad{rationalPower(f(x))} returns r where the Puiseux series ++ \spad{f(x) = g(x^r)}. laurentRep : % -> ULS ++ \spad{laurentRep(f(x))} returns \spad{g(x)} where the Puiseux series ++ \spad{f(x) = g(x^r)} is represented by \spad{[r,g(x)]}. degree: % -> RN ++ \spad{degree(f(x))} returns the degree of the leading term of the ++ Puiseux series \spad{f(x)}, which may have zero as a coefficient. laurent: % -> ULS ++ \spad{laurent(f(x))} converts the Puiseux series \spad{f(x)} to a ++ Laurent series if possible. Error: if this is not possible. laurentIfCan: % -> Union(ULS,"failed") ++ \spad{laurentIfCan(f(x))} converts the Puiseux series \spad{f(x)} ++ to a Laurent series if possible. ++ If this is not possible, "failed" is returned. add zero? x == zero? laurentRep x retract(x:%):ULS == laurent x retractIfCan(x:%):Union(ULS,"failed") == laurentIfCan x @ \section{domain UPXSCONS UnivariatePuiseuxSeriesConstructor} <>= )abbrev domain UPXSCONS UnivariatePuiseuxSeriesConstructor ++ Author: Clifton J. Williamson ++ Date Created: 9 May 1989 ++ Date Last Updated: 30 November 1994 ++ Basic Operations: ++ Related Domains: ++ Also See: ++ AMS Classifications: ++ Keywords: series, Puiseux, Laurent ++ Examples: ++ References: ++ Description: ++ This package enables one to construct a univariate Puiseux series ++ domain from a univariate Laurent series domain. Univariate ++ Puiseux series are represented by a pair \spad{[r,f(x)]}, where r is ++ a positive rational number and \spad{f(x)} is a Laurent series. ++ This pair represents the Puiseux series \spad{f(x^r)}. UnivariatePuiseuxSeriesConstructor(Coef,ULS):_ Exports == Implementation where Coef : Ring ULS : UnivariateLaurentSeriesCategory Coef I ==> Integer L ==> List NNI ==> NonNegativeInteger OUT ==> OutputForm PI ==> PositiveInteger RN ==> Fraction Integer ST ==> Stream Coef LTerm ==> Record(k:I,c:Coef) PTerm ==> Record(k:RN,c:Coef) ST2LP ==> StreamFunctions2(LTerm,PTerm) ST2PL ==> StreamFunctions2(PTerm,LTerm) Exports ==> UnivariatePuiseuxSeriesConstructorCategory(Coef,ULS) Implementation ==> add --% representation Rep := Record(expon:RN,lSeries:ULS) getExpon: % -> RN getULS : % -> ULS getExpon pxs == pxs.expon getULS pxs == pxs.lSeries --% creation and destruction puiseux(n,ls) == [n,ls] laurentRep x == getULS x rationalPower x == getExpon x degree x == getExpon(x) * degree(getULS(x)) 0 == puiseux(1,0) 1 == puiseux(1,1) monomial(c,k) == k = 0 => c :: % negative? k => puiseux(-k,monomial(c,-1)) puiseux(k,monomial(c,1)) coerce(ls:ULS) == puiseux(1,ls) coerce(r:Coef) == r :: ULS :: % coerce(i:I) == i :: Coef :: % laurentIfCan upxs == r := getExpon upxs one? denom r => multiplyExponents(getULS upxs,numer(r) :: PI) "failed" laurent upxs == (uls := laurentIfCan upxs) case "failed" => error "laurent: Puiseux series has fractional powers" uls :: ULS multExp: (RN,LTerm) -> PTerm multExp(r,lTerm) == [r * lTerm.k,lTerm.c] terms upxs == map(multExp(getExpon upxs,#1),terms getULS upxs)$ST2LP clearDen: (I,PTerm) -> LTerm clearDen(n,lTerm) == (int := retractIfCan(n * lTerm.k)@Union(I,"failed")) case "failed" => error "series: inappropriate denominator" [int :: I,lTerm.c] series(n,stream) == str := map(clearDen(n,#1),stream)$ST2PL puiseux(1/n,series str) --% normalizations rewrite:(%,PI) -> % rewrite(upxs,m) == -- rewrites a series in x**r as a series in x**(r/m) puiseux((getExpon upxs)*(1/m),multiplyExponents(getULS upxs,m)) ratGcd: (RN,RN) -> RN ratGcd(r1,r2) == -- if r1 = prod(p prime,p ** ep(r1)) and -- if r2 = prod(p prime,p ** ep(r2)), then -- ratGcd(r1,r2) = prod(p prime,p ** min(ep(r1),ep(r2))) gcd(numer r1,numer r2) / lcm(denom r1,denom r2) withNewExpon:(%,RN) -> % withNewExpon(upxs,r) == rewrite(upxs,numer(getExpon(upxs)/r) pretend PI) --% predicates upxs1 = upxs2 == r1 := getExpon upxs1; r2 := getExpon upxs2 ls1 := getULS upxs1; ls2 := getULS upxs2 (r1 = r2) => (ls1 = ls2) r := ratGcd(r1,r2) m1 := numer(getExpon(upxs1)/r) pretend PI m2 := numer(getExpon(upxs2)/r) pretend PI multiplyExponents(ls1,m1) = multiplyExponents(ls2,m2) pole? upxs == pole? getULS upxs --% arithmetic applyFcn:((ULS,ULS) -> ULS,%,%) -> % applyFcn(op,pxs1,pxs2) == r1 := getExpon pxs1; r2 := getExpon pxs2 ls1 := getULS pxs1; ls2 := getULS pxs2 (r1 = r2) => puiseux(r1,op(ls1,ls2)) r := ratGcd(r1,r2) m1 := numer(getExpon(pxs1)/r) pretend PI m2 := numer(getExpon(pxs2)/r) pretend PI puiseux(r,op(multiplyExponents(ls1,m1),multiplyExponents(ls2,m2))) pxs1 + pxs2 == applyFcn(#1 +$ULS #2,pxs1,pxs2) pxs1 - pxs2 == applyFcn(#1 -$ULS #2,pxs1,pxs2) pxs1:% * pxs2:% == applyFcn(#1 *$ULS #2,pxs1,pxs2) pxs:% ** n:NNI == puiseux(getExpon pxs,getULS(pxs)**n) recip pxs == rec := recip getULS pxs rec case "failed" => "failed" puiseux(getExpon pxs,rec :: ULS) RATALG : Boolean := Coef has Algebra(Fraction Integer) elt(upxs1:%,upxs2:%) == uls1 := laurentRep upxs1; uls2 := laurentRep upxs2 r1 := rationalPower upxs1; r2 := rationalPower upxs2 (n := retractIfCan(r1)@Union(Integer,"failed")) case Integer => puiseux(r2,uls1(uls2 ** r1)) RATALG => if zero? (coef := coefficient(uls2,deg := degree uls2)) then deg := order(uls2,deg + 1000) zero? (coef := coefficient(uls2,deg)) => error "elt: series with many leading zero coefficients" -- a fractional power of a Laurent series may not be defined: -- if f(x) = c * x**n + ..., then f(x) ** (p/q) will be defined -- only if q divides n b := lcm(denom r1,deg); c := b quo deg mon : ULS := monomial(1,c) uls2 := elt(uls2,mon) ** r1 puiseux(r2*(1/c),elt(uls1,uls2)) error "elt: rational powers not available for this coefficient domain" if Coef has "**": (Coef,Integer) -> Coef and Coef has "**": (Coef, RN) -> Coef then eval(upxs:%,a:Coef) == eval(getULS upxs,a ** getExpon(upxs)) if Coef has Field then pxs1:% / pxs2:% == applyFcn(#1 /$ULS #2,pxs1,pxs2) inv upxs == (invUpxs := recip upxs) case "failed" => error "inv: multiplicative inverse does not exist" invUpxs :: % --% values variable upxs == variable getULS upxs center upxs == center getULS upxs coefficient(upxs,rn) == one? denom(n := rn / getExpon upxs) => coefficient(getULS upxs,numer n) 0 elt(upxs:%,rn:RN) == coefficient(upxs,rn) --% other functions roundDown: RN -> I roundDown rn == -- returns the largest integer <= rn (den := denom rn) = 1 => numer rn n := (num := numer rn) quo den positive?(num) => n n - 1 roundUp: RN -> I roundUp rn == -- returns the smallest integer >= rn (den := denom rn) = 1 => numer rn n := (num := numer rn) quo den positive?(num) => n + 1 n order upxs == getExpon upxs * order getULS upxs order(upxs,r) == e := getExpon upxs ord := order(getULS upxs, n := roundDown(r / e)) ord = n => r ord * e truncate(upxs,r) == e := getExpon upxs puiseux(e,truncate(getULS upxs,roundDown(r / e))) truncate(upxs,r1,r2) == e := getExpon upxs puiseux(e,truncate(getULS upxs,roundUp(r1 / e),roundDown(r2 / e))) complete upxs == puiseux(getExpon upxs,complete getULS upxs) extend(upxs,r) == e := getExpon upxs puiseux(e,extend(getULS upxs,roundDown(r / e))) map(fcn,upxs) == puiseux(getExpon upxs,map(fcn,getULS upxs)) characteristic == characteristic$Coef -- multiplyCoefficients(f,upxs) == -- r := getExpon upxs -- puiseux(r,multiplyCoefficients(f(#1 * r),getULS upxs)) multiplyExponents(upxs:%,n:RN) == puiseux(n * getExpon(upxs),getULS upxs) multiplyExponents(upxs:%,n:PI) == puiseux(n * getExpon(upxs),getULS upxs) if Coef has "*": (Fraction Integer, Coef) -> Coef then differentiate upxs == r := getExpon upxs puiseux(r,differentiate getULS upxs) * monomial(r :: Coef,r-1) if Coef has PartialDifferentialRing(Symbol) then differentiate(upxs:%,s:Symbol) == (s = variable(upxs)) => differentiate upxs dcds := differentiate(center upxs,s) map(differentiate(#1,s),upxs) - dcds*differentiate(upxs) if Coef has Algebra Fraction Integer then coerce(r:RN) == r :: Coef :: % ratInv: RN -> Coef ratInv r == zero? r => 1 inv(r) :: Coef integrate upxs == not zero? coefficient(upxs,-1) => error "integrate: series has term of order -1" r := getExpon upxs uls := getULS upxs uls := multiplyCoefficients(ratInv(#1 * r + 1),uls) monomial(1,1) * puiseux(r,uls) if Coef has integrate: (Coef,Symbol) -> Coef and _ Coef has variables: Coef -> List Symbol then integrate(upxs:%,s:Symbol) == (s = variable(upxs)) => integrate upxs not entry?(s,variables center upxs) => map(integrate(#1,s),upxs) error "integrate: center is a function of variable of integration" if Coef has TranscendentalFunctionCategory and _ Coef has PrimitiveFunctionCategory and _ Coef has AlgebraicallyClosedFunctionSpace Integer then integrateWithOneAnswer: (Coef,Symbol) -> Coef integrateWithOneAnswer(f,s) == res := integrate(f,s)$FunctionSpaceIntegration(I,Coef) res case Coef => res :: Coef first(res :: List Coef) integrate(upxs:%,s:Symbol) == (s = variable(upxs)) => integrate upxs not entry?(s,variables center upxs) => map(integrateWithOneAnswer(#1,s),upxs) error "integrate: center is a function of variable of integration" if Coef has Field then (upxs:%) ** (q:RN) == num := numer q; den := denom q one? den => upxs ** num r := rationalPower upxs; uls := laurentRep upxs deg := degree uls if zero?(coef := coefficient(uls,deg)) then deg := order(uls,deg + 1000) zero?(coef := coefficient(uls,deg)) => error "power of series with many leading zero coefficients" ulsPow := (uls * monomial(1,-deg)$ULS) ** q puiseux(r,ulsPow) * monomial(1,deg*q*r) applyUnary: (ULS -> ULS,%) -> % applyUnary(fcn,upxs) == puiseux(rationalPower upxs,fcn laurentRep upxs) exp upxs == applyUnary(exp,upxs) log upxs == applyUnary(log,upxs) sin upxs == applyUnary(sin,upxs) cos upxs == applyUnary(cos,upxs) tan upxs == applyUnary(tan,upxs) cot upxs == applyUnary(cot,upxs) sec upxs == applyUnary(sec,upxs) csc upxs == applyUnary(csc,upxs) asin upxs == applyUnary(asin,upxs) acos upxs == applyUnary(acos,upxs) atan upxs == applyUnary(atan,upxs) acot upxs == applyUnary(acot,upxs) asec upxs == applyUnary(asec,upxs) acsc upxs == applyUnary(acsc,upxs) sinh upxs == applyUnary(sinh,upxs) cosh upxs == applyUnary(cosh,upxs) tanh upxs == applyUnary(tanh,upxs) coth upxs == applyUnary(coth,upxs) sech upxs == applyUnary(sech,upxs) csch upxs == applyUnary(csch,upxs) asinh upxs == applyUnary(asinh,upxs) acosh upxs == applyUnary(acosh,upxs) atanh upxs == applyUnary(atanh,upxs) acoth upxs == applyUnary(acoth,upxs) asech upxs == applyUnary(asech,upxs) acsch upxs == applyUnary(acsch,upxs) @ \section{domain UPXS UnivariatePuiseuxSeries} <>= )abbrev domain UPXS UnivariatePuiseuxSeries ++ Author: Clifton J. Williamson ++ Date Created: 28 January 1990 ++ Date Last Updated: June 18, 2010 ++ Basic Operations: ++ Related Domains: ++ Also See: ++ AMS Classifications: ++ Keywords: series, Puiseux ++ Examples: ++ References: ++ Description: Dense Puiseux series in one variable ++ \spadtype{UnivariatePuiseuxSeries} is a domain representing Puiseux ++ series in one variable with coefficients in an arbitrary ring. The ++ parameters of the type specify the coefficient ring, the power series ++ variable, and the center of the power series expansion. For example, ++ \spad{UnivariatePuiseuxSeries(Integer,x,3)} represents Puiseux series in ++ \spad{(x - 3)} with \spadtype{Integer} coefficients. UnivariatePuiseuxSeries(Coef,var,cen): Exports == Implementation where Coef : Ring var : Symbol cen : Coef I ==> Integer L ==> List NNI ==> NonNegativeInteger OUT ==> OutputForm RN ==> Fraction Integer ST ==> Stream Coef UTS ==> UnivariateTaylorSeries(Coef,var,cen) ULS ==> UnivariateLaurentSeries(Coef,var,cen) Exports ==> Join(UnivariatePuiseuxSeriesConstructorCategory(Coef,ULS),_ PartialDifferentialDomain(%,Variable var),_ RetractableTo UTS,CoercibleFrom Variable var) with if Coef has Algebra Fraction Integer then integrate: (%,Variable(var)) -> % ++ \spad{integrate(f(x))} returns an anti-derivative of the power ++ series \spad{f(x)} with constant coefficient 0. ++ We may integrate a series when we can divide coefficients ++ by integers. Implementation ==> UnivariatePuiseuxSeriesConstructor(Coef,ULS) add Rep := Record(expon:RN,lSeries:ULS) getExpon: % -> RN getExpon pxs == pxs.expon variable upxs == var center upxs == cen coerce(uts:UTS) == uts :: ULS :: % retractIfCan(upxs:%):Union(UTS,"failed") == (ulsIfCan := retractIfCan(upxs)@Union(ULS,"failed")) case "failed" => "failed" retractIfCan(ulsIfCan :: ULS) --retract(upxs:%):UTS == --(ulsIfCan := retractIfCan(upxs)@Union(ULS,"failed")) case "failed" => --error "retractIfCan: series has fractional exponents" --utsIfCan := retractIfCan(ulsIfCan :: ULS)@Union(UTS,"failed") --utsIfCan case "failed" => --error "retractIfCan: series has negative exponents" --utsIfCan :: UTS coerce(v:Variable(var)) == zero? cen => monomial(1,1) monomial(1,1) + monomial(cen,0) if Coef has "*": (Fraction Integer, Coef) -> Coef then differentiate(upxs:%,v:Variable(var)) == differentiate upxs if Coef has Algebra Fraction Integer then integrate(upxs:%,v:Variable(var)) == integrate upxs if Coef has coerce: Symbol -> Coef then if Coef has "**": (Coef,RN) -> Coef then roundDown: RN -> I roundDown rn == -- returns the largest integer <= rn (den := denom rn) = 1 => numer rn n := (num := numer rn) quo den positive?(num) => n n - 1 stToCoef: (ST,Coef,NNI,NNI) -> Coef stToCoef(st,term,n,n0) == (n > n0) or (empty? st) => 0 frst(st) * term ** n + stToCoef(rst st,term,n + 1,n0) approximateLaurent: (ULS,Coef,I) -> Coef approximateLaurent(x,term,n) == negative?(m := n - (e := degree x)) => 0 app := stToCoef(coefficients taylorRep x,term,0,m :: NNI) zero? e => app app * term ** (e :: RN) approximate(x,r) == e := rationalPower(x) term := ((variable(x) :: Coef) - center(x)) ** e approximateLaurent(laurentRep x,term,roundDown(r / e)) termOutput:(RN,Coef,OUT) -> OUT termOutput(k,c,vv) == -- creates a term c * vv ** k k = 0 => c :: OUT mon := k = 1 => vv vv ** (k :: OUT) c = 1 => mon c = -1 => -mon (c :: OUT) * mon showAll?:() -> Boolean -- check a global Lisp variable showAll?() == true termsToOutputForm:(RN,RN,ST,OUT) -> OUT termsToOutputForm(m,rat,uu,xxx) == l : L OUT := empty() empty? uu => 0@Coef :: OUT count : NNI := _$streamCount$Lisp n : NNI := 0 while n <= count and not empty? uu repeat if frst(uu) ~= 0 then l := concat(termOutput((n :: I) * rat + m,frst uu,xxx),l) uu := rst uu n := n + 1 if showAll?() then n := count + 1 while explicitEntries? uu and _ not eq?(uu,rst uu) repeat if frst(uu) ~= 0 then l := concat(termOutput((n :: I) * rat + m,frst uu,xxx),l) uu := rst uu n := n + 1 l := explicitlyEmpty? uu => l eq?(uu,rst uu) and frst uu = 0 => l concat(prefix("O" :: OUT,[xxx ** (((n::I) * rat + m) :: OUT)]),l) empty? l => 0@Coef :: OUT reduce("+",reverse! l) coerce(upxs:%):OUT == rat := getExpon upxs; uls := laurentRep upxs count : I := _$streamCount$Lisp uls := removeZeroes(_$streamCount$Lisp,uls) m : RN := (degree uls) * rat p := coefficients taylorRep uls xxx := zero? cen => var :: OUT paren(var :: OUT - cen :: OUT) termsToOutputForm(m,rat,p,xxx) @ \section{package UPXS2 UnivariatePuiseuxSeriesFunctions2} <>= )abbrev package UPXS2 UnivariatePuiseuxSeriesFunctions2 ++ Mapping package for univariate Puiseux series ++ Author: Scott C. Morrison ++ Date Created: 5 April 1991 ++ Date Last Updated: 5 April 1991 ++ Keywords: Puiseux series, map ++ Examples: ++ References: ++ Description: ++ Mapping package for univariate Puiseux series. ++ This package allows one to apply a function to the coefficients of ++ a univariate Puiseux series. UnivariatePuiseuxSeriesFunctions2(Coef1,Coef2,var1,var2,cen1,cen2):_ Exports == Implementation where Coef1 : Ring Coef2 : Ring var1: Symbol var2: Symbol cen1: Coef1 cen2: Coef2 UPS1 ==> UnivariatePuiseuxSeries(Coef1, var1, cen1) UPS2 ==> UnivariatePuiseuxSeries(Coef2, var2, cen2) ULSP2 ==> UnivariateLaurentSeriesFunctions2(Coef1, Coef2, var1, var2, cen1, cen2) Exports ==> with map: (Coef1 -> Coef2,UPS1) -> UPS2 ++ \spad{map(f,g(x))} applies the map f to the coefficients of the ++ Puiseux series \spad{g(x)}. Implementation ==> add map(f,ups) == puiseux(rationalPower ups, map(f, laurentRep ups)$ULSP2) @ \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}