Automatic congruences for diagonals of rational functions

computation notebook

Eric Rowland
https://ericrowland.github.io/papers.html#Automatic_congruences_for_diagonals_of_rational_functions

Introduction

This notebook contains the computations performed in the paper “Automatic congruences for diagonals of rational functions” by Eric Rowland and Reem Yassawi.

To re-run the computations, first you will need to load the IntegerSequences package.  You can load it directly from the web by evaluating the following cell.

Import["https://ericrowland.github.io/packages/IntegerSequences.m"]

You can also download it from https://ericrowland.github.io/packages.html#IntegerSequences and then evaluate the following. (If you need help, see loading a package.)

<<IntegerSequences`

To check that it loaded correctly, re-evaluate the following.

MotzkinNumber[10]

Computing automata for algebraic sequences modulo (p^α)

This function determines which residues modulo a prime power are never output by a finite automaton when fed the base-p digits of an integer:

MissingResidues[sequence:_AutomaticSequence[_],modulus_Integer]:=
With[{automaton=AutomaticSequenceAutomaton[sequence]},
Complement[Range[0,modulus-1],Replace[Cases[First[automaton],{_->state_,Except[0]}:>state], AutomatonOutputRules[automaton],{1}]]
]

version note

The computations in the paper were carried out for algebraic sequences with initial term 0.  For example:

Table[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],{n,0,10}]

AutomatonGraph[AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],2],n,Method->"Diagonal",MinimizeComplexity->False]],GraphLayout->"SpringEmbedding"]

The code now supports general algebraic sequences, in which the initial term is not required to be 0.  Such a sequence is specified using SeriesRoot.
(Additionally, the automaton is now minimized by default.)

AutomatonGraph[AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AlgebraicSequence[SeriesRoot[{Function[y,-y+1+x (y^2)],1+((O[x])^1)}]][n],2],n,Method->"Diagonal"]]]

Direct use of CatalanNumber is also supported.

AutomatonGraph[AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[CatalanNumber[n],2],n]]]

However, to maintain consistency with the published paper, the sequences in this section have initial term 0, and the automata are not minimized.

Catalan numbers (A000108)

Table[CatalanNumber[n],{n,0,10}]

Series[((1-(Sqrt[1-4 x]))/(2 x)),{x,0,10}]

Remove the constant term:

-y+1+x (y^2)/. y->y+1

(2^α)

With[{sequence=AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],2],n,Method->"Diagonal",MinimizeComplexity->False]},AutomatonGraph[AutomaticSequenceAutomaton[sequence],GraphLayout->"SpringEmbedding"]]

With[{sequence=AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],(2^2)],n,Method->"Diagonal",MinimizeComplexity->False]},AutomatonGraph[AutomaticSequenceAutomaton[sequence],GraphLayout->"SpringEmbedding"]]

With[{sequence=AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],(2^3)],n,Method->"Diagonal",MinimizeComplexity->False]},AutomatonGraph[AutomaticSequenceAutomaton[sequence],GraphLayout->"SpringEmbedding",VertexSize->.3]]

Column[AbsoluteTiming[Grid[Function[modulus,{modulus,MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False],modulus]}]/@(2^(Range[6]))]]]

AbsoluteTiming[With[{modulus=(2^7)},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

AbsoluteTiming[With[{modulus=(2^8)},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

AbsoluteTiming[With[{modulus=(2^9)},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

3 and 5

With[{sequence=AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],3],n,Method->"Diagonal",MinimizeComplexity->False]},AutomatonGraph[AutomaticSequenceAutomaton[sequence],GraphLayout->"SpringEmbedding"]]

With[{sequence=AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-y+x (((1+y))^2)],x][n],5],n,Method->"Diagonal",MinimizeComplexity->False]},AutomatonGraph[AutomaticSequenceAutomaton[sequence],GraphLayout->"SpringEmbedding"]]

Motzkin numbers (A001006)

Table[MotzkinNumber[n],{n,0,10}]

Series[((1-x-(Sqrt[1-2 x-3 (x^2)]))/(2 (x^2))),{x,0,10}]

Remove the constant term:

(x^2) (y^2)+(x-1) y+1/.y->y+1

With[{sequence=AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1+(-1+x) (1+y)+(x^2) (((1+y))^2)],x][n],8],n,Method->{"Diagonal","Degree"->"High"},MinimizeComplexity->False]},AutomatonGraph[AutomaticSequenceAutomaton[sequence],GraphLayout->"SpringEmbedding",ImageSize->500,VertexSize->.5]]

With[{modulus=(2^3)},MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,(x^2) (y^2)+(2 (x^2)+x-1) y+((x^2)+x)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False],modulus]]

AbsoluteTiming[With[{modulus=(5^2)},MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,(x^2) (y^2)+(2 (x^2)+x-1) y+((x^2)+x)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False],modulus]]]

AbsoluteTiming[With[{modulus=(13^2)},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,(x^2) (y^2)+(2 (x^2)+x-1) y+((x^2)+x)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

Riordan numbers (A005043)

Series[((1+x-(Sqrt[1-2 x-3 (x^2)]))/(2 (x+(x^2)))),{x,0,10}]

1+(-1-x) y+x (1+x) (y^2)/.y->y+1

With[{modulus=3},MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1+(-1-x) (1+y)+x (1+x) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False],modulus]]

AbsoluteTiming[With[{modulus=32},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1+(-1-x) (1+y)+x (1+x) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

directed animals (A005773)

Series[((-1+3 x-(Sqrt[1-2 x-3 (x^2)]))/(2 (-1+3 x))),{x,0,10}]

x+(1-3 x) y+(-1+3 x) (y^2)/.y->y+1

AbsoluteTiming[With[{modulus=32},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,x+(1-3 x) (1+y)+(-1+3 x) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

restricted hexagonal polyominoes (A002212)

Series[((1-x-(Sqrt[1-6 x+5 (x^2)]))/(2 x)),{x,0,10}]

1-x+(-1+x) y+x (y^2)/.y->y+1

With[{modulus=8},MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1-x+(-1+x) (1+y)+x (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]

binary trees avoiding  (A159769)

Series[((-1+2 x-2 (x^2)+(Sqrt[1-4 x+4 (x^3)]))/(2 (-2 (x^2)+(x^3)))),{x,0,10}]

-1+x+(1-2 x+2 (x^2)) y+(-2+x) (x^2) (y^2)/. y->y+1

Column[AbsoluteTiming[Grid[Function[modulus,{modulus,MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-1+x+(1-2 x+2 (x^2)) (1+y)+(-2+x) (x^2) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]}]/@(2^(Range[5]))]]]

Select[{3,7,11,15,19,23,25,27,29,31},Mod[#,4]≠3&]

AbsoluteTiming[With[{modulus=64},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,-1+x+(1-2 x+2 (x^2)) (1+y)+(-2+x) (x^2) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

Select[{3,7,9,11,13,15,19,22,23,25,27,29,31,35,37,39,43,47,51,55,57,59,61,63},!MatchQ[Mod[#,4],3]&&!MatchQ[Mod[#,32],25|29]&]

binary trees avoiding  (A159771)

Series[((1-2 x+3 (x^2)-(Sqrt[1-4 x+2 (x^2)-4 (x^3)+(x^4)]))/(4 (x^2))),{x,0,10}]

1-x+(x^2)+(-1+2 x-3 (x^2)) y+2 (x^2) (y^2)/.y->y+1

Column[AbsoluteTiming[Grid[Function[modulus,{modulus,MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1-x+(x^2)+(-1+2 x-3 (x^2)) (1+y)+2 (x^2) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False],modulus]}]/@(2^(Range[5]))]]]

AbsoluteTiming[With[{modulus=(2^6)},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1-x+(x^2)+(-1+2 x-3 (x^2)) (1+y)+2 (x^2) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

AbsoluteTiming[With[{modulus=(2^7)},MissingResidues[EchoFunction[AutomatonStateCount@*AutomaticSequenceAutomaton]@AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,1-x+(x^2)+(-1+2 x-3 (x^2)) (1+y)+2 (x^2) (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]]]

Select[{3,7,11,13,15},Mod[#,4]≠3&]

Select[{3,7,11,13,15,19,21,23,27,29,31},Mod[#,4]≠3&&Mod[#,16]≠13&]

Select[{3,7,11,13,15,19,21,23,27,29,31,35,37,39,43,45,47,51,53,55,59,61,63},Mod[#,4]≠3&&Mod[#,16]≠13&&Mod[#,32]≠21&]

Select[{3,7,11,13,15,19,21,23,27,29,31,35,37,39,43,45,47,51,53,55,59,61,63,67,71,75,77,79,83,85,87,91,93,95,99,101,103,107,109,111,115,117,119,123,125,127},Mod[#,4]≠3&&Mod[#,16]≠13&&Mod[#,32]≠21&&Mod[#,64]≠37&]

permutations avoiding 3412 and 2143 (A029759)

Series[((1-3 x)/((1-2 x) (Sqrt[1-4 x]))),{x,0,10}]

A polynomial of the correct form can be obtained by dropping the first two terms and dividing the remaining terms by 2:

Collect[(((4 x-1) (((2 x-1))^2) (y^2)+(((3 x-1))^2)/. y->1+x+2x y)/(4 x)),y,Factor]

Column[AbsoluteTiming[Grid[Function[modulus,{modulus,MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,x (1-4 x+3 (x^2)+4 (x^3))+(1+x) (((-1+2 x))^2) (-1+4 x) y+x (((-1+2 x))^2) (-1+4 x) (y^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]}]/@(2^(Range[6]))]]]

Select[{5,7,9,13,15},!MemberQ[{5,7},Mod[#,8]]&]

Select[{5,7,9,13,15,17,21,23,25,27,29,31},!MemberQ[{5,7},Mod[#,8]]&&Mod[#,16]≠9&]

Select[{5,7,9,13,15,17,21,22,23,25,27,29,31,33,37,39,41,45,47,49,51,53,55,57,59,61,63},!MemberQ[{5,7},Mod[#,8]]&&Mod[#,16]≠9&&!MemberQ[{17,27},Mod[#,32]]&]

Multiply by 2:

2{5,7}
2{9}
2{17,27}
2{22,33,51}

permutations avoiding 2143 and 1324 (A032351)

Series[((1-5 x+3 (x^2)+(x^2) (Sqrt[1-4 x]))/(1-6 x+8 (x^2)-4 (x^3))),{x,0,10}]

A polynomial of the correct form can be obtained by dropping the first three terms and dividing the remaining terms by 2:

Collect[((-1+4 x+(x^2)+2 (1-5 x+3 (x^2)) y+(-1+6 x-8 (x^2)+4 (x^3)) (y^2)/. y->1+x+2(x^2)+2(x^2)y)/(4(x^4))),y,Factor]

Column[AbsoluteTiming[Grid[Function[modulus,{modulus,MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,x (3-4 x+4 (x^2))+(-1+8 x-12 (x^2)+8 (x^3)) y+(-1+6 x-8 (x^2)+4 (x^3)) (y^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True],modulus]}]/@(2^(Range[7]))]]]

Select[{1,5},Mod[#,4]≠1&]

Select[{1,5,9,13,15},Mod[#,4]≠1&]

Select[{1,5,7,9,13,15,17,21,22,25,27,29,31},Mod[#,4]≠1&&Mod[#,16]≠15&]

Select[{1,5,7,9,13,15,17,19,21,22,23,25,27,29,30,31,33,37,38,39,41,43,45,47,49,53,54,57,59,61,63},Mod[#,4]≠1&&Mod[#,16]≠15&&!MemberQ[{7,22,27},Mod[#,32]]&]

Select[{1,5,7,9,13,15,17,19,21,22,23,25,27,29,30,31,33,37,38,39,41,43,45,46,47,49,53,54,57,59,61,62,63,65,69,70,71,73,77,79,81,83,85,86,87,89,91,93,94,95,97,101,102,103,105,107,109,111,113,115,117,118,119,121,123,125,127},Mod[#,4]≠1&&Mod[#,16]≠15&&!MemberQ[{7,22,27},Mod[#,32]]&&!MemberQ[{19,23,30,38,43},Mod[#,64]]&]

Multiply by 2:

2{1}
2{}
2{15}
2{7,22,27}
2{19,23,30,38,43}
2{46,62,70,115,119}

permutations avoiding 1342 and 2143 (A109033)

Series[((1-(Sqrt[1-8 x+16 (x^2)-8 (x^3)]))/(4 x (1-x))),{x,0,10}]

Collect[(((((4 x (1-x) y-1))^2)-(1-8 x+16 (x^2)-8 (x^3)))/(8 (-1+x) x)),y,Factor]
%/.y->y+1

Column[AbsoluteTiming[Grid[Function[modulus,{modulus,MissingResidues[AutomaticSequenceReduce[Mod[AlgebraicSequence[Function[y,x+y+2 (-1+x) x (((1+y))^2)],x][n],modulus],n,Method->"Diagonal",MinimizeComplexity->False],modulus]}]/@(2^(Range[10]))]]]

Select[{3,4,5,7},Mod[#,4]≠3&]

Select[{3,4,5,7,9,10,11,12,13,14,15},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&]

Select[{3,4,5,7,8,9,10,11,12,13,14,15,17,18,19,20,21,23,25,26,27,28,29,30,31},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&&!MemberQ[{9,10,14},Mod[#,16]]&]

Select[{3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,25,26,27,28,29,30,31,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,49,50,51,52,53,54,55,57,58,59,60,61,62,63},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&&!MemberQ[{9,10,14},Mod[#,16]]&&!MemberQ[{8,17,18},Mod[#,32]]&]

Select[{3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,49,50,51,52,53,54,55,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,89,90,91,92,93,94,95,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&&!MemberQ[{9,10,14},Mod[#,16]]&&!MemberQ[{8,17,18},Mod[#,32]]&&!MemberQ[{16,33,34,38,54},Mod[#,64]]&]

Select[{3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,49,50,51,52,53,54,55,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,217,218,219,220,221,222,223,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&&!MemberQ[{9,10,14},Mod[#,16]]&&!MemberQ[{8,17,18},Mod[#,32]]&&!MemberQ[{16,33,34,38,54},Mod[#,64]]&&!MemberQ[{24,32,65,66,70,86,120},Mod[#,128]]&]

Select[{3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,313,314,315,316,317,318,319,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,336,337,338,339,340,341,342,343,345,346,347,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,367,369,370,371,372,373,374,375,376,377,378,379,380,381,382,383,385,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420,421,422,423,424,425,426,427,428,429,430,431,432,433,434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,450,451,452,453,454,455,456,457,458,459,460,461,462,463,464,465,466,467,468,469,470,471,473,474,475,476,477,478,479,481,482,483,484,485,486,487,488,489,490,491,492,493,494,495,496,497,498,499,500,501,502,503,504,505,506,507,508,509,510,511},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&&!MemberQ[{9,10,14},Mod[#,16]]&&!MemberQ[{8,17,18},Mod[#,32]]&&!MemberQ[{16,33,34,38,54},Mod[#,64]]&&!MemberQ[{24,32,65,66,70,86,120},Mod[#,128]]&&!MemberQ[{96,129,130,134,150,176,184,240},Mod[#,256]]&]

Select[{3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,336,337,338,339,340,341,342,343,345,346,347,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,367,369,370,371,372,373,374,375,376,377,378,379,380,381,382,383,385,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420,421,422,423,424,425,426,427,428,429,430,431,432,433,434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,450,451,452,453,454,455,456,457,458,459,460,461,462,463,464,465,466,467,468,469,470,471,473,474,475,476,477,478,479,481,482,483,484,485,486,487,488,489,490,491,492,493,494,495,496,497,498,499,500,501,502,503,504,505,506,507,508,509,510,511,513,514,515,516,517,518,519,520,521,522,523,524,525,526,527,528,529,530,531,532,533,534,535,536,537,538,539,540,541,542,543,544,545,546,547,548,549,550,551,552,553,554,555,556,557,558,559,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606,607,608,609,610,611,612,613,614,615,616,617,618,619,620,621,622,623,624,625,626,627,628,629,630,631,632,633,634,635,636,637,638,639,640,641,642,643,644,645,646,647,648,649,650,651,652,653,654,655,656,657,658,659,660,661,662,663,664,665,666,667,668,669,670,671,672,673,674,675,676,677,678,679,680,681,682,683,684,685,686,687,688,689,690,691,692,693,694,695,696,697,698,699,700,701,702,703,705,706,707,708,709,710,711,712,713,714,715,716,717,718,719,720,721,722,723,724,725,726,727,728,729,730,731,732,733,734,735,736,737,738,739,740,741,742,743,744,745,746,747,748,749,750,751,752,753,754,755,756,757,758,759,760,761,762,763,764,765,766,767,768,769,770,771,772,773,774,775,776,777,778,779,780,781,782,783,784,785,786,787,788,789,790,791,792,793,794,795,796,797,798,799,800,801,802,803,804,805,806,807,808,809,810,811,812,813,814,815,816,817,818,819,820,821,822,823,825,826,827,828,829,830,831,832,833,834,835,836,837,838,839,840,841,842,843,844,845,846,847,848,849,850,851,852,853,854,855,856,857,858,859,860,861,862,863,864,865,866,867,868,869,870,871,872,873,874,875,876,877,878,879,880,881,882,883,884,885,886,887,888,889,890,891,892,893,894,895,897,898,899,900,901,902,903,904,905,906,907,908,909,910,911,912,913,914,915,916,917,918,919,920,921,922,923,924,925,926,927,928,929,930,931,932,933,934,935,936,937,938,939,940,941,942,943,944,945,946,947,948,949,950,951,952,953,954,955,956,957,958,959,960,961,962,963,964,965,966,967,968,969,970,971,972,973,974,975,976,977,978,979,980,981,982,983,984,985,986,987,988,989,990,991,993,994,995,996,997,998,999,1000,1001,1002,1003,1004,1005,1006,1007,1008,1009,1010,1011,1012,1013,1014,1015,1016,1017,1018,1019,1020,1021,1022,1023},Mod[#,4]≠3&&!MemberQ[{4,5},Mod[#,8]]&&!MemberQ[{9,10,14},Mod[#,16]]&&!MemberQ[{8,17,18},Mod[#,32]]&&!MemberQ[{16,33,34,38,54},Mod[#,64]]&&!MemberQ[{24,32,65,66,70,86,120},Mod[#,128]]&&!MemberQ[{96,129,130,134,150,176,184,240},Mod[#,256]]&&!MemberQ[{56,112,216,224,256,257,258,262,278,304,320,448},Mod[#,512]]&]

{8,17,18}-(2^4)
{16,33,34,38,54}-(2^5)
{24,32,65,66,70,86,120}-(2^6)
{96,129,130,134,150,176,184,240}-(2^7)
{56,112,216,224,256,257,258,262,278,304,320,448}-(2^8)
{48,192,312,513,514,518,534,600,640,856,880,984}-(2^9)

Grid[Table[{{(2^(α-1))+1,(2^(α-1))+2,(2^(α-1))+6,(2^(α-1))+22},(2^α)},{α,6,12}]]

Computing automata for diagonals of rational functions modulo (p^α)

Apéry numbers (A005259)

Table[AperyNumber[n],{n,0,10}]

Table[DiagonalSequence[(1/((1-x1-x2) (1-x3-x4)-x1 x2 x3 x4)),{x1,x2,x3,x4}][n],{n,0,10}]

(2^α)

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],2],n,Method->"Diagonal"]]
AutomatonGraph[%]

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(2^2)],n,Method->"Diagonal"]]
AutomatonGraph[%]

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(2^3)],n,Method->"Diagonal"]]
AutomatonGraph[%]

AbsoluteTiming[AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(2^4)],n,Method->{"Diagonal","Degree"->"High"}]]]
AutomatonGraph[%[[2]],VertexSize->.4]

AbsoluteTiming[AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(2^5)],n,Method->{"Diagonal","Degree"->"High"},Monitor->True]]]
AutomatonGraph[%[[2]],ImageSize->1000,VertexSize->.8]

(3^α)

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],3],n,Method->"Diagonal"]]
AutomatonGraph[%]

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(3^2)],n,Method->"Diagonal"]]
AutomatonGraph[%,GraphLayout->"SpringEmbedding"]

(5^α)

There are many rational expressions for which the sequence of Apéry numbers is the diagonal:

Table[DiagonalSequence[(1/((1-x1) ((1-x2) (1-x3) (1-x4) (1-x5)-x1 x2 x3))),{x1,x2,x3,x4,x5}][n],{n,0,5}]

Table[DiagonalSequence[(1/((1-x1 x2) (1-x3-x4-x1 x3 x4) (1-x5-x6-x2 x5 x6))),{x1,x2,x3,x4,x5,x6}][n],{n,0,5}]

Table[DiagonalSequence[(1/((1-x1-x2) (1-x3-x4)-x1 x2 x3 x4)),{x1,x2,x3,x4}][n],{n,0,5}]

For automaton computations, it can be advantageous to choose the expression that minimizes the number of monomials in (denominator^((p-1)(p^(α-1)))):

AbsoluteTiming[With[{denominator=(1-x1) ((1-x2) (1-x3) (1-x4) (1-x5)-x1 x2 x3),modulus=(5^2),p=5},Length[Expand[(denominator^(modulus-(modulus/p))),Modulus->modulus]]]]

AbsoluteTiming[With[{denominator=(1-x1 x2) (1-x3-x4-x1 x3 x4) (1-x5-x6-x2 x5 x6),modulus=(5^2),p=5},Length[Expand[(denominator^(modulus-(modulus/p))),Modulus->modulus]]]]

AbsoluteTiming[With[{denominator=(1-x1-x2) (1-x3-x4)-x1 x2 x3 x4,modulus=(5^2),p=5},Length[Expand[(denominator^(modulus-(modulus/p))),Modulus->modulus]]]]

By default, AutomaticSequenceReduce uses the preceding rational expression in four variables, which was discovered by Armin Straub.

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],5],n,Method->"Diagonal"]]
AutomatonGraph[%,GraphLayout->"SpringEmbedding"]

AbsoluteTiming[automaton25=AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(5^2)],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True]]]

AutomatonGraph[automaton25,ImageSize->700,VertexSize->1]

After we see 1 or 3, we are in one of these five states:

Union[Cases[First[automaton25],{_->v_,1|3}:>v]]

There are some edges pointing from one of those five states to a state that is not in those five:

Cases[First[automaton25],{Alternatives@@{2,4,6,8,10}->Except[Alternatives@@{2,4,6,8,10}],_}]

Union[Cases[First[automaton25],{Alternatives@@{2,4,6,8,10}->v:Except[Alternatives@@{2,4,6,8,10}],_}:>v]]

Whenever we see 1 or 3 from those states we go back to state 6, which outputs 0:

With[{externalstates=Union[Cases[First[automaton25],{Alternatives@@{2,4,6,8,10}->v:Except[Alternatives@@{2,4,6,8,10}],_}:>v]]},Cases[First[automaton25],{Alternatives@@externalstates->v_,1|3}:>v]]

And we don’t leave that set except to go to state 6:

Union[Cases[First[automaton25],{Alternatives@@{5,6,11,7,12}->v_,_}:>v]]

AutomatonGraph[Automaton[Cases[First[automaton25],{Alternatives@@{5,6,11,7,12}->_,_}],5]]

In other words, only these vertices can be reached after seeing 1 or 3:

Union[{2,4,6,8,10},{5,6,11,7,12}]

And all the edges labeled 1 or 3 leaving one of those nine vertices point to vertex 6:

Cases[First[automaton25],{Alternatives@@{2,4,6,8,10,5,7,11,12}->v_,1|3}:>v]

So Beuker’s conjecture holds for α==2.

Delete vertices that correspond to outputs divisible by 5:

With[{graph=AutomatonGraph[automaton25,GraphLayout->"SpringEmbedding",VertexSize->.35]},VertexDelete[graph,Cases[VertexList[graph],Interpretation[label_/;Divisible[label,5],_]]]]

Mod[(23^(Range[0,20])),25]

Clear[automaton25]

(7^α)

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],7],n,Method->"Diagonal"]]
AutomatonGraph[%,GraphLayout->"SpringEmbedding"]

AbsoluteTiming[AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AperyNumber[n],(7^2)],n,Method->"Diagonal",MinimizeComplexity->False,Monitor->True]]]
AutomatonGraph[%[[2]],GraphLayout->"SpringEmbedding",ImageSize->1000,VertexSize->.7]

ζ(2) Apéry numbers (A005258)

Table[(∑_(k=0)^n)((Binomial[n,k])^2) Binomial[n+k,k],{n,0,6}]

Table[DiagonalSequence[(1/((1-x1) (1-x2) (1-x3) (1-x4)-(1-x1) x1 x2 x3)),{x1,x2,x3,x4}][n],{n,0,6}]

AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[DiagonalSequence[(1/((1-x1) (1-x2) (1-x3) (1-x4)-(1-x1) x1 x2 x3)),{x1,x2,x3,x4}][n],(3^2)],n,MinimizeComplexity->False]]
AutomatonGraph[%,GraphLayout->"SpringEmbedding"]

Computing automata for algebraic sequences by incrementally raising the exponent

central trinomial coefficients (A002426)

The following is support code for the function AutomaticSequenceReduceList below, which computes automata for an algebraic sequence modulo p,(p^2),…,(p^α).
It is quite slow relative to AutomaticSequenceReduce, however.

ExtractFactor[polynomial_,y_,series_]:=
With[{factorlist=Select[First/@Rest[FactorList[polynomial]],Normal[#/. y->series]===0&]},
If[factorlist=={},
Print["no factors!"];$Failed
,
If[Length[factorlist]≥2,
Print["multiple factors! using last"]
];
SeriesRoot[{Function@@{Last[factorlist]/.y->#},series}]
]
]

SeriesRootPlus[SeriesRoot[{p_Function,f0_SeriesData}],SeriesRoot[{q_Function,g0_SeriesData}]]:=
Module[{factor,x,y},
factor=ExtractFactor[Resultant[p[x],q[y-x],x],y,f0+g0];
factor/;factor=!=$Failed
]
SeriesRootPlus[p_SeriesRoot,q_SeriesRoot,symbolicroots__SeriesRoot]:=
SeriesRootPlus[SeriesRootPlus[p,q],symbolicroots]
SeriesRootPlus[p_SeriesRoot]:=p

SeriesRootSubtract[SeriesRoot[{p_Function,f0_SeriesData}],SeriesRoot[{q_Function,g0_SeriesData}]]:=
Module[{factor,x,y},
factor=ExtractFactor[Resultant[p[x],q[x-y],x],y,f0-g0];
factor/;factor=!=$Failed
]

SeriesRootTimes[SeriesRoot[{p_Function,f0_SeriesData}],SeriesRoot[{q_Function,g0_SeriesData}]]:=
Module[{factor,x,y},
factor=ExtractFactor[Resultant[p[x],(x^(Exponent[q[x],x]))q[y/x],x],y,f0 g0];
factor/;factor=!=$Failed
]
SeriesRootTimes[p_SeriesRoot,q_SeriesRoot,symbolicroots__SeriesRoot]:=
SeriesRootTimes[SeriesRootTimes[p,q],symbolicroots]
SeriesRootTimes[p_SeriesRoot]:=p

DilateSeries[series:HoldPattern[SeriesData][variable_,__],exponent_]:=(Normal[series]/. variable->(variable^exponent))+((O[variable])^(exponent series[[5]]))

DilateSeriesRoot[SeriesRoot[{f_Function,series:HoldPattern[SeriesData][x_,__]},options___],exponent_Integer?Positive]:=SeriesRoot[{f/. x->(x^exponent),DilateSeries[series,exponent]},options]

AutomaticSequenceReduceList[Mod[AlgebraicSequence[seriesroot:SeriesRoot[{_Function,HoldPattern[SeriesData][x_,0,__]}]][n_],modulus_?PrimePowerQ],n_,OptionsPattern[]]:=
Module[{p,α,unprojectedcurrentcomponentseries,componentautomaton,christoloperator,coefficientrules,liftedseries,Deltaseries,y,δ,christoloutputseries,deltapolynomial,Deltapolynomial,seriesterms},
{p,α}=NumberTheory`PrimePower[modulus];
unprojectedcurrentcomponentseries=seriesroot;
componentautomaton[0]=AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AlgebraicSequence[seriesroot][n],p],n,MinimizeComplexity->OptionValue[MinimizeComplexity],Monitor->OptionValue[Monitor]]];
Do[
christoloperator=AlgebraicSequenceReduce[AutomaticSequence[componentautomaton[i-1]][n],n][[0,1]][δ]/.\[FormalX]->x;
coefficientrules=CoefficientRules[christoloperator,δ];
Do[
(* With indices, liftedseries would be liftedseries[i-1,j-1], but we only need one at a time.
 liftedseries[i-1,j-1] is congruent modulo p^(i-1) to the lifting of the ith component. *)
liftedseries=If[j==2,
unprojectedcurrentcomponentseries,
(* Compute a polynomial for liftedseries[i-1,j-1]. *)
SeriesRootSubtract[If[i==1,seriesroot,SeriesRootPlus@@Table[Deltaseries[i-k,k],{k,2,i}]],Replace[SeriesRootPlus@@Table[SeriesRootTimes[SeriesRoot[{Function@@{y,y-(p^(k-1))},(p^(k-1))+((O[x])^(Deltaseries[i-1,k][[1,2,5]]))}],Deltaseries[i-1,k]],{k,2,j-1}],SeriesRootPlus[]->SeriesRoot[{#&,((O[x])^100)}]]]
];
(* Compute a polynomial for the Christol operator evaluated at liftedseries[i-1,j-1]. *)
christoloutputseries=SeriesRootPlus@@(SeriesRootTimes[SeriesRoot[{Function@@{y,y-#2},#2+((O[x])^(#1[[1]] liftedseries[[1,2,5]]))}],DilateSeriesRoot[liftedseries,#1[[1]]]]&@@@coefficientrules);
(* Into the polynomial we just computed, plug in the Christol polynomial (in δ) times (p^(j-1)). *)
deltapolynomial=christoloutputseries[[1,1]][(p^(j-1)) christoloperator];
deltapolynomial=Cancel[(deltapolynomial/(GCD@@Last/@CoefficientRules[deltapolynomial,{x,δ}]))];
(* Compute an Ore polynomial for δ, and lift it to an algebraic (Subsuperscript[Δ, i-1, j]) with integer coefficients. *)
Deltapolynomial=OrePolynomial[deltapolynomial,{x,δ},p,Modulus->p];
(* This test is not perfect because CoefficientList[((O[x])^10),x] is the empty list, so there may be a problem for the zero sequence. *)
seriesterms=Select[δ[x]/.DeleteDuplicates[SeriesSolve[Deltapolynomial==0/.δ->δ[x],δ[x],{x,0,5},Method->"ExactAlgebraics"]],MatchQ[#,HoldPattern[SeriesData][_,_,{__Integer},__]]&];
Deltaseries[i-1,j]=SeriesRoot[{Function@@{δ,Deltapolynomial},First[seriesterms]}]
,
{j,2,α+1-i}
];
unprojectedcurrentcomponentseries=SeriesRootPlus@@Table[Deltaseries[i-j+1,j],{j,2,i+1}];
componentautomaton[i]=AutomaticSequenceAutomaton[AutomaticSequenceReduce[Mod[AlgebraicSequence[unprojectedcurrentcomponentseries][n],p],n,MinimizeComplexity->OptionValue[MinimizeComplexity],Monitor->OptionValue[Monitor]]]
,
{i,1,α-1}
];
Table[AutomaticSequence[componentautomaton[i]][n],{i,0,α-1}]
]
Options[AutomaticSequenceReduceList]={MinimizeComplexity->True,Monitor->False};

Series[(1/(Sqrt[1-2 x-3 (x^2)])),{x,0,10}]

AutomaticSequenceAutomaton/@AutomaticSequenceReduceList[Mod[AlgebraicSequence[SeriesRoot[{Function[y,(x+1) (3 x-1) (y^2)+1],Series[(1/(Sqrt[1-2 x-3 (x^2)])),{x,0,15}]}]][n],(2^2)],n,MinimizeComplexity->False]
AutomatonGraph/@%