Indefinite sum

Falling factorials

Falling factorials provide the discrete analog of the power rule from differential calculus. In infinitesimal calculus, ${\displaystyle {\frac {d}{dx}}x^{n}=nx^{n-1}}$. In the calculus of finite differences, the falling factorial

${\displaystyle (x)_{n}=x^{\underline {n}}=x(x-1)(x-2)\cdots (x-n+1)={\frac {\Gamma \left(x+1\right)}{\Gamma \left(x-n+1\right)}}}$

plays the role of ${\displaystyle x^{n}}$, and the forward difference operator satisfies

${\displaystyle \Delta (x)_{n}=n(x)_{n-1}.}$

The indefinite sum of a falling factorial is given by the discrete analog of the power rule for integration:

${\displaystyle \sum _{x}(x)_{n}={\frac {(x)_{n+1}}{n+1}}+C(x),\quad n\neq -1.}$

Equivalently, using the Gamma function:

${\displaystyle \sum _{x}{\frac {\Gamma \left(x+1\right)}{\Gamma \left(x-n+1\right)}}={\frac {\Gamma \left(x+1\right)}{\left(n+1\right)\Gamma \left(x-n\right)}}+C(x),\quad n\neq -1.}$

For the case where ${\displaystyle n=-1}$, the solution is the digamma function with a shift, ${\displaystyle \psi \left(x+1\right)+C(x)}$, which naturally extends the harmonic numbers.

Summation by parts

Indefinite summation by parts is the discrete analog of integration by parts. It is derived from the product rule for the forward difference operator.

Product rule. For two functions ${\displaystyle u(x)}$ and ${\displaystyle v(x)}$, the product rule for the forward difference is:

${\displaystyle \Delta (u(x)v(x))=u(x)\Delta v(x)+v(x+1)\Delta u(x).}$

Introducing the shift operator ${\displaystyle \mathrm {E} }$, defined by ${\displaystyle \mathrm {E} f(x)=f(x+1)}$, this can be written more compactly as:

${\displaystyle \Delta (uv)=u\Delta v+\mathrm {E} v\,\Delta u.}$

Summation by parts. Rearranging the product rule gives:

${\displaystyle u(x)\Delta v(x)=\Delta (u(x)v(x))-v(x+1)\Delta u(x).}$

Taking the indefinite sum of both sides and using the fact that ${\displaystyle \sum _{x}\Delta F(x)=F(x)+C(x)}$ (where ${\displaystyle C(x)}$ is an arbitrary 1‑periodic function) yields the formula for summation by parts:^[11][10]

${\displaystyle \sum _{x}u(x)\Delta v(x)=u(x)v(x)-\sum _{x}v(x+1)\Delta u(x)+C(x).}$

A symmetrical form, also obtained from the product rule, is:

${\displaystyle \sum _{x}f(x)\Delta g(x)+\sum _{x}g(x)\Delta f(x)=f(x)g(x)-\sum _{x}\Delta f(x)\Delta g(x)+C(x).}$

Definite summation by parts. For definite sums from ${\displaystyle a}$ to ${\displaystyle b}$, the formula becomes:

${\displaystyle \sum _{k=a}^{b}u(k)\Delta v(k)={\bigl [}u(b+1)v(b+1)-u(a)v(a){\bigr ]}-\sum _{k=a}^{b}v(k+1)\Delta u(k).}$

Complex analysis (exponential type)

Following the theory developed by Niels Erik Nørlund,^[4] the indefinite sum can be uniquely determined for analytic functions by imposing restriction on their growth in the complex plane. Specifically, by imposing minimal growth, the non-constant periodic terms can be filtered out.

Suppose ${\displaystyle f(z)}$ is analytic in a vertical strip containing the real axis, and let ${\displaystyle F(z)}$ be an analytic solution of ${\displaystyle F(z+1)-F(z)=f(z)}$ in that strip. To ensure uniqueness, require ${\displaystyle F(z)}$ to be of minimal growth, specifically to be of exponential type less than ${\displaystyle 2\pi }$ in the imaginary direction. That is, there exist constants ${\displaystyle M>0}$ and ${\displaystyle \epsilon >0}$ such that ${\displaystyle |F(z)|\leq Me^{(2\pi -\epsilon )|\Im (z)|}}$ as ${\displaystyle |\Im (z)|\to \infty }$.^[13][14]

