\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra sum.spad} \author{Stephen M. Watt, Manuel Bronstein} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package ISUMP InnerPolySum} <>= )abbrev package ISUMP InnerPolySum ++ Summation of polynomials ++ Author: SMW ++ Date Created: ??? ++ Date Last Updated: 19 April 1991 ++ Description: tools for the summation packages. InnerPolySum(E, V, R, P): Exports == Impl where E: OrderedAbelianMonoidSup V: OrderedSet R: IntegralDomain P: PolynomialCategory(R, E, V) Z ==> Integer Q ==> Fraction Z SUP ==> SparseUnivariatePolynomial Exports ==> with sum: (P, V, Segment P) -> Record(num:P, den:Z) ++ sum(p(n), n = a..b) returns \spad{p(a) + p(a+1) + ... + p(b)}. sum: (P, V) -> Record(num:P, den: Z) ++ sum(p(n), n) returns \spad{P(n)}, ++ the indefinite sum of \spad{p(n)} with respect to ++ upward difference on n, i.e. \spad{P(n+1) - P(n) = a(n)}; Impl ==> add import PolynomialNumberTheoryFunctions() import UnivariatePolynomialCommonDenominator(Z, Q, SUP Q) pmul: (P, SUP Q) -> Record(num:SUP P, den:Z) pmul(c, p) == pn := (rec := splitDenominator p).num [map(numer(#1) * c, pn)$SparseUnivariatePolynomialFunctions2(Q, P), rec.den] sum(p, v, s) == indef := sum(p, v) [eval(indef.num, v, 1 + hi s) - eval(indef.num, v, lo s), indef.den] sum(p, v) == up := univariate(p, v) lp := nil()$List(SUP P) ld := nil()$List(Z) while up ~= 0 repeat ud := degree up; uc := leadingCoefficient up up := reductum up rec := pmul(uc, 1 / (ud+1) * bernoulli(ud+1)) lp := concat(rec.num, lp) ld := concat(rec.den, ld) d := lcm ld vp := +/[(d exquo di)::Z * pi for di in ld for pi in lp] [multivariate(vp, v), d] @ \section{package GOSPER GosperSummationMethod} <>= )abbrev package GOSPER GosperSummationMethod ++ Gosper's summation algorithm ++ Author: SMW ++ Date Created: ??? ++ Date Last Updated: 19 August 1991 ++ Description: Gosper's summation algorithm. GosperSummationMethod(E, V, R, P, Q): Exports == Impl where E: OrderedAbelianMonoidSup V: OrderedSet R: IntegralDomain P: PolynomialCategory(R, E, V) Q: Join(RetractableTo Fraction Integer, Field with (coerce: P -> %; numer : % -> P; denom : % -> P)) I ==> Integer RN ==> Fraction I PQ ==> SparseMultivariatePolynomial(RN, V) RQ ==> Fraction PQ Exports ==> with GospersMethod: (Q, V, () -> V) -> Union(Q, "failed") ++ GospersMethod(b, n, new) returns a rational function ++ \spad{rf(n)} such that \spad{a(n) * rf(n)} is the indefinite ++ sum of \spad{a(n)} ++ with respect to upward difference on \spad{n}, i.e. ++ \spad{a(n+1) * rf(n+1) - a(n) * rf(n) = a(n)}, ++ where \spad{b(n) = a(n)/a(n-1)} is a rational function. ++ Returns "failed" if no such rational function \spad{rf(n)} ++ exists. ++ Note: \spad{new} is a nullary function returning a new ++ V every time. ++ The condition on \spad{a(n)} is that \spad{a(n)/a(n-1)} ++ is a rational function of \spad{n}. --++ \spad{sum(a(n), n) = rf(n) * a(n)}. Impl ==> add import PolynomialCategoryQuotientFunctions(E, V, R, P, Q) import LinearSystemMatrixPackage(RQ,Vector RQ,Vector RQ,Matrix RQ) InnerGospersMethod: (RQ, V, () -> V) -> Union(RQ, "failed") GosperPQR: (PQ, PQ, V, () -> V) -> List PQ GosperDegBd: (PQ, PQ, PQ, V, () -> V) -> I GosperF: (I, PQ, PQ, PQ, V, () -> V) -> Union(RQ, "failed") linearAndNNIntRoot: (PQ, V) -> Union(I, "failed") deg0: (PQ, V) -> I -- degree with deg 0 = -1. pCoef: (PQ, PQ) -> PQ -- pCoef(p, a*b**2) RF2QIfCan: Q -> Union(RQ, "failed") UP2QIfCan: P -> Union(PQ,"failed") RFQ2R : RQ -> Q PQ2R : PQ -> Q rat? : R -> Boolean deg0(p, v) == (zero? p => -1; degree(p, v)) rat? x == retractIfCan(x::P::Q)@Union(RN, "failed") case RN RFQ2R f == PQ2R(numer f) / PQ2R(denom f) PQ2R p == map(#1::P::Q, #1::Q, p)$PolynomialCategoryLifting( IndexedExponents V, V, RN, PQ, Q) GospersMethod(aquo, n, newV) == ((q := RF2QIfCan aquo) case "failed") or ((u := InnerGospersMethod(q::RQ, n, newV)) case "failed") => "failed" RFQ2R(u::RQ) RF2QIfCan f == (n := UP2QIfCan numer f) case "failed" => "failed" (d := UP2QIfCan denom f) case "failed" => "failed" n::PQ / d::PQ UP2QIfCan p == every?(rat?, coefficients p) => map(#1::PQ, (retractIfCan(#1::P::Q)@Union(RN, "failed"))::RN::PQ, p)$PolynomialCategoryLifting(E, V, R, P, PQ) "failed" InnerGospersMethod(aquo, n, newV) == -- 1. Define coprime polys an,anm1 such that -- an/anm1=a(n)/a(n-1) an := numer aquo anm1 := denom aquo -- 2. Define p,q,r such that -- a(n)/a(n-1) = (p(n)/p(n-1)) * (q(n)/r(n)) -- and -- gcd(q(n), r(n+j)) = 1, for all j: NNI. pqr:= GosperPQR(an, anm1, n, newV) pn := first pqr; qn := second pqr; rn := third pqr -- 3. If the sum is a rational fn, there is a poly f with -- sum(a(n), n) = q(n+1)/p(n) * a(n) * f(n). -- 4. Bound the degree of f(n). negative?(k := GosperDegBd(pn, qn, rn, n, newV)) => "failed" -- 5. Find a polynomial f of degree at most k, satisfying -- p(n) = q(n+1)*f(n) - r(n)*f(n-1) (ufn := GosperF(k, pn, qn, rn, n, newV)) case "failed" => "failed" fn := ufn::RQ -- 6. The sum is q(n-1)/p(n)*f(n) * a(n). We leave out a(n). --qnm1 := eval(qn,n,n::PQ - 1) --qnm1/pn * fn qn1 := eval(qn,n,n::PQ + 1) qn1/pn * fn GosperF(k, pn, qn, rn, n, newV) == mv := newV(); mp := mv::PQ; np := n::PQ fn: PQ := +/[mp**(i+1) * np**i for i in 0..k] fnminus1: PQ := eval(fn, n, np-1) qnplus1 := eval(qn, n, np+1) zro := qnplus1 * fn - rn * fnminus1 - pn zron := univariate(zro, n) dz := degree zron mat: Matrix RQ := zero(dz+1, (k+1)::NonNegativeInteger) vec: Vector RQ := new(dz+1, 0) while zron ~= 0 repeat cz := leadingCoefficient zron dz := degree zron zron := reductum zron mz := univariate(cz, mv) while mz ~= 0 repeat cmz := leadingCoefficient(mz)::RQ dmz := degree mz mz := reductum mz dmz = 0 => vec(dz + minIndex vec) := -cmz qsetelt!(mat, dz + minRowIndex mat, dmz + minColIndex(mat) - 1, cmz) (soln := particularSolution(mat, vec)) case "failed" => "failed" vec := soln::Vector RQ (+/[np**i * vec(i + minIndex vec) for i in 0..k])@RQ GosperPQR(an, anm1, n, newV) == np := n::PQ -- polynomial version of n -- Initial guess. pn: PQ := 1 qn: PQ := an rn: PQ := anm1 -- Find all j: NNI giving common factors to q(n) and r(n+j). j := newV() rnj := eval(rn, n, np + j::PQ) res := resultant(qn, rnj, n) fres := factor(res)$MRationalFactorize(IndexedExponents V, V, I, PQ) js := [rt::I for fe in factors fres | (rt := linearAndNNIntRoot(fe.factor,j)) case I] -- For each such j, change variables to remove the gcd. for rt in js repeat rtp:= rt::PQ -- polynomial version of rt gn := gcd(qn, eval(rn,n,np+rtp)) qn := (qn exquo gn)::PQ rn := (rn exquo eval(gn, n, np-rtp))::PQ pn := pn * */[eval(gn, n, np-i::PQ) for i in 0..rt-1] [pn, qn, rn] -- Find a degree bound for the polynomial f(n) which satisfies -- p(n) = q(n+1)*f(n) - r(n)*f(n-1). GosperDegBd(pn, qn, rn, n, newV) == np := n::PQ qnplus1 := eval(qn, n, np+1) lplus := deg0(qnplus1 + rn, n) lminus := deg0(qnplus1 - rn, n) degp := deg0(pn, n) k := degp - max(lplus-1, lminus) lplus <= lminus => k -- Find L(k), such that -- p(n) = L(k)*c[k]*n**(k + lplus - 1) + ... -- To do this, write f(n) and f(n-1) symbolically. -- f(n) = c[k]*n**k + c[k-1]*n**(k-1) +O(n**(k-2)) -- f(n-1)=c[k]*n**k + (c[k-1]-k*c[k])*n**(k-1)+O(n**(k-2)) kk := newV()::PQ ck := newV()::PQ ckm1 := newV()::PQ nkm1:= newV()::PQ nk := np*nkm1 headfn := ck*nk + ckm1*nkm1 headfnm1 := ck*nk + (ckm1-kk*ck)*nkm1 -- Then p(n) = q(n+1)*f(n) - r(n)*f(n-1) gives L(k). pk := qnplus1 * headfn - rn * headfnm1 lcpk := pCoef(pk, ck*np*nkm1) -- The degree bd is now given by k, and the root of L. k0 := linearAndNNIntRoot(lcpk, mainVariable(kk)::V) k0 case "failed" => k max(k0::I, k) pCoef(p, nom) == not monomial? nom => error "pCoef requires a monomial 2nd arg" vlist := variables nom for v in vlist while p ~= 0 repeat unom:= univariate(nom,v) pow:=degree unom nom:=leadingCoefficient unom up := univariate(p, v) p := coefficient(up, pow) p linearAndNNIntRoot(mp, v) == p := univariate(mp, v) not one? degree p => "failed" (p1 := retractIfCan(coefficient(p, 1))@Union(RN,"failed")) case "failed" or (p0 := retractIfCan(coefficient(p, 0))@Union(RN,"failed")) case "failed" => "failed" rt := -(p0::RN)/(p1::RN) negative? rt or not one? denom rt => "failed" numer rt @ \section{package SUMRF RationalFunctionSum} <>= )abbrev package SUMRF RationalFunctionSum ++ Summation of rational functions ++ Author: Manuel Bronstein ++ Date Created: ??? ++ Date Last Updated: 19 April 1991 ++ Description: Computes sums of rational functions; RationalFunctionSum(R): Exports == Impl where R: Join(IntegralDomain, RetractableTo Integer) P ==> Polynomial R RF ==> Fraction P FE ==> Expression R SE ==> Symbol Exports ==> with sum: (P, SE) -> RF ++ sum(a(n), n) returns \spad{A} which ++ is the indefinite sum of \spad{a} with respect to ++ upward difference on \spad{n}, i.e. \spad{A(n+1) - A(n) = a(n)}. sum: (RF, SE) -> Union(RF, FE) ++ sum(a(n), n) returns \spad{A} which ++ is the indefinite sum of \spad{a} with respect to ++ upward difference on \spad{n}, i.e. \spad{A(n+1) - A(n) = a(n)}. sum: (P, SegmentBinding P) -> RF ++ sum(f(n), n = a..b) returns \spad{f(a) + f(a+1) + ... f(b)}. sum: (RF, SegmentBinding RF) -> Union(RF, FE) ++ sum(f(n), n = a..b) returns \spad{f(a) + f(a+1) + ... f(b)}. Impl ==> add import RationalFunction R import GosperSummationMethod(IndexedExponents SE, SE, R, P, RF) innersum : (RF, SE) -> Union(RF, "failed") innerpolysum: (P, SE) -> RF sum(f:RF, s:SegmentBinding RF) == (indef := innersum(f, v := variable s)) case "failed" => summation(f::FE,map(#1::FE,s)$SegmentBindingFunctions2(RF,FE)) eval(indef::RF, v, 1 + hi segment s) - eval(indef::RF, v,lo segment s) sum(an:RF, n:SE) == (u := innersum(an, n)) case "failed" => summation(an::FE, n) u::RF sum(p:P, s:SegmentBinding P) == f := sum(p, v := variable s) eval(f, v, (1 + hi segment s)::RF) - eval(f,v,lo(segment s)::RF) innersum(an, n) == (r := retractIfCan(an)@Union(P, "failed")) case "failed" => an1 := eval(an, n, -1 + n::RF) (u := GospersMethod(an/an1, n, new$SE)) case "failed" => "failed" an1 * eval(u::RF, n, -1 + n::RF) sum(r::P, n) sum(p:P, n:SE) == rec := sum(p, n)$InnerPolySum(IndexedExponents SE, SE, R, P) rec.num / (rec.den :: P) @ \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}