\\ Copyright 2014, 2015, 2016, 2017 Kevin Ryde
\\
\\ recurrence-guess.gp is free software; you can redistribute it and/or
\\ modify it under the terms of the GNU General Public License as published
\\ by the Free Software Foundation; either version 3, or (at your option)
\\ any later version.
\\
\\ recurrence-guess.gp is distributed in the hope that it will be useful,
\\ but WITHOUT ANY WARRANTY; without even the implied warranty of
\\ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
\\ Public License for more details.
\\
\\ You should have received a copy of the GNU General Public License along
\\ with recurrence-guess.gp. See the file COPYING. If not, see
\\ .
\\ Usage: read("recurrence-guess.gp"); \\ loads pol-pfrac.gp too
\\
\\ recurrence_guess([0, 1, 5, 19, 65, 211, 665, 2059, 6305, 19171]);
\\ =>
\\ prints a report
\\
\\ recurrence_guess(v) takes a vector of numbers (and some other types) and
\\ attempts to find a linear recurrence which generates them. It prints a
\\ report of
\\
\\ * Recurrence factors
\\ * Characteristic polynomial
\\ * Generating function
\\ * Power expression, if possible
\\
\\ Other Modules Required:
\\
\\ pol-pfrac.gp version 2,
\\ http://user42.tuxfamily.org/pari-pol-pfrac/index.html
\\
\\ Values Length:
\\
\\ A recurrence of m terms requires minimum 2*m+1 values. The first 2*m
\\ determine the recurrence factors and 1 further value confirms it. See
\\ the help string of variable "recurrence_guess_minimum_confirm" to adjust
\\ how much confirmation is demanded.
\\
\\ If you know a maximum order m and give 2*m+1 values then something found
\\ at or below m will be correct. You might know a maximum for example if
\\ values are from a set of m mutual recurrences or a state machine of m
\\ states. For reference, a few common recurrence orders are
\\
\\ * power 2^k etc = order 1
\\ * polynomial to n^p = order p+1, eg. n^2 is order 3
\\ * sum of sequences <= sum of their orders
\\ (a union of the roots of the characteristic polynomials)
\\ * cumulative = order+1
\\
\\ When forming a complicated sum or cumulative total of recurrences etc it
\\ can be convenient to let the guess find resulting powers or generating
\\ function rather than working through formulas.
\\
\\ If you don't know a maximum m then it's always possible initial values
\\ satisfy a recurrence but more values, if you had them, would not or would
\\ need a longer recurrence. If many confirming values pass the recurrence,
\\ meaning an order m with 2*m much smaller than the vector v, then that
\\ gives some confidence, but no proof.
\\
\\ Output:
\\
\\ The output is a little rough and will change. The intention is to show
\\ interesting forms, usually to be cut and pasted into code or a document
\\ if useful.
\\
\\ Recurrence factors are shown high to low and low to high, with
\\ pseudo-expressions to try to make the direction clear. It's probably
\\ usual to write high to low for human readers, but low to high suits
\\ program code since the indices match a vector of values going low to
\\ high.
\\
\\ Roots of the characteristic polynomial are given since the values grow as
\\ a sum of powers of those. The power expression uses those roots when
\\ they're exact (see Value Types below).
\\
\\ A generating function (an "ordinary" generating function) is given in
\\ plain and partial fractions form. Generating functions are helpful for
\\ manipulating recurrences since they effectively represent all values of
\\ the sequence. The partial fractions can show when the recurrence is a
\\ sum of smaller recurrences, or recurrence plus a constant or periodic
\\ term. But it could happen that the full gf is simpler.
\\
\\ Another use for the gf partial fractions is to show what types of
\\ sequences might be summed to give values. For example 1/(1-2*x)^2 is 2^n
\\ and n*2^n. If you have some sequences and would like multiples of them
\\ to sum to the values then lindep() can find coefficients. lindep() can
\\ take a matrix of values, or vector of generating functions. lindep()
\\ won't act on t_RFRACs (as of gp 2.7.x) so for generating functions
\\ multiply by a common denominator to make plain polynomials.
\\
\\ The power expression is derived from the gf partial fractions. There's
\\ an internal check that what's shown does in fact give the values, but the
\\ mangling is slightly subtle so you might double check. The parts
\\
\\ vector_modulo([7,5,8],n)
\\
\\ mean 7,5,8 according as n==0,1,2 mod 3 respective, etc. Periodic factors
\\ like this are numerator terms of the gf partial fractions. The gf
\\ partial fractions can break down to small terms where their periodic
\\ nature is not obvious. The power expression combines terms of the same
\\ power and polynomial factor.
\\
\\ The power expression is written as index n=0 for the first value v[1].
\\ This corresponds to the generating function first term x^0. If a start
\\ at n=1 or similar is desired then try an initial dummy 0 in the input
\\ vector v. This will probably show as an "initial exception" (see below)
\\ which can be ignored.
\\
\\ The power expression always uses positive base powers, and usually reals
\\ (depending on the input values type). It's possible something like
\\ (-2)^k and negating corresponding odd periodic factors could simplify.
\\ But if a base appears in a few terms then for consistency you might
\\ prefer them all the same.
\\
\\ The power expression uses 3^floor(n/8) etc for half-power terms. If
\\ there's periodic factors then it might be possible to simplify by
\\ changing floor to ceil, or shifting to say floor((n+3)/8). But floor is
\\ printed for consistency since several terms might be similar.
\\
\\ Polynomial n^p factors use rationals, not floor() like the powers.
\\ Occasionally a result might be simplified by floor(). For example
\\
\\ recurrence_guess(vector(100,n,n--; floor(n/3)))
\\ =>
\\ 1/3*n + vector_modulo([0, -1/3, -2/3],n)
\\
\\ and you can notice the constant parts push 1/3*n down to floor(1/3*n).
\\ If a limit f(n)/n as n->infinity is wanted then plain high term may be
\\ clearer.
\\
\\ Initial Exceptions:
\\
\\ If a few initial values don't satisfy the recurrence then they become low
\\ 0s in the recurrence factors. For example
\\
\\ \\ exceptions then Fibonacci
\\ recurrence_guess([99, 99, 99, 1, 2, 3, 5, 8, 13, 21, 34]);
\\ =>
\\ v[n-5]* [0, 0, 0, 1, 1] *v[n-1] = v[n]
\\
\\ In the gf partial fractions the initial values become constant terms (not
\\ t_RFRACs), such as 98 + 99*x + 98*x^2. Those terms are differences from
\\ what the recurrence extended backwards would have given. The power
\\ expression shows initial exceptions explicitly.
\\
\\ For the values length, each initial exception requires 1 further value in
\\ addition to itself for guessing (because it lengthens the recurrence
\\ terms by 1 and each recurrence term requires 2 values). If guessing from
\\ just a few values then it might help to try omitting initial exceptions
\\ or possible exceptions.
\\
\\ High 0s in the recurrence factors are a "delay" in the recurrence.
\\ A simple example is
\\
\\ recurrence_guess([1,2,3, 10,20,30, 100,200,300]);
\\ =>
\\ v[n-3]* [10, 0, 0] *v[n-1] = v[n]
\\ vector_modulo([1, 2, 3],n) * 10^floor(n/3)
\\
\\ Each term is 10* its third previous. Simple repeated values come out
\\ similarly, for example [1,1, 2,2, 4,4, 8,8].
\\
\\ Value Types:
\\
\\ Values can be any exact Pari/GP type, so integers, rationals, and complex
\\ numbers or quads with integer or rational parts, and polynomials with any
\\ of those as coefficients. GP 2.7.x cannot mix quads of different
\\ discriminant or quads with complex, so be sure all terms are compatible.
\\
\\ The factors in the recurrence will be the same type as the values
\\ usually. Quads print as "w" in the usual way, so when cutting and
\\ pasting output be sure to change that to a desired name.
\\
\\ The gf partial fractions and the power expression also use the same type
\\ as the values. For example on rational values a gf shows say 1/(1-2*x^2)
\\ and the power expression is 2^floor(n/2), but if the values include a
\\ sqrt2=quadgen(8) then the gf can go down to 1/(1+sqrt2*x) and
\\ 1/(1-sqrt2*x), and the power expression uses sqrt2^n. Putting a factor
\\ of sqrt2 or I or similar through the values can induce the further
\\ breakdown if desired, though sometimes that can make periodic terms less
\\ clear. Maybe there should be a better way to express a base type.
\\
\\ For complex and quads it may be worth trying real and imaginary parts
\\ separately to see if they have simpler separate forms. Powers like
\\ (3+4*I)^n don't simplify. They have the same recurrence for real and
\\ imaginary parts (the two being a pair of Lucas sequences).
\\
\\ Values can be t_POL polynomial forms in some variable or multiple
\\ variables, but not 'x since that's used internally and for the generating
\\ function. The effect of t_POL values is to find a recurrence satisfying
\\ the sequence y^0 polynomial terms, and the sequence of y^1 terms, etc.
\\ This might be useful if values are something parameterized and a
\\ recurrence or power for the general case is desired. It may be worth
\\ trying terms separately to see if they are simpler individually.
\\
\\ Values which are powers of t_POL powers of a variable give that variable
\\ in the recurrence coefficients. For example
\\
\\ recurrence_guess(vector(10,n,n--; 's^n))
\\ =>
\\ v[n] = v[n-1]* s
\\
\\ Pari as of 2.9.x doesn't cannot factorize multi-variable polynomials (a
\\ polynomial with polynomial coefficients) so the characteristic polynomial
\\ is only shown unfactored and no generating function partial fractions.
\\
\\ Polynomial Sequences:
\\
\\ Values from a polynomial expression like n^2+n+1 = [1,3,7,13,21,31,43]
\\ give linear recurrences and can be found by recurrence_guess(). If you
\\ know the values are polynomial then use polinterpolate() instead as it's
\\ faster and requires only p+1 values for terms up to n^p whereas a
\\ recurrence is order p+1 so requires 2*(p+1)+1 values. If you have enough
\\ values then recurrence_guess() has the advantage of noticing powers and
\\ initial exceptions too.
\\
\\ Recurrence terms for a polynomial tend to be unenlightening. The
\\ characteristic polynomial has root 1 repeated, and similar in the
\\ generating function partial fraction denominators. For example n^2
\\ appears as characteristic polynomial (x-1)^3 and gf partial fractions
\\ 1/(1-x), 1/(1-x)^2, 1/(1-x)^3.
\\
\\ In the gf partial fractions, a term 1/(1-x)^k is values
\\ (n+1)*(n+2)*...*(n+(k-1)) / (k-1)!. This means that values like n^2 come
\\ out as a sum of various such gf parts which cancel at n powers other than
\\ n^2. Such a sum might be helpful if you are in fact working in such
\\ polynomial products.
\\
\\ Powers times a polynomial like (n+1)*7^k go similarly but the repeated
\\ root is then power 7 etc and gf partial fractions are like 1/(1-7*x)^2
\\ etc.
\\
\\ Program Use:
\\
\\ Guessing is done (and nothing printed) by
\\
\\ f = recurrence_guess_values_to_factors(v);
\\ g = recurrence_guess_values_to_gf(v);
\\
\\ returning a vector of factors or t_RFRAC generating function
\\ respectively. See the help strings for details.
\\ recurrence_guess_values_to_factors(v) is the basic operation and the gf
\\ follows from
\\
\\ g = recurrence_guess_factors_and_values_to_gf(f,v);
\\
\\ The code works by a simple matsolve() to find f
\\
\\ / v[1] v[2] v[3] \ / f[1] \ / v[4] \
\\ matsolve | v[2] v[3] v[4] | | f[2] | = | v[5] |
\\ \ v[3] v[4] v[5] / \ f[2] / \ v[6] /
\\
\\ This asks to find recurrence factors f so these f times consecutive v
\\ values give the next value, for m many such next values. This matsolve
\\ is attempted for orders m=1 upwards until finding a solution which
\\ satisfies the rest of v too. The code for this is only a few lines. The
\\ bloat is all in pretty-printing the result.
\\
\\ A small to medium recurrence is normally found and confirmed quickly.
\\ A very long v without a solution might take a long time.
\\
\\ Bill Allombert points out too that bestapprPade(Ser(v)) can give a
\\ generating function for values v. This can be faster than matsolve(),
\\ but tends to give a big expression for a non-recurrence so its result
\\ should be confirmed.
\\
\\ Compatibility:
\\
\\ Requires gp 2.7 or higher due to ".." subvectors and probably more.
\\
\\ Bugs:
\\
\\ There's a hard-coded limit on the length of periodic powers which are
\\ recognised for the print. For a given gf denominator, the period is
\\ found by seeking a power p and constant c for which den is a factor of
\\ (1-c*x^p). c becomes the base in base^floor(n/spread).
\\
\\ The best way to choose bases to show for powers is not settled. With
\\ something like (1+I)^n there's a choice of that or [periodic]*2^floor(n/2).
\\ Similarly 12th roots, but maybe 2,3,4,6,8,12 are then only exact nth
\\ roots to encounter. Maybe two power forms should be shown, one with
\\ bases corresponding to the gf partial fractions, one with nth roots of
\\ unity spun up to positive reals (and combined when same magnitudes), if
\\ that's different.
\\
\\ The power expression always starts n=0 but it could be helpful to set the
\\ start to say n=1, or bigger if values are from the middle of a sequence.
\\ That should be no harder than a gf shift and a factor on the power
\\ coefficients.
\\
\\ History:
\\
\\ Version 1 - the first version.
\\ Version 2 - fix last value missed from confirm, allow t_POL values.
\\ Version 3 - fix for printing powers in gf denominator.
\\ Version 4 - fix for initial exceptions in powers self-test.
\\ Version 5 - fix and improve some periodic powers.
\\ Version 6 - more merging of periodic constant terms,
\\ fix for self-test of some powers.
\\ Version 7 - no factorize of multi-variable characteristic and gf.
\\
\\ Home Page:
\\
\\ http://user42.tuxfamily.org/pari-recurrence-guess/index.html
\\
\\-----------------------------------------------------------------------------
\\ configuration
recurrence_guess_minimum_confirm = 1;
{
addhelp(recurrence_guess_minimum_confirm,
"recurrence_guess_minimum_confirm = integer
This is the minimum number of confirming values which recurrence_guess() etc
must see before accepting a recurrence.
An order m recurrence is determined by 2*m values.
recurrence_guess_minimum_confirm=1 then demands vector v has at least 1 more
value satisfying the recurrence. If this was not so than every even length
v would produce an order m recurrence, and with no confidence at all that it
might be right.")};
recurrence_guess_variable = 'n;
{
addhelp(recurrence_guess_variable,
"recurrence_guess_variable = 'n
The name of the index variable to show in coefficients and power expressions
printed. It should be a symbol like 'n. This is presently only used for
printouts.
(Generating functions always show x as the formal variable, and give that in
polynomials returned by the various ...to_gf().)")};
\\-----------------------------------------------------------------------------
\\ generic helpers
\\ g is a polynomial generating function ("t_RFRAC" usually).
\\ Return a vector of its first n terms, starting at coefficient of x^0.
\\ The formal variable is the given x, or default variable(g).
recurrence_guess_INTERNAL_gf_terms(g,n,var=variable(g)) =
{
\\ print("recurrence_guess_INTERNAL_gf_terms ",g," "n);
if(g==0,return(vector(n)));
my(zeros = min(n,valuation(g,var)),
v = Vec(g + O(var^n)));
if(zeros>=0, concat(vector(zeros,i,0), v),
v[-zeros+1 .. #v]);
}
\\ g is a polynomial generating function ("t_RFRAC" usually).
\\ Return the coefficient of its x^n term.
\\ The formal variable is the given x, or default variable(g).
recurrence_guess_INTERNAL_gf_term(g,n,var=variable(g)) = \
recurrence_guess_INTERNAL_gf_terms(g,n+1,var)[n+1];
\\ f is a vector of recurrence factors.
\\ v is a vector of values.
\\ Return 1 if the values all satisfy the recurrence.
recurrence_guess_INTERNAL_check(f,v) =
{
\\ pos+#f = #v last check, so last pos = #v-#f
for(pos=1,#v-#f,
my(n = sum(i=0,#f-1, v[pos+i]*f[i+1]));
if(n!=v[pos+#f], return(0)));
1;
}
\\ return a lex() style compare of a <=> b,
\\ with comparison first by abs() magnitude then arg() angle
recurrence_guess_INTERNAL_lex_by_polar(a,b) =
{
lex(apply(recurrence_guess_INTERNAL_lex_by_polar_key,a),
apply(recurrence_guess_INTERNAL_lex_by_polar_key,b));
}
recurrence_guess_INTERNAL_lex_by_polar_key(z) =
{
my(a=if(z,arg(z))); if(a<0,a+=2*Pi); \\ 0 <= a < 2*Pi
[norm(z), a, real(z), imag(z)];
}
\\ p is a polynomial in 'x.
\\ If it is a divisor of some 1 - c*x^e then return that 1 - c*x^e.
\\ If not then return 0.
\\ Positive real c is preferred, ready to be the base shown in powers.
recurrence_guess_INTERNAL_pol_periodic(p) =
{
\\ In the loop have 1+c*x^e = p*q where q is a quotient made from
\\ successive terms t/low. Those terms are chosen to make 0s except for
\\ the low 1. poldegree(c) <= poldegree(p) so memory use is capped, until
\\ the return value.
my(low=polcoeff(p,0,'x),
c=1-p/low); \\ so lowest term of c is 1, then subtract that 1
for(e=1,100,
c /= 'x;
my(t=polcoeff(c,0));
\\ print(" periodic at e="e" c="c" t="t);
if(poldegree(c,'x)<1, \\ c is a constant, good
\\ try to avoid imaginary c by squaring or cubing
for(i=2,3,
if(imag(t),
my(tp=t^i);
if(! imag(tp), t=tp; e*=i)));
\\ try for positive c by squaring
\\ have 1-(-c)*x^e, multiply 1+(-c)*x^e to be 1-(c^2)*x^(2*e)
if(recurrence_guess_INTERNAL_is_negative(t),
my(tp=t^2);
if(! recurrence_guess_INTERNAL_is_negative(tp),
t=tp; e*=2));
\\ print(" periodic final t="t" e="e);
return(1 - t*'x^e));
\\ print(" sub ", t/low * p);
c -= t/low * p);
0; \\ not found
}
\\ Return a generating function for
\\ base^floor(n/spread)
\\ * n^expon
\\ * vector_modulo(periodic_vector,n)
recurrence_guess_INTERNAL_besp_to_gf(base,expon,spread,periodic_vector) =
{
recurrence_guess_values_to_gf
(vector((2*spread*#periodic_vector*(expon+1)+1 )*4+10,
n,n--; \\ 0 upwards
base^floor(n/spread)
* n^expon
* recurrence_guess_INTERNAL_vector_modulo(periodic_vector,n)));
}
\\ x is a number or a polynomial.
\\ If it needs parens around for printing an exponent ^123 after then
\\ return a stringized "(x)". If not then return x unchanged.
recurrence_guess_INTERNAL_parens_for_power(x) =
if(recurrence_guess_INTERNAL_want_parens_for_power(x), Str("(",x,")"), x);
recurrence_guess_INTERNAL_want_parens_for_power(x) =
{
x=simplify(x); \\ reduce poly constant only to a number
if(type(x)=="t_POL",
x!='x, \\ poly term x alone no need parens
!(real(x)==0 && imag(x)==1) \\ number I or w no need parens
&& !(real(x)>=0 \\ real >=0 no need parens
&& imag(x)==0
&& denominator(real(x))==1)); \\ but rational 9/2 needs parens
}
\\ x is a number.
\\ If it needs parens around for printing with a multiply *123 after then
\\ return a stringized "(x)". If not then return x unchanged.
recurrence_guess_INTERNAL_parens_for_multiply(x) = \
if(recurrence_guess_INTERNAL_want_parens_for_multiply(x), Str("(",x,")"), x);
recurrence_guess_INTERNAL_want_parens_for_multiply(x) =
{
x=simplify(x); \\ reduce poly constant only to a number
if(type(x)=="t_POL",
my(count=0);
for(i=0,poldegree(x),
my(c=polcoeff(x,i));
count+=(c!=0);
if(count>1 || recurrence_guess_INTERNAL_want_parens_for_multiply(c),
return(1))),
real(x) && imag(x));
}
\\ v is a vector of periodic factors starting from i==0 mod #v.
\\ Return the factor for i.
recurrence_guess_INTERNAL_vector_modulo(v,i) = v[(i%#v)+1];
\\ v1 and v2 are vectors of periodic factors.
\\ Return a vector which is their sum.
\\ If the periods are different then they are replicated as necessary to
\\ their least common multiple.
recurrence_guess_INTERNAL_periodic_vector_binop(func,v1,v2) =
{
vector(lcm(#v1,#v2), i,i--;
func(recurrence_guess_INTERNAL_vector_modulo(v1,i),
recurrence_guess_INTERNAL_vector_modulo(v2,i)));
}
\\ v is a vector of periodic factors like [5,7,-3].
\\ If the factors repeat themselves like [5,7,-3, 5,7,-3] then reduce.
recurrence_guess_INTERNAL_periodic_vector_reduce(v) =
{
my(ds=Vecrev(divisors(#v))[^1]);
for(len=1,floor(#v/2),
if(#v % len, next());
for(i=len+1,#v,
if(v[i] != v[i-len],
next(2)));
return(v[1..len]));
v;
}
\\ p is a polynomial. Return the number of non-zero terms.
recurrence_guess_INTERNAL_pol_num_terms(p,var=variable(Pol(p))) = \
sum(i=0,poldegree(p,var), polcoeff(p,i,var)!=0);
\\ p is a polynomial. Return its lowest exponent non-zero term.
recurrence_guess_INTERNAL_pol_lowest_coeff(p,var=variable(Pol(p))) = \
for(i=0,poldegree(p,var), if(polcoeff(p,i,var), return(polcoeff(p,i,var)))); 0;
\\ return true if x<0
recurrence_guess_INTERNAL_is_negative(x) =
{
while(type(x)=="t_POL",
x=simplify(x);
x=pollead(x));
real(x)<0 || (real(x)==0 && imag(x)<0);
}
\\-----------------------------------------------------------------------------
\\ p is a polynomial or t_RFRAC.
\\ Return a string of it with powers ascending.
\\ if want_plus is true then ensure the string starts with a + or - sign.
recurrence_guess_INTERNAL_ascending_str(p,want_plus, \
var=variable(Pol(numerator(p)))) =
{
my(debug=0);
if(debug,print("recurrence_guess_INTERNAL_ascending_str()\n var="var" p="p));
if(p==0,return("0"));
my(den = denominator(p));
if(poldegree(den,var) != 0,
\\ t_RFRAC like (x+2)/(x^2+5)
p=numerator(p);
if(debug,print("den ",den," num ",p));
\\ notice a power in the denominator
my(den_expon=1, den_factors=iferr(factor(den),e,'cannot_factorize));
if(den_factors!='cannot_factorize
&& matsize(den_factors)[1]==1
&& den_factors[1,2] > 1,
\\ factor() is not actually equal to den, only up to a constant
\\ factor so put that into the numerator
den_expon = den_factors[1,2];
p *= den_factors[1,1]^den_expon / den;
den=den_factors[1,1];
if(debug,print(" to power den ",den," expon ",den_expon)));
\\ divide out so constant term of den is 1
my(den_low=pollead(polrecip(den),var));
den /= den_low;
p /= den_low^den_expon;
if(debug,print(" to constant1 den ",den," num ",p,
" den_low was ",den_low));
p=simplify(p);
my(p_num_terms = recurrence_guess_INTERNAL_pol_num_terms(p,var),
p_low = recurrence_guess_INTERNAL_pol_lowest_coeff(p,var));
while(type(p_low)=="t_POL",p_low=simplify(pollead(p_low)));
my(p_parens = (p_num_terms != 1
|| (real(p_low) != 0 && imag(p_low) != 0)),
p_negate = want_plus && ! p_parens && (real(p_low)<0 || imag(p_low)<0),
den_num_terms = recurrence_guess_INTERNAL_pol_num_terms(den),
den_low = pollead(polrecip(den),var),
den_parens = (den_num_terms > 1
|| (den_expon!=1 && den_low != 1)
|| (real(den_low) != 0 && imag(den_low) != 0)));
if(p_negate, p=-p);
if(debug,print(" p_parens="p_parens" str "Str(p)" type "type(p)));
return(concat([if(want_plus,if(p_negate,"- ","+ "),""),
if(p_parens,"(",""),
recurrence_guess_INTERNAL_ascending_str(p,0,var),
if(p_parens,")",""),
"/",
if(den_parens,"(",""),
recurrence_guess_INTERNAL_ascending_str(den,0,var),
if(den_parens,")",""),
if(den_expon!=1,Str("^",den_expon),"")
])));
my(str="");
my(join=0);
for(i=0,poldegree(p,var),
my(coeff=polcoeff(p,i,var));
if(coeff==0, next());
if(join, str=concat(str," "));
if(join || want_plus,
if(recurrence_guess_INTERNAL_is_negative(coeff),
str=concat(str,"- "); coeff = -coeff,
str=concat(str,"+ ")));
str = concat([str,coeff*variable(p)^i]);
join=1);
if(debug,print(" ret str "str));
str;
}
\\ f is recurrence factors which give v is a vector of values.
\\ g is a t_RFRAC generating function in 'x.
\\ gparts is a vector of g as partial fractions (so vecsum(gparts)==g).
\\ Print powers like 2^n etc for the recurrence (using gparts).
recurrence_guess_INTERNAL_show_powers(f,v,g,gparts) =
{
my(debug=0);
if(gparts=='cannot_factorize,return());
print(" as powers");
my(powers_list = List([]),
initial_exceptions=[],
glist = List([]));
for(i=1,#gparts,
my(gpart = gparts[i], \\ t_RFRAC generating function part
num=numerator(gpart), \\ t_POL
den=denominator(gpart));
if(debug, print("gpart "gpart"\n num "num"\n den "den));
\\ If den=1 then gpart is non t_RFRAC initial exceptions.
if(den==1,
initial_exceptions=v[1 .. poldegree(gpart,'x)+1];
print(" initial exceptions ",initial_exceptions);
next());
\\ If den=power like (1-2*x)^3 then reduce to den=1-2*x and set pow=3.
\\ factor() can be out by a constant factor when complex coefficients,
\\ so adjust with num = num*newden^pow/den.
my(pow=1, den_factors=factor(den));
if(matsize(den_factors)[1]==1,
pow = den_factors[1,2];
num *= den_factors[1,1]^pow / den;
den = den_factors[1,1]);
if(debug, print(" pow "pow"\n num "num"\n den "den));
\\ periodic_pol is like 1-c*x^p
my(periodic_pol=recurrence_guess_INTERNAL_pol_periodic(den));
if(periodic_pol==0,
\\ this gpart not an exact power or periodic factors of an exact
\\ power, print as the generating function
listput(glist,gpart);
next());
if(debug, print(" periodic_pol ",periodic_pol));
my(spread=poldegree(periodic_pol,'x)); \\ p high power in 1-c*x^p
if(debug, print(" spread "spread));
if(debug, print(" quotient "periodic_pol/den));
num *= (periodic_pol/den)^pow;
if(debug, print(" num with qs "num" degree "poldegree(num,'x)));
my(base=-pollead(periodic_pol,'x)); \\ c factor in 1-c*x^p
if(debug, print(" base "base));
\\ pow_poly = (n+1)*(n+2)*(n+3)*(n+4)/24 for pow=5, etc pow-1 terms
my(pow_poly=prod(i=1,pow-1, 'x/spread + i) / (pow-1)!);
if(debug, print(" pow_poly ",pow_poly," factorial ",(pow-1)!));
\\ 1/(1-c*x^p) is pow_poly at [1,0,0,...] every "spread" many terms
\\ shift up for each term of num
\\ eg. (n/3+1)*(n/3+2)*2^floor(n/3)
\\ periodic_vectors[i] is periodic vector for term n^i, formed by
\\ picking terms out of pow_poly and num factor
my(periodic_vectors=vector(poldegree(pow_poly,'x)+1, i,
vector(spread))); \\ zeros to start
for(i=0,poldegree(num,'x),
my(p= polcoeff(num,i)*subst(pow_poly,'x,'x-i) / base^floor(i/spread));
if(debug, print(" i="i" p="p" degree "poldegree(p,'x)));
for(e=0,poldegree(p,'x),
periodic_vectors[e+1][(i%spread)+1] += polcoeff(p,e,'x)));
if(debug, print(" periodic_vectors "periodic_vectors));
\\ reduce base and spread if possible
\\ power is base^floor(n/spread)
\\ try to take spread-th root of base to reduce
\\ eg. 4^floor(n/2) becomes 2^n * vector_modulo([1,1/2],n)
my(spread_divisors=divisors(spread));
forstep(s=#spread_divisors,2, -1, \\ highest to lowest divisor
my(sd=spread_divisors[s],
reduced_base);
if( \\ ispower(I,2,&b) likes to say I is square of b=0.707+0.707*I,
\\ but want to stay exact here not go to floats.
\\ Also, ispower() doesn't take quads, so protect with iferr.
iferr(ispower(base,sd,&reduced_base),e,0)
&& type(real(reduced_base))!="t_REAL"
&& type(imag(reduced_base))!="t_REAL",
my(reduced_spread = spread/sd,
\\ periodic vector adjustment
\\ base^floor(i/spread) /reduced_base^floor(i/reduced_spread)
v=vector(spread,i,i--;
1/reduced_base^floor(i/reduced_spread)));
if(debug, print(" reduce sd="sd" base "base" -> "reduced_base
" spread "spread" -> "reduced_spread
" mul v="v));
base = reduced_base;
spread = reduced_spread;
for(i=1,#periodic_vectors,
periodic_vectors[i]
= recurrence_guess_INTERNAL_periodic_vector_binop
((x,y)->x*y, periodic_vectors[i], v));
break()));
for(expon=0,#periodic_vectors-1,
listput(powers_list, [base,
expon,
spread,
periodic_vectors[expon+1]]));
);
if(debug, print(" powers_list ",powers_list));
\\ sort powers_list descending base^(1/spread), and descending expon within
powers_list = Vec(powers_list);
my(spread_lcm=lcm(apply(v->v[3],powers_list)));
my(besp_to_key=(v->[recurrence_guess_INTERNAL_lex_by_polar_key
(v[1]^(spread_lcm/v[3])), \\ base^(1/spread)
v[2] ])); \\ expon
powers_list = vecsort(powers_list,
(x,y)->lex(besp_to_key(x),besp_to_key(y)),
4); \\ descending
if(debug, print(" powers_list sort ",powers_list));
\\ merge powers_list entries of same base,expon,spread
powers_list = List(powers_list);
forstep(i=#powers_list,2, -1,
if(powers_list[i-1][1] == powers_list[i][1] \\ base \
&& powers_list[i-1][2] == powers_list[i][2] \\ expon | all same
&& (powers_list[i-1][3] == powers_list[i][3] \\ spread /
|| (powers_list[i][1]==1 && powers_list[i][2]==0)), \\ or 1^n*n^0
powers_list[i-1][4] \\ merge into previous
= recurrence_guess_INTERNAL_periodic_vector_binop
((x,y)->x+y, powers_list[i-1][4], powers_list[i][4]);
listpop(powers_list,i)));
if(debug, print("powers_list merged ",powers_list));
\\ shorten periodic vectors, and delete possible zeros
forstep(i=#powers_list,1, -1,
my(v=recurrence_guess_INTERNAL_periodic_vector_reduce(powers_list[i][4]));
if(v==[0],
listpop(powers_list,i),
powers_list[i][4]=v));
\\ print
my(func_vec = vector(#powers_list));
my(plus=" ");
for(i=1,#powers_list,
my(base = powers_list[i][1],
expon = powers_list[i][2],
spread = powers_list[i][3],
periodic = powers_list[i][4]);
print1(plus); plus = " + ";
\\ periodic part [5,7]
my(periodic_want_parens=0);
my(periodic_str=
if(\\ periodic terms
#periodic>1,
Str("vector_modulo(",periodic,","recurrence_guess_variable")"),
\\ single term == 1, skip
periodic[1]==1, "",
\\ single term != 1
periodic_want_parens=1; Str(periodic[1])));
\\ polynomial part n^2
my(poly_str="");
if(expon>0,
poly_str=Str(recurrence_guess_variable);
if(expon>1, poly_str=Str(poly_str,"^",expon)));
\\ base part 2^n
my(power_str="");
if(base != 1,
power_str=Str(recurrence_guess_INTERNAL_parens_for_power(base),"^");
if(spread==1,
power_str=Str(power_str,recurrence_guess_variable),
power_str=Str(power_str,
"floor(",recurrence_guess_variable,"/",spread,")")));
if(periodic_want_parens && (poly_str!="" || power_str!=""),
periodic_str=recurrence_guess_INTERNAL_parens_for_multiply(periodic[1]));
my(times="");
if(periodic_str!="", print1(periodic_str); times = " * ");
if(poly_str!="", print1(times,poly_str); times = " * ");
if(power_str!="", print1(times,power_str); times = " * ");
if(times=="",print1("1")); \\ if everything else was collapsed
print();
func_vec[i] = (k)->k^expon
* recurrence_guess_INTERNAL_vector_modulo(periodic,k)
* base^floor(k/spread);
);
apply(gpart->print(" + gf ",
recurrence_guess_INTERNAL_ascending_str(gpart,0,'x)),
glist);
glist = vecsum(Vec(glist));
if(debug, print("glist ",glist));
if(debug, print("initial_exceptions ",initial_exceptions));
my(func = (n)->\\ if(debug, print(" func(n="n")"));
if(n<#initial_exceptions, initial_exceptions[n+1],
recurrence_guess_INTERNAL_gf_term(glist,n)
+ sum(i=1,#func_vec, func_vec[i](n))));
my(t=vector(#v,n,n--; func(n)));
if(t != v,
error("oops, powers different from values\ngot ",t,"\nwant ",v));
for(i=1,#powers_list,
my(g=recurrence_guess_INTERNAL_besp_to_gf(powers_list[i][1],
powers_list[i][2],
powers_list[i][3],
powers_list[i][4]));
if(debug, print(" i="i" gf "g));
glist+=g);
if(debug, print("glist+besp ",glist));
glist += Polrev(initial_exceptions
- recurrence_guess_INTERNAL_gf_terms(glist,#initial_exceptions));
if(debug, print("glist+exceptions ",glist));
if(g!=glist,
error("oops, powers reconstruct generating function different\n",
"got ",glist,"\nwant ",g));
}
\\ Maybe this could have general use pretty printing partial fractions.
\\ Or maybe it would belong in pol-pfrac.gp.
recurrence_guess_INTERNAL_print_gf(g,gparts=0) =
{
my(debug=0);
print(" ",recurrence_guess_INTERNAL_ascending_str(g,0,'x));
if(gparts=='cannot_factorize,return());
if(pol_partial_fractions=='pol_partial_fractions, read("pol-pfrac.gp"));
if(!gparts, gparts=pol_partial_fractions(g));
if(debug,print("gparts="gparts));
if(g,
my(dens=List([]));
for(i=1,#gparts,
my(d=denominator(gparts[i]));
if(polcoeff(d,0), d/=polcoeff(d,0)); \\ low coeff 1
if(debug,print("d="d));
if(d==1,next());
for(j=1,#dens,
if(dens[j]%d==0, dens[j]=d; next(2)); \\ smaller power replace
if(d%dens[j]==0, next(2))); \\ bigger power
listput(dens,d)); \\ new
if(debug,print("dens="dens" "dens_exp));
my(num=g);
my(dens_exp=vector(#dens));
for(i=1,#dens,
my(d=dens[i]);
while(denominator(num)%d==0, dens_exp[i]++;num*=d));
\\ iferr() to guard against cannot factorize multi-variable poly
my(nparts=iferr(pol_partial_fractions(1/num), e,[]),
nums=List([]));
if(#nparts,
for(i=1,#nparts,
my(d=denominator(nparts[i]));
if(polcoeff(d,0), d/=polcoeff(d,0)); \\ low coeff 1
for(j=1,#nums,
if(nums[j]%d==0, nums[j]=d; next(2)); \\ smaller power replace
if(d%nums[j]==0, next(2))); \\ bigger power
listput(nums,d)); \\ new
my(nums_exp=List(vector(#nums)));
for(i=1,#nums,
my(d=nums[i]);
while(num%d==0, nums_exp[i]++;num/=d; if(d==1,break())));
if(num!=1,listinsert(nums,num,1);listinsert(nums_exp,1,1));
if(debug,print("nums="nums" "nums_exp));
if(#nums>1 || #dens>1,
print1(" = ");
my(join="");
for(i=1,#nums,
my(p=nums[i],
want_parens=if(nums_exp[i]>1,
recurrence_guess_INTERNAL_want_parens_for_power(p),
recurrence_guess_INTERNAL_want_parens_for_multiply(p)));
print1(join,
if(want_parens,"(",""),
recurrence_guess_INTERNAL_ascending_str(p,0,'x),
if(want_parens,")",""),
if(nums_exp[i]>1,Str("^",nums_exp[i]),""));
join=" * ");
print1(" / ");
join="";
if(#dens>1,print1("( "));
for(i=1,#dens,
my(p=dens[i],
want_parens=if(dens_exp[i]>1,
recurrence_guess_INTERNAL_want_parens_for_power(p),
recurrence_guess_INTERNAL_want_parens_for_multiply(p)));
print1(join,
if(want_parens,"(",""),
recurrence_guess_INTERNAL_ascending_str(p,0,'x),
if(want_parens,")",""),
if(dens_exp[i]>1,Str("^",dens_exp[i]),""));
join=" * ");
if(#dens>1,print1(" )"));
print(""))));
if(#gparts==1, print(" same in partial fractions"),
print(" = partial fractions");
for(i=1,#gparts,
my(p=gparts[i]);
print(" ",recurrence_guess_INTERNAL_ascending_str(p,i>1,'x))));
}
\\-----------------------------------------------------------------------------
\\ C for confirm so V-C for recurrence is 2*L terms,
\\ so L=floor((V-C)/2) biggest to try
recurrence_guess_values_to_factors(v) =
{
for(len=1, floor((#v-recurrence_guess_minimum_confirm)/2),
my(m = matrix(len,len,i,j, v[i+j-1]),
a = vectorv(len,i, v[i+len]));
my(f = iferr(matsolve(m,a),
e,next(), errname(e)=="e_INV"));
f = simplify(Vec(f)); \\ row vector, simplify down any t_POL terms
if(recurrence_guess_INTERNAL_check(f,v),
return(f)));
[];
}
{
addhelp(recurrence_guess_values_to_factors,
"f = recurrence_guess_values_to_factors(v)
v is a vector of values.
Return a vector f of linear recurrence coefficients for the values.
Or if a recurrence cannot be found then return an empty vector [].
The coefficients are in ascending order the same order as the values v.
So for example if an order 3 recurrence is found then
f[1]*v[1] + f[2]*v[2] + f[3]*v[3] == v[4]
")};
recurrence_guess_factors_and_values_to_gf(f,v) =
{
my(f = recurrence_guess_values_to_factors(v));
if(f == [], return(0));
my(f_pol = Pol(f),
low_pol = Polrev(v[1 .. #f]),
low_product_mod = (f_pol*low_pol) % ('x)^(#f - 1),
g = ( 'x*low_product_mod - low_pol )/('x*f_pol-1));
my(t=recurrence_guess_INTERNAL_gf_terms(g,#v,'x));
if(t != v,
error("oops, gf different from v\ngot ",t,"\nwant ",v));
g;
}
{
addhelp(recurrence_guess_factors_and_values_to_gf,
"gf = recurrence_guess_factors_and_values_to_gf(f,v)
v is a vector of values and f a vector of linear recurrence coefficients.
Return a polynomial generating function (\"t_RFRAC\") which is the linear
recurrence with initial values from v.
Must have length #v >= #f. As a self-test it is confirmed that the gf
returned does in fact give all of v, both the initial values and any extras
given (error() if not).")};
recurrence_guess_values_to_gf(v) =
{
my(f = recurrence_guess_values_to_factors(v));
if(f == [], return(0));
recurrence_guess_factors_and_values_to_gf(f,v);
}
{
addhelp(recurrence_guess_values_to_gf,
"gf = recurrence_guess_values_to_gf(v)
v is a vector of values.
Return a polynomial generating function (\"t_RFRAC\") for a linear recurrence
satisfying the values. Or return an empty vector [] if no recurrence can be
found.")};
recurrence_guess(v) =
{
my(f = recurrence_guess_values_to_factors(v));
if(f == [],
print("No linear recurrence found\n");
return());
print("Recurrence length=", #f);
print(" coefficients");
print1(" v["recurrence_guess_variable"-",#f"]* ");
if(#f==1,
print1(recurrence_guess_INTERNAL_parens_for_multiply(f[1])),
print1(f," *v["recurrence_guess_variable"-1] "));
print(" = v[",recurrence_guess_variable,"]");
if(#f>=2 && f[1]==0,
print(" leading 0 coeffs are initial exceptions before recurrence"));
print1(" v[",recurrence_guess_variable,"] = v["recurrence_guess_variable"-1]* ");
if(#f==1,
print(recurrence_guess_INTERNAL_parens_for_multiply(f[1])),
print(Vecrev(f)," *v["recurrence_guess_variable"-",#f"]"));
print();
print(" characteristic polynomial ");
my(p = Polrev(concat(-f,[1])));
print(" ",p);
my(pfactors = iferr(factor(p), e,
\\ pari 2.9 cannot factorize multi-variable polynomials
print(" (cannot factorize: ",errname(e),")");
'cannot_factorize));
if(pfactors!='cannot_factorize,
my(pfactor_strs = vector(#pfactors[,1],i,
if(pfactors[i,2]==1,
Str(pfactors[i,1]),
Str("(",pfactors[i,1],")^",pfactors[i,2]))),
width = vecmax(apply(length, pfactor_strs)));
if(#pfactors[,1] > 1 || pfactors[1,2] > 1,
print(" = factors");
for(i=1,#pfactors[,1],
my(roots=polroots(pfactors[i,1]*1.0));
printf(" %-*s %s ", width, pfactor_strs[i],
if(#roots==1,"root","roots"));
roots = vecsort(roots, recurrence_guess_INTERNAL_lex_by_polar);
for(j=1,#roots,
if(j>1,print1(", "));
if(imag(roots[j])==0, printf("%.6g",real(roots[j])),
real(roots[j])==0, printf("I*%.6g",imag(roots[j])),
printf(" %.6g+I*%.6g",real(roots[j]),imag(roots[j]))));
print()));
);
print();
print(" generating function");
my(g = recurrence_guess_factors_and_values_to_gf(f,v));
my(t=recurrence_guess_INTERNAL_gf_terms(g,#v));
if(t!=v, error("oops, gf different from values\ngot ",t,"\nwant ",v));
if(pol_partial_fractions=='pol_partial_fractions, read("pol-pfrac.gp"));
my(gparts = iferr(pol_partial_fractions(g), \\ over vecsum(0*v) maybe ?
e, 'cannot_factorize));
recurrence_guess_INTERNAL_print_gf(g,gparts);
print();
recurrence_guess_INTERNAL_show_powers(f,v,g,gparts);
print(); \\ undef return value so no result print interactively
}
\\-----------------------------------------------------------------------------