## Friday, July 2, 2010

### Symbolic infinite series

I've written a little module for symbolically evaluating hypergeometric infinite series. Technically, a series is hypergeometric if the ratio between successive terms is a rational function of the index k; more intuitively, most "nice-looking" series (such as power series of most special functions) involving powers, factorials, gamma functions, binomial coefficients (etc.) are hypergeometric.

The code currently uses SymPy, but should be possible to convert to Sage's symbolics (it mostly requires basic pattern matching, and evaluation of standard special functions). Perhaps "rewrite" is a better word than "convert" because the current code is a messy hack. It's buggy, and only handles special cases in places where a general approach could be used (such as for rational functions).

The basic idea behind the code is to perform the evaluation in two steps. Firstly, the given series is converted into canonical form. I've used the standard generalized hypergeometric series pFq. This is not always an optimal representation, but usually does the job.

In the second step, parameter transformations and table lookups are used to reduce the pFq series into elementary (or simpler) functions when possible. If this step fails, the pFq function is returned, giving a closed-form solution that can be manipulated or evaluated numerically.

The second step is the most difficult, because the transformation algorithms that need to be used for hypergeometric functions are very complicated in general (and there are all kinds of special cases). So far I've only implemented some simple transformations. For the final lookup stage, it would be straightforward to also implement expansion into Bessel functions, incomplete gamma functions, etc.

Below is a sample of some working test cases, with numerical evaluation of both the evaluated result and the original series as verification. Some of the sums don't quite simplify fully, i.e. there is an unevaluated pFq function left. The formulas were mostly generated with SymPy's latex function, so they are not always as pretty as could be.

>>> hypsum(1 / fac(k), k)E>>> _.evalf()2.71828182845905>>> mpmath.nsum(lambda k: 1 / mpmath.fac(k), [0,mpmath.inf])2.71828182845905

>>> hypsum(k**6 / fac(k), k)203*E>>> _.evalf()551.811211177186>>> mpmath.nsum(lambda k: k**6 / mpmath.fac(k), [0,mpmath.inf])551.811211177186

>>> hypsum(4 * d * k * z**(k+2) / fac(k), k)4*d*z**3*exp(z)>>> _.evalf(subs={d:'1/2',z:'1/4'})0.0401257942714919>>> mpmath.nsum(lambda k: 4 * 0.5 * k * 0.25**(k+2) / mpmath.fac(k), [0,mpmath.inf])0.0401257942714919

>>> hypsum(1 / rf(2,k), k)-1 + E>>> _.evalf()1.71828182845905>>> mpmath.nsum(lambda k: 1 / mpmath.rf(2,k), [0,mpmath.inf])1.71828182845905

>>> simplify(hypsum(k**3 * z**(k) / fac(k), k))z*exp(z) + z**3*exp(z) + 3*z**2*exp(z)>>> _.evalf(subs={d:'1/2',z:'1/4'})0.581824016936633>>> mpmath.nsum(lambda k: k**3 * 0.25**(k) / mpmath.fac(k), [0,mpmath.inf])0.581824016936633

>>> simplify(hypsum(1 / binomial(2*k,k), k))4/3 + 2*pi*3**(1/2)/27  >>> _.evalf()1.73639985871872>>> mpmath.nsum(lambda k: 1 / mpmath.binomial(2*k,k), [0,mpmath.inf])1.73639985871872

>>> hypsum(1 / binomial(3*k,k), k)2*pi*3**(1/2)*3F2([1/2, 1, 1], [1/3, 2/3], 4/27)/(3*gamma(1/3)*gamma(2/3))>>> _.evalf()1.41432204432182>>> mpmath.nsum(lambda k: 1 / mpmath.binomial(3*k,k), [0,mpmath.inf])1.41432204432182

>>> hypsum(k**4 / binomial(2*k,k), k)32/3 + 238*pi*3**(1/2)/243>>> _.evalf()15.9961018356512>>> mpmath.nsum(lambda k: k**4 / mpmath.binomial(2*k,k), [0,mpmath.inf])15.9961018356512

