\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra fs2ups.spad} \author{Clifton J. Williamson} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package FS2UPS FunctionSpaceToUnivariatePowerSeries} <>= )abbrev package FS2UPS FunctionSpaceToUnivariatePowerSeries ++ Author: Clifton J. Williamson ++ Date Created: 21 March 1989 ++ Date Last Updated: 2 December 1994 ++ Basic Operations: ++ Related Domains: ++ Also See: ++ AMS Classifications: ++ Keywords: elementary function, power series ++ Examples: ++ References: ++ Description: ++ This package converts expressions in some function space to power ++ series in a variable x with coefficients in that function space. ++ The function \spadfun{exprToUPS} converts expressions to power series ++ whose coefficients do not contain the variable x. The function ++ \spadfun{exprToGenUPS} converts functional expressions to power series ++ whose coefficients may involve functions of \spad{log(x)}. FunctionSpaceToUnivariatePowerSeries(R,FE,Expon,UPS,TRAN,x):_ Exports == Implementation where R : Join(GcdDomain,RetractableTo Integer,_ LinearlyExplicitRingOver Integer) FE : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_ FunctionSpace R) with coerce: Expon -> % Expon : OrderedRing UPS : Join(UnivariatePowerSeriesCategory(FE,Expon),Field,_ TranscendentalFunctionCategory) with differentiate: % -> % ++ differentiate(x) returns the derivative of x since we ++ need to be able to differentiate a power series integrate: % -> % ++ integrate(x) returns the integral of x since ++ we need to be able to integrate a power series TRAN : PartialTranscendentalFunctions UPS x : Symbol B ==> Boolean BOP ==> BasicOperator I ==> Integer NNI ==> NonNegativeInteger K ==> Kernel FE L ==> List RN ==> Fraction Integer S ==> String SY ==> Symbol PCL ==> PolynomialCategoryLifting(IndexedExponents K,K,R,SMP,FE) POL ==> Polynomial R SMP ==> SparseMultivariatePolynomial(R,K) SUP ==> SparseUnivariatePolynomial Polynomial R Problem ==> Record(func:String,prob:String) Result ==> Union(%series:UPS,%problem:Problem) SIGNEF ==> ElementaryFunctionSign(R,FE) Exports ==> with exprToUPS : (FE,B,S) -> Result ++ exprToUPS(fcn,posCheck?,atanFlag) converts the expression ++ \spad{fcn} to a power series. If \spad{posCheck?} is true, ++ log's of negative numbers are not allowed nor are nth roots of ++ negative numbers with n even. If \spad{posCheck?} is false, ++ these are allowed. \spad{atanFlag} determines how the case ++ \spad{atan(f(x))}, where \spad{f(x)} has a pole, will be treated. ++ The possible values of \spad{atanFlag} are \spad{"complex"}, ++ \spad{"real: two sides"}, \spad{"real: left side"}, ++ \spad{"real: right side"}, and \spad{"just do it"}. ++ If \spad{atanFlag} is \spad{"complex"}, then no series expansion ++ will be computed because, viewed as a function of a complex ++ variable, \spad{atan(f(x))} has an essential singularity. ++ Otherwise, the sign of the leading coefficient of the series ++ expansion of \spad{f(x)} determines the constant coefficient ++ in the series expansion of \spad{atan(f(x))}. If this sign cannot ++ be determined, a series expansion is computed only when ++ \spad{atanFlag} is \spad{"just do it"}. When the leading term ++ in the series expansion of \spad{f(x)} is of odd degree (or is a ++ rational degree with odd numerator), then the constant coefficient ++ in the series expansion of \spad{atan(f(x))} for values to the ++ left differs from that for values to the right. If \spad{atanFlag} ++ is \spad{"real: two sides"}, no series expansion will be computed. ++ If \spad{atanFlag} is \spad{"real: left side"} the constant ++ coefficient for values to the left will be used and if \spad{atanFlag} ++ \spad{"real: right side"} the constant coefficient for values to the ++ right will be used. ++ If there is a problem in converting the function to a power series, ++ a record containing the name of the function that caused the problem ++ and a brief description of the problem is returned. ++ When expanding the expression into a series it is assumed that ++ the series is centered at 0. For a series centered at a, the ++ user should perform the substitution \spad{x -> x + a} before calling ++ this function. exprToGenUPS : (FE,B,S) -> Result ++ exprToGenUPS(fcn,posCheck?,atanFlag) converts the expression ++ \spad{fcn} to a generalized power series. If \spad{posCheck?} ++ is true, log's of negative numbers are not allowed nor are nth roots ++ of negative numbers with n even. If \spad{posCheck?} is false, ++ these are allowed. \spad{atanFlag} determines how the case ++ \spad{atan(f(x))}, where \spad{f(x)} has a pole, will be treated. ++ The possible values of \spad{atanFlag} are \spad{"complex"}, ++ \spad{"real: two sides"}, \spad{"real: left side"}, ++ \spad{"real: right side"}, and \spad{"just do it"}. ++ If \spad{atanFlag} is \spad{"complex"}, then no series expansion ++ will be computed because, viewed as a function of a complex ++ variable, \spad{atan(f(x))} has an essential singularity. ++ Otherwise, the sign of the leading coefficient of the series ++ expansion of \spad{f(x)} determines the constant coefficient ++ in the series expansion of \spad{atan(f(x))}. If this sign cannot ++ be determined, a series expansion is computed only when ++ \spad{atanFlag} is \spad{"just do it"}. When the leading term ++ in the series expansion of \spad{f(x)} is of odd degree (or is a ++ rational degree with odd numerator), then the constant coefficient ++ in the series expansion of \spad{atan(f(x))} for values to the ++ left differs from that for values to the right. If \spad{atanFlag} ++ is \spad{"real: two sides"}, no series expansion will be computed. ++ If \spad{atanFlag} is \spad{"real: left side"} the constant ++ coefficient for values to the left will be used and if \spad{atanFlag} ++ \spad{"real: right side"} the constant coefficient for values to the ++ right will be used. ++ If there is a problem in converting the function to a power ++ series, we return a record containing the name of the function ++ that caused the problem and a brief description of the problem. ++ When expanding the expression into a series it is assumed that ++ the series is centered at 0. For a series centered at a, the ++ user should perform the substitution \spad{x -> x + a} before calling ++ this function. localAbs: FE -> FE ++ localAbs(fcn) = \spad{abs(fcn)} or \spad{sqrt(fcn**2)} depending ++ on whether or not FE has a function \spad{abs}. This should be ++ a local function, but the compiler won't allow it. Implementation ==> add ratIfCan : FE -> Union(RN,"failed") carefulNthRootIfCan : (UPS,NNI,B,B) -> Result stateProblem : (S,S) -> Result polyToUPS : SUP -> UPS listToUPS : (L FE,(FE,B,S) -> Result,B,S,UPS,(UPS,UPS) -> UPS)_ -> Result isNonTrivPower : FE -> Union(Record(val:FE,exponent:I),"failed") powerToUPS : (FE,I,B,S) -> Result kernelToUPS : (K,B,S) -> Result nthRootToUPS : (FE,NNI,B,S) -> Result logToUPS : (FE,B,S) -> Result atancotToUPS : (FE,B,S,I) -> Result applyIfCan : (UPS -> Union(UPS,"failed"),FE,S,B,S) -> Result tranToUPS : (K,FE,B,S) -> Result powToUPS : (L FE,B,S) -> Result newElem : FE -> FE smpElem : SMP -> FE k2Elem : K -> FE contOnReals? : S -> B bddOnReals? : S -> B iExprToGenUPS : (FE,B,S) -> Result opsInvolvingX : FE -> L BOP opInOpList? : (SY,L BOP) -> B exponential? : FE -> B productOfNonZeroes? : FE -> B powerToGenUPS : (FE,I,B,S) -> Result kernelToGenUPS : (K,B,S) -> Result nthRootToGenUPS : (FE,NNI,B,S) -> Result logToGenUPS : (FE,B,S) -> Result expToGenUPS : (FE,B,S) -> Result expGenUPS : (UPS,B,S) -> Result atancotToGenUPS : (FE,FE,B,S,I) -> Result genUPSApplyIfCan : (UPS -> Union(UPS,"failed"),FE,S,B,S) -> Result applyBddIfCan : (FE,UPS -> Union(UPS,"failed"),FE,S,B,S) -> Result tranToGenUPS : (K,FE,B,S) -> Result powToGenUPS : (L FE,B,S) -> Result ZEROCOUNT : I := 1000 -- number of zeroes to be removed when taking logs or nth roots ratIfCan fcn == retractIfCan(fcn)@Union(RN,"failed") carefulNthRootIfCan(ups,n,posCheck?,rightOnly?) == -- similar to 'nthRootIfCan', but it is fussy about the series -- it takes as an argument. If 'n' is EVEN and 'posCheck?' -- is truem then the leading coefficient of the series must -- be POSITIVE. In this case, if 'rightOnly?' is false, the -- order of the series must be zero. The idea is that the -- series represents a real function of a real variable, and -- we want a unique real nth root defined on a neighborhood -- of zero. n < 1 => error "nthRoot: n must be positive" deg := degree ups if (coef := coefficient(ups,deg)) = 0 then deg := order(ups,deg + ZEROCOUNT :: Expon) (coef := coefficient(ups,deg)) = 0 => error "log of series with many leading zero coefficients" -- if 'posCheck?' is true, we do not allow nth roots of negative -- numbers when n in even if even?(n :: I) then if posCheck? and ((signum := sign(coef)$SIGNEF) case I) then (signum :: I) = -1 => return stateProblem("nth root","negative leading coefficient") not rightOnly? and not zero? deg => -- nth root not unique return stateProblem("nth root","series of non-zero order") (ans := nthRootIfCan(ups,n)) case "failed" => stateProblem("nth root","no nth root") [ans :: UPS] stateProblem(function,problem) == -- records the problem which occured in converting an expression -- to a power series [[function,problem]] exprToUPS(fcn,posCheck?,atanFlag) == -- converts a functional expression to a power series --!! The following line is commented out so that expressions of --!! the form a**b will be normalized to exp(b * log(a)) even if --!! 'a' and 'b' do not involve the limiting variable 'x'. --!! - cjw 1 Dec 94 --not member?(x,variables fcn) => [monomial(fcn,0)] (poly := retractIfCan(fcn)@Union(POL,"failed")) case POL => [polyToUPS univariate(poly :: POL,x)] (sum := isPlus fcn) case L(FE) => listToUPS(sum :: L(FE),exprToUPS,posCheck?,atanFlag,0,#1 + #2) (prod := isTimes fcn) case L(FE) => listToUPS(prod :: L(FE),exprToUPS,posCheck?,atanFlag,1,#1 * #2) (expt := isNonTrivPower fcn) case Record(val:FE,exponent:I) => power := expt :: Record(val:FE,exponent:I) powerToUPS(power.val,power.exponent,posCheck?,atanFlag) (ker := retractIfCan(fcn)@Union(K,"failed")) case K => kernelToUPS(ker :: K,posCheck?,atanFlag) error "exprToUPS: neither a sum, product, power, nor kernel" polyToUPS poly == -- converts a polynomial to a power series zero? poly => 0 -- we don't start with 'ans := 0' as this may lead to an -- enormous number of leading zeroes in the power series deg := degree poly coef := leadingCoefficient(poly) :: FE ans := monomial(coef,deg :: Expon)$UPS poly := reductum poly while not zero? poly repeat deg := degree poly coef := leadingCoefficient(poly) :: FE ans := ans + monomial(coef,deg :: Expon)$UPS poly := reductum poly ans listToUPS(list,feToUPS,posCheck?,atanFlag,ans,op) == -- converts each element of a list of expressions to a power -- series and returns the sum of these series, when 'op' is + -- and 'ans' is 0, or the product of these series, when 'op' is * -- and 'ans' is 1 while not null list repeat (term := feToUPS(first list,posCheck?,atanFlag)) case %problem => return term ans := op(ans,term.%series) list := rest list [ans] isNonTrivPower fcn == -- is the function a power with exponent other than 0 or 1? (expt := isPower fcn) case "failed" => "failed" power := expt :: Record(val:FE,exponent:I) one? power.exponent => "failed" power powerToUPS(fcn,n,posCheck?,atanFlag) == -- converts an integral power to a power series (b := exprToUPS(fcn,posCheck?,atanFlag)) case %problem => b positive? n => [(b.%series) ** n] -- check lowest order coefficient when n < 0 ups := b.%series; deg := degree ups if (coef := coefficient(ups,deg)) = 0 then deg := order(ups,deg + ZEROCOUNT :: Expon) (coef := coefficient(ups,deg)) = 0 => error "inverse of series with many leading zero coefficients" [ups ** n] kernelToUPS(ker,posCheck?,atanFlag) == -- converts a kernel to a power series (sym := symbolIfCan(ker)) case Symbol => (sym :: Symbol) = x => [monomial(1,1)] [monomial(ker :: FE,0)] empty?(args := argument ker) => [monomial(ker :: FE,0)] not member?(x, variables(ker :: FE)) => [monomial(ker :: FE,0)] empty? rest args => arg := first args is?(ker,"abs" :: Symbol) => nthRootToUPS(arg*arg,2,posCheck?,atanFlag) is?(ker,"%paren" :: Symbol) => exprToUPS(arg,posCheck?,atanFlag) is?(ker,"log" :: Symbol) => logToUPS(arg,posCheck?,atanFlag) is?(ker,"exp" :: Symbol) => applyIfCan(expIfCan,arg,"exp",posCheck?,atanFlag) tranToUPS(ker,arg,posCheck?,atanFlag) is?(ker,"%power" :: Symbol) => powToUPS(args,posCheck?,atanFlag) is?(ker,"nthRoot" :: Symbol) => n := retract(second args)@I nthRootToUPS(first args,n :: NNI,posCheck?,atanFlag) stateProblem(string name operator ker,"unknown kernel") nthRootToUPS(arg,n,posCheck?,atanFlag) == -- converts an nth root to a power series -- this is not used in the limit package, so the series may -- have non-zero order, in which case nth roots may not be unique (result := exprToUPS(arg,posCheck?,atanFlag)) case %problem => result ans := carefulNthRootIfCan(result.%series,n,posCheck?,false) ans case %problem => ans [ans.%series] logToUPS(arg,posCheck?,atanFlag) == -- converts a logarithm log(f(x)) to a power series -- f(x) must have order 0 and if 'posCheck?' is true, -- then f(x) must have a non-negative leading coefficient (result := exprToUPS(arg,posCheck?,atanFlag)) case %problem => result ups := result.%series not zero? order(ups,1) => stateProblem("log","series of non-zero order") coef := coefficient(ups,0) -- if 'posCheck?' is true, we do not allow logs of negative numbers if posCheck? then if ((signum := sign(coef)$SIGNEF) case I) then (signum :: I) = -1 => return stateProblem("log","negative leading coefficient") [logIfCan(ups) :: UPS] if FE has abs: FE -> FE then localAbs fcn == abs fcn else localAbs fcn == sqrt(fcn * fcn) signOfExpression: FE -> FE signOfExpression arg == localAbs(arg)/arg atancotToUPS(arg,posCheck?,atanFlag,plusMinus) == -- converts atan(f(x)) to a power series (result := exprToUPS(arg,posCheck?,atanFlag)) case %problem => result ups := result.%series; coef := coefficient(ups,0) (ord := order(ups,0)) = 0 and coef * coef = -1 => -- series involves complex numbers return stateProblem("atan","logarithmic singularity") cc : FE := negative? ord => atanFlag = "complex" => return stateProblem("atan","essential singularity") (rn := ratIfCan(ord :: FE)) case "failed" => -- this condition usually won't occur because exponents will -- be integers or rational numbers return stateProblem("atan","branch problem") if (atanFlag = "real: two sides") and (odd? numer(rn :: RN)) then -- expansions to the left and right of zero have different -- constant coefficients return stateProblem("atan","branch problem") lc := coefficient(ups,ord) (signum := sign(lc)$SIGNEF) case "failed" => -- can't determine sign atanFlag = "just do it" => plusMinus = 1 => pi()/(2 :: FE) 0 posNegPi2 := signOfExpression(lc) * pi()/(2 :: FE) plusMinus = 1 => posNegPi2 pi()/(2 :: FE) - posNegPi2 --return stateProblem("atan","branch problem") left? : B := atanFlag = "real: left side"; n := signum :: Integer (left? and n = 1) or (not left? and n = -1) => plusMinus = 1 => -pi()/(2 :: FE) pi() plusMinus = 1 => pi()/(2 :: FE) 0 atan coef [(cc :: UPS) + integrate(plusMinus * differentiate(ups)/(1 + ups*ups))] applyIfCan(fcn,arg,fcnName,posCheck?,atanFlag) == -- converts fcn(arg) to a power series (ups := exprToUPS(arg,posCheck?,atanFlag)) case %problem => ups ans := fcn(ups.%series) ans case "failed" => stateProblem(fcnName,"essential singularity") [ans :: UPS] tranToUPS(ker,arg,posCheck?,atanFlag) == -- converts ker to a power series for certain functions -- in trig or hyperbolic trig categories is?(ker,"sin" :: SY) => applyIfCan(sinIfCan,arg,"sin",posCheck?,atanFlag) is?(ker,"cos" :: SY) => applyIfCan(cosIfCan,arg,"cos",posCheck?,atanFlag) is?(ker,"tan" :: SY) => applyIfCan(tanIfCan,arg,"tan",posCheck?,atanFlag) is?(ker,"cot" :: SY) => applyIfCan(cotIfCan,arg,"cot",posCheck?,atanFlag) is?(ker,"sec" :: SY) => applyIfCan(secIfCan,arg,"sec",posCheck?,atanFlag) is?(ker,"csc" :: SY) => applyIfCan(cscIfCan,arg,"csc",posCheck?,atanFlag) is?(ker,"asin" :: SY) => applyIfCan(asinIfCan,arg,"asin",posCheck?,atanFlag) is?(ker,"acos" :: SY) => applyIfCan(acosIfCan,arg,"acos",posCheck?,atanFlag) is?(ker,"atan" :: SY) => atancotToUPS(arg,posCheck?,atanFlag,1) is?(ker,"acot" :: SY) => atancotToUPS(arg,posCheck?,atanFlag,-1) is?(ker,"asec" :: SY) => applyIfCan(asecIfCan,arg,"asec",posCheck?,atanFlag) is?(ker,"acsc" :: SY) => applyIfCan(acscIfCan,arg,"acsc",posCheck?,atanFlag) is?(ker,"sinh" :: SY) => applyIfCan(sinhIfCan,arg,"sinh",posCheck?,atanFlag) is?(ker,"cosh" :: SY) => applyIfCan(coshIfCan,arg,"cosh",posCheck?,atanFlag) is?(ker,"tanh" :: SY) => applyIfCan(tanhIfCan,arg,"tanh",posCheck?,atanFlag) is?(ker,"coth" :: SY) => applyIfCan(cothIfCan,arg,"coth",posCheck?,atanFlag) is?(ker,"sech" :: SY) => applyIfCan(sechIfCan,arg,"sech",posCheck?,atanFlag) is?(ker,"csch" :: SY) => applyIfCan(cschIfCan,arg,"csch",posCheck?,atanFlag) is?(ker,"asinh" :: SY) => applyIfCan(asinhIfCan,arg,"asinh",posCheck?,atanFlag) is?(ker,"acosh" :: SY) => applyIfCan(acoshIfCan,arg,"acosh",posCheck?,atanFlag) is?(ker,"atanh" :: SY) => applyIfCan(atanhIfCan,arg,"atanh",posCheck?,atanFlag) is?(ker,"acoth" :: SY) => applyIfCan(acothIfCan,arg,"acoth",posCheck?,atanFlag) is?(ker,"asech" :: SY) => applyIfCan(asechIfCan,arg,"asech",posCheck?,atanFlag) is?(ker,"acsch" :: SY) => applyIfCan(acschIfCan,arg,"acsch",posCheck?,atanFlag) stateProblem(string name operator ker,"unknown kernel") powToUPS(args,posCheck?,atanFlag) == -- converts a power f(x) ** g(x) to a power series (logBase := logToUPS(first args,posCheck?,atanFlag)) case %problem => logBase (expon := exprToUPS(second args,posCheck?,atanFlag)) case %problem => expon ans := expIfCan((expon.%series) * (logBase.%series)) ans case "failed" => stateProblem("exp","essential singularity") [ans :: UPS] -- Generalized power series: power series in x, where log(x) and -- bounded functions of x are allowed to appear in the coefficients -- of the series. Used for evaluating REAL limits at x = 0. newElem f == -- rewrites a functional expression; all trig functions are -- expressed in terms of sin and cos; all hyperbolic trig -- functions are expressed in terms of exp smpElem(numer f) / smpElem(denom f) smpElem p == map(k2Elem,#1::FE,p)$PCL k2Elem k == -- rewrites a kernel; all trig functions are -- expressed in terms of sin and cos; all hyperbolic trig -- functions are expressed in terms of exp null(args := [newElem a for a in argument k]) => k::FE iez := inv(ez := exp(z := first args)) sinz := sin z; cosz := cos z is?(k,"tan" :: Symbol) => sinz / cosz is?(k,"cot" :: Symbol) => cosz / sinz is?(k,"sec" :: Symbol) => inv cosz is?(k,"csc" :: Symbol) => inv sinz is?(k,"sinh" :: Symbol) => (ez - iez) / (2 :: FE) is?(k,"cosh" :: Symbol) => (ez + iez) / (2 :: FE) is?(k,"tanh" :: Symbol) => (ez - iez) / (ez + iez) is?(k,"coth" :: Symbol) => (ez + iez) / (ez - iez) is?(k,"sech" :: Symbol) => 2 * inv(ez + iez) is?(k,"csch" :: Symbol) => 2 * inv(ez - iez) (operator k) args CONTFCNS : L S := ["sin","cos","atan","acot","exp","asinh"] -- functions which are defined and continuous at all real numbers BDDFCNS : L S := ["sin","cos","atan","acot"] -- functions which are bounded on the reals contOnReals? fcn == member?(fcn,CONTFCNS) bddOnReals? fcn == member?(fcn,BDDFCNS) exprToGenUPS(fcn,posCheck?,atanFlag) == -- converts a functional expression to a generalized power -- series; "generalized" means that log(x) and bounded functions -- of x are allowed to appear in the coefficients of the series iExprToGenUPS(newElem fcn,posCheck?,atanFlag) iExprToGenUPS(fcn,posCheck?,atanFlag) == -- converts a functional expression to a generalized power -- series without first normalizing the expression --!! The following line is commented out so that expressions of --!! the form a**b will be normalized to exp(b * log(a)) even if --!! 'a' and 'b' do not involve the limiting variable 'x'. --!! - cjw 1 Dec 94 --not member?(x,variables fcn) => [monomial(fcn,0)] (poly := retractIfCan(fcn)@Union(POL,"failed")) case POL => [polyToUPS univariate(poly :: POL,x)] (sum := isPlus fcn) case L(FE) => listToUPS(sum :: L(FE),iExprToGenUPS,posCheck?,atanFlag,0,#1 + #2) (prod := isTimes fcn) case L(FE) => listToUPS(prod :: L(FE),iExprToGenUPS,posCheck?,atanFlag,1,#1 * #2) (expt := isNonTrivPower fcn) case Record(val:FE,exponent:I) => power := expt :: Record(val:FE,exponent:I) powerToGenUPS(power.val,power.exponent,posCheck?,atanFlag) (ker := retractIfCan(fcn)@Union(K,"failed")) case K => kernelToGenUPS(ker :: K,posCheck?,atanFlag) error "exprToGenUPS: neither a sum, product, power, nor kernel" opsInvolvingX fcn == opList := [op for k in tower fcn | unary?(op := operator k) _ and member?(x,variables first argument k)] removeDuplicates opList opInOpList?(name,opList) == for op in opList repeat is?(op,name) => return true false exponential? fcn == -- is 'fcn' of the form exp(f)? (ker := retractIfCan(fcn)@Union(K,"failed")) case K => is?(ker :: K,"exp" :: Symbol) false productOfNonZeroes? fcn == -- is 'fcn' a product of non-zero terms, where 'non-zero' -- means an exponential or a function not involving x exponential? fcn => true (prod := isTimes fcn) case "failed" => false for term in (prod :: L(FE)) repeat (not exponential? term) and member?(x,variables term) => return false true powerToGenUPS(fcn,n,posCheck?,atanFlag) == -- converts an integral power to a generalized power series -- if n < 0 and the lowest order coefficient of the series -- involves x, we are careful about inverting this coefficient -- the coefficient is inverted only if -- (a) the only function involving x is 'log', or -- (b) the lowest order coefficient is a product of exponentials -- and functions not involving x (b := exprToGenUPS(fcn,posCheck?,atanFlag)) case %problem => b positive? n => [(b.%series) ** n] -- check lowest order coefficient when n < 0 ups := b.%series; deg := degree ups if (coef := coefficient(ups,deg)) = 0 then deg := order(ups,deg + ZEROCOUNT :: Expon) (coef := coefficient(ups,deg)) = 0 => error "inverse of series with many leading zero coefficients" xOpList := opsInvolvingX coef -- only function involving x is 'log' (null xOpList) => [ups ** n] (null rest xOpList and is?(first xOpList,"log" :: SY)) => [ups ** n] -- lowest order coefficient is a product of exponentials and -- functions not involving x productOfNonZeroes? coef => [ups ** n] stateProblem("inv","lowest order coefficient involves x") kernelToGenUPS(ker,posCheck?,atanFlag) == -- converts a kernel to a generalized power series (sym := symbolIfCan(ker)) case Symbol => (sym :: Symbol) = x => [monomial(1,1)] [monomial(ker :: FE,0)] empty?(args := argument ker) => [monomial(ker :: FE,0)] empty? rest args => arg := first args is?(ker,"abs" :: Symbol) => nthRootToGenUPS(arg*arg,2,posCheck?,atanFlag) is?(ker,"%paren" :: Symbol) => iExprToGenUPS(arg,posCheck?,atanFlag) is?(ker,"log" :: Symbol) => logToGenUPS(arg,posCheck?,atanFlag) is?(ker,"exp" :: Symbol) => expToGenUPS(arg,posCheck?,atanFlag) tranToGenUPS(ker,arg,posCheck?,atanFlag) is?(ker,"%power" :: Symbol) => powToGenUPS(args,posCheck?,atanFlag) is?(ker,"nthRoot" :: Symbol) => n := retract(second args)@I nthRootToGenUPS(first args,n :: NNI,posCheck?,atanFlag) stateProblem(string name operator ker,"unknown kernel") nthRootToGenUPS(arg,n,posCheck?,atanFlag) == -- convert an nth root to a power series -- used for computing right hand limits, so the series may have -- non-zero order, but may not have a negative leading coefficient -- when n is even (result := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem => result ans := carefulNthRootIfCan(result.%series,n,posCheck?,true) ans case %problem => ans [ans.%series] logToGenUPS(arg,posCheck?,atanFlag) == -- converts a logarithm log(f(x)) to a generalized power series (result := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem => result ups := result.%series; deg := degree ups if (coef := coefficient(ups,deg)) = 0 then deg := order(ups,deg + ZEROCOUNT :: Expon) (coef := coefficient(ups,deg)) = 0 => error "log of series with many leading zero coefficients" -- if 'posCheck?' is true, we do not allow logs of negative numbers if posCheck? then if ((signum := sign(coef)$SIGNEF) case I) then (signum :: I) = -1 => return stateProblem("log","negative leading coefficient") -- create logarithmic term, avoiding log's of negative rationals lt := monomial(coef,deg)$UPS; cen := center lt -- check to see if lowest order coefficient is a negative rational negRat? : Boolean := ((rat := ratIfCan coef) case RN) => negative?(rat :: RN) => true false false logTerm : FE := mon : FE := (x :: FE) - (cen :: FE) pow : FE := mon ** (deg :: FE) negRat? => log(coef * pow) term1 : FE := (deg :: FE) * log(mon) log(coef) + term1 [monomial(logTerm,0) + log(ups/lt)] expToGenUPS(arg,posCheck?,atanFlag) == -- converts an exponential exp(f(x)) to a generalized -- power series (ups := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem => ups expGenUPS(ups.%series,posCheck?,atanFlag) expGenUPS(ups,posCheck?,atanFlag) == -- computes the exponential of a generalized power series. -- If the series has order zero and the constant term a0 of the -- series involves x, the function tries to expand exp(a0) as -- a power series. negative?(deg := order(ups,1)) => stateProblem("exp","essential singularity") positive? deg => [exp ups] lc := coefficient(ups,0); xOpList := opsInvolvingX lc not opInOpList?("log" :: SY,xOpList) => [exp ups] -- try to fix exp(lc) if necessary expCoef := normalize(exp lc,x)$ElementaryFunctionStructurePackage(R,FE) opInOpList?("log" :: SY,opsInvolvingX expCoef) => stateProblem("exp","logs in constant coefficient") result := exprToGenUPS(expCoef,posCheck?,atanFlag) result case %problem => result [(result.%series) * exp(ups - monomial(lc,0))] atancotToGenUPS(fe,arg,posCheck?,atanFlag,plusMinus) == -- converts atan(f(x)) to a generalized power series (result := exprToGenUPS(arg,posCheck?,atanFlag)) case %problem => trouble := result.%problem trouble.prob = "essential singularity" => [monomial(fe,0)$UPS] result ups := result.%series; coef := coefficient(ups,0) -- series involves complex numbers (ord := order(ups,0)) = 0 and coef * coef = -1 => y := differentiate(ups)/(1 + ups*ups) yCoef := coefficient(y,-1) [monomial(log yCoef,0) + integrate(y - monomial(yCoef,-1)$UPS)] cc : FE := negative? ord => atanFlag = "complex" => return stateProblem("atan","essential singularity") (rn := ratIfCan(ord :: FE)) case "failed" => -- this condition usually won't occur because exponents will -- be integers or rational numbers return stateProblem("atan","branch problem") if (atanFlag = "real: two sides") and (odd? numer(rn :: RN)) then -- expansions to the left and right of zero have different -- constant coefficients return stateProblem("atan","branch problem") lc := coefficient(ups,ord) (signum := sign(lc)$SIGNEF) case "failed" => -- can't determine sign atanFlag = "just do it" => plusMinus = 1 => pi()/(2 :: FE) 0 posNegPi2 := signOfExpression(lc) * pi()/(2 :: FE) plusMinus = 1 => posNegPi2 pi()/(2 :: FE) - posNegPi2 --return stateProblem("atan","branch problem") left? : B := atanFlag = "real: left side"; n := signum :: Integer (left? and n = 1) or (not left? and n = -1) => plusMinus = 1 => -pi()/(2 :: FE) pi() plusMinus = 1 => pi()/(2 :: FE) 0 atan coef [(cc :: UPS) + integrate(differentiate(ups)/(1 + ups*ups))] genUPSApplyIfCan(fcn,arg,fcnName,posCheck?,atanFlag) == -- converts fcn(arg) to a generalized power series (series := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem => series ups := series.%series negative?(deg := order(ups,1)) => stateProblem(fcnName,"essential singularity") positive? deg => [fcn(ups) :: UPS] lc := coefficient(ups,0); xOpList := opsInvolvingX lc null xOpList => [fcn(ups) :: UPS] opInOpList?("log" :: SY,xOpList) => stateProblem(fcnName,"logs in constant coefficient") contOnReals? fcnName => [fcn(ups) :: UPS] stateProblem(fcnName,"x in constant coefficient") applyBddIfCan(fe,fcn,arg,fcnName,posCheck?,atanFlag) == -- converts fcn(arg) to a generalized power series, where the -- function fcn is bounded for real values -- if fcn(arg) has an essential singularity as a complex -- function, we return fcn(arg) as a monomial of degree 0 (ups := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem => trouble := ups.%problem trouble.prob = "essential singularity" => [monomial(fe,0)$UPS] ups (ans := fcn(ups.%series)) case "failed" => [monomial(fe,0)$UPS] [ans :: UPS] tranToGenUPS(ker,arg,posCheck?,atanFlag) == -- converts op(arg) to a power series for certain functions -- op in trig or hyperbolic trig categories -- N.B. when this function is called, 'k2elem' will have been -- applied, so the following functions cannot appear: -- tan, cot, sec, csc, sinh, cosh, tanh, coth, sech, csch is?(ker,"sin" :: SY) => applyBddIfCan(ker :: FE,sinIfCan,arg,"sin",posCheck?,atanFlag) is?(ker,"cos" :: SY) => applyBddIfCan(ker :: FE,cosIfCan,arg,"cos",posCheck?,atanFlag) is?(ker,"asin" :: SY) => genUPSApplyIfCan(asinIfCan,arg,"asin",posCheck?,atanFlag) is?(ker,"acos" :: SY) => genUPSApplyIfCan(acosIfCan,arg,"acos",posCheck?,atanFlag) is?(ker,"atan" :: SY) => atancotToGenUPS(ker :: FE,arg,posCheck?,atanFlag,1) is?(ker,"acot" :: SY) => atancotToGenUPS(ker :: FE,arg,posCheck?,atanFlag,-1) is?(ker,"asec" :: SY) => genUPSApplyIfCan(asecIfCan,arg,"asec",posCheck?,atanFlag) is?(ker,"acsc" :: SY) => genUPSApplyIfCan(acscIfCan,arg,"acsc",posCheck?,atanFlag) is?(ker,"asinh" :: SY) => genUPSApplyIfCan(asinhIfCan,arg,"asinh",posCheck?,atanFlag) is?(ker,"acosh" :: SY) => genUPSApplyIfCan(acoshIfCan,arg,"acosh",posCheck?,atanFlag) is?(ker,"atanh" :: SY) => genUPSApplyIfCan(atanhIfCan,arg,"atanh",posCheck?,atanFlag) is?(ker,"acoth" :: SY) => genUPSApplyIfCan(acothIfCan,arg,"acoth",posCheck?,atanFlag) is?(ker,"asech" :: SY) => genUPSApplyIfCan(asechIfCan,arg,"asech",posCheck?,atanFlag) is?(ker,"acsch" :: SY) => genUPSApplyIfCan(acschIfCan,arg,"acsch",posCheck?,atanFlag) stateProblem(string name operator ker,"unknown kernel") powToGenUPS(args,posCheck?,atanFlag) == -- converts a power f(x) ** g(x) to a generalized power series (logBase := logToGenUPS(first args,posCheck?,atanFlag)) case %problem => logBase expon := iExprToGenUPS(second args,posCheck?,atanFlag) expon case %problem => expon expGenUPS((expon.%series) * (logBase.%series),posCheck?,atanFlag) @ \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. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}