The Absent - Minded Passengers Problem : 
A Motivating Challenge Solved by Computer Algebra
[arXiv:2003.01921]
-a supplementing notebook

Carsten.Schneider@risc.jku.at

Exit[]

<<Sigma.m;

Section 3: the generating function approach

We define the sum F=(Subscript[f, n]^((k)))(w) (see equation (3.3)) discovered by Doron Zeilberger [arXiv:2001.06839]:

F= (((Underoverscript[(∑^(.)), (r=0), k])(((({{k}, {r}})))^(.)) r(Subscript[!, .]) (((1-w))^(Underscript[k-r, .])) (w^(Underscript[r, .])) (1+r w) (Subscript[(2+r w), -1-k+n]))/(n(Subscript[!, .])));

and compute the moments (Subscript[M, l])(10,2) for l=0,1,2,3:

FSpec=F/.{n->10,k->2};

MSpec0=FSpec/.w ->1

MSpec1=w*D[FSpec,w]/.w-> 1

MSpec2=w*D[w*D[FSpec,w],w]/.w->1

MSpec3=w*D[w*D[w*D[FSpec,w],w],w]/.w->1

which allows us to compute the moment (Subscript[m, 3])(10,2)

mSpec3=-MSpec0*MSpec1^3+3*MSpec1^3-3*MSpec2*MSpec1+MSpec3

The direct approach (and some functions defined for later)

Besides the harmonic numbers (built into the packages Sigma and HarmonicSums)

S[o,n]=(Underoverscript[(∑^(.)), (i=1), n])(1/(i^o))

we will use the notation

(OverBar[S])[o,n]=(Underoverscript[(∑^(.)), (i=1), (n-1)])(1/(i^o));

by defining

(OverBar[S])[A_,B_Integer]:=S[A,B-1]

Furthermore, we use

ZeroIfNegative[k_Integer]:=If[k<0,0,1]

In the following we assume k≥1. Thus we can define

ZeroIfNegative[k]:=1
ZeroIfNegative[k-1]:=1

Finally, we define

FSumEv[l_]:=Sum[ZeroIfNegative[r](F/.SumName[A_,_]->A),{r,k-l,k}]//QNormalize

which enables us to compute (Subscript[M, l])(n,k) for l=0,1,2,3,... (see equation (3.4)) of the article:

M0=FSumEv[0]/.w->1