>>> hypsum(fac(k+1)*fac(k+2)/fac(k+3)/fac(k+4), k)3F2([1, 2, 3], [4, 5], 1)/72>>> _.evalf()0.0217325998184402>>> mpmath.nsum(lambda k: mpmath.fac(k+1)*mpmath.fac(k+2)/mpmath.fac(k+3)/mpmath.fac(k+4), [0,mpmath.inf])0.0217325998184402

>>> hypsum(1/((3+k)**2*(8+6*k+k**2)), k)3F2([1, 2, 3], [4, 5], 1)/72>>> _.evalf()0.0217325998184402>>> mpmath.nsum(lambda k: 1/((3+k)**2*(8+6*k+k**2)), [0,mpmath.inf])0.0217325998184402

>>> hypsum(fac(2*k+1)/(fac(2*k+2)*(k+1)), k)pi**2/12>>> _.evalf()0.822467033424113>>> mpmath.nsum(lambda k: mpmath.fac(2*k+1)/(mpmath.fac(2*k+2)*(k+1)), [0,mpmath.inf])0.822467033424113

>>> hypsum(fac(2*k+1)/(fac(3*k+2)), k)2*pi*3**(1/2)*2F2([1, 3/2], [4/3, 5/3], 4/27)/(27*gamma(4/3)*gamma(5/3))>>> _.evalf()0.553106730441975>>> mpmath.nsum(lambda k: mpmath.fac(2*k+1)/(mpmath.fac(3*k+2)), [0,mpmath.inf])0.553106730441975

>>> hypsum(4**k*fac(k+R32)*fac(k+R52)/(fac(2*k+6)), k)3*pi/256>>> _.evalf()0.0368155389092554>>> mpmath.nsum(lambda k: 4**k*mpmath.fac(k+1.5)*mpmath.fac(k+2.5)/(mpmath.fac(2*k+6)), [0,mpmath.inf], method='e')0.0368155389092358

>>> hypsum((-1)**k * z**(2*k+1) / fac(k) / (2*k+1), k)pi**(1/2)*erf(z)/2>>> _.evalf(subs={d:'1/2',z:'1/4'})0.244887887180256>>> mpmath.nsum(lambda k: (-1)**k * 0.25**(2*k+1) / mpmath.fac(k) / (2*k+1), [0,mpmath.inf])0.244887887180256

>>> hypsum(z**(2*k+1) / (2*k+1), k)atanh(z)>>> _.evalf(subs={d:'1/2',z:'1/4'})0.255412811882995>>> mpmath.nsum(lambda k: 0.25**(2*k+1) / (2*k+1), [0,mpmath.inf])0.255412811882995

>>> hypsum(z**k / (2*k+1), k)atanh(z**(1/2))/z**(1/2)>>> _.evalf(subs={d:'1/2',z:'1/4'})1.09861228866811>>> mpmath.nsum(lambda k: 0.25**k / (2*k+1), [0,mpmath.inf])1.09861228866811

>>> hypsum((-z)**k / (k+1), k)log(1 + z)/z>>> _.evalf(subs={d:'1/2',z:'1/4'})0.892574205256839>>> mpmath.nsum(lambda k: (-0.25)**k / (k+1), [0,mpmath.inf])0.892574205256839

>>> hypsum(fac(k-R12)/((1+2*k)*fac(k))*z**(2*k), k)pi**(1/2)*asin(z)/z>>> _.evalf(subs={d:'1/2',z:'1/4'})1.79145636509746>>> mpmath.nsum(lambda k: mpmath.fac(k-0.5)/((1+2*k)*mpmath.fac(k))*0.25**(2*k), [0,mpmath.inf])1.79145636509746

>>> hypsum(z**k / fac(k+2), k)-(1 + z - exp(z))/z**2>>> _.evalf(subs={d:'1/2',z:'1/4'})0.544406667003864>>> mpmath.nsum(lambda k: 0.25**k / mpmath.fac(k+2), [0,mpmath.inf])0.544406667003864

