Moments in terms of Cumulants

This recurrence relation generates expressions for the moments in terms of the cumulants, using formulae derived in: 
G. Bunker, “Application of the Ratio Method of EXAFS Analysis to Disordered Systems”, 
Nuclear Instruments and Methods in Physics Research (1983)
If you use it, please cite it.

writefile=False;(*flag whether to write out the expressions to a file in the same directory (folder) as this notebook*)

nmax=50;(*maximum order of cumulant expression to generate*)(*each increase in nmax by 5 required about a factor of 3 in execution time*)
nextn=D[#,γ]/.{D[c[n_,γ],γ]->c[n+1,γ]}&;
Print["timing=",AbsoluteTiming[tableau=NestList[nextn,Exp[c[0,γ]],nmax];]]

tableau[[1;;5]]//ColumnForm

Length[tableau[[-1]]](*# of terms in highest order expression*)

after it’s computed, the following transformation packages up the “tableau” into a more conventional form.

moments=Drop[tableau//.{c[n_,γ]->c[n],c[0]->0},1];

moments[[1;;5]]//ColumnForm

moments[[1;;5]]/.p[n_]->(Subscript[p, n])//ColumnForm

%/.c[1]->0

this writes out the file - it can take some space

If[writefile,SetDirectory[NotebookDirectory[ ]];Export["moments.txt",moments,"Lines"]];

check that it exported OK - for some reason it reads it as a string, so use ToExpression

If[writefile,test=ToExpression@Import["moments.txt","Lines"];test==moments]

time the direct inversion from moments to cumulants (presumably similar the other way)

(*eqns=Thread[Table[(Subscript[p, n]),{n,1,nmax}]==moments[[1;;nmax]]];
vars=Table[(Subscript[c, n]),{n,1,nmax}];
AbsoluteTiming[soln=Solve[eqns,vars ];]*)

empirically the timing for inverting the system of equations seems to increase by about a factor 1.57 for each increase in n by 1
so it hits the wall badly
it’s  5.77 second for nmax=20
according to this scaling, directly solving for moments from cumulants would take months to solve for n=50.
But I don’t think it would work at all because of memory limitations.
Good thing we can calculate by direct recurrence, which takes a minute or two up to n=50

5.77*1.57^(30)/3600 hr

% day/(24 hr)