M1=w*D[FSumEv[1]//ToMathematica,w]/.w->1//ToSigma//CollectProdSum

M2=w *D[w*D[FSumEv[2]//ToMathematica,w],w]/.w->1//ToSigma//CollectProdSum

M3=w*D[w *D[w*D[FSumEv[3]//ToMathematica,w],w],w]/.w->1//ToSigma//CollectProdSum;

M4=w*D[w*D[w *D[w*D[FSumEv[4]//ToMathematica,w],w],w],w]/.w->1//CollectProdSum;

Next, we synchronize the arising objects to a common basis (i.e., all arising sums and products are algebraically among each other):

ML=SigmaReduceList[{M0,M1,M2,M3,M4},{k,n},Assumptions->k≥4];
ML=Collect[ML,ZeroIfNegative[_],CollectProdSum]

ML//Variables

Finally, we can compute the moments (Subscript[m, l])(n,k) for l=0,1,2,3,... by using the formula (3.2) of the article:

m[l_]:=Collect[Sum[Binomial[l,i](-1)^(l-i)ML[[i+1]]ML[[2]]^(l-i),{i,0,l}],ZeroIfNegative[_],CollectProdSum]

m0=m[0]

m1=m[1]

m2=m[2]

We observe that the expression

(((-1+k) k ZeroIfNegative[-2+k])/((-1+n) n))

vanishes for all k >= 1. To remove such unnecessary expressions we introduce the following functions:

ZeroIfNegative[f_,k_,v_]:=
Module[{start,res,const},
res=CollectProdSum[f/.(OverBar[S])->S,3]/.S->(OverBar[S]);
If[Head[res]===Plus,
res=Apply[List,res];
const=Apply[Plus,Select[res,FreeQ[#,(OverBar[S])[_,_]]&]];
res=Select[res,!FreeQ[#,(OverBar[S])[_,_]]&];
If[const=!=0,
res=Append[res,const];
];
If[Length[res]>1,
Return[Apply[Plus,Map[ZeroIfNegative[#,k,v]&,res]]]
]
];
start=-v-1;
While[start≥1&&(f/.k->start)===0,
start=start-1;
];
f*ZeroIfNegative[k-start-1]
]

ReduceExpr[f_]:=
Module[{res},
res=Collect[f/.S[A_,B_]->(OverBar[S])[A,B]+1/B^A,{ZeroIfNegative[_],(OverBar[S])[_,_]},Factor];
res=res/.(A_.)*ZeroIfNegative[B_]:>ZeroIfNegative[A,k,B-k];
res=Collect[res,{ZeroIfNegative[_],(OverBar[S])[_,_]},Factor];
res
]

and simplify  (Subscript[m, 2 ])(n,k) further:

m2=ReduceExpr[m[2]]

Similarly, we can compute (Subscript[m, l ])(n,k) for l=3,4

m3=ReduceExpr[m[3]]

m4=ReduceExpr[m[4]]

Given these symbolic expressions we can specialize, e.g., to the cases k=1 and k=2:

Collect[m4/.k->1,(OverBar[S])[_,_],Together]

Collect[m4/.k->2,(OverBar[S])[_,_],Together]

An efficient bypass via the exponential generating function

Finally, we follow a more efficient approach by computing the coefficients of ordinary generating function via its exponential generating function.  In the following we consider the moments up to l=7:

lMax=7;

Step 1: Calculate the exponential moments (Subscript[(OverBar[M]), l])(n,k)=(OverBar[M])[l]

(OverBar[M])[l_]:=Nest[D[#,w]&,FSumEv[l]//ToMathematica,l]/.w->1//ToSigma

(OverBar[ML])=Table[(OverBar[M])[l],{l,0,lMax}];

var=Variables[(OverBar[ML])];
subst=MapThread[Rule,{var,SigmaReduceList[var,{k,n}]}];
(OverBar[ML])=(OverBar[ML])/.subst;
(OverBar[ML])=Collect[(OverBar[ML]),ZeroIfNegative[_],CollectProdSum];

(OverBar[ML])[[4]]

Step 2: Calculate the ordinary moments (Subscript[M, l])(n,k)=M[l]

MLFast=Prepend[Table[Sum[(OverBar[ML])[[i+1]]StirlingS2[l,i],{i,1,l}],{l,1,lMax}],1];

MLFast[[4]]

Step 3: Calculate the desired moments (Subscript[m, l])(n,k)=m[l]

mFast[l_]:=Collect[Sum[Binomial[l,i](-1)^(l-i)MLFast[[i+1]]MLFast[[2]]^(l-i),{i,0,l}],ZeroIfNegative[_],CollectProdSum]

mL=Map[ReduceExpr,Table[mFast[l],{l,0,lMax}]];

Take[mL,5]-{m0,m1,m2,m3,m4}

mL[[7]]

Section 4: asymptotic expansions

The expansion for (Subscript[m, l])(n,k) for O(1/(n^4))

Finally, we compute the expansions the derived moments using the package

<<HarmonicSums.m

Here we switch back to our usual harmonic sum notation S[o,n]:

mm2=Collect[mL[[3]]/.(OverBar[S])[o_,a_]->S[o,a]-1/a^o,ZeroIfNegative[_],CollectProdSum]

and execute the HarmonicSums - command

m2Expand=SExpansion[mm2,n,3];

Finally, we carry out some simplifications and get:

m2Expand=Collect[m2Expand,{ZeroIfNegative[_],LG[_],EulerGamma,n},CollectProdSum]/.LG[n]->Log[n]+EulerGamma

The expansion for all computed moments l≤lMax

We repeat the above calculations for the full list: First, we switch back to the harmonic sum notation

mmL=Collect[mL/.(OverBar[S])[o_,a_]->S[o,a]-1/a^o,ZeroIfNegative[_],CollectProdSum];

and execute the HarmonicSums-command together with some simplificaitons:

mLExpand=SExpansion[mmL,n,3];

mLExpand=Collect[mLExpand,{ZeroIfNegative[_],LG[_],EulerGamma,n},CollectProdSum]/.LG[n]->Log[n]+EulerGamma;

and get the expansions in the following form

mLExpand[[1]]

mLExpand[[2]]

mLExpand[[3]]

mLExpand[[4]]

mLExpand[[5]]

mLExpand[[6]]

mLExpand[[7]]

mLExpand[[8]]