\documentclass[12pt]{article}
\usepackage[T1]{fontenc}
\usepackage{lmodern,amsmath,amssymb,amsthm}
\usepackage[margin=1in]{geometry}
\usepackage{needspace}
\usepackage[colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue]{hyperref}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{example}[theorem]{Example}
\newcommand{\Z}{\mathbb Z}
\newcommand{\QQ}{\mathbb Q}
\newcommand{\Fp}{\mathbb F_p}
\newcommand{\fall}[2]{(#1)_{#2}}
\DeclareMathOperator{\id}{id}
\title{Integrality of iterated averages\\and pointwise prime powers}
\author{GPT-6 Astra, edited and proofread by Darij Grinberg\footnote{Section 8 has not been proofread.}}
\date{\today}

\begin{document}

\maketitle

\begin{abstract}
Suppose that we repeatedly replace the values of an integer-valued
function on $\Z$ by averages over $p$ consecutive arguments, where
$p$ is prime. If the first $\ell$ rounds of averaging give integers,
will this remain true when we begin with the $p$th powers of the
original values? We prove that the answer is yes for $\ell\le p$.
For $p=2$ and $p=3$, the answer is yes for every $\ell$, and therefore
also when integrality is required at every stage. For larger primes,
a known counterexample first fails at round $p+1$. The proofs use
finite differences and polynomial interpolation. We also explain
connections with earlier work on divisibility of finite differences
and Stirling numbers.
\medskip

\textbf{Manifest.} This is an answer to MathOverflow question
\#68595 (which also gives a new answer to \#66959).
It was produced by GPT-6 Astra (Researcher subscription).
The main proof was found in one 3-minute shot, and expanded and
polished over several rounds of edits.
I have proofread the whole text except for Section~\ref{sec:connect}.

This text is in the public domain. --Darij
\end{abstract}

\section{Statement and finite differences}
Start with an integer-valued function on the integers. Replace each value
by the average of that value and the next $p-1$ values, where $p$ is a
prime. The new values need not be integers. If they are, we can repeat
the operation and ask whether they remain integers at each stage.
A classical olympiad problem, proposed by D. Tamarkin at the All-Russian
Mathematical Olympiad 1993 \cite{Tamarkin}, concerns this process for
$p=2$: It claims that if every stage is integer-valued, then the same
is true when we start with the squares of the original values.
In this note, we consider the case of general $p$ and taking the entries
to the $p$th power; we also ask what happens when we assume integrality
for only finitely many stages. The results revealed are weaker but no
less beautiful.

We begin by fixing the notation.
We consider the $\QQ$-vector space $\QQ^\Z$, which consists of all
functions $f : \Z \to \QQ$.
These functions can be regarded as infinite sequences
\[
\left(\ldots, f(-2), f(-1), f(0), f(1), f(2), \ldots\right)
\]
(infinite in both directions).
Such a function $f : \Z \to \QQ$ is said to be \emph{integral} if all
its values $f(n)$ are integers, i.e., if $f \in \Z^\Z$.
Sums, differences, and integer multiples of integral functions are again
integral; in other words, $\Z^\Z$ is a $\Z$-submodule of $\QQ^\Z$.

We shall consider linear (that is, $\QQ$-linear) operators on $\QQ^\Z$
such as the shift and finite-difference operators defined below.
Here linearity means that $A(af+bg)=aAf+bAg$ for rational numbers
$a,b$ and functions $f,g$.
If $A$ is such an operator and $f \in \QQ^\Z$ is a function, then we
will write $Af(n)$ for the value $(Af)(n)$.

Fix a prime $p$. Define a linear operator
$F_p \colon \QQ^\Z \to \QQ^\Z$ (called the \emph{averaging operator}) by
\[
 (F_pf)(n)=\frac1p\sum_{r=0}^{p-1}f(n+r)
 \qquad\text{for all }f \in \QQ^\Z \text{ and } n\in\Z.
\]
That is, $F_p$ acts on a function $f$ by replacing each value
$f(n)$ by the average of $p$ consecutive values (starting with
$f(n)$).
Furthermore, for a function $f\colon\Z\to\QQ$, define the function
$f^p : \Z\to\QQ$ by 
\[
 f^p(n)=f(n)^p
 \qquad\text{for all }n\in\Z.
\]

Operator powers mean repeated application: $F_p^2f=F_p(F_pf)$,
and $F_p^0f=f$. These are different from the pointwise powers $f^p$.

A function $f : \Z \to \QQ$ will be called
\emph{$\ell$-average-integral} for a given
nonnegative integer $\ell$ if $F_p^k f$ is integral for every
$0\le k\le\ell$. Note that an $\ell$-average-integral function
$f$ is thus automatically integral (since we can take $k=0$ and
obtain the integrality of $F_p^k f = F_p^0 f = f$).
Thus the parameter $\ell$ in ``$\ell$-average-integral''
records how many successive averaging steps
we require to give integers.

We say that $f : \Z \to \QQ$ is \emph{average-integral} if
$F_p^k f$ is integral for every $k\ge0$.
Of course, this is the same as being $\ell$-average-integral
for all $\ell \geq 0$ simultaneously.
For example, when $p=2$, the function
$f \colon n \mapsto 2n$ is average-integral: its $k$th
average-iterate $F_p^k f$ is $n \mapsto 2n+k$.

We will prove the following two properties of
$\ell$-average-integrality:

\begin{theorem}\label{thm:main}
Let $0\le\ell\le p$. If $f\colon\Z\to\Z$ is
$\ell$-average-integral, then $f^p$ is $\ell$-average-integral.
\end{theorem}

\begin{theorem}\label{thm:small-primes}
Let $\ell\ge0$ be arbitrary.
\begin{enumerate}
\item For $p=2$, the $\ell$-average-integral functions form a ring
under pointwise addition and pointwise multiplication.
\item For $p=3$, if $f$ is $\ell$-average-integral, then $f^3$ is
$\ell$-average-integral.
\end{enumerate}
\end{theorem}

\begin{corollary}\label{cor:inf-depth23}
The set of all average-integral functions is closed under
pointwise products for $p=2$ and under pointwise cubes
(i.e., under the map $f \mapsto f^3$) for $p=3$.
\end{corollary}

Theorem~\ref{thm:main} answers a question asked in~\cite{MO-finite}.
Corollary~\ref{cor:inf-depth23} was
already discussed in~\cite{MO-original}. For primes $p>3$,
the bound in Theorem~\ref{thm:main} cannot be increased:
Proposition~\ref{prop:sharp}
gives a function whose $p$th power first fails at round $p+1$.

Our proofs use several ideas suggested by Zarathustra Brady (``zeb''
on MathOverflow) in his answer to \cite{MO-original}.
They have two major ingredients in common.
First we translate integrality of averages
into divisibility of finite differences. Then we use polynomial
interpolation to study these divisibilities on a finite interval.
This proves Theorem~\ref{thm:main} in Section~\ref{sec:interpolation}.
To handle arbitrarily many rounds for $p=2,3$
(and thus prove Theorem~\ref{thm:small-primes} and
Corollary~\ref{cor:inf-depth23}), we will need a fuller
description of the polynomials that remain integral under averaging.

Our main tool is the \emph{finite difference operator}
$\Delta : \QQ^\Z \to \QQ^\Z$ defined by
\[
\Delta f(n)=f(n+1)-f(n) \qquad \text{ for all } f \in \QQ^\Z \text{ and } n \in \Z.
\]
Unlike averaging, taking a finite difference always preserves
integrality; i.e., we have $\Delta f \in \Z^\Z$ for each $f \in \Z^\Z$.
We will express the condition on repeated averages as a
condition on the divisibility of repeated differences.

Define the \emph{shift operator} $E : \QQ^\Z \to \QQ^\Z$ by
\[
Ef(n)=f(n+1) \qquad \text{ for all } f \in \QQ^\Z \text{ and } n \in \Z.
\]
Thus, $\Delta=E-\id$. All three operators $F_p$, $E$, and $\Delta$ are
$\QQ$-linear endomorphisms of $\QQ^{\Z}$. The latter two preserve
$\Z^{\Z}$.
The first few differences of a function $f : \Z \to \QQ$ are
\[
 \Delta f(n)=f(n+1)-f(n),\qquad
 \Delta^2f(n)=f(n+2)-2f(n+1)+f(n).
\]

Set $q=p-1$.
All divisibility statements about functions below are pointwise:
for example, $\Delta^jf\in p^k\Z^\Z$ means that
$p^k$ divides $\Delta^jf(n)$ for every integer $n$.
The next lemma says that -- as far as integrality is concerned
-- one averaging step ``corresponds'' to a block
of $q=p-1$ difference steps. For instance, when $p=3$ and $f$ is
integral, two integral averaging steps are equivalent to the two conditions
$3\mid\Delta^2f(n)$ and $9\mid\Delta^4f(n)$ for every $n$.

\begin{lemma}[Difference criterion]\label{lem:difference}
For a function $f\colon\Z\to\QQ$ and any nonnegative integer $\ell\ge0$,
the following conditions are equivalent:
\begin{enumerate}
\item The function $f$ is $\ell$-average-integral.
\item We have $\Delta^{kq}f\in p^k\Z^{\Z}$ for all $0\le k\le\ell$.
\item We have $\Delta^j f\in p^{\lfloor j/q\rfloor}\Z^{\Z}$
for all $0\le j\le\ell q$.
\end{enumerate}
\end{lemma}

\begin{proof}
We assume that $f$ is integral, since all three of the ostensibly
equivalent conditions imply this\footnote{Indeed, the second
condition, applied to $k = 0$, says precisely that $f \in \Z^\Z$,
that is, $f$ is integral. The third condition, applied to
$j = 0$, likewise yields the same conclusion. And of course,
the first condition yields integrality of $f$ as well, as we
already know.}.

We first compare one averaging step with $q$ difference steps. If
$A(X)=\sum_i a_iX^i$ is a polynomial, then $A(E)$ means the operator
$\sum_i a_iE^i$. Since $(E^if)(n)=f(n+i)$, we can rewrite the
definition of $F_p$ as
\begin{align}
 F_p = \dfrac{1}{p} (\id+E+\cdots+E^{p-1}).
 \label{eq:Fp=}
\end{align}
Hence,
\begin{align}
 pF_p=\id+E+\cdots+E^{p-1}.
 \label{eq:pFp=}
\end{align}
Moreover, from $\Delta = E - \id$ and $q = p-1$, we obtain
$\Delta^q=(E-\id)^{p-1}$.
These formulas explain why the following polynomial congruence is useful:
\begin{equation}\label{eq:averaging-congruence}
 1+X+\cdots+X^{p-1}\equiv(X-1)^{p-1}\pmod p.
\end{equation}
To see this congruence, reduce modulo $p$ and use
$X^p-1=(X-1)^p$. Dividing both sides by $X-1$ in the polynomial ring
$\mathbb F_p[X]$ gives the congruence \eqref{eq:averaging-congruence}.

That is, in $\Z[X]$, we have
\[
 1+X+\cdots+X^{p-1}
 = (X-1)^{p-1} + pB(X)
 \qquad \text{ for some } B(X) \in\Z[X].
\]
Now substitute the shift operator $E$ for $X$ in
this equality. This gives the operator identity
\[
 pF_p=\Delta^q+pB(E)
\]
(since $pF_p=\id+E+\cdots+E^{p-1}$ and $\Delta^q=(E-\id)^{p-1}$).
Putting
\begin{align}
D=\Delta^q/p
 \label{eq:D=}
\end{align}
and dividing by $p$, we obtain
\begin{align}
 F_p=D+B(E),
 \label{eq:Fp=D+BE}
\end{align}
so that $D = F_p-B(E)$ and $B(E) = F_p-D$.
The important point is that $B(E)$ preserves $\Z^\Z$:
it is an integer linear combination of shifts. For a given integral
function $f$, the identity \eqref{eq:Fp=D+BE} therefore shows that
$F_pf$ is integral if and only if $Df$ is integral.
For example, when $p=3$, we have $B(X)=X$ and
\[
 F_3f(n)=f(n+1)+\frac{f(n+2)-2f(n+1)+f(n)}{3}.
\]
The numerator in the fraction is $\Delta^2f(n)$, making the
one-step divisibility condition particularly transparent.

The same comparison works for every iterate. All three operators
$F_p, D, B(E)$ are
polynomials in $E$, so they commute, and the binomial theorem applies.
If $F_p^if$ is integral for $0\le i\le k$, then
(from $D = F_p-B(E)$) we see that
\[
 D^kf=(F_p-B(E))^kf
 =\sum_{i=0}^k\binom ki(-B(E))^{k-i}F_p^if
\]
is integral as well: each summand is obtained from an
integral function (viz., $F_p^if$) by applications of the
operator $B(E)$ (which preserves $\Z^\Z$) and integer linear combinations.
Conversely, suppose that $D^if$ is integral for every $0\le i\le k$.
Using $F_p=D+B(E)$ and the commutativity of $D$ and $B(E)$, we obtain
\[
 F_p^kf=(D+B(E))^kf
 =\sum_{i=0}^k\binom ki B(E)^{k-i}D^if.
\]
Each $D^if$ is integral by assumption, and applying $B(E)$ any
number of times preserves integrality. Multiplication by the integer
$\binom ki$ also preserves integrality, as does summing the resulting
functions. Thus $F_p^kf$ is integral.
Altogether, we thus have shown that $F_p^kf$ is integral for every
$0\le k \le \ell$ if and only if $D^kf$ is integral for every
$0\le k \le \ell$.
Since $D^kf=\Delta^{kq}f/p^k$, this proves the equivalence
of conditions 1 and 2 in the lemma.

It remains to show that conditions 2 and 3 are equivalent.
Clearly, condition 3 implies condition 2 (by applying it to $j = kq$).
Conversely, consider a difference order $j$ that is not necessarily a
multiple of $q$. Dividing $j$ by $q$ with remainder, we write $j$
as $j=kq+r$, where $0\le r<q$; then
$k=\lfloor j/q\rfloor\le\ell$. Condition 2 thus says that
$\Delta^{kq}f=p^kg$ for an integral function $g$.
Taking $r$ more differences then gives
$\Delta^jf=p^k\Delta^rg$ (since $j=kq+r$), which is
still divisible by $p^k$ at every integer.
That is, $\Delta^j f \in p^k \Z^\Z = p^{\lfloor j/q\rfloor}\Z^\Z$.
This yields condition 3.
\end{proof}

\section{The polynomial congruence}\label{sec:congruence}
Our next step is a statement about polynomials. We will later use it
on a polynomial that agrees with $f$ on a finite interval. A small
amount of notation lets us keep track of powers of $p$ without being
concerned about other prime factors.

Let $R=\Z_{(p)}$, the subring of $\QQ$ consisting of fractions (of
the form $r/s$ with $r, s \in \Z$) whose
denominators are not divisible by $p$. For example, if $p=3$, then
$1/2$ belongs to $R$, whereas $1/3$ does not. A statement such as
$a\in p^kR$ (for a given nonnegative integer $k$)
means that, after cancellation, the numerator of $a$ has
at least $k$ factors of $p$ and its denominator has none. When $a$ is
an integer, this is just the usual divisibility by $p^k$.

For a nonzero rational number $a$, write $v_p(a)$ for the exponent of
$p$ in its numerator minus the exponent in its denominator, and set
$v_p(0)=+\infty$. Thus $a\in p^kR$ is equivalent to $v_p(a)\ge k$.
We will use Legendre's formula (also known as de Polignac's formula)
\begin{equation}\label{eq:legendre}
 v_p(t!)=\sum_{i\ge1}\left\lfloor\frac{t}{p^i}\right\rfloor
 \qquad\text{for all }t\ge0.
\end{equation}
Indeed, among $1,\ldots,t$, each multiple of $p$ contributes one factor
of $p$, each multiple of $p^2$ contributes a further factor, and so on.
The sum counts every factor exactly once.

We use the same operators $E,\Delta,F_p$ on polynomials as on functions:
\[
 (EH)(x)=H(x+1),\qquad
 (\Delta H)(x)=H(x+1)-H(x),\qquad
 (F_pH)(x)=\frac1p\sum_{r=0}^{p-1}H(x+r).
\]
These definitions are compatible with evaluation at integers.
More precisely, for any $H \in \QQ[x]$,
let $\operatorname{ev}(H) \in \QQ^\Z$ be the function
$n\mapsto H(n)$. Then
\begin{equation}\label{eq:evaluation-commutes}
 \operatorname{ev}(AH)=A\operatorname{ev}(H)
 \qquad\text{for }A=E,\Delta,F_p.
\end{equation}
For example, both sides for $A=\Delta$, evaluated at $n$, equal
$H(n+1)-H(n)$. Thus it makes no difference whether we first take
a polynomial difference and then evaluate, or first evaluate and
then take the difference of the resulting function. The same holds
for any number of iterations. We will use this compatibility when
passing between polynomials and functions.

For finite differences, a convenient $R$-module basis of the
polynomial ring $R[x]$ is given by the \emph{falling factorials}
\[
 \fall{x}{j}=x(x-1)\cdots(x-j+1)\qquad \text{ for all } j \geq 0.
\]
The first few are $\fall{x}{0} = 1$, $\fall{x}{1} = x$,
$\fall{x}{2} = x(x-1)$, and $\fall{x}{3} = x(x-1)(x-2)$.
Each has leading coefficient $1$, with one polynomial of each degree.
We can therefore expand any polynomial in $R[x]$ uniquely in this
basis: subtract a multiple of the falling factorial of highest degree,
then repeat with the remainder. For example,
\[
 x^2=\fall{x}{2}+\fall{x}{1},\qquad
 x^3=\fall{x}{3}+3\fall{x}{2}+\fall{x}{1}.
\]
So each polynomial $P \in R[x]$ can be written uniquely as
$\sum_{j \geq 0} c_j\fall{x}{j}$ with $c_j \in R$ (such that $c_j = 0$
for all but finitely many $j$). The coefficients $c_j$ are called
the \emph{falling-factorial coefficients} of $P$.
This holds equally well for any commutative ring in place of $R$.

The falling factorials also satisfy $\Delta\fall{x}{0}=0$ and
\begin{equation}\label{eq:falling-difference}
 \Delta\fall{x}{j}=j\fall{x}{j-1}\qquad\text{for all }j\ge1.
\end{equation}
This follows by factoring $\fall{x}{j-1}$ out of
$\fall{x+1}{j}-\fall{x}{j}$; the remaining factor is
$(x+1)-(x-j+1)=j$.
Iterating~\eqref{eq:falling-difference} and setting $x=0$ shows that, if
$H(x)=\sum_{j\ge0}c_j\fall{x}{j}$, then
\begin{equation}\label{eq:newton-coefficient}
 \Delta^m H(0)=m!c_m\qquad\text{for all }m\ge0.
\end{equation}
Indeed, terms $c_j\fall{x}{j}$ with $j<m$ have disappeared
after $m$ differences, and those
with $j>m$ still have a factor $x$,
so only the term for $j=m$ remains when we evaluate at $x=0$.

We first isolate a fact about falling-factorial coefficients over
$\Fp$, the field of residues modulo $p$. The falling factorials form
a basis of $\Fp[x]$ by the same leading-coefficient argument as above.
For $P\in\Fp[x]$, write $[\fall{x}{i}]P$ for the $\fall{x}{i}$-coefficient of
$P$ in this basis. We set $[\fall{x}{i}]P=0$ when $i<0$;
this is a convention for coefficients and does not define falling
factorials with negative indices.

Multiplication by $x^p-x$ has a particularly simple effect on this
basis: it shifts every index up by $p$.

\begin{lemma}\label{lem:falling-shift}
Put $y=x^p-x\in\Fp[x]$. Then:
\begin{enumerate}
\item We have
\begin{equation}\label{eq:block-factorial}
 y\fall{x}{i}=\fall{x}{p+i}
 \qquad\text{for all }i\ge0.
\end{equation}
\item More generally, we have
\begin{equation}\label{eq:block-factorials}
 y^a\fall{x}{j}=\fall{x}{ap+j}
 \qquad\text{for all }a,j\ge0.
\end{equation}
\item For every $P\in\Fp[x]$, we have
\begin{equation}\label{eq:falling-coefficient-shift}
 [\fall{x}{i}](yP)=[\fall{x}{i-p}]P
 \qquad\text{for all integers }i.
\end{equation}
\item More generally, for every $P\in\Fp[x]$, we have
\begin{equation}\label{eq:falling-coefficient-shifts}
 [\fall{x}{i}](y^aP)=[\fall{x}{i-ap}]P
 \qquad\text{for all }a\ge0\text{ and all integers }i.
\end{equation}
\end{enumerate}
\end{lemma}
\begin{proof}
We first recall the classical finite-field property
\[
 x(x-1)\cdots(x-p+1)=x^p-x,
\]
which holds because
both sides are monic polynomials of degree $p$ whose roots
are precisely the elements of $\Fp$.
This identity is also known as Lagrange's congruence;
see~\cite[Section~10]{Gessel}, where it is used in an umbral proof
of Touchard's congruence for Bell numbers.
Now, we study the product
$\fall{x}{p+i} = x(x-1)\cdots(x-p-i+1)$. Its first $p$ factors
have product
\[
 x(x-1)\cdots(x-p+1)=x^p-x=y.
\]
Its remaining $i$ factors are
$x-p,x-p-1,\ldots,x-p-i+1$, which (because we are in characteristic
$p$) are precisely the same as $x,x-1,\ldots,x-i+1$; thus their
product is $\fall{x}{i}$. So we obtain
$y\fall{x}{i}=\fall{x}{p+i}$. This proves \eqref{eq:block-factorial}
(including the case $i=0$, when the remaining product is empty).
Applying \eqref{eq:block-factorial} successively $a$ times yields
\eqref{eq:block-factorials}.

To prove \eqref{eq:falling-coefficient-shift},
write $P=\sum_{j\ge0}c_j\fall{x}{j}$. Part (1) gives
\[
 yP=\sum_{j\ge0}c_j\fall{x}{p+j}.
\]
Thus the $\fall{x}{i}$-coefficient is $c_{i-p}$ if $i\ge p$, and
is zero otherwise. This is exactly \eqref{eq:falling-coefficient-shift},
with our convention for negative indices.
To prove \eqref{eq:falling-coefficient-shifts}, apply
\eqref{eq:falling-coefficient-shift} successively $a$ times.
\end{proof}

We now use this shift to describe which coefficients of a $p$th
power can occur at indices less than $p^2$.

\begin{lemma}\label{lem:frobenius-support}
Let $H\in\Fp[x]$, and let $Q$ be the remainder when $H$ is divided
by $x^p-x$. Suppose that $\deg Q\le d$, where $0\le d<p$, and write
\[
 H(x)^p=\sum_{j\ge0}\gamma_j\fall{x}{j}
 \qquad \text{ with } \gamma_j \in \Fp
\]
(so that $\gamma_j$ is the $\fall{x}{j}$-coefficient of $H(x)^p$
for each $j\geq 0$).
Then
\[
 \gamma_{ap+b}=0
 \qquad\text{for all }0\le a,b<p\text{ with }a+b>d.
\]
Thus a nonzero $\fall{x}{j}$-coefficient at an index $j<p^2$ can occur only when
the quotient and remainder of $j$ upon division by $p$ sum to at most $d$.
\end{lemma}
\begin{proof}
Put $y=x^p-x$ and write $H=Q+yT$. We use the \emph{Frobenius
endomorphism} of $\Fp[x]$, namely the map $U\mapsto U^p$.
This map preserves sums and products: in particular,
$(U+V)^p=U^p+V^p$, by the binomial theorem in characteristic $p$.
Also, it fixes every coefficient in $\Fp$, so $Q(x)^p=Q(x^p)$.
Thus, from $H = Q+yT$, we obtain
\begin{equation}\label{eq:frobenius}
 H(x)^p= Q(x)^p + y^p T(x)^p = Q(x^p) + y^p T(x)^p
 = Q(x+y)+y^pT(x)^p \qquad
\end{equation}
(since $x^p = x+y$).
We examine the falling-factorial coefficients of the two terms
separately; each time we will show that the $\fall{x}{ap+b}$-coefficient
is $0$ whenever $0\leq a,b<p$ with $a+b>d$.

For the second term, apply the coefficient-shift identity
\eqref{eq:falling-coefficient-shifts} with $a=p$ and $P=T(x)^p$.
It gives
\[
 [\fall{x}{i}]\bigl(y^pT(x)^p\bigr)
   =[\fall{x}{i-p^2}]\bigl(T(x)^p\bigr)=0
 \qquad\text{for all }0\le i<p^2.
\]
The last equality holds because $i-p^2<0$. Thus the second term
contributes nothing at any of the indices under consideration.

It remains to examine $Q(x+y)$. Expand it first as a polynomial in
two formal variables $x$ and $y$, before making the substitution
$y=x^p-x$. Because $\deg Q\le d$, every term is a multiple of
$x^by^a$ with $a+b\le d$. Express $x^b$ as a linear combination
of $\fall{x}{t}$ for $0\le t\le b$. Multiplication by $y^a$ then
turns these into $\fall{x}{ap+t}$, by~\eqref{eq:block-factorials}.
Thus, $Q(x+y)$ is expanded as an $\Fp$-linear combination of
falling factorials $\fall{x}{ap+t}$ whose indices $ap+t$ satisfy
\[
 %ap+t\qquad\text{with }
 a,t\ge0\text{ and }a+t\le d.
\]
Since $d<p$, we have $a,t<p$, so $a$ and $t$ are precisely the
quotient and remainder when the index is divided by $p$.
This proves the claimed vanishing.
\end{proof}

The next lemma extracts the divisibility we need from just one
condition: the sum of $P$ over $p$ consecutive arguments is divisible
by $p$.

\begin{lemma}\label{lem:polynomial}
Suppose $P\in R[x]$ satisfies
\begin{equation}\label{eq:sum-condition}
 \sum_{r=0}^{p-1}P(r)\in pR.
\end{equation}
Then
\[
 \Delta^{k(p-1)}(P^p)(0)\in p^kR
 \qquad\text{for all }1\le k\le p.
\]
\end{lemma}
\begin{proof}
Write $P(x)^p=\sum_{j\ge0}c_j\fall{x}{j}$ with $c_j\in R$.
For $m=k(p-1)$ and $1\le k\le p$, we have $m<p^2$ and
thus, by~\eqref{eq:legendre}, we have
\begin{equation}\label{eq:factorial-valuation}
 v_p(m!)=\left\lfloor\frac mp\right\rfloor=k-1.
\end{equation}
(For the last equality, note that $m/p=k-k/p$. Since
$1\le k\le p$, we have $0<k/p\le1$, and hence
$k-1\le m/p<k$. Thus $\lfloor m/p\rfloor=k-1$.)
By~\eqref{eq:newton-coefficient}, the difference we want to study is
\[
 \Delta^m(P^p)(0)=m!c_m.
\]
We need this product to be divisible by $p^k$. Formula
\eqref{eq:factorial-valuation} says that $m!$ already supplies exactly
$k-1$ factors of $p$. Thus we need just one more factor from $c_m$:
our remaining task is to show $c_m\in pR$ for $m=k(p-1)$.

This is a question modulo $p$. Fractions in $R$ can be reduced modulo
$p$ by inverting their denominators in $\Fp$ (there is a canonical ring
isomorphism $R/pR \cong \Z/p\Z = \Fp$), so we may now work in $\Fp[x]$.
Here $\Fp$ is the field of residues modulo $p$, and $\overline P$
denotes the reduction of $P$. The coefficients in the falling-factorial
expansion also reduce modulo $p$: reducing the identity
$P(x)^p=\sum_{j\ge0}c_j\fall{x}{j}$ modulo $p$
gives $\overline P^{\,p} = \sum_{j\ge0} \overline{c_j} \fall{x}{j}$,
so that the falling-factorial coefficients of $\overline P^{\,p}$ are
the $\overline{c_j}$.
That is, $\overline{c_j} = [\fall{x}{j}](\overline P^{\,p})$ for each
$j\geq 0$.

We first replace $\overline P$ by a polynomial of degree less than $p$
that takes the same values on $\Fp$. Put $y=x^p-x$ and divide
$\overline P$ by this monic polynomial $y$ to obtain
\[
 \overline P=Q+yT,\qquad \text{ where } Q,T\in\Fp[x] \text{ with } \deg Q<p.
\]
Fermat's little theorem says that $y(r)=0$ for every $r\in\Fp$.
Hence $Q(r)=\overline P(r)$ for every $r \in \Fp$,
and thus our hypothesis \eqref{eq:sum-condition} rewrites as
\begin{equation}
 \sum_{r\in\Fp}Q(r)=0.
 \label{eq:sum-condition-Q}
\end{equation}
This forces the coefficient of $x^{p-1}$ in $Q$ to vanish.
The reason is the famous power-sum identity
\begin{equation}\label{eq:residue-sums}
 \sum_{r\in\Fp}r^b=
 \begin{cases}
 0,&\text{if } 0\le b<p-1,\\
 -1,&\text{if } b=p-1
 \end{cases}
 \qquad \text{ for every } 0 \leq b < p,
\end{equation}
which we prove in Appendix~\ref{app:power-sums}.
Applying \eqref{eq:residue-sums} term by term to $Q$ shows that
$\sum_{r\in\Fp} Q(r)$ is the negative of the $x^{p-1}$-coefficient
of $Q$ (since $\deg Q < p$).
Thus, \eqref{eq:sum-condition-Q} tells us that the latter
coefficient is $0$. In other words,
\[
 \deg Q\le p-2.
\]

We can therefore apply Lemma~\ref{lem:frobenius-support} to
$H=\overline P$ and $d=p-2$. It says that the falling-factorial
coefficient $[\fall{x}{ap+b}]\left(\overline P^{\,p}\right)$
at an index $ap+b$ vanishes whenever $0 \leq a,b < p$ and
$a+b>p-2$.
Our desired indices $k(p-1)$ (for $1 \leq k \leq p$) have
exactly this property, since they can be written as
\[
 k(p-1)=ap+b
 \qquad\text{for } a = k-1 \text{ and } b = p-k
\]
(and these values $a$ and $b$ satisfy $0\le a<p$ and $0\le b<p$ and
$a + b = (k-1) + (p-k) = p-1 > p-2$). Consequently,
$[\fall{x}{k(p-1)}]\left(\overline P^{\,p}\right)$ vanishes.
This says exactly that $c_{k(p-1)}\in pR$, supplying the one factor
of $p$ still needed to complete the proof.
\end{proof}

\begin{remark}
For example, when $p=5$ and $\deg Q\le3$, the possible indices less
than $25$ in the falling-factorial expansion of $\overline P^{\,5}$ are
\[
 0,1,2,3,\quad 5,6,7,\quad 10,11,\quad 15.
\]
None is one of $4,8,12,16,20$, the indices needed for the first five
averaging depths. The factorials at these latter indices have respective
$5$-adic valuations $0,1,2,3,4$; the lemma supplies one further factor of $5$.
\end{remark}

\section{From local interpolation to arbitrary functions}\label{sec:interpolation}

Recall that our interest lies in functions $f:\Z\to\QQ$, not
in polynomials.
However, a value of the $m$-th finite difference $\Delta^m f$
depends on only $m+1$ values of $f$.
Thus, even though our function $f$ need not itself be a polynomial, we can
replace it by a polynomial separately on each interval of the required
length. To be precise, expanding
\[
\Delta^m = (E-\id)^m = \sum_{r=0}^m(-1)^{m-r}\binom mrE^r
\]
gives
\begin{equation}\label{eq:finite-difference-values}
 \Delta^mf(n)=\sum_{r=0}^m(-1)^{m-r}\binom mr f(n+r).
\end{equation}
This formula makes precise how an interval of agreement shrinks
when we take differences. Suppose two functions $f,g:\Z\to\QQ$
agree at every integer in an interval $[a,b]$, where $a\le b$ are
integers. Then
\begin{equation}\label{eq:local-difference-agreement}
 \Delta^j f(t)=\Delta^j g(t)
 \qquad\text{for all }0\le j\le b-a
 \text{ and }a\le t\le b-j.
\end{equation}
Indeed, the formula for $\Delta^j f(t)$ uses only the values at
$t,t+1,\ldots,t+j$. All these arguments lie in $[a,b]$, so each
value of $f$ in the sum can be replaced by the corresponding value
of $g$. Thus each further application of $\Delta$
loses the rightmost point of the interval of agreement, and
$j$ iterations of $\Delta$ lose its last $j$ points.

In particular, suppose a polynomial $P \in \QQ[x]$ satisfies $P(r)=f(n+r)$
for all $0\le r\le m$ (where $n$ and $m\geq 0$ are two fixed integers).
Applying the same observation to the functions
$r\mapsto P(r)$ and $r\mapsto f(n+r)$ gives
\begin{equation}\label{eq:interpolation-differences}
 (\Delta^jP)(r)=(\Delta^jf)(n+r)
 \qquad\text{for all }0\le j\le m\text{ and }0\le r\le m-j.
\end{equation}
Here polynomial differences agree with differences of the evaluated
function by~\eqref{eq:evaluation-commutes}. For $j=m$, only $r=0$
remains: $\Delta^mP(0)=\Delta^mf(n)$. The same reasoning applies
to $P^p$ and $f^p$, since taking pointwise powers preserves the
original equalities of values.

Newton interpolation constructs the polynomial we need directly
from these differences. Recall that for each $j \geq 0$,
the polynomial $\binom{x}{j}=\fall{x}{j}/j! \in \QQ[x]$ is
integer-valued\footnote{Indeed, its values at nonnegative
integers are the standard binomial coefficients, whereas its
values at negative integers can be reduced to the nonnegative
ones via the upper negation formula
$\binom{-a}{j}=(-1)^j\binom{a+j-1}{j}$.}.
For any function $f : \Z \to \QQ$, the polynomial
\begin{equation}
 P(x)=\sum_{j=0}^m\Delta^jf(n)\binom{x}{j}
 \label{eq:newtonP}
\end{equation}
(known as the \emph{$m$-th Newton interpolant} of $f$
based at $n$)
agrees with $f(n+x)$ at the integers $x=0,\ldots,m$.
For example, to match three values $u=f(n)$, $v=f(n+1)$, and
$w=f(n+2)$, the formula gives
\[
 P(x)=u+(v-u)x+(w-2v+u)\frac{x(x-1)}2.
\]
Substitution of $x=0,1,2$ gives $u,v,w$, respectively. In
the general case, the formula
\begin{equation}
 P(r) = f(n+r) \qquad \text{ for all } 0 \leq r \leq m
 \label{eq:newtonPagree}
\end{equation}
is verified by the same binomial theorem we used for the
operators: since $E=\id+\Delta$, for $0\le r\le m$ we have
\[
 f(n+r)=(E^rf)(n)
 =((\id+\Delta)^r f)(n)
 =\sum_{j=0}^r\binom rj\Delta^jf(n)=P(r)
\]
(here, all addends with $j>r$ in $P(r)$ vanish because $\binom rj=0$).
We apply this formula using the divisibility supplied by
Lemma~\ref{lem:difference}.

\begin{proof}[Proof of Theorem~\ref{thm:main}]
The case $\ell=0$ says only that $f^p$ is integral, which is clear
since $f$ is integral.
Suppose $1\le\ell\le p$.
Set $q = p-1$.
By Lemma~\ref{lem:difference} (equivalence
of the first two conditions, applied to $f^p$ instead of $f$),
our goal is to show that $\Delta^{kq}(f^p) \in p^k \Z^\Z$ for all
$0\le k \le \ell$.
In other words, we must prove that
\begin{equation}
 p^k\mid\Delta^{kq}(f^p)(n)
 \qquad\text{for all }n\in\Z\text{ and }0\le k\le\ell.
 \label{pf:thm:main:goal}
\end{equation}
Fix such $n$ and $k$, and put $m=kq$; thus,
$m = kq \leq \ell q\le pq = p(p-1) < p^2$.
We assume that $k\geq 1$, since the case $k=0$ is obvious.
Define the polynomial
\begin{equation}\label{eq:interpolant}
 P(x)=\sum_{j=0}^{m}\Delta^j f(n)\binom{x}{j}
     =\sum_{j=0}^{m}\frac{\Delta^j f(n)}{j!}\fall{x}{j}.
\end{equation}
Newton interpolation -- i.e., the formula
\eqref{eq:newtonPagree} -- gives
\begin{equation}\label{eq:interpolation-values}
 P(r)=f(n+r)\qquad\text{for all }0\le r\le m.
\end{equation}

To apply Lemma~\ref{lem:polynomial}, we must check its two
hypotheses. First, the coefficients of $P$ must belong to $R$.
This is not automatic from the interpolation formula: division by
$j!$ might put $p$ in a denominator. However,
the difference bounds prevent
this. Indeed, Lemma~\ref{lem:difference} (this time applied
to $f$ itself) yields that each $j \leq m$ satisfies
$\Delta^j f \in p^{\lfloor j/(p-1)\rfloor}\Z^\Z$
(since $j \leq m \leq \ell q$), so that
$v_p(\Delta^j f(n)) \ge \lfloor j/(p-1)\rfloor$. Since
$j\le m < p^2$, the formula \eqref{eq:legendre} yields
\[
 v_p(j!)=\left\lfloor\frac jp\right\rfloor
 \le\left\lfloor\frac{j}{p-1}\right\rfloor.
\]
The numerator $\Delta^jf(n)$ in \eqref{eq:interpolant}
therefore contains at least as many
factors of $p$ as the denominator $j!$. After cancellation, the
$\fall{x}{j}$-coefficient therefore belongs to $R$. Since each
falling factorial has integer coefficients, this proves $P\in R[x]$.

Second, the sum of $P$ over $0,\ldots,p-1$ must belong to $pR$.
Indeed, since $k \geq 1$, we have $m = kq \ge q = p-1$, so that each
$r \in \left\{0,1,\ldots,p-1\right\}$ satisfies $r \leq m$
and therefore $P(r) = f(n+r)$ by~\eqref{eq:interpolation-values}.
Hence,
\[
 \sum_{r=0}^{p-1}P(r)
 =\sum_{r=0}^{p-1}f(n+r)=p(F_pf)(n)\in p\Z
\]
(since $1 \leq \ell$ and the $\ell$-average-integrality of $f$
entail that $F_pf$ is integral, so $(F_pf)(n) \in \Z$).
Thus, $\sum_{r=0}^{p-1}P(r) \in p\Z \subseteq pR$.
Lemma~\ref{lem:polynomial} consequently yields
$\Delta^{k(p-1)}(P^p)(0) \in p^k R$.
In view of $m = kq = k(p-1)$, this rewrites as
$\Delta^{m}(P^p)(0) \in p^k R$.
But \eqref{eq:interpolation-differences}
(applied to $P^p$, $f^p$, $m$ and $0$
instead of $P$, $f$, $j$ and $r$) yields
$\Delta^m(P^p)(0) = \Delta^m(f^p)(n)$.
Hence,
\[
 \Delta^m(f^p)(n) = \Delta^m(P^p)(0) \in p^kR.
\]
The left side is an integer, so it is therefore
divisible by $p^k$. In other words,
$p^k \mid \Delta^m(f^p)(n) = \Delta^{kq}(f^p)(n)$
(since $m = kq$).
This proves \eqref{pf:thm:main:goal}.
As $n$ and $k$ were arbitrary,
Lemma~\ref{lem:difference} proves that $f^p$ is
$\ell$-average-integral.
Theorem~\ref{thm:main} is thus proved.
\end{proof}

\section{Polynomials and repeated averaging}
The finite-depth theorem is now proved. To study larger depths, it
helps to answer a more specific question: which polynomials remain
integer-valued after arbitrarily many averaging steps?

A monomial $x^t$ need not have this property, but multiplying it by a
suitable power of $p$ will make it work. The exponent that appears is
\[
 d(t)=\left\lfloor\frac{t}{p-1}\right\rfloor-v_p(t!).
\]
There is a reason for this expression. The difference criterion
(Lemma~\ref{lem:difference}) asks
for $\lfloor t/(p-1)\rfloor$ factors of $p$ in a $t$th difference,
while the $t$th difference of $x^t$ is $t!$.
Thus $d(t)$ counts the
additional factors needed at this highest difference order.
The next lemma shows that these factors also suffice for \textbf{every}
order. We then obtain the polynomial criterion from zeb's
answer~\cite{zeb}, and use it both for $p=2,3$ and for the counterexample
when $p>3$.

\begin{lemma}\label{lem:monomials}
For $t\ge0$, put
\[
 d(t)=\left\lfloor\frac{t}{p-1}\right\rfloor-v_p(t!).
\]
Then:
\begin{enumerate}
\item[(a)] We have $d(t)\ge0$ for all $t \geq 0$.
\item[(b)] We have $d(pt)=d(t)$ for all $t\ge0$.
\item[(c)] We have $d(t) \geq 1$ for all $t > 0$ satisfying $p-1 \mid t$.
\item[(d)] For any $t \geq 0$, the function $n\mapsto p^{d(t)}n^t$ is average-integral.
\end{enumerate}
\end{lemma}
\begin{proof}
For $t>0$, Legendre's formula~\eqref{eq:legendre} gives
\[
 v_p(t!)=\sum_{i\ge1}\left\lfloor\frac{t}{p^i}\right\rfloor
 <\sum_{i\ge1}\frac{t}{p^i}=\frac{t}{p-1}.
\]
This implies part (c) (since $\left\lfloor\frac{t}{p-1}\right\rfloor
= \frac{t}{p-1}$ when $p-1\mid t$) and also
part (a) (since $v_p(t!)$ is an integer, so that
$v_p(t!) < \frac{t}{p-1}$ entails
$v_p(t!) \le \left\lfloor\frac{t}{p-1}\right\rfloor$).
% \[
% d(t)\ge0,
% \]
% and even
% \begin{align}
% d(t)\ge1 \qquad \text{ when } t>0\text{ and }p-1\mid t.
% \label{eq:dtge1}
% \end{align}
Also $d(0)=0$.
The same formula~\eqref{eq:legendre} yields
\[
 v_p((pt)!)=t+v_p(t!),\qquad
 \left\lfloor\frac{pt}{p-1}\right\rfloor
 =\left\lfloor t + \frac{t}{p-1}\right\rfloor
 =t+\left\lfloor\frac{t}{p-1}\right\rfloor,
\]
and hence $d(pt) = d(t)$ for all $t \geq 0$;
thus, part (b) is proved.

It remains to prove part (d).
First, we shall work with the polynomials
$p^{d(t)} x^t$ for $t \ge 0$; then we will evaluate
them at integers.

Define an $R$-submodule $M_p$ of $R[x]$ by
\begin{align}
 M_p=\left\{\sum_{t=0}^N c_tx^t:
       N\ge0,\ c_t\in p^{d(t)}R\text{ for all }0\le t\le N\right\}.
\label{eq:Mp=}
\end{align}
In words, a polynomial belongs to $M_p$ when its $x^t$-coefficient
for each $t$ contains at least $d(t)$ factors of $p$.
Equivalently, the elements of $M_p$ are the finite $R$-linear combinations
of $p^{d(t)}x^t$. We will prove that averaging
(i.e., the operator $F_p$ on polynomials) preserves $M_p$.
It suffices to check that $F_p(p^{d(t)}x^t) \in M_p$ for all $t \geq 0$,
since averaging is linear.

For $t=0$ the assertion is immediate: averaging fixes constants.
For $t>0$, start with the definition of $F_p$. In each summand,
expand $(x+r)^t$ by the binomial theorem. This gives
\begin{align*}
 F_p(p^{d(t)}x^t)
 &=\frac{p^{d(t)}}p\sum_{r=0}^{p-1}(x+r)^t\\
 &=p^{d(t)-1}\sum_{r=0}^{p-1}\sum_{j=0}^t
       \binom tj x^j r^{t-j}\\
 &=\sum_{j=0}^t p^{d(t)-1}\binom tj
       \left(\sum_{r=0}^{p-1}r^{t-j}\right)x^j.
\end{align*}
Write $S_u=\sum_{r=0}^{p-1}r^u$. The term with $j=t$ has
$S_0=p$ and $\binom tt=1$, so it is $p^{d(t)}x^t$.
Separating this term, we obtain
\[
 F_p(p^{d(t)}x^t)
 =p^{d(t)}x^t+
   \sum_{j=0}^{t-1}p^{d(t)-1}\binom tj S_{t-j}x^j.
\]
The leading term of this polynomial is $p^{d(t)}x^t$, which
is already in $M_p$. Fix a lower degree $j<t$
and write $u=t-j>0$. Then, the coefficient of $x^j$ in
$F_p(p^{d(t)}x^t)$ is
\begin{align}
 C_j=p^{d(t)-1}\binom tj S_u.
 \label{eq:Cj=}
\end{align}
We must show $C_j \in p^{d(j)}R$, that is, $v_p(C_j)\ge d(j)$.
Two estimates give this.

First, the binomial coefficient accounts for the factorials in $d$.
Since $t=j+u$, the factorial formula for $\binom tj$ gives
\[
 v_p\binom tj=v_p(t!)-v_p(j!)-v_p(u!).
\]
Adding $d(t)=\lfloor t/(p-1)\rfloor-v_p(t!)$ cancels
$v_p(t!)$, leaving the exact equality
\[
 d(t)+v_p\binom tj
 =\left\lfloor\frac{t}{p-1}\right\rfloor-v_p(j!)-v_p(u!).
\]
Since $t=j+u$, we also have
\[
 \left\lfloor\frac{t}{p-1}\right\rfloor
 = \left\lfloor\frac{j+u}{p-1}\right\rfloor
 \ge \left\lfloor\frac{j}{p-1}\right\rfloor
     +\left\lfloor\frac{u}{p-1}\right\rfloor,
\]
since all $x,y \in \mathbb R$ satisfy $\left\lfloor x+y\right\rfloor
\ge \left\lfloor x\right\rfloor + \left\lfloor y\right\rfloor$.
Subtracting $v_p(j!)+v_p(u!)$ from both sides yields
\begin{equation}\label{eq:d-binomial-bound}
 d(t)+v_p\binom tj\ge d(j)+d(u).
\end{equation}
Second, we claim that
\begin{equation}\label{eq:d-Su-bound}
 d(u)+v_p(S_u)\ge1.
\end{equation}
Indeed, if $p-1$ divides $u$, then we already proved $d(u)\ge1$
in Lemma~\ref{lem:monomials}(c) (since $u>0$), and so
\eqref{eq:d-Su-bound} follows.
If $p-1$ does not divide $u$, then $S_u$ is divisible by $p$
(an easy consequence of \eqref{eq:residue-sums}\footnote{In
more detail: Assume that $p-1$ does not divide $u$. Then,
$u$ has a nonzero remainder $w$ upon division by $p-1$.
Thus, $0 < w < p-1$ and $u \equiv w \mod p-1$. Hence,
Fermat's Little Theorem yields the congruence
$r^u \equiv r^w \mod p$
for each integer $r$ not divisible by $p$; but the
same congruence holds for all integers $r$ divisible by $p$
because $u$ and $w$ are positive.
Summing this congruence over all $r \in \{0,1,\ldots,p-1\}$,
we find $S_u \equiv S_w \mod p$.
But \eqref{eq:residue-sums} shows that the sum $S_w$ is
divisible by $p$, since $0 < w < p-1$; thus,
$S_u \equiv S_w \equiv 0 \mod p$.}), thus $v_p(S_u)\ge 1$,
so that \eqref{eq:d-Su-bound} follows again.
In either case, \eqref{eq:d-Su-bound} is proved.

Now, using \eqref{eq:Cj=}, \eqref{eq:d-binomial-bound}
and \eqref{eq:d-Su-bound} in sequence, we find
\[
 v_p(C_j)
 =d(t)-1+v_p\binom tj+v_p(S_u)
 \ge d(j)+d(u)-1+v_p(S_u)
 \ge d(j),
\]
so $C_j \in p^{d(j)}R$.
That is, the $x^j$-coefficient of $F_p(p^{d(t)}x^t)$ is
divisible by $p^{d(j)}$ in $R$.
Since this holds for all $j$, we thus have
$F_p(p^{d(t)}x^t) \in M_p$.
Since this holds for all $t \geq 0$, we thus obtain
$F_pM_p\subseteq M_p$.

Now, fix $t\ge0$ and let
$f:\Z\to\Z$ be the function given by $f(n)=p^{d(t)}n^t$.
We must show that $f$ is average-integral,
i.e., that $F_p^k f$ is integral for every $k\ge0$.
Since $p^{d(t)}x^t\in M_p$
and $F_pM_p\subseteq M_p$, induction on $k$ shows that
\[
F_p^k(p^{d(t)}x^t)\in M_p\subseteq R[x]
\qquad \text{ for every }k\ge0.
\]
Compatibility with evaluation, \eqref{eq:evaluation-commutes},
therefore gives
\begin{align}
 (F_p^k f)(n)\in R
 \qquad\text{for all }k\ge0\text{ and }n\in\Z
\label{eq:FpknfinR}
\end{align}
(since $f = \operatorname{ev}(p^{d(t)}x^t)$
and thus $F_p^k f = \operatorname{ev}(F_p^k(p^{d(t)}x^t))$).
On the other hand, $f$ is integral and each application of $F_p$
adds $p$ values and divides their sum by $p$. Inductively, this
shows that
$p^k(F_p^k f)(n)$ is an integer. Thus $(F_p^k f)(n)$ belongs to
$\Z[1/p]$, the ring of rational numbers whose reduced denominators
are powers of $p$.
Combining this with \eqref{eq:FpknfinR}, we see that
$(F_p^k f)(n) \in R \cap \Z[1/p] = \Z$
(indeed, a rational number in $\Z[1/p]$ can have only $p$ in its
denominator, whereas a rational number in $R$ cannot have $p$
there; thus, a rational number in $R \cap \Z[1/p]$ must have
denominator $1$).
That is, $(F_p^k f)(n)$ is an integer for every $k\ge0$ and
$n\in\Z$. This proves that $f$ is average-integral.
Thus, part (d) is proved.
\end{proof}


Let us next record the analogue of Lemma~\ref{lem:difference} for
$R$-valued functions. As before, $p$ is our fixed prime and $q=p-1$.

\begin{lemma}[Difference criterion over $R$]\label{lem:difference-R}
For a function $f\colon\Z\to\QQ$ and a nonnegative integer $\ell$,
the following conditions are equivalent:
\begin{enumerate}
\item We have $F_p^k f\in R^\Z$ for all $0\le k\le\ell$.
\item We have $\Delta^{kq}f\in p^kR^\Z$ for all $0\le k\le\ell$.
\item We have $\Delta^j f\in p^{\lfloor j/q\rfloor}R^\Z$
for all $0\le j\le\ell q$.
\end{enumerate}
\end{lemma}
\begin{proof}
This is an analogue of Lemma~\ref{lem:difference} with some $\Z$'s
replaced by $R$'s; the exact same proof applies, mutatis mutandis.
(The only thing we use about $R$ is that $R$ is a subring of $\QQ$.)
% Use the operators $D=\Delta^q/p$ and $B(E)$ from
% \eqref{eq:D=} and~\eqref{eq:Fp=D+BE}. Since $B(E)$ is an integer
% linear combination of shifts, it preserves $R^\Z$. It also commutes
% with $D$ and $F_p$, since all three are polynomials in $E$.
% Consequently, the binomial expansions give
% \[
 % D^kf=\sum_{i=0}^k\binom ki(-B(E))^{k-i}F_p^if,
 % \qquad
 % F_p^kf=\sum_{i=0}^k\binom ki B(E)^{k-i}D^if.
% \]
% If all the functions $F_p^if$ for $0\le i\le k$ are $R$-valued,
% then every summand in the first expansion is $R$-valued, and hence
% so is $D^kf$. The second expansion gives the converse implication.
% Since $D^kf=\Delta^{kq}f/p^k$, this proves (1)$\Leftrightarrow$(2).
%
% To deduce (3) from (2), write $j=kq+s$ with $0\le s<q$.
% Then $k=\lfloor j/q\rfloor$ and $\Delta^{kq}f=p^kg$ for some
% $g\in R^\Z$. Taking $s$ further differences gives
% $\Delta^jf=p^k\Delta^sg\in p^kR^\Z$, since differences preserve
% $R$-valued functions. Conversely, taking $j=kq$ in (3) gives (2).
\end{proof}

We can now describe all the polynomials $P$ for which every iterated
average $F_p^k P$ takes values in $R$ at all integers.
The description has several useful forms: one bounds the
ordinary monomial coefficients, and the other bounds the coefficients
in the binomial (or, equivalently, falling-factorial) basis.
We will use the $R$-submodule
\begin{align}
 M_p=\left\{\sum_{t=0}^N c_tx^t:
       N\ge0,\ c_t\in p^{d(t)}R\text{ for all }0\le t\le N\right\}
\label{eq:Mp=2}
\end{align}
of $R[x]$ (which already appeared in the proof of
Lemma~\ref{lem:monomials}).
%defined in \eqref{eq:Mp=}.
% In the notation below, membership in $M_p$
% means simply that the coefficient of $x^j$ belongs to $p^{d(j)}R$
% for every $j$.

We recall that each polynomial $P(x) \in \QQ[x]$ can be written as a
$\QQ$-linear combination $\sum_j b_j \fall{x}{j}$ of falling
factorials. Since $\fall{x}{j} = j! \binom{x}{j}$, we can rewrite
this as $\sum_j b_j j! \binom{x}{j}$; so every polynomial
$P(x) \in \QQ[x]$ can be written as a
$\QQ$-linear combination $\sum_j a_j \binom{x}{j}$ of
binomial coefficients. Again, the coefficients $a_j$ here are
unique.

\Needspace{12\baselineskip}
\begin{proposition}[Polynomial criterion]\label{prop:criterion}
For $P\in\QQ[x]$, the following conditions are equivalent:
\begin{enumerate}
\item We have $F_p^kP(n)\in R$ for every $k\ge0$ and $n\in\Z$.
\item Writing $P(x)=\sum_j c_jx^j$, we have $c_j\in p^{d(j)}R$
for every $j$; in other words, $P\in M_p$.
\item Writing $P(x)=\sum_j b_j\fall{x}{j}$, we have
$b_j\in p^{d(j)}R$ for every $j$.
\item Writing $P(x)=\sum_j a_j\binom{x}{j}$, we have
$a_j\in p^{\lfloor j/(p-1)\rfloor}R$ for every $j$.
\item We have $\Delta^{r(p-1)} P(n) \in p^r R$ for each $n \in \Z$ and $r \geq 0$.
\item We have $\Delta^j P(n) \in p^{\lfloor j/(p-1) \rfloor}R$ for each $n \in \Z$ and $j \geq 0$.
\end{enumerate}
If $P$ is integer-valued on $\Z$, these conditions are also equivalent
to average-integrality of the function $n\mapsto P(n)$.
\end{proposition}

\begin{proof}
We first compare the conditions involving values and differences.
Put $f(n)=P(n)$ for all $n\in\Z$.
Apply Lemma~\ref{lem:difference-R} to $f$,
using the compatibility with evaluation in~\eqref{eq:evaluation-commutes}.
Allowing $\ell$ to range over all nonnegative integers proves
(1)$\Leftrightarrow$(5)$\Leftrightarrow$(6).

Next we prove (6)$\Rightarrow$(4)$\Rightarrow$(5).
Since $\binom{x}{j}=\fall{x}{j}/j!$, the expansion
$P=\sum_j a_j\binom{x}{j}$ can also be written as
\[
 P=\sum_j\frac{a_j}{j!}\fall{x}{j}.
\]
Thus \eqref{eq:newton-coefficient} gives
$\Delta^jP(0)=j!(a_j/j!)=a_j$.
Applying condition (6) at $n=0$, we thus obtain
$a_j\in p^{\lfloor j/(p-1)\rfloor}R$, which is (4).

To prove (4)$\Rightarrow$(5), divide~\eqref{eq:falling-difference} by $j!$
to obtain the normalized difference identity
\[
 \Delta\binom{x}{j}=\binom{x}{j-1}\quad\text{for all }j\ge1,
 \qquad \Delta\binom{x}{0}=0.
\]
So each iteration of $\Delta$, when applied to a binomial-coefficient
polynomial $\binom{x}{j}$, lowers the binomial index $j$ by $1$,
until the term vanishes.
So $\Delta^k \binom{x}{j} = \binom{x}{j-k}$ for all $j,k\geq 0$,
where we interpret the right-hand side as $0$ if $j-k<0$.
Thus, under condition (4), for every $r\ge0$ we have
\begin{equation}\label{eq:divided-newton}
 \Delta^{r(p-1)}P(x)
 =\sum_{j\ge r(p-1)}a_j
       \binom{x}{j-r(p-1)}.
\end{equation}
For each index in this sum, $j\ge r(p-1)$, so
$\lfloor j/(p-1)\rfloor\ge r$. Condition (4) therefore ensures
$a_j \in p^rR$. Also, the binomial polynomials $\dbinom{x}{k}$
take integer values on all of $\Z$ (as we have already seen).
Consequently, \eqref{eq:divided-newton} shows that
$\Delta^{r(p-1)}P$ takes values in $p^r R$ for every
$r$, which is exactly condition (5).
This completes the equivalence of conditions
(1), (4), (5), and (6).

The equivalence of (3) and (4) is also easy, since the respective
coefficients are related by $a_j = j! b_j$
(because $\fall{x}{j} = j! \binom{x}{j}$), and thus we have
$a_j \in p^{\lfloor j/(p-1)\rfloor} R$ if and only if
$b_j \in p^{d(j)}R$.

Next we prove (2)$\Rightarrow$(1). Under (2), write
\[
 P(x)=\sum_{j=0}^N u_j p^{d(j)}x^j,\qquad \text{ with } u_j\in R.
\]
For each $j$, put $f_j(n)=p^{d(j)}n^j$. Lemma~\ref{lem:monomials}(d)
says that $f_j$ is average-integral, so that
$F_p^k f_j$ is integral for every $k\ge0$.
Since $\operatorname{ev} P=\sum_{j=0}^N u_j f_j$ and $F_p^k$ is $\QQ$-linear
(and by \eqref{eq:evaluation-commutes}),
\[
 (F_p^k P)(n)=\sum_{j=0}^N \underbrace{u_j}_{\in R}
                \underbrace{(F_p^k f_j)(n)}_{\in \Z} \in R
 \qquad\text{for all }k\ge0\text{ and }n\in\Z.
\]
This is condition (1).

For the converse, we take a bird's-eye view.
Do not fix $P$.
As we recall, $M_p$ is the set of all polynomials
$P \in \QQ[x]$ satisfying condition (2).
Let $M'_p$ be the set of all polynomials
$P \in \QQ[x]$ satisfying condition (1).
Both $M_p$ and $M'_p$ are $R$-submodules of $\QQ[x]$
(since each $F_p^k$ is $R$-linear).
Since we have proved (2)$\Rightarrow$(1), we know
that $M_p \subseteq M'_p$. We shall now show that
$M'_p \subseteq M_p$.

Indeed, let $P \in M'_p$. We must show that $P \in M_p$.
We shall show this by strong induction on $\deg P$;
the base case $P = 0$ is obvious.
Suppose $P\ne0$ and let $c_tx^t$ be its leading term.
The polynomial $P$ satisfies condition (1) (since
$P \in M'_p$), and thus also condition (6)
(since (1) and (6) are equivalent),
which at $j=t$ and $n=0$ gives
\[
 \Delta^tP(0)\in p^{\lfloor t/(p-1)\rfloor}R.
\]
Since $\fall{x}{t}$ is monic of degree $t$, the polynomial
$P-c_t\fall{x}{t}$ has degree less than $t$, and thus
is a linear combination of $\fall{x}{0},\ \fall{x}{1},\ \ldots,
\fall{x}{t-1}$. Thus the
$\fall{x}{t}$-coefficient of $P$ is $c_t$, and so
\eqref{eq:newton-coefficient} gives $\Delta^tP(0)=t!c_t$.
Therefore,
\[
 t!c_t=\Delta^tP(0)\in p^{\lfloor t/(p-1)\rfloor}R,
 \qquad\text{hence}\qquad c_t\in p^{d(t)}R.
\]
Therefore, $c_t x^t \in M_p \subseteq M'_p$.
Hence, $P - c_t x^t \in M'_p$ as well
(because $P$ and $c_t x^t$ belong to $M'_p$, and since
$M'_p$ is an $R$-module). By the induction hypothesis, this
entails $P - c_t x^t \in M_p$, since $P - c_t x^t$ has
smaller degree than $P$.
Adding this to $c_t x^t \in M_p$, we find $P \in M_p$,
since $M_p$ is an $R$-module.
So the induction is complete, and $M'_p \subseteq M_p$
is proved; that is, condition (1) implies condition (2).


% Write $c_t=p^{d(t)}u_t$ with $u_t\in R$, and define a function
% $h : \Z \to \QQ$ by $h(n)=c_tn^t$.
% By Lemma~\ref{lem:monomials}(d) and linearity,
% $(F_p^k h)(n)=u_t(F_p^k f_t)(n)\in R$ for every $k\ge0$ and $n\in\Z$,
% where $f_t(n)=p^{d(t)}n^t$.
% Thus the leading term $c_tx^t$ of $P$ has property (1). Subtracting
% it from $P$
% leaves a polynomial of smaller degree, whose associated function is
% $f-h$. Since $F_p^k(f-h)=F_p^k f-F_p^k h$, this polynomial still
% has property (1). Repeating the
% argument proves the required divisibility for every coefficient,
% which is (2). The zero polynomial presents no exception.

Finally, suppose $P$ is integer-valued, so $f$ is integral.
For every $k\ge0$ and $n\in\Z$, the number $p^k(F_p^k f)(n)$
is an integer, so $(F_p^k f)(n)$ is a rational number whose
denominator is a power of $p$. Under condition (1), the number
$(F_p^k f)(n)$ also belongs to $R$,
so the denominator argument at the end of the proof of
Lemma~\ref{lem:monomials}(d) shows that
it is an integer. Hence $f$ is average-integral. Conversely, if
$f$ is average-integral, then $(F_p^k f)(n)\in\Z\subseteq R$
for every $k,n$, which is exactly condition (1).
\end{proof}

For example, when $p=3$ we have $d(0)=d(1)=0$, $d(2)=1$, and
$d(3)=0$. Thus a polynomial of degree at most $3$ in $R[x]$
has all its averages valued in $R$ exactly when its $x^2$-coefficient
is divisible by $3$. In particular, $x$, $3x^2$, and $x^3$ are
average-integral on $\Z$, but $x^2$ is not. Indeed,
\[
 F_3(x^2)=x^2+2x+\frac53.
\]
This example illustrates why the coefficient conditions are useful:
we can decide about all averaging stages without calculating them.

The criterion also lets us return to arbitrary functions without
losing control at large depths. Although we match the function only
on a finite interval, the interpolating polynomial will have integral
averages at \emph{every} depth. This useful fact is the content of
the following lemma.

\begin{lemma}[Local lifting]\label{lem:lifting}
Let $f \in \QQ^\Z$ be $\ell$-average-integral, let $n\in\Z$, and put
$m=\ell(p-1)$. There is an integer-valued polynomial $P\in M_p$ such that
\[
 P(r)=f(n+r)\qquad\text{for all }0\le r\le m.
\]
In fact, the Newton interpolant
\begin{align}
P(x)=\sum_{j=0}^{m}\Delta^j f(n)\binom{x}{j}
\label{eq:newtonP2}
\end{align}
has these properties,
and the function $r\mapsto P(r)$ is average-integral (at every depth).
\end{lemma}
\begin{proof}
Define $P$ by \eqref{eq:newtonP2}.
The matching values $P(r) = f(n+r)$ follow from \eqref{eq:newtonPagree}.
Also, each coefficient $\Delta^jf(n)$ is an integer
(since $f$ is integral), and the
binomial polynomials $\binom{x}{j}$ take integer values on all of $\Z$,
as noted above. Hence $P$ is
integer-valued on all of $\Z$, not just on the interval where it
matches $f$.

For each $j\le m$, Lemma~\ref{lem:difference} yields
$\Delta^j f \in p^{\lfloor j/(p-1)\rfloor} \Z^\Z$, so that
the integer
$\Delta^jf(n)$ is divisible by $p^{\lfloor j/(p-1)\rfloor}$;
in particular, $\Delta^jf(n) \in p^{\lfloor j/(p-1)\rfloor} R$.
This accounts for all the addends in \eqref{eq:newtonP2}.
% For $j>m$, the binomial-basis coefficient of $P$ is zero, so it
% satisfies the same condition.
Thus, the integer-valued polynomial $P \in \QQ[x]$ meets
condition (4) of Proposition~\ref{prop:criterion}
(with $a_j = \Delta^jf(n)$ for all $j \leq m$, and $a_j = 0$
for all $j > m$). Therefore, that proposition gives
$P\in M_p$ and average-integrality of the function
$r \mapsto P(r)$.
% Here average-integrality still refers to averaging over $p$ consecutive
% values, for the same fixed prime $p$. To see explicitly why the
% conclusion concerns ordinary integers, fix $k\ge0$. The polynomial
% criterion gives $(F_p^k g)(s)\in R$ for every $s\in\Z$,
% where $g(s)=P(s)$. Since $g$ is integral,
% $p^k(F_p^k g)(s)\in\Z$ as well. Thus $(F_p^k g)(s)$ belongs
% to both $R$ and $\Z[1/p]$, and is therefore an integer.
% We have proved $F_p^k g\in\Z^\Z$ for every $k\ge0$.
% No other prime needs a separate check: applying $F_p^k$ has introduced
% no denominator factors other than $p$.
\end{proof}

\section{Every depth for \texorpdfstring{$p=2$ and $p=3$}{p=2 and p=3}}
We now have the ingredients to prove Theorem~\ref{thm:small-primes}.
The key distinction between $p=2$ and $p=3$ will be
visible in the coefficients of a product and a cube.

\begin{proof}[Proof of Theorem~\ref{thm:small-primes}]
We first work with polynomials. Proposition~\ref{prop:criterion} reduces the question
to coefficient divisibility: we must show that $M_2$ is closed under
products and that $M_3$ is closed under cubes. Afterwards, local
interpolation will transfer the result to arbitrary functions.

For $p=2$, the definition of the function $d$ (see Lemma~\ref{lem:monomials})
simplifies to $d(j)=j-v_2(j!)$; therefore
\[
 d(i+j)=d(i)+d(j)-v_2\binom{i+j}{i}
 \le d(i)+d(j)
 \qquad\text{for all }i,j\ge0.
\]
The product of the two basic polynomials $2^{d(i)}x^i$ and
$2^{d(j)}x^j$ has coefficient $2^{d(i)+d(j)}$ at degree $i+j$.
The inequality we just proved says that this contains at least the required
$d(i+j)$ factors of $2$. The product therefore belongs to $M_2$.
More explicitly, if $P=\sum_i a_i2^{d(i)}x^i$ and
$Q=\sum_j b_j2^{d(j)}x^j$, with $a_i,b_j\in R$, then
\[
 PQ=\sum_{i,j}a_i b_j\bigl(2^{d(i)}x^i\bigr)
                         \bigl(2^{d(j)}x^j\bigr).
\]
Each product in parentheses belongs to $M_2$ by the preceding
calculation. Since $M_2$ is closed under finite $R$-linear
combinations, $PQ\in M_2$. Closure under addition and additive inverses follows from
the definition, and the constant polynomial $1$ belongs to $M_2$
because $d(0) = 0$.
Thus $M_2$ is a subring of $\QQ[x]$.

For $p=3$, products need not stay in $M_3$. However, when we cube a
sum, every mixed term has a coefficient divisible by $3$. We will
see that this extra factor is enough.

Here we have $d(j)=\lfloor j/2\rfloor-v_3(j!)$. Thus, for all
$i,j,k\ge0$, we have
\begin{align*}
 d(i+j+k)-d(i)-d(j)-d(k)
 &=\left\lfloor\frac{i+j+k}{2}\right\rfloor
   -\left\lfloor\frac i2\right\rfloor
   -\left\lfloor\frac j2\right\rfloor
   -\left\lfloor\frac k2\right\rfloor
   -v_3\binom{i+j+k}{i,j,k}.
\end{align*}
Here the multinomial coefficient is an integer, and the difference of
floors is at most $1$. Indeed, write $i=2a+\epsilon$,
$j=2b+\eta$, and $k=2c+\theta$, where each of
$\epsilon,\eta,\theta$ is $0$ or $1$. The difference of floors then
becomes $\lfloor(\epsilon+\eta+\theta)/2\rfloor\le1$.
Also $v_3\binom{i+j+k}{i,j,k}\ge0$, since this multinomial
coefficient is an integer. Consequently,
\begin{equation}\label{eq:triple-bound}
 d(i+j+k)\le d(i)+d(j)+d(k)+1
 \qquad\text{for all }i,j,k\ge0.
\end{equation}
Write $P=\sum_j c_jx^j\in M_3$, so that $v_3(c_j)\ge d(j)$
for every $j$. Multiplying three copies of $P$ gives
$P^3=\sum_{i,j,k}c_ic_jc_kx^{i+j+k}$.
Group the ordered triples of indices according to which indices coincide.
There is one ordering when all three are equal, three when exactly two
are equal, and six when all are distinct. Hence
\[
 P^3=\sum_j c_j^3x^{3j}
     +3\sum_{i\ne j}c_i^2c_jx^{2i+j}
     +6\sum_{i<j<k}c_ic_jc_kx^{i+j+k}.
\]
All sums are finite. First consider the pure terms
$c_j^3x^{3j}$. The $3$-valuations $v_3(c_j^3) = 3 v_3(c_j)$
of their coefficients are at least $3d(j)$
(since $v_3(c_j) \geq d(j)$),
whereas membership in $M_3$ requires only $d(3j)=d(j)$,
by Lemma~\ref{lem:monomials}(b). Since $d(j)\ge0$, we have
$3d(j)\ge d(j)$.
These terms therefore belong to $M_3$.

Each mixed term is
$3c_i^2c_jx^{2i+j}$ for some $i\ne j$, or
$6c_ic_jc_kx^{i+j+k}$ for some $i<j<k$.
In the first case the $3$-valuation of its coefficient is at least
$1+2d(i)+d(j)\ge d(2i+j)$, by~\eqref{eq:triple-bound} with indices
$i,i,j$. In the second case it is at least
$1+d(i)+d(j)+d(k)\ge d(i+j+k)$, again by~\eqref{eq:triple-bound}.
Thus every mixed term also belongs
to $M_3$, and so does their sum with the pure terms. This proves
$P^3\in M_3$.

We now transfer these polynomial statements to functions. Fix a
depth $\ell$ and an integer $n$. In the $p=2$ case, let $f,g
\in \QQ^\Z$ be $\ell$-average-integral. Lemma~\ref{lem:lifting}
(applied here with $m=\ell$, since $p-1=1$)
supplies integer-valued polynomials $P,Q\in M_2$ with
\begin{align}
 P(s)=f(n+s),\qquad Q(s)=g(n+s)
 \qquad\text{for all }0\le s\le\ell.
 \label{eq:Ps=fn+s}
\end{align}
Their product $PQ$ belongs to $M_2$ (since $M_2$ is a ring)
and is integer-valued.
Put $h(s)=P(s)Q(s)$ for all $s\in\Z$.
Proposition~\ref{prop:criterion} (specifically, its last
sentence) gives $F_2^r h\in\Z^\Z$ for every $r\ge0$, while
\eqref{eq:Ps=fn+s} shows that
\begin{align}
 h(s)=(fg)(n+s) \qquad \text{for all }0\le s\le\ell.
 \label{eq:hs=fgn+s}
\end{align}

Consider an arbitrary function $w : \Z \to \QQ$.
Any value of $F_p^r w$ is determined by only
$r(p-1)+1$ consecutive values of $w$. This can
be seen by writing each application of $F_p$ as a sum of shifts.
For $r\ge0$, the equality \eqref{eq:Fp=} gives
\[
 F_p^r=\frac1{p^r}(\id+E+\cdots+E^{p-1})^r
 =\frac1{p^r}\sum_{a_1,\ldots,a_r\in\left\{0,1,\ldots,p-1\right\}}E^{a_1+\cdots+a_r}
\]
by distributivity
(where for $r = 0$, the sum consists of a single addend equal
to $E^0$, corresponding to the empty $0$-tuple).
Here $a_i$ records which shift was chosen from the $i$th factor.
Applying this identity to a function $w$ and evaluating at $n$ gives
\begin{equation}\label{eq:iterated-average-values}
 (F_p^r w)(n)=\frac1{p^r}
 \sum_{a_1,\ldots,a_r\in\left\{0,1,\ldots,p-1\right\}} w(n+a_1+\cdots+a_r)
\end{equation}
for all $r \geq 0$.
In the $p=2$ case, all the inner sums $a_1+\cdots+a_r$
lie between $0$ and $r$.
Thus, for $0\le r\le\ell$, if we compute $(F_2^r(fg))(n)$
and $(F_2^r h)(0)$ using \eqref{eq:iterated-average-values},
then each summand in the formula for
$(F_2^r(fg))(n)$ equals the corresponding summand for $(F_2^r h)(0)$,
by~\eqref{eq:hs=fgn+s}.
Consequently,
\[
 (F_2^r(fg))(n)=(F_2^r h)(0)\in\Z
 \qquad \left(\text{since } F_2^r h \in \Z^\Z \right).
\]
Since $n$ was arbitrary, this shows that $F_2^r(fg) \in \Z^\Z$.
Since $r \leq \ell $ was arbitrary, this means that
$fg$ is $\ell$-average-integral.
So the set of all $\ell$-average-integral functions is closed
under multiplication.
Linearity of $F_2$ also gives closure under sums and additive inverses,
and the constant function $1$ is $\ell$-average-integral
because $F_2$ fixes it. Hence, the set of all
$\ell$-average-integral functions for $p=2$ is a ring
under pointwise addition and multiplication.

For $p=3$, let us again fix an $\ell$-average-integral
function $f \in \QQ^\Z$, and
apply the same argument to an integer-valued
polynomial $P\in M_3$ that satisfies
\begin{align}
 P(s)=f(n+s)
 \qquad\text{for all }0\le s\le 2\ell.
 \label{eq:Ps=fn+s2l}
\end{align}
(Such a $P$ again exists by Lemma~\ref{lem:lifting}.)
As we saw above, $P \in M_3$ entails $P^3 \in M_3$.
Put $h(s)=P(s)^3$.
Since $P^3$ is integer-valued and belongs to $M_3$,
Proposition~\ref{prop:criterion} (again, its last
sentence, applied to $P^3$) shows that the function
$h$ is average-integral, i.e., we have
$F_3^r h\in\Z^\Z$ for all $r\ge0$.
For $0\le r\le\ell$, the sum $a_1+\cdots+a_r$ in
\eqref{eq:iterated-average-values} lies between $0$ and $2r\le2\ell$.
At all these arguments, we have $h(s)=P(s)^3=f(n+s)^3$
by \eqref{eq:Ps=fn+s2l}, and hence
\[
 (F_3^r(f^3))(n)=(F_3^r h)(0)\in\Z.
\]
Since $n$ was arbitrary, this shows that
$F_3^r(f^3)\in\Z^\Z$, and therefore
$f^3$ is $\ell$-average-integral.
\end{proof}

\begin{proof}[Proof of Corollary~\ref{cor:inf-depth23}]
Suppose $f,g$ are average-integral for $p=2$.
For any $k\ge0$, apply the $p=2$ case of Theorem~\ref{thm:small-primes}
with $\ell=k$ to see that $fg$ is $k$-average-integral;
thus, $F_2^k(fg)\in\Z^\Z$. This shows that $fg$ is average-integral.

Likewise, if $f$ is average-integral for $p=3$, then
Theorem~\ref{thm:small-primes} for $p=3$ and
$\ell=k$ gives $F_3^k(f^3)\in\Z^\Z$ for every $k\ge0$.
Thus $f^3$ is average-integral.
\end{proof}

\begin{remark}
The polynomial strategy for $p=2,3$ was already sketched by Gjergji Zaimi
in his June 16, 2011 answer to~\cite{MO-original}. The proof above supplies
the full coefficient estimates and the local lifting argument.
The distinction between products and cubes is necessary: for $p=3$,
$f(n)=n$ is average-integral, since $F_3^rf(n)=n+r$, whereas
\[
 (F_3(f^2))(n)=n^2+2n+\frac53
\]
is never integral.
\end{remark}

\section{Sharpness for primes greater than three}\label{sec:sharp}
We now verify zeb's counterexample from~\cite{zeb} and use
Theorem~\ref{thm:main} to identify its first failing depth.

\begin{proposition}[zeb's example, with exact failing depth]\label{prop:sharp}
Let $p>3$ be prime and put $f(n)=1+n^{p+1}$. Then $f$ is
average-integral, while $f^p$ is $p$-average-integral but is not
$(p+1)$-average-integral. In particular, $F_p^{p+1}(f^p)$ is the first
nonintegral iterate.
\end{proposition}
For a first look at the obstruction, take $p=5$. The relevant
expansion is
\[
 (1+x^6)^5=1+5x^6+10x^{12}+10x^{18}+5x^{24}+x^{30}.
\]
The polynomial criterion requires two factors of $5$ at degree $24$,
since $d(24)=6-4=2$, but the coefficient is only $5$.
All other terms satisfy the criterion. This already shows that
some iterate fails. The proof below also determines the first
failing iterate, by using a falling factorial of degree $24$.

\begin{proof}
Since $p>3$, Legendre's formula gives $v_p((p+1)!)=1$ and
$\lfloor(p+1)/(p-1)\rfloor=1$. Thus $d(p+1)=0$, so
Lemma~\ref{lem:monomials}(d)
shows that $n\mapsto n^{p+1}$ is average-integral.
Since $F_p$ fixes the constant function $1$ and is linear,
$F_p^k f=1+F_p^k(n\mapsto n^{p+1})$ is integral for every $k\ge0$.
Thus $f$ is average-integral. We now examine $f^p$.
The idea is to separate its polynomial into
terms that are already known to be average-integral and one term
that fails the difference test. The binomial theorem first gives
\[
 (1+x^{p+1})^p=\sum_{r=0}^p\binom pr(x^{p+1})^r
             =\sum_{r=0}^p\binom pr x^{r(p+1)}.
\]
The terms with $r=0$ and $r=p$ are $1$ and $x^{p^2+p}$.
The term with $r=p-1$ is $px^{(p-1)(p+1)}$.
Set $N=p^2-1=(p-1)(p+1)$ and separate these three terms to obtain
\[
 (1+x^{p+1})^p
 =1+\sum_{r=1}^{p-2}\binom pr x^{r(p+1)}
   +p x^N+x^{p^2+p}.
\]
For $1\le r\le p-2$, we have $r(p+1)<p^2$, so
\[
 v_p\bigl((r(p+1))!\bigr)=r,\qquad
 d(r(p+1))=\left\lfloor\frac{2r}{p-1}\right\rfloor\le1.
\]
Also $v_p((p^2+p)!)=p+2=\lfloor(p^2+p)/(p-1)\rfloor$, and hence
$d(p^2+p)=0$. For $1\le r\le p-2$, the coefficient $\binom pr$ is divisible
by $p$, which supplies the at most one factor demanded by
$d(r(p+1))\le1$. The last term needs no additional factor, since
$d(p^2+p)=0$. By Lemma~\ref{lem:monomials}(d), all the displayed
summands except possibly $px^N$ are therefore average-integral.

To test the remaining term $px^N$, we replace $x^N$ by the falling
factorial $\fall{x}{N}$, whose $N$th difference is especially simple.
We must check that the lower-degree terms introduced by this
replacement are harmless.

For $0\le j<N$, write $j=ap+b$ with $0\le b<p$.
Then $0\le a\le p-1$ and $a+b\le2p-3$: equality $a=b=p-1$ is
excluded by $j<N$. Legendre's formula gives $v_p(j!)=a$, whence
\[
 d(j)=\left\lfloor\frac{ap+b}{p-1}\right\rfloor-a
     =\left\lfloor\frac{a+b}{p-1}\right\rfloor\le1.
\]
Since $x^N-\fall{x}{N}$ is an integer polynomial
of degree less than $N$, every coefficient of $p(x^N-\fall{x}{N})$ is divisible by $p$.
The inequalities $d(j)\le1$ show that each of its monomials satisfies
the polynomial criterion. Thus $p(x^N-\fall{x}{N})$ is
average-integral. We conclude that
\begin{align}
 (1+x^{p+1})^p=p\fall{x}{N}+A(x),
 \label{eq:1+xp+1p=pxN+A}
\end{align}
where $A$ defines an average-integral function $a(n)=A(n)$.
Put $h(n)=p\fall{n}{N}$. Evaluating the identity
\eqref{eq:1+xp+1p=pxN+A} at integers
gives $f^p=h+a$, so linearity gives
\[
 F_p^k h=F_p^k(f^p)-F_p^k a
 \qquad\text{for all }k\ge0.
\]
Since $F_p^k a$ is integral for every $k$, if $f^p$ were
$(p+1)$-average-integral, then $h$ would be $(p+1)$-average-integral
as well. Lemma~\ref{lem:difference}, with
$N=(p+1)(p-1)$, would force
$\Delta^N h$ to belong to
$p^{\left\lfloor N/(p-1) \right\rfloor} \Z^\Z
= p^{p+1} \Z^\Z$, so that all values of the polynomial
$\Delta^N\bigl(p\fall{x}{N}\bigr)$ at integers
would be multiples of $p^{p+1}$.
In other words, all values of the polynomial
$\Delta^N \fall{x}{N}$ would be multiples of $p^p$.
But this is false, since iterating
\eqref{eq:falling-difference} $N$ times shows that
$\Delta^N \fall{x}{N}$ is the constant polynomial $N!$,
which is not a multiple of $p^p$
(indeed, Legendre's formula shows that $v_p(N!)=p-1$).
This contradiction shows that $f^p$ cannot be
$(p+1)$-average-integral.
Finally, Theorem~\ref{thm:main} shows that $f^p$ is $p$-average-integral,
so failure first occurs at iteration $p+1$.
\end{proof}

\begin{remark}
The proof of Lemma~\ref{lem:polynomial} also helps explain why the
depth $p+1$ is different from the depths $1,\ldots,p$.
There we wrote $P(x)^p=\sum_j c_j\fall{x}{j}$ and used
$\Delta^m(P^p)(0)=m!c_m$.
For $m=k(p-1)$ with $1\le k\le p$, the factorial $m!$ supplies
exactly $k-1$ factors of $p$. The coefficient argument modulo $p$
supplies the remaining factor by showing that $c_m\in pR$.
At $k=p+1$, however, we have $m=p^2-1$ and still
$v_p(m!)=p-1$. To obtain $\Delta^m(P^p)(0)\in p^{p+1}R$,
we would now need $c_m\in p^2R$, rather than just $c_m\in pR$.
Thus the argument modulo $p$ no longer suffices.
Proposition~\ref{prop:sharp} shows that, for every $p>3$, this
stronger divisibility can actually fail.
\end{remark}

\section{Non-primes}\label{sec:composite}

The assumption that the averaging length $p$ is prime is essential
for Theorem~\ref{thm:main}. Here are some examples of what can
go wrong when we average over composite-length windows.

For an integer $m\ge2$, define an operator $F_m : \QQ^\Z \to \QQ^\Z$
by
\[
 F_m f(n)=\frac{f(n)+f(n+1)+\cdots+f(n+m-1)}{m}
\]
(generalizing the $F_p$ we have studied above).
An $m$-periodic function $f$ exhibits very simple behavior
under averaging:
every block of $m$ consecutive values has the same sum, so $F_m f$
is constant. If that sum is zero, then $F_m f=0$, and hence $F_m^k f=0$
for every $k\ge1$.

\begin{example}
For $m=4$, let $f$ be the $4$-periodic function whose values in one period
are $1,-1,0,0$. Then $F_4 f=0$, whereas
\[
 F_4(f^4)=\frac{1+1+0+0}{4}=\frac12.
\]
Thus, even when $F_4^k f$ is integral for every $k\ge0$, its fourth
power can fail at the very first averaging step. The same example,
with $m-2$ zeros in place of the two zeros, works for every even
$m>2$: it gives $F_m f=0$ and $F_m(f^m)=2/m$.
\end{example}

There is a similar obstruction whenever $m$ is divisible by the
square of a prime:

\begin{example}
Suppose that $r$ is prime and $r^2\mid m$.
Take an $m$-periodic function $f$ with $r$ entries equal to $1$, one entry
equal to $-r$, and all remaining entries equal to zero in each period.
There is room for these entries because $r+1\le r^2\le m$.
Their sum is zero, so again $F_m^k f=0$ for every $k\ge1$.
On the other hand,
\[
 F_m(f^m)=\frac{r+(-r)^m}{m}.
\]
The numerator is congruent to $r$ modulo $r^2$, and therefore is not
divisible by $m$. Hence $f^m$ fails to be integer-valued after one
averaging step. In particular, this rules out an analogue of the
power-preservation theorem for every composite prime power.
\end{example}

Odd squarefree lengths can fail as well:

\begin{example}
For $m=15$, take the
$15$-periodic function with values $1,1,-2$ followed by twelve zeros
in each period. Again the sum is zero, so $F_{15}^k f=0$ for every
$k\ge1$. But
\[
 F_{15}(f^{15})=\frac{1+1+(-2)^{15}}{15}
 =-\frac{10922}{5}\notin\Z.
\]
Thus the power-preservation assertion fails at the first step even
for $m=3\cdot5$.

More generally, suppose $m\ge3$ is odd and some integer $a$ satisfies
$a^m\not\equiv a\pmod m$. We may choose $0\le a\le m-1$,
since the congruence depends only on the residue of $a$ modulo $m$.
Take an $m$-periodic function with $a$ entries equal to $1$, one
entry equal to $-a$, and $m-a-1$ entries equal to zero in each period.
Its sum is zero, whereas the sum of its $m$th powers is
$a-a^m$, which is not divisible by $m$.
This gives the same kind of counterexample.
\end{example}

Conversely, if
\[
 a^m\equiv a\pmod m\qquad\text{for all }a\in\Z
\]
(that is, $m$ is either prime or a Carmichael number),
then every integer-valued function $f$ satisfies
\[
 \sum_{r=0}^{m-1}f(n+r)^m
 \equiv\sum_{r=0}^{m-1}f(n+r)\pmod m
 \qquad\text{for all }n\in\Z.
\]
Consequently, integrality of $F_m f$ implies integrality of
$F_m(f^m)$. Thus, for odd $m$, the displayed congruence for all $a$
is exactly the condition for power preservation at the first step.
This does not imply preservation at all averaging steps.
We can, however, obtain integrality at some larger depths,
provided we increase the exponent along with the depth.

\begin{proposition}\label{prop:composite-growing-powers}
Let $m\ge2$ be an integer such that
\[
 a^m\equiv a\pmod m\qquad\text{for all }a\in\Z.
\]
Let $f\colon\Z\to\Z$ satisfy $F_m f\in\Z^\Z$.
Then, for every integer $t\ge0$, we have
\[
 F_m^k(f^{m^t})\in\Z^\Z
 \qquad\text{for all }0\le k\le t+1.
\]
\end{proposition}
\begin{proof}
Subtracting two consecutive averages gives
\[
 f(n+m)-f(n)
 =m\bigl((F_m f)(n+1)-(F_m f)(n)\bigr).
\]
Thus
\begin{align}
f(n+m)\equiv f(n)\pmod m \qquad \text{ for every } n\in\Z
\label{eq:fn+m}
\end{align}
(since $F_m f\in\Z^\Z$).

We first check how raising to the $m$th power improves this
congruence. Suppose $a$ and $b$ are two integers such that
$a\equiv b\pmod{m^r}$, where $r\ge1$,
and write $a=b+m^r u$ with $u\in\Z$. Expanding
$a^m = (b + m^r u)^m$ using the binomial theorem gives
\[
 a^m-b^m=m^{r+1}u b^{m-1}
 +\sum_{i=2}^m\binom mi m^{ri}u^i b^{m-i}.
\]
Every term is divisible by $m^{r+1}$, since $ri\ge r+1$
when $i\ge2$. Thus, $a^m \equiv b^m \pmod{m^{r+1}}$.
Applying this observation $t$ times yields
\[
 a^{m^t}\equiv b^{m^t}\pmod{m^{r+t}}.
\]
Thus, from \eqref{eq:fn+m}, we obtain
\[
 f(n+m)^{m^t}\equiv f(n)^{m^t}\pmod{m^{t+1}}
 \qquad\text{for all }n\in\Z.
\]

Let $b : \Z \to \Z$ be the $m$-periodic function agreeing
with $f^{m^t}$ at
$0,1,\ldots,m-1$. The last congruence shows that
\[
 f^{m^t}=b+m^{t+1}c
\]
for an integer-valued function $c$ on $\Z$.
Repeatedly using the hypothesis $a^m\equiv a\pmod m$ gives
$a^{m^t}\equiv a\pmod m$ for all $a \in \Z$. Consequently,
\[
 \sum_{s=0}^{m-1}b(s)
 =\sum_{s=0}^{m-1}f(s)^{m^t}
 \equiv\sum_{s=0}^{m-1}f(s)
 =m(F_m f)(0)\equiv0\pmod m.
\]
Since $b$ is $m$-periodic, $F_m b$ is therefore an integer constant,
and every further average of $b$ is the same constant.
On the other hand, $m^k F_m^k c$ is integer-valued: expanding
the $k$ averaging operations expresses it as a sum of integer
values of $c$. Hence $m^{t+1}F_m^k c$ is integer-valued for
$0\le k\le t+1$. Applying $F_m^k$ to
$f^{m^t}=b+m^{t+1}c$ proves the assertion.
\end{proof}

\begin{example}[The averaging length $561$]\label{exa:561}
The composite integer $m=561=3\cdot11\cdot17$ is known to be
a Carmichael number; thus it satisfies the
hypothesis of Proposition~\ref{prop:composite-growing-powers}.
% Indeed, $2$, $10$, and $16$ all divide $560$, so Fermat's little
% theorem gives $a^{561}\equiv a$ modulo each of $3$, $11$, and $17$
% for every integer $a$. Since these primes are distinct, the same
% congruence holds modulo $561$.
Hence, integrality of $F_{561}f$ guarantees
(by Proposition~\ref{prop:composite-growing-powers}) that
the first two averages
of $f^{561}$, the first three averages of $f^{561^2}$, and, more
generally, the first $t+1$ averages of $f^{561^t}$ are integral.

Keeping the exponent equal to $561$, however, does not give
preservation at every depth. Take $f(n)=n$. Then
\[
 (F_{561}^k f)(n)=n+280k
 \qquad\text{for all }k\ge0\text{ and }n\in\Z,
\]
so every average of $f$ is integer-valued. Nevertheless,
a mechanical computation (see Appendix~\ref{app:python}) shows that
\begin{equation}\label{eq:561-obstruction}
 \sum_{a_1,\ldots,a_5=0}^{560}(a_1+\cdots+a_5)^{561}
 \equiv162\pmod{243}.
\end{equation}
The left-hand side is $561^5(F_{561}^5(f^{561}))(0)$
(by \eqref{eq:iterated-average-values}).
Since $243=3^5$ divides $561^5$, the congruence shows that
$(F_{561}^5(f^{561}))(0)$ is not an integer.
\end{example}

\section{Connections with the literature}\label{sec:connect}
The averaging problem also belongs to a longer history of divisibility
conditions on finite differences. We describe some connections here,
including two results from which parts of our conclusions can be
deduced. This section is independent of the proofs above. The
paragraphs on $p$-adic functions provide additional context and may
be skipped without affecting the rest of the note.

\subsection{An exact antecedent for the binary product theorem}
Sun~\cite[Definition 2.1]{Sun} calls a function
$u\colon\Z_{\ge0}\to\Z_{(p)}$ \emph{$p$-regular} if
\[
 \sum_{r=0}^{j}(-1)^r\binom jr u(r)\in p^j\Z_{(p)}
 \qquad\text{for all }j\ge0.
\]
These sums equal $(-1)^j\Delta^ju(0)$. His Product Theorem,
\cite[Theorem 2.3, p.~200]{Sun}, states that products of $p$-regular
functions are $p$-regular. For $p=2$, it gives our infinite-depth product assertion. To see the
connection, let $f,g$ be average-integral functions and fix $n\in\Z$.
The functions $u(r)=f(n+r)$ and $v(r)=g(n+r)$ on the nonnegative
integers are $2$-regular by Lemma~\ref{lem:difference}, since $q=1$.
Sun's theorem gives $2^j\mid\Delta^j(fg)(n)$ in $\Z_{(2)}$ for
every $j$. These differences are ordinary integers, so the same
divisibility holds in $\Z$. Since $n$ was arbitrary, the difference
criterion proves that $fg$ is average-integral.
For $p=3$ our bound is $\lfloor j/2\rfloor$, not $j$, so Sun's theorem
does not cover our class.

There is a related formulation using $p$-adic integers. Recall that
$\Z_p$ records compatible residues modulo $p,p^2,p^3,\ldots$; in its
notion of convergence, $a_j\to0$ means that for each $k$, all
sufficiently late terms are divisible by $p^k$. This ring is different
from $R=\Z_{(p)}$, whose elements are rational numbers.

Kellner~\cite[Definition 3.1 and Theorem 6.1]{Kellner} studies the
$p$-adic Kummer space $\mathcal K_{p,2}$ defined by
$\Delta^j u(x)\in p^j\Z_p$ and proves that it is a ring. At $p=2$,
this is exactly the continuous $p$-adic version of the present class.

\subsection{Mahler expansions and the exact coefficient condition}
Newton interpolation has an infinite-series analogue over $\Z_p$.
Mahler's interpolation theorem~\cite{Mahler} identifies continuous
functions $u\colon\Z_p\to\Z_p$ with uniformly convergent series
\[
 u(x)=\sum_{j\ge0}a_j\binom{x}{j},\qquad
 a_j=\Delta^ju(0)\longrightarrow0.
\]
This extension procedure also allows us to recover Sun's earlier
product theorem from Kellner's ring theorem, for every prime $p$.
Indeed, let $u\colon\Z_{\ge0}\to R$ be $p$-regular and put
$a_j=\Delta^ju(0)$. Embed $R$ in $\Z_p$.
Since $a_j\in p^jR$, the series
\[
 U(x)=\sum_{j\ge0}a_j\binom{x}{j}
\]
converges uniformly on $\Z_p$. Newton interpolation gives
$U(n)=u(n)$ for every nonnegative integer $n$.
The identity $\Delta\binom{x}{j}=\binom{x}{j-1}$ for $j\ge1$
then gives
\[
 \Delta^k U(x)=\sum_{j\ge k}a_j\binom{x}{j-k}
 \in p^k\Z_p
 \qquad\text{for all }k\ge0\text{ and }x\in\Z_p.
\]
Here each binomial polynomial takes values in $\Z_p$, and the
series still converges after division by $p^k$.
Thus $U\in\mathcal K_{p,2}$.
If $v$ is another $p$-regular function, its extension $V$ likewise
belongs to $\mathcal K_{p,2}$. Kellner's ring theorem gives
$UV\in\mathcal K_{p,2}$, so $\Delta^k(uv)(0)\in p^k\Z_p$.
This finite difference also belongs to $R$. Since
$R\cap p^k\Z_p=p^kR$, it follows that $uv$ is $p$-regular,
as asserted by Sun's theorem.

In this language the continuous $p$-adic version of average-integrality is
precisely
\begin{equation}\label{eq:mahler-bound}
 a_j\in p^{\lfloor j/(p-1)\rfloor}\Z_p\qquad\text{for all }j\ge0.
\end{equation}
For necessity, apply the difference criterion and evaluate the
$j$th difference at zero. For sufficiency, note that
$\Delta\binom{x}{j}=\binom{x}{j-1}$ for $j\ge1$, by
\eqref{eq:falling-difference} after division by $j!$.
Taking $r(p-1)$ differences
therefore shifts the binomial indices by $r(p-1)$, just as in
\eqref{eq:divided-newton}. Every remaining coefficient is divisible
by $p^r$, and division by $p^r$ still leaves a convergent series with
values in $\Z_p$. The difference criterion now gives the desired
integrality of the averages.

Every average-integral function has such a continuous extension,
including agreement at negative integers. Here is a direct justification
of the latter point. The bounds on $a_j=\Delta^jf(0)$ make the Mahler
series converge uniformly and agree with $f$ at nonnegative integers.
For a fixed negative integer $n$, put
\[
 A_j=\Delta^jf(n),\qquad
 P_N(x)=\sum_{j=0}^N A_j\binom{x-n}{j}.
\]
The same bounds give a uniform limit $u_n$ on $\Z_p$, and Newton
interpolation gives $u_n(x)=f(x)$ for every integer $x\ge n$.
Thus $u_n$ and the extension constructed at zero agree on the dense set
$\Z_{\ge0}$, hence everywhere; in particular they agree at $n$.
This proves the asserted extension. Conversely, the restriction of a
function satisfying~\eqref{eq:mahler-bound} is average-integral whenever
its values on $\Z$ are ordinary integers.

\subsection{A second proof using Stirling numbers}
The Stirling number $S(N,m)$ of the second kind counts the partitions
of an $N$-element set into $m$ nonempty blocks. These integers also
express ordinary powers in the falling-factorial basis:
\[
 x^N=\sum_{m=0}^N S(N,m)\fall{x}{m}.
\]
For a nonnegative integer $x$, this identity counts maps from an
$N$-element set to an $x$-element set: first partition the domain
into its nonempty fibers, then choose distinct images for the blocks.
It holds as a polynomial identity because it holds at every
nonnegative integer $x$. Thus a divisibility theorem about $S(N,m)$
can be used in place of our elementary coefficient argument.

In \cite[Lemma 2.1]{Adelberg}, Adelberg proves a lower bound
on $v_p(S(N,m))$. In terms of factorial valuations, his bound is
\begin{equation}\label{eq:adelberg}
 v_p(S(N,m))\ge v_p(N!)-v_p(m!)
 -\left\lfloor\frac{N-m}{p-1}\right\rfloor
\end{equation}
for all $N\ge m\ge1$.
This bound yields another proof of Theorem~\ref{thm:main}, as follows.

Let $f$ be $\ell$-average-integral, where $1\le\ell\le p$,
and fix $n\in\Z$. Lemma~\ref{lem:lifting} supplies an
integer-valued polynomial $P\in M_p$ with
$P(r)=f(n+r)$ for $0\le r\le\ell(p-1)$.
Write $P(x)=\sum_j c_jx^j$. By definition of $M_p$,
$c_j\in p^{d(j)}R$. Since $d(j)\ge0$, reducing modulo $p$ therefore gives
\[
 \overline P(x)=\sum_{j:\,d(j)=0}\overline{c_j}x^j.
\]
In characteristic $p$, the intermediate binomial coefficients
$\binom pa$, for $1\le a<p$, vanish. Hence $(U+V)^p=U^p+V^p$;
applying this identity repeatedly to the finite sum above gives
\[
 \overline P(x)^p
 =\sum_{j:\,d(j)=0}\overline{c_j}^{\,p}x^{pj}
 =\sum_{j:\,d(j)=0}\overline{c_j}x^{pj}.
\]
The last equality uses $c^p=c$ in $\Fp$.
Every exponent $N=pj$ in this sum satisfies $d(N)=d(j)=0$,
by Lemma~\ref{lem:monomials}(b).

As in the elementary proof, we want the falling-factorial coefficient
at $m=k(p-1)$ to be divisible by $p$. Stirling numbers describe exactly
this change of basis. We will check that every monomial surviving
modulo $p$ has the required coefficient divisible by $p$.
For $m=k(p-1)$ with $1\le k\le p$, we have $v_p(m!)=k-1$.
Therefore the right side of~\eqref{eq:adelberg} simplifies to
\[
 v_p(N!)-(k-1)
 -\left(\left\lfloor\frac{N}{p-1}\right\rfloor-k\right)
 =1-d(N).
\]
Taking $N=pj$ gives $p\mid S(pj,m)$ whenever $pj\ge m$ and $d(j)=0$;
for $pj<m$ the Stirling number is zero. It follows that
\[
 [\fall{x}{m}]\bigl(\overline P^{\,p}\bigr)
 =\sum_{j:\,d(j)=0}\overline{c_j}\,\overline{S(pj,m)}=0.
\]
Thus, if $b_m$ is the $\fall{x}{m}$-coefficient of $P^p$, then
$b_m\in pR$. By~\eqref{eq:newton-coefficient},
\[
 \Delta^m(P^p)(0)=m!b_m\in p^kR,
\]
since $v_p(m!)=k-1$.
For $1\le k\le\ell$, we have $m=k(p-1)\le\ell(p-1)$.
The polynomials $P^p$ and the shifted function $r\mapsto f(n+r)^p$
therefore agree at $r=0,\ldots,m$. Applying
\eqref{eq:interpolation-differences} gives
\[
 \Delta^{k(p-1)}(f^p)(n)=\Delta^m(P^p)(0)\in p^kR.
\]
The left side is an integer, so it is divisible by $p^k$.
Since $n$ and $k$ were arbitrary, Lemma~\ref{lem:difference} gives
$F_p^k(f^p)\in\Z^\Z$ for $0\le k\le\ell$.
The case $\ell=0$ follows directly from integrality of $f$.

Thus the finite-depth conclusion can be deduced from this published
valuation bound and the polynomial criterion already present in the
2011 MathOverflow discussion. The proof in Section~\ref{sec:congruence} gives a direct
finite-field argument for the particular vanishing needed here.

% \subsection{Scope of the search}
% The original binary square theorem is Tamarkin's 1993 olympiad problem;
% an archived statement and induction proof are available in~\cite{Tamarkin}.
% The 2011 discussion~\cite{MO-original} explicitly records the truth of
% the ternary case, Zaimi's polynomial approach, and zeb's polynomial
% criterion and counterexample. These are antecedents, not new claims of
% this note. The search for this note (September 27, 2026), using the
% averaging formulation as well as $p$-regular functions, Kummer spaces,
% Mahler coefficients, and Stirling-number valuations, did not locate a
% journal statement of the exact ternary averaging theorem or of the
% finite-depth assertion in the form of Theorem~\ref{thm:main}. This is a
% report of what was found, not a claim that no such publication exists.

\appendix
\section{Power sums over a prime field}\label{app:power-sums}
We prove~\eqref{eq:residue-sums}, working throughout in $\Fp$.

When $b=0$, our convention $r^0=1$ makes the sum a sum of $p$ ones,
which is zero in $\Fp$. When $b=p-1$, Fermat's little theorem gives
$r^{p-1}=1$ for every nonzero $r$, while $0^{p-1}=0$. Hence the sum
is $p-1=-1$ in $\Fp$.

It remains to consider $1\le b<p-1$. Choose a nonzero residue $a$
with $a^b\ne1$. Such a residue exists because the polynomial
$X^b-1$ has at most
$b$ roots, fewer than the $p-1$ nonzero residues. Write
$S=\sum_{r\in\Fp}r^b$. Multiplication by $a$ permutes $\Fp$, so
\[
 S=\sum_{r\in\Fp}(ar)^b=a^b\sum_{r\in\Fp}r^b=a^bS.
\]
Since $a^b-1\ne0$ in the field $\Fp$, the equality forces $S=0$.

\section{Verification of Example~\ref{exa:561}}\label{app:python}
The following exact computation in Python verifies
\eqref{eq:561-obstruction}. After $r$ passes through the loop,
\texttt{counts[s]} is, modulo $243$, the number of $r$-tuples
with entries in $\{0,1,\ldots,560\}$ whose sum is congruent to
$s$ modulo $243$. Appending an entry $a$ gives the recurrence
used in the loop. The final sum weights each residue $s$ by
$s^{561}$ modulo $243$.
\begin{verbatim}
modulus = 243
counts = [1] + [0] * (modulus - 1)
for r in range(5):
	counts = [
		sum(counts[(s - a) % modulus] for a in range(561)) % modulus
		for s in range(modulus)
	]
remainder = sum(counts[s] * pow(s, 561, modulus)
			for s in range(modulus)) % modulus
assert remainder == 162
\end{verbatim}
This establishes failure of average-integrality
by the fifth step; it does not assert
that the fifth step is the first failure.

\begin{thebibliography}{9}
\bibitem{MO-finite}
Darij Grinberg,
\emph{Generalization of Tamarkin's ARO 1993, final round, problem 10/8:
part II}, MathOverflow question 68595, June 23, 2011.
\url{https://mathoverflow.net/q/68595}.

\bibitem{MO-original}
Darij Grinberg,
\emph{Generalization of Tamarkin's ARO 1993, final round, problem 10/8:
still a conjecture?}, MathOverflow question 66959, June 5, 2011.
\url{https://mathoverflow.net/q/66959}.

\bibitem{zeb}
zeb,
Answer to MathOverflow question 66959, June 16, 2011
(counterexample and the weighted-monomial criterion).
\url{https://mathoverflow.net/a/68000}.

\bibitem{Sun}
Zhi-Hong Sun,
\emph{Congruences concerning Bernoulli numbers and Bernoulli polynomials},
Discrete Applied Mathematics \textbf{105} (2000), 193--223.
\href{https://doi.org/10.1016/S0166-218X(00)00184-0}
{doi:10.1016/S0166-218X(00)00184-0}.

\bibitem{Kellner}
Bernd C. Kellner,
\emph{Classification of $p$-adic functions satisfying Kummer type congruences},
preprint, 2009, arXiv:0909.0743v2.
\url{https://arxiv.org/abs/0909.0743}.

\bibitem{Mahler}
Kurt Mahler,
\emph{An interpolation series for continuous functions of a $p$-adic variable},
Journal f\"ur die reine und angewandte Mathematik \textbf{199} (1958), 23--34.
\href{https://doi.org/10.1515/crll.1958.199.23}
{doi:10.1515/crll.1958.199.23}.

\bibitem{Adelberg}
Arnold Adelberg,
\emph{The $p$-adic analysis of Stirling numbers via higher order Bernoulli numbers},
International Journal of Number Theory \textbf{14} (2018), no.~10, 2767--2779.
Preprint: \url{https://arxiv.org/abs/1805.00995}; Lemma~2.1.

\bibitem{Gessel}
Ira M. Gessel,
\emph{Applications of the classical umbral calculus},
Algebra Universalis \textbf{49} (2003), no.~4, 397--434.
\href{https://doi.org/10.1007/s00012-003-1813-5}
{doi:10.1007/s00012-003-1813-5}.
Preprint: \url{https://arxiv.org/abs/math/0108121}.

\bibitem{Tamarkin}
D. Tamarkin,
Problem 8 for grade 10, final round of the 1993 All-Russian Mathematical
Olympiad. Archived statement and solution, Problems.ru, problem 109520.
\url{https://problems.ru/view_problem_details_new.php?id=109520}.

\end{thebibliography}
\end{document}