View source for OEIS/Holonomic Project

This project tries to define as many existing [https://oeis.org OEIS] sequences as possible by '''holonomic recurrences'''. 

:Whereas the usual linear recurrence has constant coefficients, a holonomic recurrence has coefficients that are polynomials in <code>n</code> (the sequence index). Holonomic recurrences can be derived from linear differential equations, and are applicable to a wide class of series - see [https://en.wikipedia.org/wiki/Holonomic_function Wikipedia: Holonomic function].
===Example===
 [https://oeis.org/A000957 A000957] Fine's sequence (or Fine numbers): 
   number of relations of valence >= 1 on an n-set; 
   also number of ordered rooted trees with n edges having root of even degree.
 0, 1, 0, 1, 2, 6, 18, 57, 186, 622, 2120, 7338, 25724, 91144, 325878, ...
 G.f.: (1-sqrt(1-4*x))/(3-sqrt(1-4*x))
 D-finite with recurrence: 2*n*a(n) +(12-7*n)*a(n-1) +2*(3-2*n)*a(n-2)=0. - R. J. Mathar
===Types of recurrences===
There is a hierarchy of sequence types which are included in this project, with increasing complexity and computational effort:
* constant
* periodic
** finite (followed by infinitely many zeroes), polynomials
** continued fractions of sqrt(n)
* recurrent with constant coefficients (cf. [https://oeis.org/wiki/Index_to_OEIS:_Section_Rec OEIS Index entries for linear recurrences with constant coefficients])
** defined by a rational generating function (fraction of two polynomials in x), among them sequences for Coexter groups and (lattice) coordination sequences
* holonomic recurrent (D-finite, hypergeometric functions, algebraic equations)
====Initial terms====
For all types the recurrence frequently starts only after some initial terms which do not participate in the recurrence equation.
====Exponential generating function====
The project also handles e.g.f.s, in which case the terms are multiplied by <code>n!</code> during the generation.of the sequence.
===Embedding in jOEIS===
[https://github.com/archmageirvine/joeis jOEIS] is a project which aims to implement all OEIS sequences in Java. Currently (Feb. 2021) there are more than 104,000 sequences implemented. They are compiled in a single jar file, and any sequence can be called repeatedly in a uniform way to produce the next terms.

In jOEIS there is a series of basic Java classes which generate sequences of the holonomic type. They can all be found under [https://github.com/archmageirvine/joeis/tree/master/src/irvine/oeis https://github.com/archmageirvine/joeis/tree/master/src/irvine/oeis]:
 BlockMultAddSequence
 ContinuedFractionOfSqrtSequence
 CoordinationSequence
 CoxeterSequence
 FiniteSequence
 GeneratingFunctionSequence
 HolonomicRecurrence
 LatticeCoordinationSequence
 LinearRecurrence
 PaddingSequence
 PeriodicSequence
 PrependSequence
Currently the following numbers of sequences are covered by this project:
     36 block
   2911 conti
    946 coord
   2251 coxet
   4601 finit
   8520 gener
   5902 holon
    340 latti
  20035 linea
    336 paddi
   1202 perio
    700 prepe
 ------------
  47780 total
===Program for holonomic recurrences===
Recurrences are attractive because it is rather simple to implement them in any programming language. A CAS may have special functions for them (e.g. LinearRecurrence and RecurrenceTable in Mathematica). The algorithm is straightforward, at it needs only fixed space even in the holonomic case.
In jOEIS, [https://github.com/archmageirvine/joeis/tree/master/src/irvine/oeis/HolonomicRecurrence.java HolonomicRecurrence.java] uses a ring buffer of size <code>k</code> for the recurrence variables <code>a(n), a(n-1), ... a(n-k+1)</code>, and a second array where the powers of the index <code>n</code> are maintained. The length of this array corresponds with the degree of the equation (the highest power of <code>n</code>).
The program is not too big, and could be rewritten in other languages, for example Python.
===Parameterization===
The program uses these parameters:
* recurrence equation
* initial terms
* index of the first term to be returned (OEIS offset1)
====Equation====
Given that the recurrence equation is written in the form of an ''annihilator'':
 p[0]*1 + p[1]*a[n-k+1] + p[2]*a[n-k+2] + ... + p[k-1]*a[n-k+k-1] + p[k]*a[n] = 0
then the corresponding parameter is represented by the array
 [p[0], p[1], p[2], ... p[k]]
where each polynomial p[i] is represented by an array [c0, c1, c2 ... cm] that defines the exponents of <code>n</code>. 
The initial example equation
 2*n*a(n) +(12-7*n)*a(n-1) +2*(3-2*n)*a(n-2)=0
is represented by:
 [[0],[6,-4],[12,-7],[0,2]]
=====Constant term=====
The first polynomial allows for a constant term, but normally the homogenous form is used with <code>p[0] = [0]</code>. 
=====Index shifting=====
Usually the last recurrence variable is <code>a(n)</code>, but it is also possible to ''shift'' the recurrence index, for example to have
  2*a(n) - 2*a(n+1) - a(n+2) = 0
Such a shifting is specified by an additional parameter <code>dist</code> which is 0 by default. The example would get <code>dist=2</code>.
In the future, all equations will be recomputed such that the constant term as well as the shift become always zero.
====Initial terms====
For an ''order'' (number of variables - 1) <code>k</code> a recurrence normally needs <code>k</code> initial terms to be computable. With the current implementation it is possible 
* that there are more initial terms, in which case the recurrence only applies to the last <code>k</code> terms,
* that there are fewer initial terms, in which case a suitable number of zeroes is implied before the first term.
=====''Black startable'' sequences=====
As a special possibility a recurrence can be started "out of nothing", that is the initial terms can be omitted altogether. The implementation then implies a leading sequence term one, and a suitable number of zeroes before. We call this property ''black startable'' (an expression used for power stations), and there is a surprising number of sequences that have it.
===Database table===
The main result of the project is a database table with suitable parameter sets for the OEIS sequences that are covered. The table is create with the following MySQL statements:
 --  Table for holonomic recurrences in the OEIS
 DROP    TABLE  IF EXISTS holref;
 CREATE  TABLE            holref
    ( aseqno   VARCHAR(10) NOT NULL  -- A-number
    , callcode VARCHAR(32)           -- "holos", generate calls for HolonomicRecurrence
    , offset1  INT                   -- offset
    , matrix   VARCHAR(8192)         -- coefficients of polynomials in n
    , init     VARCHAR(4096)         -- initial terms for the recurrence
    , dist     INT                   -- optional shift of the index n (0)
    , gftype   INT                   -- 0 for o.g.f, 1 for e.g.f
    , keyword  VARCHAR(128)          -- list of keywords for source sequence type etc.
    , PRIMARY KEY(aseqno, callcode)
    );
 COMMIT;
The tab-separated data loaded into the database table can be found under '''[http://www.teherba.org/OEIS-mat/holref.txt holref.txt]'''.

=== Parameters of hypergeometric series ===
There is a number of OEIS sequences that can be generated by a recurrence
 f(n) * a(n) = g(n) * a(n-1)
For example, formulae like
 [https://oeis.org/A217364 A217364] a(n) = 2^n*binomial(5*n, n)/(4*n+1).
 [https://oeis.org/A274136 A274136] a(n) = (n+1)*(2*n+2)!/(n+2)
 [https://oeis.org/A275655 A275655] a(n) = binomial(6*n,3*n)*binomial(2*n,n).
are of this type.
In these cases, the parameters for <code>HolonomicRecurrence</code> can be computed by the following Maple procedure:
 with(StringTools):
 holo := proc(as,f); local b, num, den, fd, str, filename;
   b := simplify(simplify(-f(n)/f(n-1), GAMMA));
   printf("***** b = %a\n", b);
   num := seq(coeff(numer(b),n,j),j=0..degree(numer(b)));
   den := seq(coeff(denom(b),n,j),j=0..degree(denom(b)));
   str := cat("[[0],", DeleteSpace(convert([num],string)), ",", DeleteSpace(convert([den],string)), "]");
   printf("%a\tholos\t0\t\%s\t[%a]\t0\t\n", as, str, f(0));
 end;
We obtain, for example:
 holo(A217364, proc(n) 2^n*binomial(5*n, n)/(4*n+1) end);
 ***** b = (-3125*n^4+6250*n^3-4375*n^2+1250*n-120)/(128*n^4-64*n^3-8*n^2+4*n)
 A217364 holos   0       [[0],[-120,1250,-4375,6250,-3125],[0,4,-8,-64,128]]     [1]    0
Now we follow the Maple file above by an empty line and a ''pattern specification'' line:
 holo($(PARM0),proc(n) $(PARM3) end proc);
Then it forms a pattern file for the Perl program [http://github.com/gfis/OEIS-mat/common/callmaple.pl callmaple.pl] which merges the pattern file with a list of recurrence expressions in "seq4" format, and the output is in the same format with all recurrence expressions replaced by the parameters suitable for <code>HolonomicRecurrence</code>.