Computer Algebra for
Concrete Mathematics

Part 4: The Holonomic Paradigm

Peter.Paule@risc.jku.at

June 24, 2014, JKU

Some Symbolic Packages

The following packages of my RISC Combinatorics Group are used:

NOTE. First you have to pick them up at the web page given above...

<<RISC`fastZeil`

Fast Zeilberger Package version 3.60
written by Peter Paule, Markus Schorn, and Axel Riese
Copyright 1995-2009, Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

<<RISC`GeneratingFunctions`

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

InputForms: Binomials

(SubStar[(({{n_}, {k_}}))]):=Binomial[n,k]

(SubStar[(({{a}, {3}}))])

The Holonomic Paradigm

Holonomic Sequences

Definition (holonomic sequences).

(Subscript[((Subscript[f, n])), n≥0])  is   holonomic :⟺ 
∃ polynomials  (Subscript[p, i])(x),not all 0, such that

(Subscript[p, d])(n)(Subscript[f, n+d])+⋯+(Subscript[p, 0])(n)(Subscript[f, n])=0

for all n≥0 such that (Subscript[p, d])(n)≠ 0.

NOTE.  The recurrence (resp. the recurrence operator) has order d; its degree is
              the maximum of the degrees of the (Subscript[p, d])(n),..., (Subscript[p, 0])(n) .

Examples

((SubStar[(({{n+1}, {k}}))])/(SubStar[(({{n}, {k}}))])) /.k->7//FullSimplify

■  (Subscript[f, n,k]):=(SubStar[(({{n}, {k}}))]) :
(n+k+1)(Subscript[f, n+1,k])-(n+1)(Subscript[f, n,k])=0,
and (k+1)(Subscript[f, n,k+1])-(n-k)(Subscript[f, n,k])=0
■ (Subscript[f, n]):=(∑_(k=0)^n)(SubStar[(({{n-k}, {k}}))]) ((Subscript[f, n+2])-(Subscript[f, n+1])-(Subscript[f, n])=0)
■ Chebyshev: (Subscript[T, n])(x)
            ((Subscript[T, n+2])(x)-2x (Subscript[T, n+1])(x)+(Subscript[T, n])(x)=0)
■ Bessel: (Subscript[J, n])(x) 
      (x (Subscript[J, n+2])(x)-2n (Subscript[J, n+1])(x)+x (Subscript[J, n])(x)=0)

?*ChebyshevT*

?*BesselJ*

\[SelectionPlaceholder] Binary trees

t[n_]:=(1/(n+1))(SubStar[(({{2n}, {n}}))])

((t[n+1])/(t[n]))//FullSimplify

4-(6/(2+n))

Hypergeometric Sequences and Gosper Algorithm

Definition: hypergeometric sequences

(Subscript[((Subscript[t, n])), n≥0])  is   hypergeometric :⟺ 
∃ non-zero polynomials  (Subscript[p, i])(x) such that

(Subscript[p, 1])(n)(Subscript[t, n+1])+(Subscript[p, 0])(n)(Subscript[t, n])=0

for all n≥0 such that (Subscript[p, 1])(n)≠ 0.

\[SelectionPlaceholder] Geometric sequences

t[n_]:=(w^n)

((t[n+1])/(t[n]))//FullSimplify

Gosper' s Algorithm

Input: hypergeometric sequence (Subscript[((Subscript[t, n])), n≥0])

Output: hypergeometric sequence (Subscript[((Subscript[s, n])), n≥0]) such that

(Subscript[t, n]) = (Subscript[s, n+1])-(Subscript[s, n])

for all n≥0 sufficiently large; if no such solution exists,
       the algorithm detects this.

NOTE.  Gosper’s algorithm is implemented in the RISC package `fastZeil`:

?Gosper

Unprotect[Show];Gosper[(w^k),{k,0,n}]

Prove[]

