Some D-finite and Some Possibly D-finite Sequences in the OEIS

Manuel Kauers and Christoph Koutschan

The paper has been published in the Journal in Integer Sequences, and can be found at  https://cs.uwaterloo.ca/journals/JIS/ or at https://arxiv.org/abs/2303.02793

(* Specify directory where the data files are stored. *)
SetDirectory[NotebookDirectory[]];
(* A procedure to load and display relevant information about a sequence. *)
InitializeSeq[num_]:=
Module[{str,ord},
str=StringSplit/@Import["b"<>ToString[num]<>".txt","Lines"];
data=ToExpression/@Last/@str;
Print["Sequence: ",StringTake[ToString[Take[data,10]],{2,-2}]<>", ..."];
Print["Length: ",Length[data]];
Print["Offset: ",offset=ToExpression[str[[1,1]]]];
Print["Recurrence: ",rec=Get["rec"<>ToString[num]<>".m"]];
ord=Max[Cases[rec,a[n+nn_.]:>nn,Infinity]];
Print["Check: ",MatchQ[
Table[rec,{n,offset,offset+Length[data]-ord-1}]/.a[n_]:>data[[n-offset+1]],{(0)..}]];
];
(* Load RISC packages, see http://www3.risc.jku.at/research/combinat/software/ergosum/ *)
<<RISC`GeneratingFunctions`;
<<RISC`Guess`;
<<RISC`HolonomicFunctions`;

Package GeneratingFunctions version 0.8 written by Christian Mallinger
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Guess Package version 0.52
written by Manuel Kauers
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

1  Introduction

1.1  A237684

InitializeSeq[237684]

1.1  A039836

InitializeSeq["039836"]

1.2  A187990

InitializeSeq[187990]

rec1=Collect[rec/((n-26)(n-27)),a[_],Factor]

Table[rec1,{n,1,49}]/.a[n_]:>data[[n]]

RSolve[{rec1==0,a[1]==117},a[n],n]

Table[%[[1,1,2]],{n,27}]

Take[data,27]

ListPlot[Differences[Take[data,{15,35}]],Joined->True]

(* Assuming that n denotes the maximal possible value, substitute n -> n+4 to get OEIS sequence. *)
Enumerate187990[n_]:=
Sort/@DeleteCases[Join[
Flatten[Table[{-x1,-x2,-x3,x3,x2,x1},{x1,1,n},{x2,x1,n},{x3,x2,n}],2],
Flatten[Table[{-x1,-x2,x2-1,x2-1,x1-1,x1-1},{x1,1,n},{x2,1,x1}],1],
Flatten[Table[If[x1≠x2+1,{-x1,-x2,-x2,x2+1,x1-1,x1-1},{}],{x1,1,n},{x2,1,n-1}],1],
Flatten[Table[{-x1,-x1,-x2,-x2,x2+1,x1+1},{x1,1,n-1},{x2,1,x1}],1],
Flatten[Table[If[x1≠x2+1,{-x1,-x2,x2,x1-2,x1-2,x1-1},{}],{x1,2,n},{x2,1,n}],1],
Flatten[Table[If[x1≠x2,{-x1,-x1+1,-x1+1,-x2,x2,x1+1},{}],{x1,2,n-1},{x2,1,n}],1],
Table[{-x1-3,x1,x1,x1,x1,x1+2},{x1,0,n-3}],
Table[{-x1-3,x1,x1,x1+1,x1+1,x1+1},{x1,0,n-3}],
Table[{-x1-4,x1,x1,x1+1,x1+2,x1+3},{x1,0,n-4}],
Table[{-x1-2,-x1,-x1,-x1,-x1,x1+3},{x1,1,n-3}],
Table[{-x1-1,-x1-1,-x1-1,-x1,-x1,x1+3},{x1,1,n-3}],
Table[{-x1-3,-x1-2,-x1-1,-x1,-x1,x1+4},{x1,1,n-4}]
],{}];
Enumerate187990[5]

Length[%]

(* Count all instances produced above. *)
{Sum[1,{i,n},{j,i,n},{k,j,n}],
Sum[1,{i,n},{j,i}],
(n-1+Sum[n-2,{i,2,n}]),
Sum[1,{i,n-1},{j,i}],
Sum[n-1,{i,2,n}],
Sum[n-1,{i,2,n-1}],
n-3+2*(n-2),
n-4+2*(n-3)}

Together[Total[%]]

(* Indeed, we get the same polynomial, that solved the guessed recurrence. *)
Together[%/.n->n+4]

3  Transfer Matrix Method

Note: all transfer matrices here are transposed, because Mathematica computes M.v more efficiently than v.M.

3.1  A177317

InitializeSeq[177317]

(* All solutions for n=2. *)
Select[Permutations[Join@@Table[Range[5],{2}]],Max[Abs[Differences[#]]]≤1&]

Length[%]

(* Transfer matrix *)
MatrixForm[tmat=Table[If[Abs[i-j]≤1,x[i],0],{i,5},{j,5}]/.x[5]->1]

vec=svec=Append[Array[x,4],1];
t0=AbsoluteTime[];
Table[
vec=Expand[tmat.vec];
If[Mod[n5,5]===0,
cf=Coefficient[Total[vec],(Times@@Array[x,4])^(n5/5)];
Print["n=", n5/5," (timing = ",-t0+(t0=AbsoluteTime[]),"): ",cf];
];
,{n5,2,60}];

ByteCount[vec]

(* The rational generating function. *)
mat1=IdentityMatrix[5]-t*tmat;
rat=Total[Expand[Table[
(-1)^(i+j)*Det[Delete[#,i]&/@Delete[mat1,j]],{i,5},{j,5}].svec]]/Det[mat1]

(* Sanity check (by series expansion). *)
ser=Normal[Series[rat,{x[1],0,4},{x[2],0,4},{x[3],0,4},{x[4],0,4},{t,0,19}]];
Table[Coefficient[ser,(x[1] x[2] x[3] x[4])^i*t^(5i-1)],{i,4}]

Timing[ct1=CreativeTelescoping[rat/((x[1] x[2] x[3] x[4])^(n+1) t^(5n)),Der[x[1]],{Der[x[2]],Der[x[3]],Der[x[4]],Der[t],S[n]}][[1]];
#[ct1]&/@{ByteCount,UnderTheStaircase}
]

Timing[
ct2=CreativeTelescoping[ct1,Der[x[2]]][[1]];
#[ct2]&/@{ByteCount,UnderTheStaircase}
]

Timing[
ct3=FindCreativeTelescoping[ct2,Der[x[3]]][[1]];
#[ct3]&/@{ByteCount,UnderTheStaircase}
]

Timing[
ct4=FindCreativeTelescoping[ct3,Der[x[4]]][[1]];
#[ct4]&/@{ByteCount,UnderTheStaircase}
]

Timing[
ct5=FindCreativeTelescoping[ct4,Der[t]][[1]];
#[ct5]&/@{ByteCount,UnderTheStaircase}
]

ct5

(* Compare with the guessed recurrence: *)
Together[ApplyOreOperator[ct5,a[n]]/rec]

3.2  A199250

InitializeSeq[199250]

(* Produce more terms with the guessed recurrence, and perform two checks:
values agree with existing data, all new values are integers. *)
re2l=RE2L[Prepend[Table[a[i]==data[[i]],{i,22}],rec==0],a[n],{1,1000}];
{Take[re2l,Length[data]]===data,And@@(IntegerQ/@re2l)}

(* Guess minimal-order recurrence. *)
mrec=Collect[GuessMinRE[Take[re2l,98],a[n],StartPoint->1],a[_],Factor]

(* Guess recurrence for even indices (ae[n] = a[2n]). *)
rece=Collect[Numerator[Together[
GuessUnivRE[Take[re2l,{2,200,2}],ae[n],Order->4,Degree->8,StartPoint->1][[1]]]],ae[_],Factor]

(* Guess recurrence for odd indices (ao[n] = a[2n+1]). *)
reco=Collect[Numerator[Together[
GuessUnivRE[Take[re2l,{1,200,2}],ao[n],Order->4,Degree->10,StartPoint->0][[1]]]],ao[_],Factor]

(* These are all possible rows: *)
rows=DeleteCases[Tuples[Range[0,3],2],{a_,a_}]

(* Construct the transfer matrix. *)
tmat=Table[With[{r1=rows[[j]],r2=rows[[i]]},
If[r1[[1]]===r2[[1]]||r1[[2]]===r2[[2]],0,x^Count[r2,0]*y^Count[r2,1]*z^Count[r2,2]]
],{i,12},{j,12}];
MatrixForm[tmat]

(* Compute values with the transfer matrix approach. *)
vec=Prepend[Table[0,{11}],x y];
vvec=First[DeleteCases[#,0]]&/@tmat;
tmat1=tmat/.{x->1,y->1,z->1};
Timing[
data1=Table[
vec=Expand[vvec*(tmat1.vec)];
val=Table[n/2,{4}];
If[OddQ[n],val+={1,1,-1,-1}/2];
val=(Times@@({x,y,z,1}^#))&/@Permutations[val];
val=Coefficient[Total[vec],val];
Total[val/2],{n,2,100}];
]

(* Sanity check. *)
Take[data1,55]-Rest[data]

mat1=IdentityMatrix[12]-t*tmat;
rat=Factor[Sum[t x y*(-1)^(i+1)*Det[Delete[#,i]&/@Rest[mat1]],{i,12}]/Det[mat1]]

(* Test coefficient extraction formula. *)
ser=Normal[Series[rat,{t,0,10},{x,0,5},{y,0,5},{z,0,5}]];
Table[With[{n=nn/2},1/2*Coefficient[ser,(x y z)^n t^(2n)]],{nn,2,10,2}]===Take[data,{2,10,2}]
Table[With[{n=(nn-1)/2},1/2*Coefficient[(x+y+x y+z+x z+y z)*ser,(x y z)^(n+1) t^(2n+1)]],{nn,3,9,2}]===Take[data,{3,9,2}]

Timing[ct1=CreativeTelescoping[rat/(x y z)^(n+1)/t^(2n+1),Der[t],{Der[x],Der[y],Der[z],S[n]}][[1]];]

Timing[ct2=FindCreativeTelescoping[ct1,Der[x]][[1]];]

Timing[ct3=FindCreativeTelescoping[ct2,Der[y]][[1]];]

Timing[ct4=FindCreativeTelescoping[ct3,Der[z]][[1]];]

Exponent[opeven=ct4[[1]],{S[n],n}]

(* Sanity check. *)
test=ApplyOreOperator[opeven,f[n]];
Table[test,{n,0,100}]/.f[n_]:>re2l[[2n]]

(* This operator is a left multiple of the guessed recurrence. *)
OreReduce[opeven,{ToOrePolynomial[rece,ae[n]]}]

(* For the odd case, modify the rational function accordingly. *)
ratodd=Together[Total[(rat/#)&/@((t x y z)t^(2n+1)(x y z)^n*{x,y,z,x y,x z,y z})]]

Timing[ct1=CreativeTelescoping[ratodd,Der[t],{Der[x],Der[y],Der[z],S[n]}][[1]];]

Timing[ct2=FindCreativeTelescoping[ct1,Der[x]][[1]];]

Timing[ct3=FindCreativeTelescoping[ct2,Der[y]][[1]];]

Timing[ct4=FindCreativeTelescoping[ct3,Der[z]][[1]];]

Exponent[opodd=ct4[[1]],{S[n],n}]

(* Sanity check. *)
test=ApplyOreOperator[opodd,f[n]];
Table[test,{n,0,100}]/.f[n_]:>re2l[[2n+1]]

(* This operator is a left multiple of the guessed recurrence. *)
OreReduce[opodd,{ToOrePolynomial[reco,ao[n]]}]

(* The two operators can be combined to an operator annihilating the whole sequence. *)
op=First[DFinitePlus[DFiniteSubstitute[{opeven},{n->n/2}],DFiniteSubstitute[{opodd},{n->(n-1)/2}]]];
Exponent[op,{S[n],n}]

(* Verify that it is compatible with the guessed order-8 recurrence. *)
OreReduce[op,ToOrePolynomial[{mrec},a[n]]]

3.3  A250556

InitializeSeq[250556]

(* The C-finite recurrence can be guessed by omitting the first 11 terms. *)
GuessMinRE[Take[data,{12,47}],f[n]]

(* Produce more terms with the guessed recurrence, and perform two checks:
values agree with existing data, all new values are integers. *)
re2l=RE2L[Prepend[Table[a[i]==data[[i]],{i,22}],rec==0],a[n],{1,1000}];
{Take[re2l,Length[data]]===data,And@@(IntegerQ/@re2l)}

(* Minimal-order recurrence *)
mrec=Collect[Numerator[Together[GuessMinRE[Take[re2l,{12,61}],a[n],StartPoint->11]]],a[_],Factor]

(* Naive implementation for enumeration. *)
D2[a_]:=Differences[Differences[a]];
MyTest[v_]:=MemberQ[v.#&/@Tuples[{1,-1},Length[v]],0];
Do[Print[n,": ",Timing[Length[Select[Tuples[{0,1,2,3},n+2],MyTest[D2[#]]&]]]],{n,7}]

(* Derive a global bound for the partial signed sums. *)
(* Example for bnd=18 *)
D2[{3,0,3,0,3,1,3,0,1}]

(* Example for bnd=19 *)
D2[{1,3,0,2,0,3,0,3,0,3,1}]

(* Construct the transfer matrix, using the optimal bound B=19. *)
bnd=19;
Timing[
(* Possible start configurations (after one step). *)
s1=Flatten[Table[{{a2,a3},Intersection[{1,-1}*(a1-2a2+a3),Range[0,bnd]]}
,{a1,0,3},{a2,0,3},{a3,0,3}],2];
(* Transition function. *)
NextState[{{a1_Integer,a2_Integer},set_List},k_Integer]:=
{{a2,k},Intersection[Abs[Flatten[Outer[Plus,set,{-1,1}*(a1-2a2+k)]]],Range[0,bnd]]};
(* Construct set of reachable states. *)
states=FixedPoint[Function[s,Union[s,Union@@(Table[NextState[#,k],{k,0,3}]&/@s)]],Union[s1]];
(* start vector v_{init} and end vector v_{final}. *)
svec=Table[0,{Length[states]}];
(svec[[Position[states,#][[1,1]]]]+=1)&/@s1;
evec=If[MemberQ[#2,0],1,0]&@@@states;
(* Build the transfer matrix. *)
tmat=row0=SparseArray[{Table[0,{Length[states]}]}];
Function[s,
row=row0;
Do[row[[1,Position[states,NextState[s,k]][[1,1]]]]=1,{k,0,3}];
tmat=Join[tmat,row];
]/@states;
tmat=Transpose[Rest[tmat]];
{Dimensions[tmat],Count[Normal[tmat],1,2],ByteCount/@{tmat,Normal[tmat]}}
]

dim=Dimensions[tmat][[1]]

(* Compute values using the transfer matrix. *)
vec=svec;Table[vec=tmat.vec;vec.evec,{10}]

(* Computation is very fast. Values agree with those produced by the guessed recurrence. *)
Timing[vec=svec;Rest[re2l]===Table[vec=tmat.vec;vec.evec,{999}]]

(* Number of nonzero entries in v_init and v_final. *)
Count[Normal[#],1|2]&/@{svec,evec}

(* Positions of these nonzero entries. *)
epos=Flatten[Position[evec,1]];
spos=Flatten[Position[svec,1|2]]

(* Estimated time when applying the determinant formula naively: 60*720*3s = 36h. *)
With[{i=227,j=1392},Timing[
Factor[(-1)^(i+j)*Det[Delete[#,i]&/@Delete[IdentityMatrix[dim]-t*tmat,j]]]
]]

(* Do one column at a time. *)
mat1=IdentityMatrix[dim]-t*tmat;
Timing[
nums=Function[i,Det[ReplacePart[mat1,i->evec]]]/@spos;
]

(* The rational generating function. *)
rat=t*Together[(svec[[spos]].nums)/Det[mat1]]

(* Sanity check. *)
Series[rat,{t,0,46}]-data.t^Range[47]

Appendix: show that B=19

(* Find all pairs of multisets (S1,S2) with values in {2,...,6} such that
 Total[S1]=Total[S2]=b and there are no nontrivial subsets T1 of S1 and
T2 of S2 with Total[T1]=Total[T2]. *)
Test0[b_]:=Union[Join@@(
Test1[b,Append[#,6],Complement[Range[2,5],#]]&/@Subsets[Range[2,5],{0,3}])];
Test1[b_,s1_,s2_]:=
Module[{sol1,sol2},
sol1=Solve[Prepend[x[#]≥0&/@s1,Total[#*x[#]&/@s1]==b],x/@s1,Integers];
sol2=Solve[Prepend[x[#]≥0&/@s2,Total[#*x[#]&/@s2]==b],x/@s2,Integers];
Return[Join@@Flatten[Outer[Test2,sol1,sol2,1],1]];
];
Test2[sol1_,sol2_]:=
Module[{s1=sol1,s2=sol2,ps1,ps2},
{s1,s2}=Flatten/@({s1,s2}/.HoldPattern[x[i_]->k_]:>Table[i,{k}]);
{ps1,ps2}=Union[Total/@Subsets[#,{1,Length[#]-1}]]&/@{s1,s2};
Return[If[Intersection[ps1,ps2]==={},{Sort[{s1,s2}]},{}]];
];

Test0[18]

Join@@Table[Test0[b],{b,19,30}]

(* Test which arrays can be built from the given second differences,
and whose partial sums are at least b. *)
PartialSums[v_]:=FoldList[Plus,First[v],Rest[v]];
Diff2Array[d2_List]:=
Module[{n2=Length[d2]+2,sol},
sol=Solve[Join[Thread[D2[Array[x,n2]]==d2],Table[0≤x[i]≤3,{i,n2}]],Integers];
If[sol=!={},Array[x,n2]/.sol,{}]];
Test3[b_,s1_,s2_]:=
Module[{perms,svect,diffs},
perms=Abs[Select[Permutations[Join[s1,-s2]],Max[Abs[PartialSums[#]]]≥b&]];
svect=(-1)^Range[Length[Join[s1,s2]]];
diffs=Join[#,-#]&[svect*#&/@perms];
Return[Union@@(Diff2Array/@diffs)];
];

Test3[18,{6,6,6},{4,4,5,5}]

Test3[19,{1,1,6,6,6},{5,5,5,5}]

Test3[19,{2,6,6,6},{5,5,5,5}]

Test3[19,{5,5,5,5},{4,4,4,4,4}]

Test3[19,{6,6,6,6},{4,5,5,5,5}]

Test3[19,{6,6,6,6,6},{5,5,5,5,5,5}]

3.4  A264947

The OEIS lists only 20 terms of this sequence. We computed 80 terms (on compute-purley, crashed probably because out-of-memory), but were not able to guess a recurrence with them.

(* Compute values with the transfer matrix method. *)
cols=Tuples[{0,1,2,3},4];
tmat=Table[If[Or@@MapThread[SameQ,cols[[{i,j}]]],0,1],{i,256},{j,256}];
vec=vvec=((x^Count[#,0]*y^Count[#,1]*z^Count[#,2])&/@cols);
Do[
Print[n,": ",Timing[
vec=Expand[vvec*(tmat.vec)];
Coefficient[Total[vec],(x y z)^n]/24
]],{n,2,10}]

(* We continue this computation on RICAM's compute server radon1. *)
(* Terms n=2,...,53 were computed on grantley processors. *)
(* Terms n=54,...,80 on purley processors. *)
(* Memory consumption at n=80 was more than 215 GB. *)
str="2:{0.069961,60}
3:{0.637427,3201}
4:{2.19072,184740}
5:{7.41572,11375145}
6:{12.8547,730983420}
7:{29.5854,48402531561}
8:{42.3991,3282992503164}
9:{83.2498,226854309720993}
10:{110.189,15915758107113276}
11:{186.79,1130694005695927761}
12:{230.252,81177583723495750340}
13:{354.157,5880587303767912833417}
14:{420.631,429300706847441007321756}
15:{613.332,31551853305004056089812729}
16:{714.901,2332682592774096715790438556}
17:{976.644,173364378813531387237639561417}
18:{1112.93,12944676911197784353139378818620}
19:{1443.15,970602277004604478741978783654137}
20:{1608.48,73052030408618314498132204294989636}
21:{2025.26,5517104052146769136190406619946748625}
22:{2219.31,417970269693293204826940794292755282876}
23:{2729.9,31755527174421370850740379081100981773217}
24:{2970.76,2418975934412719050064221587794842339432332}
25:{3594.49,184710921342968669134163175096236855105044281}
26: {3886.21, 14135887348401610659310525755660141025805396412}
27: {4612.65, 1084057674533418597309164076479502490092692413401}
28: {4926.76, 83294687297695996487296133362481319140025771576036}
29: {5812.04, 6411520652304194600222219852762237278437530027109633}
30: {6187.07, 494347612655302255756667794797862393706004006990449116}
31: {7212.1, 38175540463975303292909694032596913469684571657360319457}
32: {7641.88, 2952414235489779837568180412502990267046488842542798111932}
33: {8759.51, 228649620771671501516499860320220154028314817910323524247169}
34: {9268.76, 17730862612235836591255750514436505927262075267548129408360380}
35: {10586.9, 1376649781393386397202658291679150213747668160928012174182043073}
36: {11110.7, 107009411222489145387956320384330468059442500701300626830913088676}
37: {12585.4, 8327188593517364995020327803326626009226948946693602804486316239977}
38: {13175.2, 648676069852905410796568057694024558259206275558011650142953096934972}
39: {15001.6, 50581001404316899018968316681587103308256422699912397306712929275865673}
40: {15508.6, 3947805659698685988683929734716872359238515550819491407967142406746868508}
41: {17342.9, 308399147222261800676330778658313601219073933946560673207166426335652031041}
42: {18083.1, 24112447409196765040878502756573179517492816198672275640090583636080034701756}
43: {20125.5, 1886787835319799938379818105928076723411184655104047209600302507158968998421489}
44: {20949., 147755130878557648268499773809712039244342634235501019196169297111769580735381956}
45: {23201., 11579370093128825450825933404080247878556065927068626000807392982817372618932202793}
46: {24201.8, 908105385301233478422538418034223400562994578658360475063321050637135075476054023516}
47: {26825.1, 71266181824726823921324040654867167193476523124997584531075712029317752760708333643257}
48: {28162.4, 5596473279872653409465262836355974132939541418081163631059581045467246292915567223014188}
49: {31032.8, 439761900249949927702660613606859404111446506457518493714468286067207732192957942294411465}
50: {32360.8, 34576592592762435155379942361008198747574980142123925550748615338429692566552624268045657916}
51: {35625.3, 2720182511711922680539482335341321276207325058505358344875059030057886710720822649964987650873}
52: {36998.4, 214119179002310180277790942534261104480999332069957117683550314081407443802463255672863982869060}
53: {40546.8, 16863422525092775259252335452842925622415206628251884617823220642863303695917030005288928030168529}
54: {31914.8, 1328801018696525777602445344173083632435553041457883803581803860918220856542739915902745466980311804}
55: {35993.5, 104758718179385737040026753049752525322347714180154345197928332394270658804300441247363834496259033153}
56: {37668.5, 8262828447336975530272573022369359792924359207406249293637243738718731372563460177389371732247647022060}
57: {39518.8, 652031199651779504888570623509874000388057141565552251584213310527666888959569975204911306830850688058457}
58: {39763.6, 51475694632693150107072891343755372168964819023491604515873799535111461308380490581426802331324769881443772}
59: {45254.8, 4065591279329931975158235402983295762669836107682393646134658539394692483489560663912807432021880226860638521}
60: {43636.8, 321237817485629068522445220415087399267431674153841598407716075255629895954797870222360788950526924517817179364}
61: {51463.2, 25392482346456698670646502854863139433321194310353726542797843251102949564306389470269197438622890523485826653281}
62: {49365.6, 2007953123100431854794314540037277408576704345488472126168170164181844048850732419071052334349195204253850332170652}
63: {54508.6, 158842436473528707512143792478680194148582530808323362946293065887278051769919703297171455138765693933499384187423009}
64: {55498.3, 12570106624527408113699615647995077654475021071876168273040747672144021267503343411793251591939503628864654208191115068}
65: {60498.1, 995098243893926458432903060576503230957066209999728231404752148278766913289813555976041858163555628471860362517815080993}
66: {54860.7, 78803023470047286558686199017340702508568762207108290549219331990771185871303618722560069077300055215480706234168728622524}
67: {66118.2, 6242596674229713086395357023615891509393326780357638438529182729730831346773767457402457146544400041445780138407520204679713}
68: {66369.7, 494685175538530733759081009215086333916650877173709150372053040636928548315747879920573161026398793698953496217726153364192292}
69: {70763.5, 39212963913823620993090949297287060246036209842940053638210655052758868077437501127122288972578558017793246272111934555532195529}
70: {70749.7, 3109307705134351325206283317455059077470968790483650476580296799186296334222927974359895757333615331575088952930618800132008315772}
71: {78957.2, 246619415658832597738416791502561390072125769600195798796215632558487202771110585021525326346036935805950800163408177257142535721737}
72: {80284.5, 19566663520666839840882651733210307855863950492933725157905822358646622480537959531427242680748867679477352618536455611436161412837500}
73: {86086.6, 1552847506947478482992935002289033140191687797989293116953567941285570920516507894067948193685887334569937922171070122349127475119818753}
74: {88913.2, 123270758766115761772260128591410880610874536574768815842069292551032514386675747704698864912873205168595481506871186636814260381915562684}
75: {93262.9, 9788302789344638394740988630216846588153937349783925069367345065408852770059346761784920672668456252188351450260757054760109098898363576113}
76: {100431., 777441554111519740898680565476153729510551339316083109439116794734881344678382582465198276256701446196196172088864748856495178757995718880324}
77: {101555., 61764398059971578808332533818681508065912573772099647446045251500133874859941707723264465334444502595531258397976849316160669503430010687766761}
78: {113436., 4908129096776091045320075917502905956325325264117748476732078466402607695408199236964804195727997948580993469737529180096064602461891862206108188}
79: {110954., 390120111786328113538562455616786929332168155592273598475569710202503605372602188863559933755059072571477476112028654125245706976834159866392851033}
80: {123981., 31015779471397020829939886469410531327177774799037091316989277835091394621986947957747347675373127854226194498480349096934751135478191634525159660796}";
data=Prepend[ToExpression[Last[StringSplit[#,{",","}"}]]]&/@StringSplit[str,"\n"],1];timings=Prepend[ToExpression[StringSplit[#,{",","{"}][[2]]]&/@StringSplit[str,"\n"],0];

(* Total CPU time (in days). *)
Total[timings]/3600/24

4  Lattice Walks

4.1  A265234

InitializeSeq[265234]

(* Compute values using a transfer matrix. *)
rows=DeleteCases[Tuples[{0,1,2,3},4],{___,a_,a_,___}];
tmat=Table[1,{i,Length[rows]},{j,Length[rows]}];
vec=vvec=((x^Count[#,0]*y^Count[#,1]*z^Count[#,2])&/@rows);
Prepend[Table[vec=Expand[vvec*(tmat.vec)];Coefficient[Total[vec],(x y z)^n]/24,{n,2,6}],1]

(* Define the step-set polynomial. *)
ssp=Total[(x^Count[#,0]*y^Count[#,1]*z^Count[#,2])&/@rows]

rat=1/(1-t*ssp)

(* Sanity check by series expansion. *)
ser=Normal[Series[rat,{t,0,6}]];
Table[Coefficient[ser,(t x y z)^n]/24,{n,6}]

(* Transformed rational function for diagonal extraction by residue. *)
rat1=Together[(rat/.{t->t/x,x->x/y,y->y/z})/(x y z)]

Timing[
ct1=CreativeTelescoping[rat1,Der[x],Der/@{t,y,z}][[1]];
#[ct1]&/@{ByteCount,UnderTheStaircase}
]

Timing[
ct2=CreativeTelescoping[ct1,Der[y]][[1]];
#[ct2]&/@{ByteCount,UnderTheStaircase}
]

Timing[
ct3=CreativeTelescoping[ct2,Der[z]][[1]];
#[ct3]&/@{ByteCount,UnderTheStaircase}
]

(* This is the differential equation satisfied by the diagonal. *)
Factor[ct3]

(* Convert it to a recurrence for its Taylor coefficients. *)
recop=CreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[ct3,S[n]-1],OreAlgebra[Der[t],S[n]]],Annihilator[t^(-n-1),{Der[t],S[n]}]],Der[t]][[1,1]]

(* Compare with the guessed recurrence. *)
Together[rec/ApplyOreOperator[recop,a[n]]]

View it as the extraction of a coefficient from powers of the step-set polynomial:

Timing[Table[Coefficient[ssp^n,(x y z)^n]/24,{n,10}]]

Timing[poly=1;Table[poly=Expand[poly*ssp];Coefficient[poly,(x y z)^n]/24,{n,1,10}]]

Timing[poly=1;tab=Table[poly=Expand[poly*ssp];Coefficient[poly,(x y z)^n]/24,{n,1,56}];]

4.2  A172572

InitializeSeq[172572]

(* We get a simpler recurrence by dividing out a hypergeometric factor. *)
Factor[DFiniteTimesHyper[ToOrePolynomial[rec,a[n]],1/Binomial[3n,n]]]

9-fold multinomial sum, reduced to 6-fold by Chu-Vandermonde

(* The entries of the multinomial> *)
expr=Array[c,15]/.Solve[Array[c,15].Permutations[{1,1,0,0,0,0}]==Table[n,{6}]][[1]];
vars=Union[Cases[expr,c[_],Infinity]];
TableForm[expr]

(* Sanity check. *)
Clear[A172572];
A172572[n_]:=A172572[n]=
Sum@@Prepend[Table[{vars[[i]],0,n},{i,9}],Multinomial@@expr];
Do[Print[Timing[A172572[n]]],{n,3}]

(* Adapt the ranges of summation for faster computation. *)
Clear[A172572];
A172572[n_]:=A172572[n]=
Sum[Multinomial[c[1],c[2],c[3],n-c[1]-c[2]-c[3]-c[4],c[4],c[6],c[7],n-c[1]-c[6]-c[7]-c[8],c[8],c[10],-n+c[1]+c[2]+c[3]+c[6]+c[7]+c[10],n-c[2]-c[6]-c[10]-c[11],2 n-c[1]-c[2]-c[3]-c[4]-c[6]-c[7]-c[8]-c[10]-c[11],c[11],-n+c[1]+c[2]+c[4]+c[6]+c[8]+c[11]],
{c[1],0,n},
{c[2],0,n-c[1]},
{c[3],0,n-c[1]-c[2]},
{c[4],0,n-c[1]-c[2]-c[3]},
{c[6],0,Min[n-c[1],n-c[2],2n-c[1]-c[2]-c[3]-c[4]]},
{c[7],0,Min[n-c[1]-c[6],2 n-c[1]-c[2]-c[3]-c[4]-c[6]]},
{c[8],0,Min[n-c[1]-c[6]-c[7],2 n-c[1]-c[2]-c[3]-c[4]-c[6]-c[7]]},
{c[10],Max[0,n-c[1]-c[2]-c[3]-c[6]-c[7]],Min[n-c[2]-c[6],2 n-c[1]-c[2]-c[3]-c[4]-c[6]-c[7]-c[8]]},
{c[11],Max[0,n-c[1]-c[2]-c[4]-c[6]-c[8]],Min[n-c[2]-c[6]-c[10],2 n-c[1]-c[2]-c[3]-c[4]-c[6]-c[7]-c[8]-c[10]]}];

Do[Print[n,": ",Timing[A172572[n]]],{n,14}]

Take[data,10]-Table[A172572[n],{n,10}]

(* We can apply the Chu-Vandermonde identity. *)
ChuVandermonde=sum[1/k!/(m-k)!/(r-k)!/(n-r+k)!,{k,0,r}]==(m+n)!/r!/(m+n-r)!/m!/n!;

expr1=FunctionExpand[Multinomial@@expr]/.Gamma[n_]->(n-1)!;
expr2=Select[expr1,Not[FreeQ[#,c[11]]]&]

test=ChuVandermonde/.{k->c[11],r->n-c[2]-c[6]-c[10],m->2 n-c[1]-c[2]-c[3]-c[4]-c[6]-c[7]-c[8]-c[10],n->c[1]+c[4]+c[8]-c[10]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]];
expr2=Select[expr1,Not[FreeQ[#,c[8]]]&]

test=ChuVandermonde/.{k->c[8],r->n-c[1]-c[6]-c[7],n->n+c[4]-c[6]-c[7]-c[10],m->2 n-c[1]-c[2]-c[3]-c[4]-c[6]-c[7]-c[10]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]];
expr2=Select[expr1,Not[FreeQ[#,c[4]]]&]

test=ChuVandermonde/.{k->c[4],r->n-c[1]-c[2]-c[3],n->2 n-c[1]-c[2]-c[3]-c[6]-c[7]-c[10],m->2 n-c[1]-c[2]-c[3]-c[6]-c[7]-c[10]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]]

(* Hence, we end up with a 6-fold hypergeometric sum. *)
Clear[A172572];
A172572[n_]:=A172572[n]=
Sum[((3 n)! (4 n-2 c[1]-2 c[2]-2 c[3]-2 c[6]-2 c[7]-2 c[10])!)/(c[1]! c[2]! (n-c[1]-c[2]-c[3])! c[3]! c[6]! (n-c[1]-c[6]-c[7])! c[7]! (n-c[2]-c[6]-c[10])! (n-c[3]-c[7]-c[10])! ((((2 n-c[1]-c[2]-c[3]-c[6]-c[7]-c[10])!))^2) c[10]! (-n+c[1]+c[2]+c[3]+c[6]+c[7]+c[10])!),
{c[1],0,n},
{c[2],0,n-c[1]},
{c[3],0,n-c[1]-c[2]},
{c[6],0,Min[n-c[1],n-c[2]]},
{c[7],0,Min[n-c[1]-c[6],n-c[3]]},
{c[10],Max[0,n-c[1]-c[2]-c[3]-c[6]-c[7]],Min[n-c[2]-c[6],n-c[3]-c[7]]}];

Do[Print[n,": ",Timing[A172572[n]]],{n,20}]

(* By exploiting symmetries, we can make the enumeration more efficient. *)
Clear[A172572];
A172572[n_]:=A172572[n]=
Sum[Switch[Length[Union[c/@{1,2,3}]],1,1,2,3,3,6]*
Sum[((3 n)! (4 n-2 c[1]-2 c[2]-2 c[3]-2 c[6]-2 c[7]-2 c[10])!)/(c[1]! c[2]! (n-c[1]-c[2]-c[3])! c[3]! c[6]! (n-c[1]-c[6]-c[7])! c[7]! (n-c[2]-c[6]-c[10])! (n-c[3]-c[7]-c[10])! ((((2 n-c[1]-c[2]-c[3]-c[6]-c[7]-c[10])!))^2) c[10]! (-n+c[1]+c[2]+c[3]+c[6]+c[7]+c[10])!),
{c[6],0,Min[n-c[1],n-c[2]]},
{c[7],0,Min[n-c[1]-c[6],n-c[3]]},
{c[10],Max[0,n-c[1]-c[2]-c[3]-c[6]-c[7]],Min[n-c[2]-c[6],n-c[3]-c[7]]}],
{c[1],0,n},
{c[2],c[1],n-c[1]},
{c[3],c[2],n-c[1]-c[2]}];

Timing[Do[Print[n,": ",Timing[A172572[n]]],{n,33}]]

Table[A172572[n],{n,33}]-data

Timing[Do[Print[n,": ",Timing[A172572[n]]],{n,34,44}]]

Interpretation as 6-dimensional walks

steps=Permutations[{1,1,0,0,0,0}];
Timing[With[{n=20},
vals=Table@@Prepend[Table[{c[i],0,If[i===1,n,c[i-1]]},{i,6}],0];
vals[[1,1,1,1,1,1]]=1;
Do[
Do[
vals[[n1+1,n2+1,n3+1,n4+1,n5+1,nsum-n1-n2-n3-n4-n5+1]]=
Total[Function[s,
pos=Reverse[Sort[{n1,n2,n3,n4,n5,nsum-n1-n2-n3-n4-n5}-s]];
If[Min[pos]<0,0,vals[[##]]&@@(pos+1)]
]/@steps],
{n1,0,Min[n,nsum]},
{n2,0,Min[n1,nsum-n1]},
{n3,0,Min[n2,nsum-n1-n2]},
{n4,0,Min[n3,nsum-n1-n2-n3]},
{n5,Ceiling[(nsum-n1-n2-n3-n4)/2],Min[n4,nsum-n1-n2-n3-n4]}];
,{nsum,2,6n,2}];
Table[vals[[i+1,i+1,i+1,i+1,i+1,i+1]],{i,n}]-Take[data,n]
]]

(* This variant is slightly slower, but uses a bit less memory. *)
Clear[Walks6D];
Walks6D[n_Integer,steps_List]:=Walks6D[n,steps]=
Module[{vals={{{{{1}}}}},seq={},n0,n1},
Do[
vals=Table[
n0=Ceiling[(nsum-n1-n2-n3-n4)/2];
Join[
Table[0,{n0}],
Table[Total[Function[s,
pos=Reverse[Sort[{n1,n2,n3,n4,n5,nsum-n1-n2-n3-n4-n5}-s]];
If[Min[pos]<0,0,vals[[##]]&@@Most[pos+1]]
]/@steps],
{n5,n0,Min[n4,nsum-n1-n2-n3-n4]}]],
{n1,0,Min[n,nsum]},
{n2,0,Min[n1,nsum-n1]},
{n3,0,Min[n2,nsum-n1-n2]},
{n4,0,Min[n3,nsum-n1-n2-n3]}];
If[IntegerQ[n1=nsum/6+1],AppendTo[seq,vals[[n1,n1,n1,n1,n1]]]];
,{nsum,2,6n,2}];
Return[seq];
];

Walks6D[33,Permutations[{1,1,0,0,0,0}]]-data

Timing[MaxMemoryUsed[Walks6D[44,Permutations[{1,1,0,0,0,0}]];]]

4.2  A172671

InitializeSeq[172671]

(* Simpler recurrence by hypergeometric transformation (due to E. Kaltofen). *)
Factor[DFiniteTimesHyper[ToOrePolynomial[rec,a[n]],(n!)^3/(3n)!]]

15-fold multinomial sum, reduced to 11-fold by Chu-Vandermonde

expr=Array[c,21]/.Solve[Array[c,21].Join[Permutations[{2,0,0,0,0,0}],Permutations[{1,1,0,0,0,0}]]==Table[n,{6}]][[1]];
vars=Union[Cases[expr,c[_],Infinity]];
TableForm[expr]

(* We can apply the Chu-Vandermonde identity. *)
ChuVandermonde=sum[1/k!/(m-k)!/(r-k)!/(n-r+k)!,{k,0,r}]==(m+n)!/r!/(m+n-r)!/m!/n!;

expr1=FunctionExpand[Multinomial@@expr]/.Gamma[n_]->(n-1)!;
expr2=Select[expr1,Not[FreeQ[#,c[17]]]&]

test=ChuVandermonde/.{k->c[17],r->n-2 c[3]-c[8]-c[12]-c[16],m->2 n-c[1]-c[2]-c[3]-c[4]-c[5]+c[6]-c[7]-c[8]-c[9]-c[10]-c[12]-c[13]-c[14]-c[16],n->c[1]+c[2]-c[3]-c[4]+c[5]-c[6]+c[7]+c[10]+c[14]-c[16]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]];

expr2=Select[expr1,Not[FreeQ[#,c[14]]]&]

test=ChuVandermonde/.{k->c[14],r->n-2 c[2]-c[7]-c[12]-c[13],m->2 n-c[1]-c[2]-c[3]-c[4]-c[5]+c[6]-c[7]-c[8]-c[9]-c[10]-c[12]-c[13]-c[16],n->n+c[1]-c[2]-c[3]-c[4]+c[5]-c[6]+c[10]-c[12]-c[13]-c[16]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]];

expr2=Select[expr1,Not[FreeQ[#,c[10]]]&]

test=ChuVandermonde/.{k->c[10],r->n-2 c[1]-c[7]-c[8]-c[9],m->2 n-c[1]-c[2]-c[3]-c[4]-c[5]+c[6]-c[7]-c[8]-c[9]-c[12]-c[13]-c[16],n->2 n-c[1]-c[2]-c[3]-c[4]+c[5]-c[6]-c[7]-c[8]-c[9]-c[12]-c[13]-c[16]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]];

expr2=Select[expr1,Not[FreeQ[#,c[6]]]&]

test=ChuVandermonde/.{k->c[6],r->2 n-c[1]-c[2]-c[3]-c[4]+c[5]-c[7]-c[8]-c[9]-c[12]-c[13]-c[16],m->-n+c[1]+c[2]+c[3]+c[4]-c[5]+c[7]+c[8]+c[9]+c[12]+c[13]+c[16],n->4 n-2 c[1]-2 c[2]-2 c[3]-2 c[4]-2 c[7]-2 c[8]-2 c[9]-2 c[12]-2 c[13]-2 c[16]}

expr2===test[[1,1]]

expr1=expr1/expr2*test[[2]];

expr1

Clear[A172671];
A172671[n_]:=A172671[n]=
Sum[((3 n)! (3 n-c[1]-c[2]-c[3]-c[4]-c[5]-c[7]-c[8]-c[9]-c[12]-c[13]-c[16])!)/(c[1]! c[2]! c[3]! c[4]! (n-2 c[5])! c[5]! c[7]! c[8]! (n-2 c[1]-c[7]-c[8]-c[9])! c[9]! c[12]! (n-2 c[2]-c[7]-c[12]-c[13])! c[13]! (n-2 c[3]-c[8]-c[12]-c[16])! (n-2 c[4]-c[9]-c[13]-c[16])! (2 n-c[1]-c[2]-c[3]-c[4]+c[5]-c[7]-c[8]-c[9]-c[12]-c[13]-c[16])! c[16]! (-n+c[1]+c[2]+c[3]+c[4]-c[5]+c[7]+c[8]+c[9]+c[12]+c[13]+c[16])!),
{c[1],0,n/2},
{c[5],0,n/2},
{c[2],0,n/2},
{c[3],0,n/2},
{c[4],0,n/2},
{c[7],0,Min[n-2c[1],n-2 c[2]]},
{c[8],0,Min[n-2 c[1]-c[7],n-2 c[3]]},
{c[9],0,Min[n-2 c[1]-c[7]-c[8],n-2 c[4]]},
{c[12],0,Min[n-2 c[2]-c[7],n-2 c[3]-c[8]]},
{c[13],0,Min[n-2 c[2]-c[7]-c[12],n-2 c[4]-c[9]]},
{c[16],Max[0,n-c[1]-c[2]-c[3]-c[4]+c[5]-c[7]-c[8]-c[9]-c[12]-c[13]],Min[n-2 c[3]-c[8]-c[12],n-2 c[4]-c[9]-c[13]]}]

Do[Print[n,": ",Timing[A172671[n]]],{n,8}]

Table[A172671[n],{n,8}]-Take[data,8]

6-dimensional walks

(* Using the same procedure Walks6D from A172572 *)
Timing[MaxMemoryUsed[
test=Walks6D[33,Join[Permutations[{1,1,0,0,0,0}],Permutations[{2,0,0,0,0,0}]]]
]]

test===data

Timing[MaxMemoryUsed[
test=Walks6D[40,Join[Permutations[{1,1,0,0,0,0}],Permutations[{2,0,0,0,0,0}]]]
]]

4.3  A188818

InitializeSeq[188818]

(* Naive transfer matrix implementation. *)
A188818[n_]:=
Module[{states,tmat},
states=Tuples[{0,1},n];
tmat=SparseArray[Outer[If[MemberQ[Rest[#1]-Most[#2],1]||MemberQ[Most[#1]-Rest[#2],1],0,1]&,states,states,1]];
vec=Table[1,{2^n}];
Do[vec=tmat.vec,{n-1}];
Return[Total[vec]];
];
Do[Print[Timing[A188818[n]]],{n,10}]

(* Enumerate all even solutions for n=5. *)
pos=Select[Tuples[Range[5],2],EvenQ[Total[#]]&];
sol=DeleteCases[If[#2≤#5≤#8&&#1≤#4≤#7≤#10≤#13&&#6≤#9≤#12&&#6≤#4≤#2&&#11≤#9≤#7≤#5≤#3&&#12≤#10≤#8,{##},{}]&@@@Tuples[{0,1},Length[pos]],{}];
(* This is in agreement with e_5=40. *)
Length[sol]

tab=Table["",{5},{5}];
MatrixForm[Reverse[Transpose[ReplacePart[tab,MapThread[Rule,{pos,Range[13]}]]]]]

(* All even solutions for n=5. *)
MatrixForm[Reverse[Transpose[ReplacePart[tab,MapThread[Rule,{pos,#}]]]]]&/@sol

Lattice Walks

(* Krattenthaler: Catalan paths (a,b) -> (c,d) with x≥y. *)
myL[{a_,b_},{c_,d_}]:=Binomial[c+d-a-b,c-a]-Binomial[c+d-a-b,c-b+1];
(* myD: Dyck paths from (1,y1) -> (n,y2) with y≥1, and where 1+y1 = n+y2 (mod 2). *)
(* If y1 is odd, we map (x,y) to ((x+y-2)/2,(x-y)/2) and apply myL: *)
Simplify[myL@@({(#1+#2-2)/2,(#1-#2)/2}&@@@{{x1,y1},{x2,y2}})]
(* If y1 is even, we map (x,y) to ((x+y-1)/2,(x-y+1)/2) and apply myL: *)
Simplify[myL@@({(#1+#2-1)/2,(#1-#2+1)/2}&@@@{{x1,y1},{x2,y2}})]
(* Since both expressions are equal, we define: *)
myD[y1_,y2_,n_]:=Binomial[n-1,(n-y1+y2-1)/2]-Binomial[n-1,(n+y1+y2-1)/2];

(* We start with even paths (on positions (x,y) with x+y even). *)

(* Sum over all Dyck paths with y1=1,3,...,n+2, y2=1,3,...,n+2 (if n odd)
or y2=2,4,...,n+2 (if n even). *)
(* No path can touch both lower and upper boundary. *)
(* Mirror those whose start and end points are s.t. it can potentially touch the upper boundary. *)
Table[Sum[If[y1+y2≤n+1,myD[y1,y2,n],myD[n+3-y1,n+3-y2,n]],{y1,1,n+2,2},{y2,n+2,1,-2}],{n,15}]

(* Combine those paths which appear twice. *)
Table[Sum[If[y1+y2≤n+1,2,1]*myD[y1,y2,n],{y1,1,n+2,2},{y2,n+3-y1,1,-2}],{n,15}]

(* Substitute y1=2k+1, y2=n+2-2l. *)
Table[
Sum[If[l≥k+1,2,1]*myD[2k+1,n-2l+2,n],{k,0,Floor[(n+1)/2]},{l,k,Floor[(n+1)/2]}],{n,15}]

(* Eliminate the case distinction by separating the case l=k. *)
Table[
Sum[2*Sum[myD[2k+1,n-2l+2,n],{l,k+1,Floor[(n+1)/2]}]+myD[2k+1,n-2k+2,n],{k,0,Floor[(n+1)/2]}],{n,15}]

(* Insert the definition in terms of binomial coefficients. *)
Table[
Sum[2*Sum[Binomial[n-1,n-k-l]-Binomial[n-1,n+k-l+1],{l,k+1,Floor[(n+1)/2]}]+Binomial[n-1,n-2 k],{k,0,Floor[(n+1)/2]}],{n,15}]

(* Simplify: Sum[Binomial[n-1,n-2 k], {k,0,(n+1)/2}] = 2^(n-2) for n≥2. *)
Table[
2*Sum[Binomial[n-1,n-k-l]-Binomial[n-1,n+k-l+1],{k,0,Floor[(n+1)/2]},{l,k+1,Floor[(n+1)/2]}]+2^(n-2),{n,2,15}]

(* Now the odd paths (on positions (x,y) with x+y odd). *)

(* Sum over all Dyck paths with y1=2,4,...,n+1, y2=2,4,...,n+1 (if n odd)
or y1=2,4,...,n+2, y2=1,3,...,n+1 (if n even). *)
Table[Sum[If[y1+y2≤n+1,myD[y1,y2,n],myD[n+3-y1,n+3-y2,n]],{y1,2,n+2,2},{y2,n+1,1,-2}],{n,15}]

(* Combine those paths which appear twice. *)
Table[Sum[If[y1+y2≤n+1,2,1]*myD[y1,y2,n],{y1,2,n+2,2},{y2,n+3-y1,1,-2}],{n,15}]

(* Substitute y1=2k+2, y2=n+1-2l. *)
Table[Sum[If[l≥k+1,2,1]*myD[2k+2,n-2l+1,n],{k,0,Floor[n/2]},{l,k,Floor[n/2]}],{n,15}]

(* Eliminate the case distinction by separating the case l=k. *)
Table[Sum[2*Sum[myD[2k+2,n-2l+1,n],{l,k+1,Floor[n/2]}]+myD[2k+2,n-2k+1,n],{k,0,Floor[n/2]}],{n,15}]

(* Insert the definition in terms of binomial coefficients. *)
Table[Sum[2*Sum[Binomial[n-1,n-k-l-1]-Binomial[n-1,n+k-l+1],{l,k+1,Floor[n/2]}]+Binomial[n-1,n-2 k-1],{k,0,Floor[n/2]}],{n,15}]

(* Simplify: Sum[Binomial[n-1,n-2 k-1], {k,0,n/2}] = 2^(n-2) for n≥2. *)
Table[2*Sum[Binomial[n-1,n-k-l-1]-Binomial[n-1,n+k-l+1],{k,0,Floor[n/2]},{l,k+1,Floor[n/2]}]+2^(n-2),{n,2,15}]

(* Putting everything together. *)
Table[
(2^(n-2)+2*Sum[Binomial[n-1,n-k-l]-Binomial[n-1,n+k-l+1],{k,0,Floor[(n+1)/2]},{l,k+1,Floor[(n+1)/2]}])*
(2^(n-2)+2*Sum[Binomial[n-1,n-k-l-1]-Binomial[n-1,n+k-l+1],{k,0,Floor[n/2]},{l,k+1,Floor[n/2]}]),
{n,2,32}]===Rest[data]

Creative telescoping  for (Subscript[o, n])

Sum[Binomial[n-1,n-k-l-1]-Binomial[n-1,n+k-l+1],{k,0,Floor[n/2]},{l,k+1,Floor[n/2]}]

(* n even *)
Sum[Binomial[2m-1,2m-k-l-1]-Binomial[2m-1,2m+k-l+1],{k,0,m},{l,k+1,m}]

Timing[annoe1=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m-1,2m-k-l-1],{l,k+1,m}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

Timing[annoe2=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m-1,2m+k-l+1],{l,k+1,m}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

annoe=DFinitePlus[annoe1,annoe2]

(* n odd, n = 2m+1 *)
Sum[Binomial[2m,2m-k-l]-Binomial[2m,2m+k-l+2],{k,0,m},{l,k+1,m}]

Timing[annoo1=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m,2m-k-l],{l,k+1,m}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

Timing[annoo2=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m,2m+k-l+2],{l,k+1,m}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

annoo=DFinitePlus[annoo1,annoo2]

anno=DFinitePlus[DFiniteSubstitute[annoe,{m->n/2},Algebra->OreAlgebra[S[n]]],DFiniteSubstitute[annoo,{m->(n-1)/2},Algebra->OreAlgebra[S[n]]]];
#[anno]&/@{ByteCount,Support}

anno

annog=NormalizeCoefficients/@ToOrePolynomial[GuessUnivRE[Table[Sum[Binomial[n-1,n-k-l-1]-Binomial[n-1,n+k-l+1],{k,0,Floor[n/2]},{l,k+1,Floor[n/2]}],{n,0,60}],f[n],Order->3,Degree->4],f[n]]

OreReduce[anno[[1]],annog]

Creative telescoping  for (Subscript[e, n])

Sum[Binomial[n-1,n-k-l]-Binomial[n-1,n+k-l+1],{k,0,Floor[(n+1)/2]},{l,k+1,Floor[(n+1)/2]}]

(* n even, n = 2m *)
Sum[Binomial[2m-1,2m-k-l]-Binomial[2m-1,2m+k-l+1],{k,0,m},{l,k+1,m}]

Timing[annee1=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m-1,2m-k-l],{l,k+1,m}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

Timing[annee2=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m-1,2m+k-l+1],{l,k+1,m}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

annee=DFinitePlus[annee1,annee2]

(* n odd, n = 2m+1 *)
Sum[Binomial[2m,2m-k-l+1]-Binomial[2m,2m+k-l+2],{k,0,m+1},{l,k+1,m+1}]

Timing[anneo1=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m,2m-k-l+1],{l,k+1,m+1}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

Timing[anneo2=NormalizeCoefficients/@(First[FindCreativeTelescoping[DFiniteTimes[Annihilator[Sum[Binomial[2m,2m+k-l+2],{l,k+1,m+1}],{S[k],S[m]}],Annihilator[Gamma[k+1+e]/Gamma[k+1],{S[k],S[m]}]],S[k]-1]]/.e->0)]

anneo=DFinitePlus[anneo1,anneo2]

anne=DFinitePlus[DFiniteSubstitute[annee,{m->n/2},Algebra->OreAlgebra[S[n]]],DFiniteSubstitute[anneo,{m->(n-1)/2},Algebra->OreAlgebra[S[n]]]];
#[anne]&/@{ByteCount,Support}

anne

anneg=NormalizeCoefficients/@ToOrePolynomial[GuessUnivRE[Table[Sum[Binomial[n-1,n-k-l]-Binomial[n-1,n+k-l+1],{k,0,Floor[(n+1)/2]},{l,k+1,Floor[(n+1)/2]}],{n,0,60}],f[n],Order->3,Degree->4],f[n]]

OreReduce[anne[[1]],anneg]

Final recurrence

GBEqual[anno,DFinitePlus[anno,Annihilator[2^(n-2),S[n]]]]

GBEqual[anne,DFinitePlus[anne,Annihilator[2^(n-2),S[n]]]]

Timing[ann=DFiniteTimes[anno,anne];]

#[ann]&/@{ByteCount,Support}

OreReduce[ann[[1]],ToOrePolynomial[rec,a[n]]]

4.4  A306322

InitializeSeq[306322]

(* Nara[n,k] is the Narayana number N(n+k-1,k). *)
Nara[i_,j_]:=1/(i+j-1)*Binomial[i+j-1,i]*Binomial[i+j-1,i-1];
Table[Nara[n,n],{n,10}]

(* First expression. *)
Table[2*(Sum[Nara[i,j]-Binomial[i+j-2,i-1],{i,n},{j,n}]-Binomial[2n,n]+1+Sum[Nara[j-i,n],{i,n-1},{j,i+1,n-1}])-(Nara[n,n]-2*Binomial[2n,n]+1),{n,40}]===Rest[data]

(* Simplified expression. *)
Table[2*Sum[Sum[Nara[i,j],{i,1,n}]+(n-j-1)*Nara[j,n],{j,1,n}]-2*Binomial[2n,n]+Nara[n,n]+3,{n,40}]===Rest[data]

ann1=Annihilator[Sum[Nara[i,j],{i,1,n},{j,1,n}],S[n]];
Exponent[ann1[[1]],{S[n],n}]

ann2=Annihilator[Sum[(n-j-1)*Nara[j,n],{j,1,n}],S[n]]

ann3=Annihilator[-2*Binomial[2n,n]+Nara[n,n]+3,S[n]]

ann=DFinitePlus[ann1,ann2,ann3];
Exponent[ann[[1]],{S[n],n}]

OreReduce[ann,ToOrePolynomial[{rec},a[n]]]

5  Further Examples

5.1  A195806

InitializeSeq[195806]

(* Also a C-finite recurrence can be found. *)
crec=GuessMinRE[data,a[n]]

(* It is compatible with the guessed minimal-order recurrence. *)
OreReduce@@(ToOrePolynomial[#,a[n]]&/@{crec,{rec}})

(* Compute more terms with the guessed recurrence. *)
re2l=RE2L[Prepend[Table[a[i]==data[[i]],{i,4}],rec==0],a[n],{1,1000}];
{Take[re2l,Length[data]]===data,And@@(IntegerQ/@re2l)}

(* Interpolate the constituents of the quasi-polynomial. *)
Table[Expand[InterpolatingPolynomial[Table[{n,re2l[[n]]},{n,s,200,6}],n]],{s,6}]

(* Sanity check. *)
A195806[n_]:=Switch[Mod[n,6],
1,(2051+5912 n+5688 n^2+3900 n^3+2275 n^4+780 n^5+130 n^6)/1296,
2,(1904+5464 n+5592 n^2+3900 n^3+2275 n^4+780 n^5+130 n^6)/1296,
3,(1539+5400 n+5688 n^2+3900 n^3+2275 n^4+780 n^5+130 n^6)/1296,
4,(1808+5912 n+5688 n^2+3900 n^3+2275 n^4+780 n^5+130 n^6)/1296,
5,((1+n)^2 (2147+1170 n+1105 n^2+520 n^3+130 n^4))/1296,
0,(1296+5400 n+5688 n^2+3900 n^3+2275 n^4+780 n^5+130 n^6)/1296];
Table[A195806[n],{n,1000}]===re2l

(* Different representation *)
(a[6n+#1]==Factor[#2/.n->6n+#1])&@@@Partition[Rest[List@@DownValues[A195806][[1,2]]],2]

%/.n->1

TableForm[arr=Table[Symbol["a"<>ToString[i]<>ToString[j]],{i,5},{j,i}]/.{a11->0,a51->0,a55->0}]

TableForm[eqns=Table[Total[arr[[i]]]==Sum[arr[[5-i+j,6-i]],{j,i}]==Sum[arr[[5-i+j,j]],{j,i}],{i,2,5}]]

ineq=Join[Thread[0≤#],Thread[#≤n]]&[Variables[List@@@eqns]]/.First[Solve[eqns]]

(* Such systems of linear inequalities can (in principle) be treated by the Omega package. *)
<<RISC`Omega`

Omega Package version 2.49
written by Axel Riese
(in cooperation with George E. Andrews and Peter Paule)
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

CrudeGF=OSum[q^n*z^(a21+a22+a31+a32+a33+a41),ineq,λ]

Timing[OR[CrudeGF]]

Mathematica crashed during the elimination of the 7th slack variable.

OR[CrudeGF[[1]],(Subscript[λ, 12])]

5.1  A216940

InitializeSeq[216940]

(* Compute more terms using the guessed recurrence relation. *)
re2l=RE2L[{a[1]==data[[1]],rec==0},a[n],{1,1000}];
{Take[re2l,Length[data]]===data,And@@(IntegerQ/@re2l)}

(* It seems that the n-th sequence term is given by a degree-37 polynomial. *)
Factor[InterpolatingPolynomial[Take[re2l,100],n]]

(* The hexagonal array with side length 4. *)
TableForm[hex=Table[If[Abs[3-i]≤j≤12-Abs[3-i]&&OddQ[i+j],Symbol["a"<>ToString[i]<>ToString[(j-Abs[3-i])/2]],""],{i,0,6},{j,0,12}]]

(* We only consider E and SW because SE is implied by these *)
ineq=Select[Flatten[Join[
Table[hex[[i,j]]≤hex[[i,j+2]],{i,7},{j,11}],
Table[hex[[i,j+1]]≤hex[[i+1,j]],{i,6},{j,12}]
]],FreeQ[#,""]&]

Length[%]

(* Again, one can try to solve this system of inequalities with Omega. *)
<<RISC`Omega`

Omega Package version 2.49
written by Axel Riese
(in cooperation with George E. Andrews and Peter Paule)
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

CrudeGF=OSum[q^a63*z^Total[Most[Variables[List@@@ineq]]],ineq,λ]

Timing[OR[CrudeGF]]

(* Unfortunately, this computation didn't finish in reasonable time. *)

5.2  A194478

InitializeSeq[194478]

Binomial[n+1,2]//FunctionExpand

(* The closed form. *)
TraditionalForm[cf=HoldForm[Binomial[Binomial[n+1,2],6]-(3 Sum[Binomial[i,j]Binomial[x-i,6-j],{j,3,6},{i,1,n}]/.x->n(n+1)/2)+(3 Sum[Binomial[i,3]Binomial[j,3],{i,1,n},{j,1,i-1}]+3 Sum[Binomial[i,3]Binomial[j,3],{i,1,n},{j,1,n-i}]+3 Sum[Binomial[i-1,2]Binomial[j-1,2](Binomial[n+1,2]-i-j+1)+Binomial[i-1,3]Binomial[j-1,2]+Binomial[i-1,2]Binomial[j-1,3],{i,1,n},{j,n-i+1,n}]+3 Sum[Binomial[i-1,3]Binomial[j-1,3],{i,1,n},{j,n-i+1,n}])- Sum[(i-2)(j-2)(Sum[l-2,{l,n-Min[i,j]+1,2n-(i+j)}]+Sum[l-2,{l,2n+2-(i+j),n}]),{i,3,n},{j,n-i+1,n}]]]

Expand[FunctionExpand[ReleaseHold[cf]]]

Simplify[%,Assumptions->n>5]/.Floor[(n+1)/2]->(n+(1-(-1)^n)/2)/2

Expand[%]/.(-1)^a_:>(-1)^(a/.n->0)*(-1)^(Mod[Coefficient[a,n],2]*n)/.a_*(-1)^(n+1):>(-a)*(-1)^n

Collect[%,(-1)^n,Factor]

(* Sanity check *)
Table[%,{n,31}]===data

6  Conjectures

6.1  A215570

InitializeSeq[215570]

(* Simpler recurrence for transformed sequence. *)
Factor[DFiniteTimesHyper[ToOrePolynomial[rec,a[n]],(n!)^3*((n+1)!)^2/(5n)!]]

(* Code from OEIS, slightly simplified. *)
(* Computed 51 terms in about 2.5h on radon1 (memory usage 60GB). *)
Clear[mya,myb];
Timing[
myb[l_]:=myb[l]=
Module[{m=Length[l],n=Total[l],g},
g=l.Range[m]-(m+1)(n-1)/2;If[n<2,1,Sum[If[l[[i]]>0,myb[ReplacePart[l,i->l[[i]]-1]],0],{i,1,Min[g,m]}]]];
mya[k_]:=myb[Array[k&,5]];
Table[mya[n],{n,0,12}]
]

(* Code from OEIS, a bit more simplified. *)
Clear[mya,myb];
Timing[
myb[l_]:=myb[l]=
With[{n=Total[l]},
If[n<2,1,Sum[If[l[[i]]>0,myb[ReplacePart[l,i->l[[i]]-1]],0],{i,1,Min[l.Range[5]-3*(n-1),5]}]]];
mya[k_]:=myb[Array[k&,5]];
Table[mya[n],{n,0,12}]
]

Length[DownValues[myb]]

(* Transfer matrix method is slower, but more memory-efficient. *)
(* Computed 51 terms in 21h (memory about 2 GB) *)
Timing[With[{n=12},
seq={1};
tmat=Table[If[Abs[j-i]≤2,x[j-i+3],0],{i,3n+1},{j,3n+1}]/.x[5]->1;
vec=Prepend[Table[0,{3n}],1];
Do[
vec=PolynomialMod[Expand[tmat.vec],Array[x,4]^(n+1)];
If[Mod[i,5]===0,AppendTo[seq,Coefficient[Total[vec],(Times@@Array[x,4])^(i/5)]]];
,{i,5n}];
seq
]]

(* The additional three terms that we found. *)
data=Join[data,{109077635701149385673307834354313488656611420675978722550093654767156150247798646452158745836073766066836510193294858821035415499345237177075859552482287257000000,320518781996956464254305073090595172642716936717597758539843844836406214884389287005904094666281078845230748750512417310938341310184558733989330751809185813089240000,942993553387261719839432368142529421305313544889687234217568261481733850928388393453200584613854491160136599571045529542095321859090899240127873178680468809691523430400}];

(* Test: they also satisfy the guessed recurrence. *)
Table[rec,{n,44,47}]/.a[n_]:>data[[n+1]]

6.2  A339987

InitializeSeq[339987]

(* Simpler recurrence for transformed sequence. *)
Factor[DFiniteTimesHyper[ToOrePolynomial[rec,a[n]],(n+1)/Pochhammer[5/2,n-2]]]

6.3  A269021

InitializeSeq[269021]

(* Simpler recurrence for transformed sequence. *)
Factor[DFiniteTimesHyper[ToOrePolynomial[rec,a[n]],1/((2n)!)^2]]

6.4  A181198

InitializeSeq[181198]

(* Compute more terms of the sequence. For each n, timing and number of variables is displayed. *)
NextPartitions[n1_,n2_,n3_,n4_]:=
If[n1<n,f[n1+1,n2,n3,n4],0]+
If[n2<n1-1||n2===n-1,f[n1,n2+1,n3,n4],0]+
If[n3<n2-1||n3===n-1===n2-1,f[n1,n2,n3+1,n4],0]+
If[n4<n3-1,f[n1,n2,n3,n4+1],0];
Timing[
seq181198=Table[
pp=f[1,0,0,0];
Print[n,": ",Timing[Max[Table[pp=Expand[pp/.f[ns__]:>NextPartitions[ns]];Length[pp],{4n-2}]]]];
pp/.f[n,n,n,n-1]->1,{n,100}];
Take[seq181198,10]
]

(* The additional values also satisfy the guessed recurrence. *)
Union[Table[rec,{n,98}]/.a[n_]:>seq181198[[n]]]

(* Mathematica code from OEIS A181196 (bivariate sequence T(n,k), our A181198 is T(4,k)). *)
Clear[myb];
myb[l_List]:=myb[l]=
With[{n=Length[l]},If[Union[l]=={0},1,Sum[If[(i==1||l[[i-1]]≤l[[i]])&&l[[i]]>If[i==n,0,l[[i+1]]],myb[ReplacePart[l,i->l[[i]]-1]],0],{i,1,n}]]];myT[n_,k_]:=myb[Array[n&,k]];Timing[Table[myT[4,n],{n,30}]===Take[seq181198,30]]

Length[DownValues[myb]]

(* Conjectured closed form. *)
Table[(-64)^n*(n-1)*Pochhammer[-1/2,2n]*Pochhammer[1/2,n]/(4(3n)!)*(-1+3*Sum[(-4)^k (7 k^2-1)/((k-1) k (k+1)^2 (2 k-1)^2 (2 k+1)^3)*Binomial[3k,2k]*Binomial[k+1/2,k],{k,2,n-1}]),{n,2,100}]===Rest[seq181198]

6.4  A181199

InitializeSeq[181199]

(* Compute more terms of the sequence. For each n, timing and number of variables is displayed. *)
(* Computing the first 100 terms on radon1-grantley took 49660s. *)
NextPartitions[n1_,n2_,n3_,n4_,n5_]:=
If[n1<n,f[n1+1,n2,n3,n4,n5],0]+
If[n2<n1-1||n2===n-1,f[n1,n2+1,n3,n4,n5],0]+
If[n3<n2-1||n3===n-1===n2-1,f[n1,n2,n3+1,n4,n5],0]+
If[n4<n3-1||n4===n-1===n3-1,f[n1,n2,n3,n4+1,n5],0]+
If[n5<n4-1,f[n1,n2,n3,n4,n5+1],0];
Timing[
seq181199=Table[
pp=f[1,0,0,0,0];
Print[n,": ",Timing[Max[Table[pp=Expand[pp/.f[ns__]:>NextPartitions[ns]];Length[pp],{5n-2}]]]];
pp/.f[n,n,n,n,n-1]->1,{n,60}];
Take[seq181199,10]
]

(* Conjectured closed form. *)
Table[
1-27/4*Sum[(-1)^k*(25216 k^8+9888 k^7-14496 k^6+11208 k^5+23832 k^4+7383 k^3-1522 k^2-939 k-90)/((k+1) (2 k-1) (2 k+1) (3 k+1) (3 k+2) (3 k+3) (4 k-1) (4 k+1) (4 k+3))*(5k)!/(3k)!/(k!)^2*Sum[(-1)^i*((3 i+1) (3 i+2) (4 i+3) (137855872 i^11+860969696 i^10+2047036856 i^9+2032587274 i^8-24192441 i^7-1894061166 i^6-1671661480 i^5-524330624 i^4+36004789 i^3+62751860 i^2+13865604 i+927360))/((i+1)^2 (i+2)^2 (2 i-1) (2 i+1) (2 i+3) (25216 i^8+9888 i^7-14496 i^6+11208 i^5+23832 i^4+7383 i^3-1522 i^2-939 i-90) (25216 i^8+211616 i^7+760768 i^6+1543976 i^5+1973632 i^4+1683047 i^3+971955 i^2+353502 i+60480))*(3i)!/(i!)^3,{i,1,k-1}],{k,1,n-1}],{n,1,60}]===seq181199

6.5  A181280

InitializeSeq[181280]

(* Compute more terms using the guessed recurrence. *)
re2l=RE2L[Prepend[Table[a[i]==data[[i]],{i,10}],rec==0],a[n],{1,1000}];
{Take[re2l,Length[data]]===data,And@@(IntegerQ/@re2l)}

(* Find the minimal-order recurrence. *)
mrec=GuessUnivRE[Take[re2l,{4,51}],a[n],Order->6,Degree->5,StartPoint->4];
mrec=Collect[Numerator[Together[mrec[[1]]]],a[_],Expand]

(* The originally guessed recurrence is a left multiple of the minimal-order one. *)
OreReduce@@(ToOrePolynomial[#,a[n]]&/@{rec,{mrec}})

(* Conjectured closed form. *)
Table[(1/3)*2^(2*n-11)*(6*n^2-219*n+820)-(1/9)*2^(n-5)*(3*n+32)-(113/3)*(-1)^n*2^(3*n-14)+2^(4*n-9)-(1/3)*(-1)^n*2^(2*n-11)*(13*n-164)+(1/9)*2^(3*n-14)*(288*n-3473),{n,4,1000}]===Drop[re2l,3]

6.6  A253217

InitializeSeq[253217]

(* Conjectured differential equation. *)
deq=-2 (-144-9472 x+59516 (x^2)+293099 (x^3)-6118786 (x^4)+20605760 (x^5)+8447840 (x^6)-8859568 (x^7)+576032 (x^8)-3834112 (x^9)-672768 (x^10)+43008 (x^11)) f[x]-2 x (136+8840 x-56698 (x^2)-266263 (x^3)+5562098 (x^4)-18100720 (x^5)-9643156 (x^6)-2227000 (x^7)+37415840 (x^8)-41803520 (x^9)-10291200 (x^10)+1118208 (x^11)) f'[x]-(x^2) (-128-8208 x+50520 (x^2)+156939 (x^3)-5838066 (x^4)+25114104 (x^5)+8254632 (x^6)-85240800 (x^7)+269354016 (x^8)-71598976 (x^9)-35864064 (x^10)+4171776 (x^11)) f''[x]-(x^3) (40+3032 x-13122 (x^2)-215313 (x^3)+1094157 (x^4)-3361026 (x^5)+3620772 (x^6)-71338632 (x^7)+134615328 (x^8)-4944256 (x^9)-15910400 (x^10)+1720320 (x^11)) Derivative[3][f][x]-(-4+x) (-1+x) (x^4) (1+2 x) (-1+4 x) (-1+16 x) (-4-202 x+257 (x^2)+4672 (x^3)-26492 (x^4)-7520 (x^5)+1344 (x^6)) Derivative[4][f][x];

(* Sanity check. *)
PolynomialMod[deq/.f->Function@@{x,Sum[data[[n]]*x^n,{n,37}]},x^38]

(* Factorization of the differential operator, computed by Maple. *)
{L1,L2,L3}=Factor[ToOrePolynomial[{(-172032*x^15+1790464*x^14-1744000*x^13-14355936*x^12+14262792*x^11+2806428*x^10-3103506*x^9+490737*x^8+40371*x^7-15834*x^6+500*x^5+16*x^4)**(Der[x]^2+((7225344*x^11-68551168*x^10-38001536*x^9+655096032*x^8-257000136*x^7-120067572*x^6+26515998*x^5+181023*x^4-1528359*x^3-207288*x^2+53464*x+3248)/(7*x+2)/(172032*x^11-1790464*x^10+1744000*x^9+14355936*x^8-14262792*x^7-2806428*x^6+3103506*x^5-490737*x^4-40371*x^3+15834*x^2-500*x-16))**Der[x]+(52684800*x^12-473285120*x^11-1237651072*x^10+3741714912*x^9+709591920*x^8-282597864*x^7+133579476*x^6+114181086*x^5-23095209*x^4-8060712*x^3+114296*x^2+1168*x+384)/x/(7*x+2)^2/(172032*x^11-1790464*x^10+1744000*x^9+14355936*x^8-14262792*x^7-2806428*x^6+3103506*x^5-490737*x^4-40371*x^3+15834*x^2-500*x-16)),Der[x]+(168*x^5-356*x^4-133*x^3+18*x^2+56*x+4)/(7*x+2)/(x-4)/x/(x-1)/(2*x+1)/(4*x-1),Der[x]+(8*x^3+3*x-2)/(x-1)/x/(2*x+1)/(4*x-1)}]]

L1**L2**L3

(* Check correctness of factorization. *)
Together[ApplyOreOperator[%,f[x]]/deq]

6.7  A098926

InitializeSeq["098926"]

re2l=RE2L[Prepend[Table[f[i-1]==data[[i]],{i,LeadingExponent[op][[1]]}],ApplyOreOperator[op,f[n]]==0],f[n],1000];
And@@(IntegerQ/@re2l)

(* Generating function: ODE and closed form. *)
f1x=(x^2-x-2)/(x(x-1))*Exp[(x+1)/(x(x-1))];
f2y=(y^5-3y^4+2y^3-2y^2-y+1)/(y(y+1)^4(y-2)^2)*Exp[-(2y^2+2)/(y(y-1)(y+1))];
f3z=(z^2(z-2)(z^8-2z^7-12z^6+28z^5-10z^4-22z^3+4z^2+4z+1))/((z-1)^2(z^5-3z^4+2z^3-2z^2-z+1)^2)*Exp[(z-1)/(z(z+1))];
ode=
DFiniteTimes[
Annihilator[f1x,Der[x]],
{DFiniteTimes[
Annihilator[f2y,Der[y]],
{Annihilator[f3z,Der[z]][[1]]**Der[z]/.z->y}
][[1]]**Der[y]/.y->x}
][[1]];
Factor[ode]

(* Sanity check. *)
PolynomialMod[ApplyOreOperator[ode,Sum[data[[n]]*x^(n+2),{n,34}]],x^36]

6.8  A164735

InitializeSeq[164735]

(* Recurrence is only valid for n ≥ 3. *)
Table[rec,{n,3,58}]/.a[n_]:>data[[n]]

(* Compute more terms, using the guessed recurrence (n=1,2 have to be omitted). *)
re2l=RE2L[Prepend[Table[a[i]==data[[i]],{i,3,14}],rec==0],a[n],{1,1000}];
{Take[re2l,{3,Length[data]}]===Drop[data,2],And@@(IntegerQ/@re2l)}

(* With period 18, we get first order recurrence, with polynomial solutions. *)
Collect[GuessMinRE[Take[re2l,{8,400,18}],a[n]],a[_],Factor]

RSolve[%==0,a[n],n]

(* Compute all 18 polynomials by interpolation. *)
pols=Table[Factor[40*InterpolatingPolynomial[Take[re2l,{s,800,18}],k]/.k->k+1],{s,18}];
pols=Prepend[Most[pols],Factor[Last[pols]/.k->k-1]]

(* Sanity check. *)
Table[{kk,r}=QuotientRemainder[n,18];pols[[r+1]]/40/.k->kk,{n,3,70}]===Drop[data,2]

(* Define the Kaprekar map. *)
Kaprekar[n_Integer]:=#1-#2&@@(FromDigits/@({Reverse[#],#}&[Sort[IntegerDigits[n]]]));
NestList[Kaprekar,86526432,3]

(* There are (at least) two different patterns that constitute Kaprekar 3-cycles *)
Clear[myZ1,myZ2]
myZ1[k_/;k≥0,x1_/;x1≥0,x3_/;x3≥0,x5_/;x5≥1,x7_/;x7≥1,x9_/;x9≥1]:=
FromDigits[Flatten[{
MapThread[Table[#1,{#2}]&,{Reverse[Range[9]],Riffle[{x9,x7,x5,x3,x1},k]}],0,
MapThread[Table[#1,{#2}]&,{Reverse[Range[0,10]],Riffle[{0,x1+1,x3,x5,x7,x9-1},k]}],1
}]];
myZ2[x1_/;x1≥0,x3_/;x3≥1,x5_/;x5≥0]:=
FromDigits[Flatten[{
6,Table[5,{x5}],4,Table[3,{x3}],Table[1,{x1}],
0,8,Table[8,{x1}],Table[6,{x3}],5,Table[4,{x5}],4}]];
myZ1[3,1,0,1,1,1]

(* Create all 3-cycles according to our conjecture (that there are no other patterns). *)
Clear[My3Cycles];
My3Cycles[n_]:=
Module[{s},
res=Flatten[Table[
s=(n-9k-2)/2;
Table[
If[x1+1==x5==x7,{},myZ1[k,x1,x3,x5,x7,s-x1-x3-x5-x7]]
,{x1,0,s-3}
,{x3,0,s-3-x1}
,{x5,1,s-2-x1-x3}
,{x7,1,s-1-x1-x3-x5}]
,{k,If[OddQ[n],1,0],Floor[(n-8)/9],2}]];
If[EvenQ[n],
res=Join[res,
Flatten[Table[myZ2[x1,x3,n/2-3-x1-x3],{x1,0,n/2-3},{x3,1,n/2-3-x1}]]];
];
Return[Union[Sort[NestList[Kaprekar,#,2]]&/@res]];
];
My3Cycles[12]

Table[Length[My3Cycles[n]],{n,70}]===data

(* Now count these 3-cycles. *)
Table[
s1=Sum[1,{k,If[OddQ[n],1,0],Floor[(n-8)/9],2}
,{x1,0,(n-9k-8)/2}
,{x3,0,(n-9k-8)/2-x1}
,{x5,1,(n-9k-6)/2-x1-x3}
,{x7,1,(n-9k-4)/2-x1-x3-x5}];
s2=Sum[1,{k,If[OddQ[n],1,0],Floor[(n-8)/9],2}
,{x1,0,(n-9k-8)/2}
,{x3,0,(n-9k-8)/2-3x1}];
s3=Sum[1,{x1,0,n/2-3},{x3,1,n/2-3-x1}];
1/3*(s1-s2)+(1-Mod[n,2])*s3,{n,70}]===data

(* Derive closed-form expressions for even n (mod 18). *)
expr=1/3*(sum[1,{k,0,Floor[(n-8)/18]}
,{x1,0,(n-18k-8)/2}
,{x3,0,(n-18k-8)/2-x1}
,{x5,1,(n-18k-6)/2-x1-x3}
,{x7,1,(n-18k-4)/2-x1-x3-x5}]
-sum[1,{k,0,Floor[(n-8)/18]}
,{x1,0,Floor[(n-18k-8)/6]}
,{x3,0,(n-18k-8)/2-3x1}])+
sum[1,{x1,0,n/2-3},{x3,1,n/2-3-x1}];
expr=Table[expr/.n->18m+2l,{l,0,8}]/.e_Floor:>FullSimplify[e,Element[m,Integers]]/.sum->Sum;
Factor[40*expr/.m->k]//TableForm

(* They agree with the interpolated polynomials from above. *)
%===Take[pols,{1,18,2}]

(* Derive closed-form expressions for odd n (mod 18). *)
expr=1/3*(sum[1,{k,0,Floor[(n-17)/18]}
,{x1,0,(n-18k-17)/2}
,{x3,0,(n-18k-17)/2-x1}
,{x5,1,(n-18k-15)/2-x1-x3}
,{x7,1,(n-18k-13)/2-x1-x3-x5}]-
sum[1,{k,0,Floor[(n-17)/18]}
,{x1,0,Floor[(n-18k-17)/6]}
,{x3,0,(n-18k-17)/2-3x1}]);
expr=Table[expr/.n->18m+2l+1,{l,0,8}]/.e_Floor:>FullSimplify[e,Element[m,Integers]]/.sum->Sum;
Factor[40*expr/.m->k]//TableForm

(* They agree with the interpolated polynomials from above. *)
%===Take[pols,{2,18,2}]