>>> hypsum((3*z**2-5)**(-4*k+3), k)-(5 - 3*z**2)**3/(1 - 1/(5 - 3*z**2)**4)>>> _.evalf(subs={d:'1/2',z:'1/4'})-111.666432272586>>> mpmath.nsum(lambda k: (3*0.25**2-5)**(-4*k+3), [0,mpmath.inf])-111.666432272586

>>> hypsum(z**(2*k+1) / fac(2*k+1), k)sinh(z)>>> _.evalf(subs={d:'1/2',z:'1/4'})0.252612316808168>>> mpmath.nsum(lambda k: 0.25**(2*k+1) / mpmath.fac(2*k+1), [0,mpmath.inf])0.252612316808168

>>> hypsum((-1)**k * z**(2*k+1) / fac(2*k+1), k)sin(z)>>> _.evalf(subs={d:'1/2',z:'1/4'})0.247403959254523>>> mpmath.nsum(lambda k: (-1)**k * 0.25**(2*k+1) / mpmath.fac(2*k+1), [0,mpmath.inf])0.247403959254523

>>> hypsum((-1)**k * d**k * (1-z)**(2*k+1) / fac(2*k), k)(1 - z)*cos(d**(1/2)*(1 - z))>>> _.evalf(subs={d:'1/2',z:'1/4'})0.646980115568008>>> mpmath.nsum(lambda k: (-1)**k * 0.5**k * (1-0.25)**(2*k+1) / mpmath.fac(2*k), [0,mpmath.inf])0.646980115568008

>>> hypsum(k * z**(2*k) / fac(2*k+1), k)cosh(z)/2 - sinh(z)/(2*z)>>> _.evalf(subs={d:'1/2',z:'1/4'})0.0104819163234500>>> mpmath.nsum(lambda k: k * 0.25**(2*k) / mpmath.fac(2*k+1), [0,mpmath.inf])0.01048191632345

>>> hypsum(gamma(k-R12)**2/(fac(k)**2), k)16>>> _.evalf()16.0000000000000>>> mpmath.nsum(lambda k: mpmath.gamma(k-0.5)**2/(mpmath.fac(k)**2), [0,mpmath.inf])16.0

>>> hypsum(1/(k+1)**2, k)pi**2/6>>> _.evalf()1.64493406684823>>> mpmath.nsum(lambda k: 1/(k+1)**2, [0,mpmath.inf])1.64493406684823

>>> hypsum(k/(k+1)**3, k)-zeta(3) + pi**2/6>>> _.evalf()0.442877163688632>>> mpmath.nsum(lambda k: k/(k+1)**3, [0,mpmath.inf])0.442877163688632

>>> simplify(hypsum((3+k)/(3+3*k)*z**k/fac(k), k))-(2 - 2*exp(z) - z*exp(z))/(3*z)>>> _.evalf(subs={d:'1/2',z:'1/4'})1.18540958339656>>> mpmath.nsum(lambda k: (3+k)/(3+3*k)*0.25**k/mpmath.fac(k), [0,mpmath.inf])1.18540958339656

>>> hypsum(k**2 * z**(2*k+1) / fac(2*k+9), k)z**3*(-1F2([2], [6, 13/2], z**2/4) + 2*1F2([3], [6, 13/2], z**2/4))/39916800>>> _.evalf(subs={d:'1/2',z:'1/4'})3.92066920165299e-10>>> mpmath.nsum(lambda k: k**2 * 0.25**(2*k+1) / mpmath.fac(2*k+9), [0,mpmath.inf])3.92066920165299e-10

>>> simplify(hypsum((k+1)*(k+2)*(k+3)*(1+2*k)*z**k, k))(6 + 42*z)/(1 - 5*z + 10*z**2 - 10*z**3 + 5*z**4 - z**5)>>> _.evalf(subs={d:'1/2',z:'1/4'})69.5308641975309>>> mpmath.nsum(lambda k: (k+1)*(k+2)*(k+3)*(1+2*k)*0.25**k, [0,mpmath.inf])69.5308641975309

