This is a spot of Pari/GP code to guess a linear recurrence from a vector of numbers (or vector of polynomials for guessing with further parameters). The result is pretty printed. For example,
read("recurrence-guess.gp"); recurrence_guess([1, 3, 9, 25, 65, 161, 385, 897, 2049, 4609]); => Recurrence length=3 coefficients v[n-3]* [4, -8, 5] *v[n-1] = v[n] v[n] = v[n-1]* [5, -8, 4] *v[n-3] characteristic polynomial x^3 - 5*x^2 + 8*x - 4 = factors (x - 2)^2 roots 2.00000 x - 1 roots 1.00000 generating function (1 - 2*x + 2*x^2)/(1 - 5*x + 8*x^2 - 4*x^3) = (1 - 2*x + 2*x^2) / ( (1 - x) * (1 - 2*x)^2 ) = partial fractions 1/(1 - x) - 1/(1 - 2*x) + 1/(1 - 2*x)^2 as powers n * 2^n + 1 OEIS %H <a href="/index/Rec#order_03">Index entries for linear recurrences with constant coefficients</a>, signature (5,-8,4). %F a(n) = 5*a(n-1) - 8*a(n-2) + 4*a(n-3).
The guess is found by a simple
matsolve(). Linear recurrences
include powers, polynomials, and polynomials times powers. Values given can
themselves be GP polynomials for parameterization, or (very) limited symbolic
calculation, or bivariate gfs guessed on one variable.
recurrence-guess.gp(52k, and sig)
recurrence-guess-15.tar.gz(53k, and sig)
tar file includes some self-tests, and the following
examples/polmod.gp script illustrating linear recurrence
t_POLMOD, which is efficient and compact but a
To install so
recurrence_guess() is always available
recurrence-guess.gp in say your
~/gp directory (which is in the GP
then in file
Give a full path (possibly starting
~/) if installed somewhere else.
Similar code can be found in
bestapprPade(Ser(vec))gives a generating function. On a long recurrence, sometimes
lindep()seems much faster than
matsolve()(would intend to use that if so). The nice output is the tedious part. Of course "nice" is a matter of personal preference and the output is still quite mechanical.
This page Copyright 2016, 2017, 2018, 2019, 2020, 2021 Kevin Ryde, except for the GPLv3 logo which is Copyright Free Software Foundation and used here in accordance with its terms.
(Back to the sitemap, or the PARI/GP section there)