EXAM QUESTION: Explain why this is delivering a proof of the `Gosper output'!

Holonomic Power Series (Functions)

Definition (holonomic power series, resp. functions).

y(x)=(∑_(n=0)^∞)(Subscript[f, n])(x^n)    is   holonomic  :⟺
∃ polynomials  (Subscript[p, i])(x),not all 0, such that

(Subscript[p, d])(x)(y^((d)))(x)+ ⋯+(Subscript[p, 0])(x)y(x)=0

Examples

■ exp(x) (y'-y=0), ■ sin(x) &  cos(x) ((y^(,,))-y=0),
■ log(1+x) ((1+x)y'=0),
■ Chebyshev  (Subscript[T, n])(x): 
(1-(x^2))(y^(,,))-x y'+(n^2)y=0,
■ Bessel  (Subscript[J, n])(x):
                (x^2)(y^(,,))+x y'+((x^2)-(n^2))y=0,
■ (∑_(n=0)^∞)(Subscript[f, n])(x^n)  where (Subscript[((Subscript[f, n])), n≥0]) is hypergeometric

\[SelectionPlaceholder] Example:binary trees (again)

t[n_]:=(1/(n+1))(SubStar[(({{2n}, {n}}))])

((t[n+1])/(t[n]))//FullSimplify //Factor

BinTreeDE=RE2DE [{(2+n) a[n+1]-2(1+2n)a[n]==0,a[0]==1},a[n],y[x]]

DSolve[{-1-(-1+2 x) y[x]-(-x+4 (x^2)) y'[x]==0,y[0]==1},{y[x]},{x}]

KEY CONNECTION: holo sequences ⟷ holo power series (resp. functions)

Theorem:

y(x)=(∑_(n=0)^∞)(Subscript[f, n])(x^n)    is   holonomic  ⟺ (Subscript[((Subscript[f, n])), n≥0])  is   holonomic

Example (binary trees again)

y(x)=(∑_(n=0)^∞)(1/(n+1))(SubStar[(({{2n}, {n}}))])(x^n)=((1-(Sqrt[1-4 x]))/(2 x))

Recall

BinTreeDE

DE2RE[BinTreeDE,y[x],a[n]]

KEY MECHANISM: algorithmic closure properties

Closure properties

Theorem: Given holonomic power series (resp. functions)
(Subscript[y, 1])(x)= (∑_(n≥0))(Subscript[f, n])(x^n), (Subscript[y, 2])(x)= (∑_(n≥0))(Subscript[g, n])(x^n) then the following objects
are also holonomic:

(Subscript[y, 1])(x)+(Subscript[y, 2])(x), (Subscript[y, 1])(x)(Subscript[y, 2])(x),(∑_(n≥0))(Subscript[f, n])(Subscript[g, n])(x^n),D (Subscript[y, 1])(x),
∫(Subscript[y, 1])(x)dx, (Subscript[y, 1])((Subscript[y, 2])(x)) [(Subscript[y, 2])(x) alg. &  (Subscript[y, 2])(0)=0],
(∑_(n≥0))(x^n) (∑_(k=0)^n)(Subscript[f, k])

?*AE*

Examples

■ holonomic:((e^x)/(Sqrt[1-4x])) ,(e^(x/(Sqrt[1-4x])));  (∑_(n≥0))(x^n) ((((∑_(k=0)^n)((SubStar[(({{n}, {k}}))])^3)))^2)  ,...

NOTE: There are also non-holonomic objects, 
e.g, (x/((e^x)-1)) (Bernoulli) or (Sqrt[n]),(n^n),...

■ (e^(x/(Sqrt[1-4x]))) is holonomic:

NOTE. f(x)=(e^x) is holonomic: f'(x)-f(x)=0.

NOTE: f(x)=(x/(Sqrt[1-4x])) is algebraic, and thus holonomic (because of closure properties):

(1-4x)f(((x))^2)-(x^2)=0.

AE2DE[(1-4x)((f[x])^2)-(x^2)==0,f[x]]

?ACompose

ACompose[f'[x]-f[x]==0,(1-4x)((f[x])^2)-(x^2)==0,f[x]]

Proving with Holonomic Closure Properties

A 21 st Century Proof of Pythagoras

sinDE = {y''[x]+y[x]==0,y[0]==0,y'[0]==1}

cosDE = {y''[x]+y[x]==0,y[0]==1,y'[0]==0}

?DECauchy

DECauchy[sinDE,sinDE,y[x]]

DECauchy[cosDE,cosDE,y[x]]

DEPlus[%7,%8,y[x]]

DEPlus[%10,y[x]==-1,y[x]]

?DEPlus

The IQ Problem

I.Q. Tests
Setzen Sie die Reihe fort: 1  1  2  3  5  8  13  21  ?
Es gibt zwei Lösungsmoeglichkeiten, eine leichte und eine schwierigere. Versuchen Sie, ob Sie beide finden koennen.
[Aufgabe 13, Denksport I für Superintelligente; ``Check Your Own I.Q.’’, Hans J. Eysenck, 1966]

BOOK SOLUTIONS

“34. (Leicht. Jede Zahl wird durch Subtraktion der folgenden von der naechstfolgenden Zahl gebildet: 2-1 = 1 usw. bis 34-21 = 13, also ist die fehlende Zahl 34.

Schwierig. Das Quadrat jeder Zahl unterscheidet sich um 1 vom Produkt der Zahlen rechts und links von ihr: (1^2)=1, 2×1=2; (2^2)=4, 1×3=3; usw. bis (21^2)=441, 13×34=442.)"

{1,1,2,3,5,8,13,21,34}

Note. This corresponds to Cassini’s Identity (1680)”:

F[n+1]F[n-1]-((F[n])^2) = (((-1))^(n+1)),         n⩾1.

We prove this automatically with the GeneratingFunctions package!

The Fibonacci sequence (Leonardo Fibonacci, 1202)

GuessRE[{1,1,2,3,5,8},F[k]]

re={-F[k]-F[1+k]+F[2+k]==0,F[0]==1,F[1]==1};

RSolve[re,F[k],k]

F[k]=(1/(Sqrt[5])) ( (((((1+(Sqrt[5]))/2)))^(k+1)) - (((((1-(Sqrt[5]))/2)))^(k+1)) )

?RSolve

RSolve[- F[k]-  F[1+k]+t F[2+k]==0,F[k],k]

Cassini Identity: Computer Proof 1

Die ``schwierige'' Loesung entspricht der Identitaet: 
Fuer n≥0 (mit F[-1]=0) gilt

F[n+1] F[n-1] -((F[n])^2) = (((-1))^(n+1)).

Preprocessing:  sequence ->  defining recurrence

Clear[F]

Sequence  (Subscript[(F[n]), n≥0])  :

re = {-F[n]-F[1+n]+F[2+n]==0,F[0]==1,F[1]==1}

?RESubsequence

Sequence  (Subscript[(F[n+1]), n≥0])  :

rePlus = RESubsequence[re, F[n], n+1]

Sequence  (Subscript[(F[n-1]), n≥0])  :

reMinus = RESubsequence[re, F[n], n-1]

Sequence  (Subscript[((((-1))^n)), n≥0])  :

reAlternate = {F[n+1]== -F[n],F[0]==1}

Hadamard Product of Sequences , resp.  Recurrences

?REHadamard

Sequence  (Subscript[(F[n+1]F[n-1]), n≥0])  :

rePlusMinus = REHadamard[rePlus,reMinus,F[n]]

Sequence  (Subscript[(((F[n])^2)), n≥0])  :

reSquared = REHadamard[re,re,F[n]]

Sequence  (Subscript[(-((F[n])^2)), n≥0])  :

reSquaredMinus = REHadamard[reSquared,{F[n]==-1},F[n]]

Proof Endspurt

Sequence  (Subscript[(F[n+1]F[n-1]-((F[n])^2)), n≥0])  :

reDiff = REPlus[rePlusMinus,reSquaredMinus,F[n]]

Sequence  (Subscript[(F[n+1]F[n-1]-((F[n])^2)+(((-1))^n)), n≥0])  :

reOut = REPlus[reDiff,reAlternate,F[n]]

The recurrence has order 3. Since all initial values are 0, 
ALL values of the sequence must be  0.                                  Q.E.D.

Cassini Identity: Computer Proof 2

Computer proof of Cassini identity

Recall

Clear[F];re

Proof.

f=DefineS[re,F[k]]

Shift[f,1]*Shift[f,-1]-f^2==DefineS[{a[n+1]==-a[n],a[0]==-1},a[n]]

More on Symbolic Summation

Consider the sum

(∑_(k=0)^(n-2))k k!

(-1)!

For which values of n is this answer correct?

Alternatively, we can use the procedure from the RISC package `fastZeil`:

Unprotect[Show];Gosper[k k!,{k,0,n-2}]

Note that -in contrast to Mathematica-  the output of `Gosper` includes 
additional information about the domain of admissible n!

Execute

Prove[]

EXAM QUESTION: Explain why this is delivering a proof of the `Gosper output'.

\[SelectionPlaceholder] Another Gosper example: binary trees

t[n_]:=(1/(n+1))(SubStar[(({{2n}, {n}}))])

Gosper[t[k] t[n-k-1],{k,0,n-1}]

An Automatic Proof of the Binomial Theorem

We prove the  binomial theorem in the form

(∑_(k=0)^n)(({{n}, {k}})) (x^k)(y^(n-k))=(((x+y))^n)

Set

g[n_] := (∑_(k=0)^n)(SubStar[(({{n}, {k}}))]) (x^k)(y^(n-k))

{g[0],g[1],g[2],g[3],g[4],g[5],g[6],g[7],g[8]} //Factor

Proof.

Unprotect[Show];Zb[(SubStar[(({{n}, {k}}))])(x^k)(y^(n-k)),{k,0,n},n,1]

Prove[]

NOTE. What we use here is Zeilberger' s extension of Gosper' s algorithm. 
            This is needed because a direct application of Gosper fails:

Gosper[(SubStar[(({{n}, {k}}))])(x^k)(y^(n-k)),{k,0,n}]