>>> hypsum(k**3 * (-z)**k / (k+1), k)-z*(2*(1 - z)/(1 + z)**3 - ((2 + 2*z)/(1 + z) - (2 + 2*z)*log(1 + z)/z)/(z*(1 + z)) - 2/(1 + z)**2)/2>>> _.evalf(subs={d:'1/2',z:'1/4'})-0.0285742052568390>>> mpmath.nsum(lambda k: k**3 * (-0.25)**k / (k+1), [0,mpmath.inf])-0.028574205256839

>>> hypsum(E**k / fac(k), k)exp(E)>>> _.evalf()15.1542622414793>>> mpmath.nsum(lambda k: mpmath.e**k / mpmath.fac(k), [0,mpmath.inf])15.1542622414793

>>> hypsum(cos(k) / fac(k), k)exp(exp(I))/2 + exp(exp(-I))/2>>> _.evalf()1.14383564379164 + .0e-20*I>>> mpmath.nsum(lambda k: mpmath.cos(k) / mpmath.fac(k), [0,mpmath.inf])1.14383564379164

>>> hypsum(cos(k) / (2*k+1), k)atanh(exp(-I/2))*exp(I/2)/2 + atanh(exp(I/2))*exp(-I/2)/2>>> _.evalf()0.975556628913311>>> mpmath.nsum(lambda k: mpmath.cos(k) / (2*k+1), [0,mpmath.inf])0.975556628913311

>>> simplify(hypsum(cos(k) * k / fac(k), k))exp(I + exp(I))/2 + exp(-I + exp(-I))/2>>> _.evalf()-0.458967373729452 + .0e-20*I>>> mpmath.nsum(lambda k: mpmath.cos(k) * k / mpmath.fac(k), [0,mpmath.inf])-0.458967373729452

>>> hypsum((-1)**k / (2*k+1), k)pi/4>>> _.evalf()0.785398163397448>>> mpmath.nsum(lambda k: (-1)**k / (2*k+1), [0,mpmath.inf])0.785398163397448

>>> hypsum((-1)**k / (2*k+6), k)-1/4 + log(2)/2>>> _.evalf()0.0965735902799727>>> mpmath.nsum(lambda k: (-1)**k / (2*k+6), [0,mpmath.inf])0.0965735902799727

>>> hypsum((-1)**k / (2*k+1)**2, k)Catalan>>> _.evalf()0.915965594177219>>> mpmath.nsum(lambda k: (-1)**k / (2*k+1)**2, [0,mpmath.inf])0.915965594177219

>>> hypsum((-1)**k / (2*k+2)**2, k)pi**2/48>>> _.evalf()0.205616758356028>>> mpmath.nsum(lambda k: (-1)**k / (2*k+2)**2, [0,mpmath.inf])0.205616758356028

>>> hypsum((-1)**k * (k+1) / (2*k+3)**2, k)3F2([3/2, 3/2, 2], [5/2, 5/2], -1)/9>>> _.evalf()0.0652837153898853>>> mpmath.nsum(lambda k: (-1)**k * (k+1) / (2*k+3)**2, [0,mpmath.inf])0.0652837153898854

>>> hypsum(1/(k**2+2*k+4), k)I*3**(1/2)*(-polygamma(0, 1 + I*3**(1/2)) + polygamma(0, 1 - I*3**(1/2)))/6>>> _.evalf()0.740267076581851>>> mpmath.nsum(lambda k: 1/(k**2+2*k+4), [0,mpmath.inf])0.740267076581851

This work was possible thanks to the support of NSF grant DMS-0757627, which is gratefully acknowledged.

Arkapravo said...

Very nice !

Aaron Meurer said...

Cool. So do you plan on backporting the symbolic solver to SymPy?

Fredrik Johansson said...

Aaron: it's possible. The basic functionality anyway, it's quite simple; if I don't do it, probably someone else can do it quite easily.