🚀 go-pugleaf

Here is a quick answer to a question of "fruskey@sol.UVic.CA (Frank  Ruskey)"
in article <1992Sep24.180403.1548@sol.UVic.CA>.
It was prepared by Bruno Salvy <Bruno.Salvy@inria.fr> and myself.

> Consider
> 
>   A(k)  =  sum( ((k)_j)^2 / (k^2-1)_j,  j = 1..k)                      (1)
> 
>              where  (k)_j = k! / (k-j)! is the "falling factorial"
> 
>              and  (k^2-1)_j = (k^2-1)! / (k^2-1-j)!
> 
> And that this can be expressed as 
>   
>   A(k)  =  F( 1, -k, -k; 1-k^2; -1) - 1                                (2)
> 
>   where  F  is the generalized hypergeometric function.
> 
> QUESTION:  We believe that asymptotically, 
> 
>    A(k) ~ (1/2)(k * Pi)^(1/2) - 1/6.  


We believe that actually  ["n" replacing "k"] the answer is

#                                       1/2                  1/2
#               1/2  1/2          23  Pi        47     457 Pi
#         1/2 Pi    n    - 1/6 + ---- ----- - ----- + ---- ----- + ...
#                                 96    1/2   270 n   9216   3/2
#                                      n                    n
#

Here is a brief indication of the computation. A Maple programme of
about 9 lines is given in the appendix. The rough steps are as follows.

1. Compute the log of the general term of the sum. We use 
  2*log(1-x/n)-log(1-x/n^2);
then expand it in series w.r.t. x for x->0. Then just replace
x by r(r+1)/2, x^2 by r(r+1)(2r+1)/6, etc. This summation gives the
expansion of the log of the general r-th term in the sum.

It is easy to prove that the terms for r > n^0.5001 are exponentially small.

2. Take out exp(-r^2/n) and expand asymptotically. In this way, we get a 
valid expansion of the general term.

         2                           3        2        4        3        2
        r               1/3 r - 1/3 r  + 1/2 r    1/6 r  - 1/2 r  - 1/2 r
 exp(- ----) (1 - r/n + ----------------------- + ------------------------+...)
         n                          2                         3
                                   n                         n

3. Estimate the asymptotic expansion of sums of terms of the form
exp(-r^2/n)*r^a/n^b. This is achieved by means of Euler Maclaurin
summation which (miraculously?) works on this example.
Euler Maclaurin summation takes care of the replacement of sums by integrals.
For instance:

> eulermac(exp(-r^2/n)*r,r,2); 
                                                               2
                                                     2        r
                      2                 2           r  exp(- ----)
                     r                 r                       n
      - 1/2 n exp(- ----) - 1/2 exp(- ----) r - 1/6 --------------
                      n                 n                  n

[Some further hacks are needed to handle erf's and exp's at infinity!]
Mathematically, you only add and subtract exponentially small terms by
replacing upper limits by +infty back and forth, so that it's no problem!

Finally, you get the result....

	Philippe Flajolet & Bruno Salvy INRIA (Rocquencourt, France)


APPENDIX: Maple programme. 
Warning check accuracy by increasing "Order". Order=10 gives only
results till O(n^(-2)), which means 5 significant terms.

[No attempt was made at readability of Maple source]


Order:=10; # Tune!
2*log(1-x/n)-log(1-x/n^2);
s:=series(",x);
s2:=map(normal,asympt(convert([seq(op(2*i-1,s)*
       sum(x^op(2*i,s),x=1..r),i=1..Order-1)],`+`),n));
exp(-r^2/n)*asympt(exp(s2+r^2/n),n);
eulmac:=readlib(eulermac)(eval(subs(O=<0>,")),r);
eulmac:=subs(exp(r^2/n)=1/exp(-r^2/n),eulmac);
# The result
res.Order:=asympt(eval(subs(O=<0>,erf=<1>,exp=<0>,eulmac)-subs(O=<0>,
        r=0,eulmac)),n);

# With Order=10: This is correct till O(1/n^2):
res10:=
91/2*Pi**(1/2)*n**(1/2)-1/6+23/96/n**(1/2)*Pi**(1/2)-47/270/n+457/9216/n**(3/2
)*Pi**(1/2)+79/1260/n**2-5837/30720/n**(5/2)*Pi**(1/2)-37/315/n**3+96347/552960
/n**(7/2)*Pi**(1/2)+32041/90720/n**4+40009/414720/n**(9/2)*Pi**(1/2)-51463/
9979200/n**5-80119/725760/n**(11/2)*Pi**(1/2)-8361533/29937600/n**6-24137/
302400/n**(13/2)*Pi**(1/2)-341273/359251200/n**7+1/120/n**(15/2)*Pi**(1/2)+
348713/314344800/n**8+1978139/2095632000/n**9+O(1/n**10):

#                                           1/2                  1/2
#                   1/2  1/2          23  Pi        47     457 Pi
#             1/2 Pi    n    - 1/6 + ---- ----- - ----- + ---- ----- + ...
#                                     96    1/2   270 n   9216   3/2
#                                          n                    n
#