Let ${\displaystyle F_{1}(z)}$ and ${\displaystyle F_{2}(z)}$ be two analytic solutions satisfying this growth condition. Their difference ${\displaystyle C(z)=F_{1}(z)-F_{2}(z)}$ is then analytic, 1‑periodic (i.e., ${\displaystyle C(z+1)=C(z)}$), and inherits the same exponential type less than ${\displaystyle 2\pi }$.

Nørlund uses a fundamental result in complex analysis (related to Carlson's theorem, the Phragmén–Lindelöf principle, and the Paley–Wiener theorem) which states that a non‑constant periodic entire function must have exponential type at least ${\displaystyle 2\pi }$.^[4] This follows from its Fourier series expansion: if ${\displaystyle C(z)}$ is non‑constant, its Fourier series contains a term ${\displaystyle a_{n}e^{2\pi inz}}$ with ${\displaystyle n\neq 0}$, which has type ${\displaystyle 2\pi |n|\geq 2\pi }$. Since ${\displaystyle C(z)}$ has type strictly less than ${\displaystyle 2\pi }$, it cannot contain any such term and therefore must be constant.

The exponential type less than ${\displaystyle 2\pi }$ in the imaginary direction on ${\displaystyle f}$ condition is sufficient but not strictly necessary. Nørlund's general definition of the principal solution is the analytic solution ${\displaystyle F}$ having the minimal possible exponential type for the given ${\displaystyle f}$.^[4] If ${\displaystyle f}$ has exponential type ${\displaystyle k}$ in imaginary direction, then the principal solution ${\displaystyle F(z)}$ will also have type ${\displaystyle k}$ in that strip, provided it converges. For example, ${\displaystyle f(z)=\sin(7z)}$ has exponential type ${\displaystyle 7}$; its principal solution exists and has type ${\displaystyle 7}$, even though ${\displaystyle 7>2\pi }$.

When ${\displaystyle f}$ has exponential type exactly ${\displaystyle 2\pi n}$ for some non‑zero integer ${\displaystyle n}$ in every strip where it is analytic (e.g. ${\displaystyle f(z)=\sin(2\pi nz)}$ has type ${\displaystyle 2\pi n}$; its antidifference contains ${\displaystyle \sin(\pi n)=0}$ in the denominator) the principal solution fails to exist (or is undefined everywhere) because it resonates with the kernel of the difference operator. In all other cases (i.e., when ${\displaystyle f}$ is meromorphic and its exponential type in some vertical strip is not an integer multiple of ${\displaystyle 2\pi }$) the principal solution exists and is uniquely determined by minimal exponential type.

Real analysis (higher‑order convexity)

In real analysis, the uniqueness condition can be given using higher‑order convexity, generalizing the Bohr-Mollerup theorem. For an integer ${\displaystyle p\geq 0}$, a function is called ${\displaystyle p}$-convex if its divided differences of order ${\displaystyle p}$ are non‑negative, and ${\displaystyle p}$-concave if those divided differences are non-positive. A function is called eventually ${\displaystyle p}$-convex (resp. eventually ${\displaystyle p}$-concave) if there exists ${\displaystyle M>0}$ such that it is ${\displaystyle p}$-convex (resp. ${\displaystyle p}$-concave) on the interval ${\displaystyle (M,\infty )}$.

Marichal and Zenaïdi proved the following uniqueness theorem, their method requiring the solution to be eventually ${\displaystyle p}$-convex or ${\displaystyle p}$-concave.^[15][16]

Theorem. Let ${\displaystyle p\geq 0}$ be an integer and let ${\displaystyle g:\mathbb {R} _{+}\to \mathbb {R} }$ satisfy ${\displaystyle \lim _{n\to \infty }\Delta ^{p}g(n)=0}$. If ${\displaystyle f:\mathbb {R} _{+}\to \mathbb {R} }$ is an eventually ${\displaystyle p}$-convex or eventually ${\displaystyle p}$-concave solution of ${\displaystyle \Delta f=g}$, then ${\displaystyle f}$ is uniquely determined up to an additive constant. Moreover, for any ${\displaystyle x>0}$,

${\displaystyle f(x)=f(1)+\lim _{n\to \infty }\left(\sum _{k=1}^{n-1}g(k)-\sum _{k=0}^{n-1}g(x+k)+\sum _{j=1}^{p}{\binom {x}{j}}\Delta ^{\,j-1}g(n)\right),}$

and the convergence is uniform on bounded subsets of ${\displaystyle \mathbb {R} _{+}}$.

Müller–Schleicher axiomatic method

In their paper How to Add a Noninteger Number of Terms^[5], Müller and Schleicher introduced an axiomatic approach to fractional summation with a real or complex number of terms. Their method extends the classical discrete sum

${\displaystyle \sum _{k=1}^{x}f(k)}$

to non-integer and complex upper limits ${\displaystyle x}$. The definition is built upon six natural axioms:

1. Continued Summation: ${\displaystyle \sum _{\nu =x}^{y}f(\nu )+\sum _{\nu =y+1}^{z}f(\nu )=\sum _{\nu =x}^{z}f(\nu )}$.
2. Translation Invariance: ${\displaystyle \sum _{\nu =x+s}^{y+s}f(\nu )=\sum _{\nu =x}^{y}f(\nu +s)}$.
3. Linearity: ${\displaystyle \sum _{\nu =x}^{y}(\lambda f(\nu )+\mu g(\nu ))=\lambda \sum _{\nu =x}^{y}f(\nu )+\mu \sum _{\nu =x}^{y}g(\nu )}$.
4. Empty Sum Condition: ${\displaystyle \sum _{\nu =1}^{1}f(\nu )=f(1)}$ (equivalent to the empty sum condition).
5. Holomorphy for Monomials: for each ${\displaystyle d\in \mathbb {N} }$, ${\displaystyle z\mapsto \sum _{\nu =1}^{z}\nu ^{d}}$ is holomorphic in ${\displaystyle \mathbb {C} }$.
6. Right-Shift Continuity: if ${\displaystyle f(z+n)\to 0}$ pointwise as ${\displaystyle n\to +\infty }$, then ${\displaystyle \sum _{\nu =x}^{y}f(\nu +n)\to 0}$; more generally, if ${\displaystyle f(z+n)}$ can be approximated by polynomials ${\displaystyle p_{n}(z+n)}$ of fixed degree with ${\displaystyle |f(z+n)-p_{n}(z+n)|\to 0}$, then:

${\displaystyle \left|\sum _{\nu =x}^{y}f(\nu +n)-\sum _{\nu =x}^{y}p_{n}(\nu +n)\right|\to 0}$.

Axioms S1–S4 force the sum to align with the ordinary finite sum when the limits are integers. Axiom S5 forces monomials to behave the same way under the generalization of fractional sums. Axiom S6 is the crucial axiom which allows one to "step back" the asymptotic region to determine the fractional sum in a finite interval. The exact conditions for the method to work are, as stated in the Definition 1.2 of the paper:

Let ${\displaystyle U\subset \mathbb {C} }$ and ${\displaystyle \sigma \in \mathbb {N} \cup \{-\infty \}}$. A function ${\displaystyle f:U\rightarrow \mathbb {C} }$ will be called fractional summable of degree ${\displaystyle \sigma }$ if the following conditions are satisfied:

• ${\displaystyle x+1\in U}$ for all ${\displaystyle x\in U;}$
• there exists a sequence of polynomials ${\displaystyle (p_{n})_{n\in \mathbb {N} }}$ of fixed degree ${\displaystyle \sigma }$ such that for all ${\displaystyle x\in U}$

${\displaystyle |f(n+x)-p_{n}(n+x)|\rightarrow 0}$ as ${\displaystyle n\rightarrow +\infty }$

• for every ${\displaystyle x,y+1\in U,}$ the limit

${\displaystyle \lim _{n\to \infty }\left(\sum _{\nu =n+x}^{n+y}p_{n}(\nu )+\sum _{\nu =1}^{n}{\bigl (}f(\nu +x-1)-f(\nu +y){\bigr )}\right),}$

exists.

In the simplest case when ${\displaystyle f(t)\to 0}$ as ${\displaystyle t\to \infty }$ (i.e., the approximating polynomials are zero), this reduces to:

${\displaystyle \nabla ^{-1}f(x)=\sum _{k=1}^{x}f(k)=\sum _{n=1}^{\infty }{\bigl (}f(n)-f(n+x){\bigr )}+C}$

Odd functions

If ${\displaystyle f(z)}$ is an odd function (${\displaystyle f(-z)=-f(z)}$), the unique analytic solution ${\displaystyle F(z)}$ satisfies:^[17]

${\displaystyle F(z)=F(-1-z).}$

This represents a point symmetry about ${\displaystyle z=-1/2}$.

Even functions

If ${\displaystyle f(z)}$ is an even function (${\displaystyle f(-z)=f(z)}$), the unique analytic solution ${\displaystyle F(z)}$ satisfies:^[17]

${\displaystyle F(z)+F(-1-z)=F(-1)}$.

Newton series

For an entire function of exponential type less than ${\displaystyle \ln \left(2\right)}$^[20] the inverse forward difference operator, ${\displaystyle \Delta ^{-1}f(x)}$, can be expressed by its Newton series expansion: ^[21][22]

${\displaystyle \sum _{x}f(x)=\sum _{k=1}^{\infty }{\binom {x}{k}}\Delta ^{k-1}f\left(0\right)+C(x)=\sum _{k=1}^{\infty }{\frac {\Delta ^{k-1}f(0)}{k!}}(x)_{k}+C(x).}$

${\displaystyle (x)_{k}={\frac {\Gamma (x+1)}{\Gamma (x-k+1)}}}$ is the falling factorial.

Bernoulli‑operator series expansion

Formally, the inverse forward difference operator can be expressed in terms of the derivative operator ${\displaystyle D={\frac {d}{dx}}}$ using the exponential generating function of the Bernoulli numbers:^[23][24][25]

$${\displaystyle \Delta ^{-1}={\frac {1}{e^{D}-1}}=\sum _{v=0}^{\infty }{\frac {B_{v}}{v!}}D^{\,v-1},}$$

where ${\displaystyle B_{v}}$ are the Bernoulli numbers defined by the generating function ${\displaystyle {\frac {t}{e^{t}-1}}=\sum _{v=0}^{\infty }B_{v}{\frac {t^{v}}{v!}}}$. Under this convention ${\displaystyle B_{1}=-{\tfrac {1}{2}}}$.

If ${\displaystyle f}$ is a polynomial, only finitely many terms of the series are non-zero as the finite difference of a monomial is a polynomial of one degree lower (following by induction, finitely many terms are required). For ${\displaystyle f(x)=x^{n}}$ one obtains the antidifference:^[24]

$${\displaystyle \sum _{x}x^{n}={\frac {B_{n+1}(x)}{n+1}}+C(x),}$$

where ${\displaystyle B_{n}(x)}$ are the Bernoulli polynomials of the first order.^[24]

If ${\displaystyle f}$ admits a Maclaurin series expansion ${\displaystyle f(x)=\sum _{n=0}^{\infty }{\frac {f^{(n)}(0)}{n!}}x^{n}}$, the antidifference of monomials in the series expansion yields the formal series:^[25]

$${\displaystyle \sum _{x}f(x)=\sum _{n=1}^{\infty }{\frac {f^{(n-1)}(0)}{n!}}B_{n}(x)+C(x).}$$

For non‑polynomials this expansion is generally asymptotic.

Relation to the inverse backward difference

If one instead expands the inverse backward difference operator, ${\displaystyle \nabla ^{-1}={\frac {e^{D}}{e^{D}-1}}}$ (which extends ${\displaystyle \sum _{k=1}^{x}f(k)}$), it admits to the same expansion, but with ${\displaystyle B_{1}=+{\tfrac {1}{2}}}$ in place of ${\displaystyle B_{1}=-{\tfrac {1}{2}}}$.

Euler–Maclaurin formula

The Euler–Maclaurin formula extends ${\displaystyle \nabla ^{-1}f(x)=\sum _{k=1}^{x}f(k)}$:^[6][13] $${\displaystyle {\begin{aligned}\nabla ^{-1}f(x)&=\int _{1}^{x}f(t)dt+{\frac {f(1)+f(x)}{2}}\\&\quad +\sum _{k=1}^{p}{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(x)-f^{(2k-1)}(1)\right)+R_{p}+C(x)\end{aligned}}}$$ where ${\displaystyle B_{2k}}$ are the even Bernoulli numbers, ${\displaystyle p}$ is an arbitrary positive integer, and ${\displaystyle R_{p}}$ is the remainder term given by:

${\displaystyle R_{p}=(-1)^{p+1}\int _{1}^{x}f^{(p)}(t){\frac {P_{p}(t)}{p!}}\,dt,}$

with ${\displaystyle P_{p}(t)=B_{p}(t-\lfloor t\rfloor )}$ being the periodized Bernoulli function related to the Bernoulli polynomials.

Laplace summation (Gregory summation formula)

Laplace's summation formula, closely related to the Gregory summation formula, can be seen as the discrete counterpart to the Euler–Maclaurin formula. The inverse forward difference ${\displaystyle \Delta ^{-1}f(x)}$:^{[26][27][12][28]}

${\displaystyle \sum _{x}f(x)=\int _{0}^{x}f(t)dt-\sum _{k=1}^{\infty }{\frac {c_{k}}{k!}}\Delta ^{k-1}f(x)+C(x)}$

where ${\displaystyle c_{k}=\int _{0}^{1}(x)_{k}dx}$ are the Cauchy numbers of the first kind.

${\displaystyle (x)_{k}={\frac {\Gamma (x+1)}{\Gamma (x-k+1)}}}$ is the falling factorial.

Truncating the series after ${\displaystyle n}$ terms leaves a remainder that can be expressed as an integral of ${\displaystyle f^{(n)}}$ times a periodic Bernoulli polynomial.^[12][28] In the notation of Charles Jordan, Gregory's formula is:^[12]

$${\displaystyle {\begin{aligned}\sum _{x=a}^{z}f(x)&=\int _{a}^{z}f(x)\,dx\\&\quad -\sum _{m=1}^{n}b_{m}{\bigl [}\Delta ^{m-1}f(z)-\Delta ^{m-1}f(a){\bigr ]}\\&\quad -b_{n}\,(z-a)\,\Delta ^{n}f(\xi ),\quad a<\xi <z,\end{aligned}}}$$ where the coefficients ${\displaystyle b_{m}}$ are the Bernoulli numbers of the second kind. Note the argument is without a shift, aligning with the inverse backward difference.

Abel–Plana formula

The indefinite sum ${\displaystyle \nabla ^{-1}f(x)=\sum _{k=1}^{x}f(k)}$ can be analytically continued by applying the standard Abel-Plana formula to the finite sum ${\displaystyle \sum _{k=1}^{n}f(k)}$ and then analytically continuing the integer limit ${\displaystyle n}$ to the variable ${\displaystyle x}$. This yields the formula:^[7] $${\displaystyle {\begin{aligned}\nabla ^{-1}f(x)&=\int _{1}^{x}f(t)dt+{\frac {f(1)+f(x)}{2}}\\&\quad +i\int _{0}^{\infty }{\frac {\left(f(x-it)-f(1-it)\right)-\left(f(x+it)-f(1+it)\right)}{e^{2\pi t}-1}}dt+C(x)\end{aligned}}}$$

This analytic continuation is valid when the conditions for the original formula are met. The sufficient conditions are:^[13][14]

1. Analyticity: ${\displaystyle f(z)}$ must be analytic in the closed vertical strip between ${\displaystyle \Re (z)=1}$ and ${\displaystyle \Re (z)=\Re (x)}$. The formula provides the analytic solution up to, but not beyond, the nearest singularities of ${\displaystyle f}$ to the line ${\displaystyle \Re (z)=1}$.
2. Growth: ${\displaystyle f(z)}$ must be of exponential type less than ${\displaystyle 2\pi }$ in this strip, satisfying ${\displaystyle |f(z)|\leq Me^{(2\pi -\epsilon )|\Im (z)|}}$ for some ${\displaystyle M>0}$, ${\displaystyle \epsilon >0}$ as ${\displaystyle |\Im (z)|\to \infty }$.