Examples of Use of MgfunFr\303\251d\303\251ric ChyzakTo illustrate the course \342\200\234Holonomic summation and integration\342\200\235 during ISSFQFT2012 (July 2012).Set upThe package to be used is part of the Algolib library, which can be downloaded from http://algo.inria.fr/libraries/. After installing it, one needs to let Maple know about it:libname := "/home/chyzak", libname:Here we use our package Mgfun (F. Chyzak + contributions by Shaoshi Chen, Cyril Germa, Lucien Pech, and Ziming Li). An homologue package for Mathematica was written by Christoph Koutschan.with(Mgfun);infolevel[:-CreativeTelescoping] := 5:An example in P. Paule's coursef := (-1)^k / 2^k * binomial(n, k) * binomial(2*k, k);seq(i = add(eval(f, {n = i, k = j}), j = 0 .. i), i = 0..10);Creative telescoping in view of summation over LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= so as to describe a sum parametrized by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= is done like this:ct := creative_telescoping(f, n::shift, k::shift);This result can be interpreted in terms of the summand as:eval(ct[1][1], _F = unapply(_f(n, k), n)) = Delta[k] ( ct[1][2] );This cannot be evaluated for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1GLDYlUSJuRidGL0YyRjk= or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1GLDYlUSJuRidGL0YyRjUtRjY2LVEiK0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZmbkY1LUkjbW5HRiQ2JFEiMUYnRjlGOQ== or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1GLDYlUSJuRidGL0YyRjUtRjY2LVEiK0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZmbkY1LUkjbW5HRiQ2JFEiMkYnRjlGOQ== !So we can only sum until LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1GLDYlUSJuRidGL0YyRjUtRjY2LVEqJnVtaW51czA7RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORmZuRjUtSSNtbkdGJDYkUSIxRidGOUY5 and have to adjust the relation:ct[1][1] = (-n-1) * _f(n, n) + (n+2) * add(_f(n+2, n+i), i = 0..2) + eval(ct[1][2], k = n) - eval(ct[1][2], k = 0);eval(%, _f = unapply(f, n, k));Now, the adjusted inhomogeneous part turns out to be zero.normal(rhs(%), expanded);So, the left-hand component of the pair obtained by creative telescoping is a recurrence satisfied by the sum, and we can solve:sol := rsolve(ct[1][1] = 0, _F(n));Check a few values:seq(eval(sol/_F(0), n = 2 * p), p = 0..5);Home work in M. Schlosser's courseOur goal, (1.2) in Gasper and Schlosser's 2004 paper:f := (c - a * (a + t))^beta * (c - (a + 1) * (a + t))^(beta - 1) / (c - (a + t)^2)^(2 * beta) * t^beta * (1 - t)^(beta - 1);GAMMA(beta)^2/2/GAMMA(2*beta) = (c - (a + 1)^2) * Int(f, t = 0 .. 1);Encoding of the integrand by 1st-order linear functional equations (we just hide large coefficients under dots):dfinite_expr_to_sys(f, _f(beta::shift, t::diff)): collect(%, {_f, diff}, p -> `...`);Mgfun does most of it by itself (20 seconds):ct := creative_telescoping(f, beta::shift, t::diff):One typically doesn't want to see the result, as it can be large, but here it is for this example:ct;The inhomogeneous part of the creative-telescoping equation behaves nicely:Qf := eval(ct[1][2], _f(beta, t) = f):eval(Qf, t = 1);eval(Qf, t = 0);Thus, we directly get a homogeneous recurrence for the integral:rec := collect(ct[1][1], _F, factor);At this point, we can check that the known closed-form for the integral satisfies this recurrence:eval(rec, _F = unapply(GAMMA(beta)^2 / 2 / GAMMA(2 * beta) / (c - (a + 1)^2), beta));simplify(%, GAMMA);Therefore, checking enough initial conditions would prove the integral identity.On the other hand, we can fake that we don't know the closed-form evaluation of the integal and use a symbolic solver to get a closed form:sol := rsolve(rec, _F(beta));Last example in C. Raab's coursef := GegenbauerC(m, mu, x) * GegenbauerC(n, nu, x) * (1 - x^2)^(nu - 1/2);ct := creative_telescoping(f, [m::shift, n::shift], x::diff):collect(map2(op, 1, ct), _F, factor);These are the recurrences also found by Clemens Raab in his lecture.The sequences LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkobWZlbmNlZEdGJDYlLUYjNictSSVtc3ViR0YkNiUtSSNtaUdGJDYmUSJhRicvJSVib2xkR1EmZmFsc2VGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictRiM2Ji1GNDYmUSJuRidGN0Y6Rj1GN0Y6Rj0vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYnLyUlc2l6ZUdRIzEyRidGNy9GPlEnbm9ybWFsRidGN0ZORkhGS0Y3Rk4= and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkobWZlbmNlZEdGJDYlLUYjNictSSVtc3ViR0YkNiUtSSNtaUdGJDYmUSJiRicvJSVib2xkR1EmZmFsc2VGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictRiM2Ji1GNDYmUSJuRidGN0Y6Rj1GN0Y6Rj0vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYnLyUlc2l6ZUdRIzEyRidGNy9GPlEnbm9ybWFsRidGN0ZORkhGS0Y3Rk4= of Ap\303\251ry's proof of irrationality of \316\266(3)Ap\303\251ry's proof of irrationality of \316\266(3) relies on showing that the sequences with general terms LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiYkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn below satisfy the same second-order recurrence. This session shows how to get this recurrence.a[n] = Sum(binomial(n, k)^2 * binomial(n + k, k)^2, k = 0..n);b[n] = Sum(binomial(n, k)^2 * binomial(n + k, k)^2 * (Sum(1/m^3, m = 1..n) + Sum((-1)^(m+1) / (2 * m^3 * binomial(n, m) * binomial(n + m, m)), m = 1..k)), k = 0..n);with(Mgfun);Inner indefinite binomial sum: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEjRERGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSShtZmVuY2VkR0YkNiQtRiM2Ji1GLDYlUSJuRidGL0YyLUkjbW9HRiQ2LVEiLEYnL0YzUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0ZFLyUqc3ltbWV0cmljR0ZFLyUobGFyZ2VvcEdGRS8lLm1vdmFibGVsaW1pdHNHRkUvJSdhY2NlbnRHRkUvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GLDYlUSJrRidGL0YyRkFGQS1GPjYtUSJ+RidGQUZDL0ZHRkVGSEZKRkxGTkZQRlIvRlZGVC1GPjYtUSomY29sb25lcTtGJ0ZBRkNGaG5GSEZKRkxGTkZQL0ZTUSwwLjI3Nzc3NzhlbUYnL0ZWRl5vLUkrbXVuZGVyb3ZlckdGJDYnLUY+Ni1RJiZTdW07RidGQS9GRFEmdW5zZXRGJy9GR0Znby9GSUYxL0ZLRmdvL0ZNRjEvRk9GMS9GUUZnb0ZSL0ZWUSwwLjE2NjY2NjdlbUYnLUYjNiYtRiw2JVEibUYnRi9GMi1GPjYtUSI9RidGQUZDRmhuRkhGSkZMRk5GUEZdb0Zfby1JI21uR0YkNiRRIjFGJ0ZBRkEtRiM2J0ZYRi8vJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0YxRjJGUC8lLGFjY2VudHVuZGVyR0ZFLUkmbWZyYWNHRiQ2KC1GIzYoLUklbXN1cEdGJDYlLUY2NiQtRiM2KS1GPjYtUSomdW1pbnVzMDtGJ0ZBRkNGaG5GSEZKRkxGTkZQL0ZTUSwwLjIyMjIyMjJlbUYnL0ZWRmVyRmhwRi8vRl9xUStbMCwxNjAsODBdRidGYXEvJTZzZWxlY3Rpb24tcGxhY2Vob2xkZXJHRjFGMkZBLUYjNidGYnAtRj42LVEiK0YnRkFGQ0ZobkZIRkpGTEZORlBGZHJGZnJGaHBGL0YyLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0YvRmdyRmFxRmlyRjItRiM2LC1GaXA2JFEiMkYnRkFGZW4tRltyNiVGYnAtRiM2JS1GaXA2JFEiM0YnRkFGL0YyRmBzLUY2NiYtRmZxNigtRiM2KEY6Ri9GZ3JGYXFGaXJGMi1GIzYnRmJwRi9GXnFGYXFGMi8lLmxpbmV0aGlja25lc3NHRmJzLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRlt1LyUpYmV2ZWxsZWRHRkVGQS9JK21zZW1hbnRpY3NHRiRRKWJpbm9taWFsRidGYHUtRiw2I1EhRictRjY2Ji1GZnE2KC1GIzYqRjpGXXNGYnBGL0ZnckZhcUZpckYyRmV0Rmd0Rml0Rlx1Rl51RkFGYHVGYHVGL0ZecUZhcUYyL0ZodEZbcUZpdEZcdUZedUZBCreative telescoping on the summand LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiZEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYpLUYvNiVRIm5GJ0YyRjUtSSNtb0dGJDYtUSIsRicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYvNiVRImtGJ0YyRjVGPS1GLzYlUSJtRidGMkY1RjJGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRkE=:d_nkm := (-1)^(m+1) / (2 * m^3 * binomial(n, m) * binomial(n + m, m));ct_d := creative_telescoping(d_nkm, [n::shift, k::shift], m::shift);The meaning of this is given by the following two relations, where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5MTY7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYlLUYvNiVRIm1GJy9GM1EldHJ1ZUYnL0Y2USdpdGFsaWNGJ0Y9Rj8vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y1 means a forward finite difference operator w.r.t. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEibUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=:eval(map(e -> e[1] = Delta[m](e[2]), eval(ct_d, _F = unapply(_f(n, k, m), n, k))), _f = d);(in other words:)eval(%, {d(n + 1, k, m) = Delta[n](d(n, k, m)) + d(n, k, m), d(n, k + 1, m) = Delta[k](d(n, k, m)) + d(n, k, m)});Deal with left-hand relationCheck it:eval(ct_d[1][1], _F = unapply(d_nkm, n, k));eval(ct_d[1][2], _f(n, k, m) = d_nkm);normal(%% - (eval(%, m = m + 1) - %), expanded);Exploit left-hand relation to get relation on sum:eval(ct_d[1][2], _f(n, k, m) = d_nkm);eval(%, m = k + 1) - eval(%, m = 1);inhom1 := map(normal@expand, %, expanded);dfinite_expr_to_sys(inhom1, u(n::shift, k::shift));sys1 := eval(%, u = unapply(DD(n + 1, k) - DD(n, k), n, k));sys1;Deal with right-hand relationinhom2 := factor(normal(expand(eval(d_nkm, m = k + 1)), expanded));dfinite_expr_to_sys(inhom2, u(n::shift, k::shift));sys2 := eval(%, u = unapply(DD(n, k + 1) - DD(n, k), n, k));sys2;Finally, gather both systems:sys_for_DD := collect(sys1 union sys2, DD, expand);sys_for_DD;Second inner indefinite binomial sum: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEkaXozRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkobWZlbmNlZEdGJDYkLUYjNiYtRiw2JVEibkYnRi9GMi1JI21vR0YkNi1RIixGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGRS8lKnN5bW1ldHJpY0dGRS8lKGxhcmdlb3BHRkUvJS5tb3ZhYmxlbGltaXRzR0ZFLyUnYWNjZW50R0ZFLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRiw2JVEia0YnRi9GMkZBRkEtRj42LVEifkYnRkFGQy9GR0ZFRkhGSkZMRk5GUEZSL0ZWRlQtRj42LVEqJmNvbG9uZXE7RidGQUZDRmhuRkhGSkZMRk5GUC9GU1EsMC4yNzc3Nzc4ZW1GJy9GVkZeby1JK211bmRlcm92ZXJHRiQ2Jy1GPjYtUSYmU3VtO0YnRkEvRkRRJnVuc2V0RicvRkdGZ28vRklGMS9GS0Znby9GTUYxL0ZPRjEvRlFGZ29GUi9GVlEsMC4xNjY2NjY3ZW1GJy1GIzYmLUYsNiVRIm1GJ0YvRjItRj42LVEiPUYnRkFGQ0ZobkZIRkpGTEZORlBGXW9GX28tSSNtbkdGJDYkUSIxRidGQUZBLUYjNidGOkYvLyUrZm9yZWdyb3VuZEdRLFsyMDAsMCwyMDBdRicvJSxwbGFjZWhvbGRlckdGMUYyRlAvJSxhY2NlbnR1bmRlckdGRS1JJm1mcmFjR0YkNigtRiM2KEZocEYvL0ZfcVErWzAsMTYwLDgwXUYnRmFxLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR0YxRjItRiM2Jy1JJW1zdXBHRiQ2JUZicC1GIzYlLUZpcDYkUSIzRidGQUYvRjIvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRi9GXnFGYXFGMi8lLmxpbmV0aGlja25lc3NHRltxLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl9zLyUpYmV2ZWxsZWRHRkVGQQ==We use a similar process, but simplified because of the simpler summand: some of the operations can be done \342\200\234by hand.\342\200\235dfinite_expr_to_sys(1/(n + 1)^3, u(n::shift, k::shift));sys_for_iz3 := eval(%, u = unapply(iz3(n + 1, k) - iz3(n, k), n, k)) union {iz3(n, k + 1) - iz3(n, k)};sys_for_iz3;cofactor := binomial(n, k)^2 * binomial(n + k, k)^2;System for the summand of the outer sum: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYvLUkjbWlHRiQ2JlEpY29mYWN0b3JGJy8lJWJvbGRHUSZmYWxzZUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RJyZzZG90O0YnRi8vRjZRJ25vcm1hbEYnLyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTi1JKG1mZW5jZWRHRiQ2JS1GIzYrLUYsNiZRI0RERidGL0YyRjUtRlI2JS1GIzYpLUYsNiZRIm5GJ0YvRjJGNS1GOTYuUSIsRidGL0Y8Rj4vRkFGNEZCRkRGRkZIRkpGTC9GUFEsMC4zMzMzMzMzZW1GJy1GLDYmUSJrRidGL0YyRjUvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnRi9GPEYvRjwtRjk2LlEiK0YnRi9GPEY+RkBGQkZERkZGSEZKL0ZNUSwwLjIyMjIyMjJlbUYnL0ZQRl1wLUYsNiZRJGl6M0YnRi9GMkY1RllGY29GZm9GL0Y8Ri9GPC1GOTYuUSJ+RidGL0Y8Rj5GQEZCRkRGRkZIRkpGTEZPLUY5Ni5RIj1GJ0YvRjxGPkZARkJGREZGRkhGSi9GTVEsMC4yNzc3Nzc4ZW1GJy9GUEZpcEZicC1JJW1zdXBHRiQ2JS1GUjYoLUkmbWZyYWNHRiQ2KC1GIzYqRmduRmNvRi9GMi8lK2ZvcmVncm91bmRHUStbMCwxNjAsODBdRicvJSxwbGFjZWhvbGRlckdGNC8lNnNlbGVjdGlvbi1wbGFjZWhvbGRlckdGNEY1LUYjNilGYG9GY29GL0YyL0ZmcVEsWzIwMCwwLDIwMF1GJ0ZocUY1LyUubGluZXRoaWNrbmVzc0dRIjBGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zlci8lKWJldmVsbGVkR0YxRmNvRi9GPC9JK21zZW1hbnRpY3NHRiRRKWJpbm9taWFsRidGanItRiM2Ji1JI21uR0YkNiVRIjJGJ0YvRjxGL0YyRjUvJTFzdXBlcnNjcmlwdHNoaWZ0R0Zici1GXHE2JS1GUjYoLUZhcTYoLUYjNixGZ25GaW9GYG9GY29GL0YyRmVxRmhxRmpxRjVGXHJGYHJGY3JGZnJGaHJGY29GL0Y8RmpyRmpyRl1zRmNzLUZSNiUtRiM2Ky1JK211bmRlcm92ZXJHRiQ2Jy1GOTYxUSYmU3VtO0YnRmNvRmZvRi8vJTBmb250X3N0eWxlX25hbWVHUSgyRH5NYXRoRidGPC9GP1EmdW5zZXRGJy9GQUZbdS9GQ0Y0L0ZFRlt1L0ZHRjQvRklGNC9GS0ZbdUZML0ZQUSwwLjE2NjY2NjdlbUYnLUYjNiotRiw2KVEibUYnRmNvRmZvRi9GMkZndEY1LUY5NjFGZ3BGY29GZm9GL0ZndEY8Rj5GQEZCRkRGRkZIRkpGaHBGanAtRmBzNihRIjFGJ0Zjb0Zmb0YvRmd0RjxGY29GZm9GL0ZndEY8LUYjNigtRiw2KUZpbkZjb0Zmb0YvRjJGZ3RGNUZjb0Zmb0YvRmd0RjxGSi8lLGFjY2VudHVuZGVyR0YxLUZhcTYoLUYjNihGW3ZGY29GZm9GL0ZndEY8LUYjNigtRlxxNiVGZnUtRiM2KC1GYHM2KFEiM0YnRmNvRmZvRi9GZ3RGPEZjb0Zmb0YvRmd0RjxGY3NGY29GZm9GL0ZndEY8L0ZhckZddkZjckZmckZockZpby1GYnQ2J0ZkdEZkdS1GIzYoLUYsNilGYm9GY29GZm9GL0YyRmd0RjVGY29GZm9GL0ZndEY8RkpGYnYtRmFxNigtRiM2KC1GXHE2JS1GUjYoLUYjNiktRjk2MVEqJnVtaW51czA7RidGY29GZm9GL0ZndEY8Rj5GQEZCRkRGRkZIRkpGXHBGXnBGW3ZGY29GZm9GL0ZndEY8RmNvRmZvRi9GZ3RGPC1GIzYqRmZ1LUY5NjFGW3BGY29GZm9GL0ZndEY8Rj5GQEZCRkRGRkZIRkpGXHBGXnBGW3ZGY29GZm9GL0ZndEY8RmNzRmNvRmZvRi9GZ3RGPC1GIzYtLUZgczYoRmJzRmNvRmZvRi9GZ3RGPC1GOTYxRmRwRmNvRmZvRi9GZ3RGPEY+RkBGQkZERkZGSEZKRkxGT0Zqdi1GUjYqLUZhcTYoRl52LUYjNihGZnVGY29GZm9GL0ZndEY8RmByRmNyRmZyRmhyRmNvRmZvRi9GZ3RGPEZqckZqci1GLDYjUSFGJy1GUjYqLUZhcTYoLUYjNipGYHZGZ3hGZnVGY29GZm9GL0ZndEY8RmN5RmByRmNyRmZyRmhyRmNvRmZvRi9GZ3RGPEZqckZqckZjb0Zmb0YvRmd0RjxGYXdGY3JGZnJGaHJGY29GZm9GL0Y8Ri9GPEZjb0Zmb0YvRjw=Now, it would be great just to put the two systems side by side in a call, but the routine requires the same name for both systems; so the following makes an error:`sys+sys`(sys_for_DD, sys_for_iz3);Therefore we help Mgfun to add the two sums:sys_for_2_indef_sums := `sys+sys`(eval(sys_for_DD, DD = u), eval(sys_for_iz3, iz3 = u));Next, we prepare a system for the left-hand factor of the product:sys_for_cofactor := dfinite_expr_to_sys(cofactor, u(n::shift, k::shift));And we do the multiplication:sys_for_summand := `sys*sys`(sys_for_cofactor, sys_for_2_indef_sums);sys_for_summand;Fourth-order recurrence for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYlLUkjbWlHRiQ2JlEiYkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiYtRi82JlEibkYnRjJGNUY4RjJGNUY4LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnRjIvRjlRJ25vcm1hbEYnHere, we use the command creative_telescoping by providing a system (40 seconds):ct_for_double_sum := creative_telescoping(LFSol(sys_for_summand), n::shift, k::shift):Read off the recurrence for the sum:rec_for_double_sum := collect(ct_for_double_sum[1][1], _F, factor);rec_for_double_sum;Second-order recurrence for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYlLUkjbWlHRiQ2JlEiYUYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiYtRi82JlEibkYnRjJGNUY4RjJGNUY4LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnRjIvRjlRJ25vcm1hbEYnct_for_single_sum := creative_telescoping(cofactor, n::shift, k::shift);rec_for_single_sum := collect(ct_for_single_sum[1][1], _F, factor);rec_for_single_sum;Prove that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYlLUkjbWlHRiQ2JlEiYkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiYtRi82JlEibkYnRjJGNUY4RjJGNUY4LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnRjIvRjlRJ25vcm1hbEYn also satisfies the second-order recurrenceObserve that a solution of the 4th-order recurrence is fully determined by its first four values, as its leading coefficient never vanishes on natural integers:factor(coeff(rec_for_double_sum, _F(n + 4)));Therefore, we can introduce the sequence starting by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiYkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMEYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdGPS1JI21vR0YkNi1RIixGJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRkgvJSpzeW1tZXRyaWNHRkgvJShsYXJnZW9wR0ZILyUubW92YWJsZWxpbWl0c0dGSC8lJ2FjY2VudEdGSC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVGLi1GIzYlLUY7NiRRIjFGJ0Y+RjJGNUZARkItRiw2JUYuLUYjNiUtRjs2JFEiMkYnRj5GMkY1RkBGQi1GLDYlRi4tRiM2JS1GOzYkUSIzRidGPkYyRjVGQEY+ and prolonged by unrolling the second-order recurrence. If we show that it satisfies the fourth-order recurrence, it will be exactly the sequence LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYkLUYjNiQtSSVtc3ViR0YkNiUtSSNtaUdGJDYlUSJiRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiUtRjQ2JVEibkYnRjdGOkY3RjovJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9GO1Enbm9ybWFsRidGRUZF. This will show that the latter also satisfies the second-order recurrence._F(n + 2) = solve(rec_for_single_sum, _F(n + 2));collect(eval(%, n = n - 2), _F, factor);simpl := subs(a = %, n -> a);eval(rec_for_double_sum, simpl(n + 4));eval(%, simpl(n + 3));normal(eval(%, simpl(n + 2)));Computing LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYzLUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2L1ErJkludGVncmFsO0YnLyUnZmFtaWx5R1EwVGltZXN+TmV3flJvbWFuRicvJSVib2xkR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJnVuc2V0RicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdRJXRydWVGJy8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRkIvJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTS1GIzYpLUkjbW5HRiQ2JVEiMEYnRjVGOEYyRjUvJSdpdGFsaWNHRkIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0ZCL0Y5USdpdGFsaWNGJy1GIzYpLUYvNi5RKCZpbmZpbjtGJ0Y1RjgvRjxGNy9GP0Y3L0ZBRjcvRkRGNy9GRkY3L0ZIRjcvRkpGN0ZLRk5GMkY1RlZGWEZlbkZnbi8lMXN1cGVyc2NyaXB0c2hpZnRHRlUvJS9zdWJzY3JpcHRzaGlmdEdGVS1JJW1zdXBHRiQ2JS1JI21pR0YkNilRInhGJ0YyLyUlc2l6ZUdRIzEyRidGNUZWLyUwZm9udF9zdHlsZV9uYW1lR1EoMkR+TWF0aEYnRmduLUYjNiktRl1wNilRImFGJ0YyRmBwRjVGVkZjcEZnbkYyRmBwRjVGVkZjcEZnbkZlby1GLzYtUTEmSW52aXNpYmxlVGltZXM7RidGOEZeb0Zfb0Zgb0Zhb0Zib0Zjb0Zkb0ZLRk4tRmpvNiUtRi82MVEvJkV4cG9uZW50aWFsRTtGJ0YyRmBwRjVGY3BGOEZeb0Zfb0Zgb0Zhb0Zib0Zjb0Zkb0ZLL0ZPUSwwLjExMTExMTFlbUYnLUYjNistRi82MVEqJnVtaW51czA7RidGMkZgcEY1RmNwRjhGXm9GX29GYG9GYW9GYm9GY29GZG8vRkxRLDAuMjIyMjIyMmVtRicvRk9GW3ItRl1wNilRI2J4RidGMkZgcEY1RlZGY3BGZ25GMkZgcEY1L0ZZUSpbMCwwLDI1NV1GJy8lKXJlYWRvbmx5R0ZCL0ZkcFEqMkR+T3V0cHV0RidGOEZlb0ZbcS1JJW1zdWJHRiQ2JS1GXXA2KVEiS0YnRjJGYHBGNUZWRmNwRmduLUYjNiktRl1wNilRIm1GJ0YyRmBwRjVGVkZjcEZnbkYyRmBwRjVGVkZjcEZnbkZnby1JKG1mZW5jZWRHRiQ2KC1GIzYqLUZdcDYpUSNjeEYnRjJGYHBGNUZWRmNwRmduRjJGYHBGNUZgckZickZkckY4RjJGYHBGNUZjcEY4RltxLUZncjYlRmlyLUYjNiktRl1wNilRIm5GJ0YyRmBwRjVGVkZjcEZnbkYyRmBwRjVGVkZjcEZnbkZnby1GYnM2KC1GIzYqLUZdcDYpUSNkeEYnRjJGYHBGNUZWRmNwRmduRjJGYHBGNUZgckZickZkckY4RjJGYHBGNUZjcEY4LUYvNi5RIn5GJ0Y1RjhGXm9GX29GYG9GYW9GYm9GY29GZG9GS0ZOLUYvNi9RMCZEaWZmZXJlbnRpYWxEO0YnRjJGNUY4RjtGPi9GQUY9RkMvRkZGPUZHRklGS0ZOLUZdcDYmRl9wRjVGVkZnbkYyRmBwRjVGOA==_F(a) = Int(_f(a, x), x = 0..infinity);for:f := x^a*exp(-b*x)*BesselK(m,c*x)*BesselK(n,d*x);under the knowledge:(close to 0)MultiSeries:-series(BesselK(nu, x), x = 0, 3) assuming nu > 5;(close to infinity)MultiSeries:-asympt(BesselK(nu, x), x, 1);The following will be used to display results more nicely.pnice := proc(expr, $) local l := [a, b, c, d, m, n]; lcoeff(expr, l) * convert(map((v, e) -> v^degree(e, v), l, expr), `*`) + `...` end proc:rnice := proc(expr, $) local de := denom(expr); pnice(numer(expr)) / `if`(de = 1, 1, pnice(de)) end proc:nice := proc(expr, $) collect(expr, {_F, _f, diff}, rnice) end proc:\342\200\234Creative telescoping\342\200\235 (15 seconds):res := creative_telescoping(f, [a::shift], x::diff):P_applied_to_F := res[1][1]:cPF := collect(P_applied_to_F, _F, nice);Q_applied_to_f := res[1][2]:cQf := collect(Q_applied_to_f, _f, nice);Meaning:eval(cPF, _F = unapply(_f(a, x), a)) = Diff( cQf, x );(now in terms of example)eval(eval(P_applied_to_F, _F = unapply(_f(a, x), a)) = diff( Q_applied_to_f, x ), _f = unapply(f, a, x)):simplify( lhs(%) / rhs(%) );After integration:cPF = Limit(cQf, x = infinity) - Limit(cQf, x = 0);Inhomogeneous term is exponentiallyl small at infinity:MultiSeries:-asympt(f, x, 2) assuming b > 0 and c > 0 and d > 0;Inhomogeneous term is 0 at 0 for large enough LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= (a few seconds):MultiSeries:-series(eval(Q_applied_to_f, _f = unapply(f, a, x)) / x^(a-m-n), x = 0, 5) assuming m > 5 and n > 5:x^(a-m-n) * collect(convert(%, polynom), x, e -> `...`);Internally, what has been obtained first is a system:map(nice, dfinite_expr_to_sys(f, _f(a::shift, x::diff)));ODE for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzY0LUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2L1ErJkludGVncmFsO0YnLyUnZmFtaWx5R1EwVGltZXN+TmV3flJvbWFuRicvJSVib2xkR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJnVuc2V0RicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdRJXRydWVGJy8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRkIvJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTS1GIzYpLUkjbW5HRiQ2JVEiMEYnRjVGOEYyRjUvJSdpdGFsaWNHRkIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0ZCL0Y5USdpdGFsaWNGJy1GIzYpLUYvNi9RKCZpbmZpbjtGJ0YyRjVGOC9GPEY3L0Y/RjcvRkFGNy9GREY3L0ZGRjcvRkhGNy9GSkY3RktGTkYyRjVGVkZYRmVuRmduLyUxc3VwZXJzY3JpcHRzaGlmdEdGVS8lL3N1YnNjcmlwdHNoaWZ0R0ZVLUYsNidGLi1GIzYmRlJGNUZWRmduRmluRmVvRmdvLUklbXN1YkdGJDYlLUkjbWlHRiQ2J1EiSkYnRjJGNUZWRmduLUYjNictRlM2JlEiMUYnRjJGNUY4RjJGNUZWRmduRmdvLUkobWZlbmNlZEdGJDYmLUYjNictRmFwNidRInhGJ0YyRjVGVkZnbkYyLyUlc2l6ZUdRIzEyRidGNUY4RjJGNUY4Rl1wLUZqcDYmLUYjNictRmFwNidRInlGJ0YyRjVGVkZnbkYyRmFxRjVGOEYyRjVGOC1GXnA2JUZgcC1GIzYnLUZTNiZRIjJGJ0YyRjVGOEYyRjVGVkZnbkZnby1GanA2Ji1GIzYoLUZhcDYnUSJjRidGMkY1RlZGZ24tSSZtc3FydEdGJDYjLUYjNiYtRmFwNidRI3h5RidGMkY1RlZGZ25GMkY1RjhGMkZhcUY1RjhGMkY1RjgtRi82L1EifkYnRjJGNUY4Rl5vRl9vRmBvRmFvRmJvRmNvRmRvRktGTi1GLzYvUTAmRGlmZmVyZW50aWFsRDtGJ0YyRjVGOEY7Rj4vRkFGPUZDL0ZGRj1GR0ZJRktGTi1GYXA2JkZqcUY1RlZGZ24tRi82LkZjc0Y1RjhGXm9GX29GYG9GYW9GYm9GY29GZG9GS0ZORmRzLUZhcDYmRmBxRjVGVkZnbkYyRmFxRjVGOA==Interested in:_F2(c) = Int(Int(_f2(c, x, y), y = 0..infinity), x = 0..infinity);f2 := BesselJ(1, x) * BesselJ(1, y) * BesselJ(2, c * sqrt(x*y)) / exp(x + y);res2 := creative_telescoping(f2, c::diff, [x::diff, y::diff]):The local identity obtained by creative telescoping is the following:subs(_F(c) = _f(c, x, y), res2[1][1]) = Diff(res2[1][2], x) + Diff(res2[1][3], y);It can be integrated over LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= aan LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=, which leads to boundary terms that are 0. Therefore, the first component of the triple is a differential equation satisfied by the parametrized double integral:res2[1][1] = 0;