# This script is based on a Maple session by Bruno Salvy (March 4,
# 2003), later updated and simplified: no didactic information here,
# only what is used in the Coq proof.

# Convention: This script has one section for each sequence to be
# considered, which start by '### Sequence x_{...} = <definition>'.
# When this sequence is explicit or is obtained by a closure
# computation of type plus or times, the section stores the
# corresponding operators into a file 'ann_x.v'; when it is obtained
# by creative telescoping, say x is the sum of y over some domain, the
# pairs (P, Q) are stored into the generated file for y, 'ann_y.v',
# before renormalizing the Ps and storing them into the file for x,
# 'ann_x.v'.

# Summary of the file (listing the closures):
#   * c_{n,k} = explicit hypergeometric term
#   * a_n = \sum_{k=0}^n c_{n,k}
#   * d_{n,k,m} = explicit hypergeometric term
#   * s_{n,k} = \sum_{m=0}^k d_{n,m}
#   * z_n = explicit special sequence
#   * u_{n,k} = s_{n,k} + z_n
#   * v_{n,k} = u_{n,k} * c_{n,k}
#   * b_n = \sum_{k=0}^n v_{n,k}


# We use Algolib, which is distributed in the git repository,
libname := ".", libname:
# as well as converters from Maple to SSR notation.
read "maple2ssr.mpl":

# Help functions.
rec_order := proc(rec, u, x, $)
  max(
    map(p -> op(p) - x,
    map2(select, has, indets(rec, specfunc(anything, u)), x))
  )
end proc:
to_skew_poly := proc(rec, u, xs, Sxs, $)
  local T := table();
  local fcalls := indets(rec, specfunc(anything, u));
  local i, shifted_u, x;
  for i to nops(xs) do T[xs[i]] := Sxs[i] end do;
  fcalls := map(
    shifted_u ->
      shifted_u = mul(T[x]^(op(select(has, shifted_u, x)) - x), x = xs),
    fcalls);
  eval(rec, fcalls)
end proc:


### Define sequence to prepare evaluations.
c_nk := binomial(n,k)^2 * binomial(n+k,k)^2:
a_n = Sum(eval(c_nk, k = l), l = 0..n):
d_nkm := (-1)^(m+1) / (2 * m^3 * binomial(n,m) * binomial(n+m,m)):
s_nk := Sum(d_nkm, m = 1..k):
z_n := Sum(1/l^3, l = 1..n):
u_nk := s_nk + z_n:
v_nk := u_nk * c_nk:
b_n := Sum(eval(v_nk, k = l), l = 0..n):


### Sequence c_{n,k} = explicit hypergeometric term.

extra1 := [n], [Sn], c:
extra2 := [n, k], [Sn, Sk], c:
extraA := extra2, mode = tdeg(Sn, Sk):
extraP := extra1, mode = "P":
extraQ := extra2, mode = "Q":

maple2ssr:-pp_skew_poly(
  Sn - factor(normal(eval(c_nk, n = n + 1) / c_nk, expanded)),
  extraA
):
maple2ssr:-write(%, append = false, file = "ann_c.v");

maple2ssr:-pp_skew_poly(
  Sk - factor(normal(eval(c_nk, k = k + 1) / c_nk, expanded)),
  extraA
):
maple2ssr:-write(%, file = "ann_c.v");


### Sequence a_n = \sum_{k=0}^n c_{n,k}.

# This summation is over natural boundaries, as shown by inhom = 0 below.
CT := Mgfun:-creative_telescoping(c_nk, n::shift, k::shift)[1]:
ct := eval(CT, {_F = unapply(c_nk, n), _f = unapply(c_nk, n, k)}):
inhom := normal(limit(ct[2], k = n+3)) - eval(ct[2], k = 0);

maple2ssr:-pp_skew_poly(to_skew_poly(CT[1], _F, [n], [Sn]), extraP):
maple2ssr:-write(%, file = "ann_c.v");

maple2ssr:-pp_skew_poly(to_skew_poly(CT[2], _f, [n, k], [Sn, Sk]), extraQ):
maple2ssr:-write(%, skip_line = false, file = "ann_c.v");

old_module := "annotated_recs_c":
maple2ssr:-pp_in_terms_of([1], [n], a, old_module, ["P_flat"], "Sn2"):
maple2ssr:-write(%, append = false, skip_line = false, file = "ann_a.v");


### Sequence d_{n,k,m} = explicit hypergeometric term.

extra2 := [n, k], [Sn, Sk], d:
extra3 := [n, k, m], [Sn, Sk, Sm], d:
extraA := extra3, mode = tdeg(Sn, Sk, Sm):
extraP := proc(i, $) extra2, mode = "P" || i end proc:
extraQ := proc(i, $) extra3, mode = "Q" || i end proc:

maple2ssr:-pp_skew_poly(
  Sn - factor(normal(eval(d_nkm, n = n + 1) / d_nkm, expanded)),
  extraA
):
maple2ssr:-write(%, append = false, file = "ann_d.v");

maple2ssr:-pp_skew_poly(
  Sk - factor(normal(eval(d_nkm, k = k + 1) / d_nkm, expanded)),
  extraA
):
maple2ssr:-write(%, file = "ann_d.v");

maple2ssr:-pp_skew_poly(
  Sm - factor(normal(eval(d_nkm, m = m + 1) / d_nkm, expanded)),
  extraA
):
maple2ssr:-write(%, file = "ann_d.v");


### Sequence s_{n,k} = \sum_{m=0}^k d_{n,m}.

# This summation is not over natural boundaries as it is an indefinite one.
CT := Mgfun:-creative_telescoping(d_nkm, [n::shift, k::shift], m::shift):
# CT[1][1] expresses \Delta_n; CT[2][1] expresses \Delta_k.
# CT[2][2] is 0.

# Operators on sum from CT[1].
eval(CT[1][2], _f(n, k, m) = d_nkm):
eval(%, m = k + 1) - eval(%, m = 1):
inh := map(normal@expand, %, expanded):
convert(Mgfun:-dfinite_expr_to_sys(inh, f(n::shift, k::shift)), list):
# Each recurrence Ps_[i] is certified by Qs_[i].
Ps_ := collect(
        eval(%, f = unapply(d(n + 1, k, m) - d(n, k, m), n, k)),
        d, normal):
Qs_ := collect(
        eval(%%, f = unapply(eval(CT[1][2], _f = d), n, k)),
        d, normal):

# Operators on sum from CT[2].
inh := factor(normal(expand(eval(d_nkm, m = k + 1)), expanded)):
convert(Mgfun:-dfinite_expr_to_sys(inh, f(n::shift, k::shift)), list):
# Each recurrence Ps__[i] is certified by Qs__[i].
Ps__ := collect(
        eval(%, f = unapply(d(n, k + 1, m) - d(n, k, m), n, k)),
        d, normal):
Qs__ := collect(
        eval(%%, f = unapply(0, n, k)),
        d, normal):

Ps := map(op, [Ps_, Ps__]):
Qs := map(op, [Qs_, Qs__]):
N := nops(Ps):

for i to N do
  eval(Qs[i], d = unapply(d_nkm, n, k, m));
  eval(%, m = m + 1) - %;
  eval(Ps[i], d = unapply(d_nkm, n, k, m));
  print(i = normal(%% - %, expanded))
end do:

Ps := map(to_skew_poly, Ps, d, [n, k], [Sn, Sk]):
Qs := map(to_skew_poly, Qs, d, [n, k, m], [Sn, Sk, Sm]):

for i to N do
  maple2ssr:-pp_skew_poly(Ps[i], extraP(i));
  maple2ssr:-write(%, file = "ann_d.v");

  maple2ssr:-pp_skew_poly(Qs[i], extraQ(i));
  maple2ssr:-write(%, skip_line = evalb(i < N), file = "ann_d.v")
end do:

# Normalize the set of Ps for s_{n,k} by a Gröbner-basis computation.
Alg := Ore_algebra:-shift_algebra([Sn, n], [Sk, k]):
MOrd := Groebner:-MonomialOrder(Alg, tdeg(Sn, Sk)):
bPs := Groebner:-Basis(Ps, MOrd):
Groebner:-LeadingMonomial(bPs, MOrd);
bPs_names := map(maple2ssr:-pp_name_of_monomial, %, [Sn, Sk]);
Groebner:-LeadingMonomial(Ps, MOrd);
Ps_names := [seq("P" || i || "_flat", i = 1..nops(Ps))];

# Transformation rules: GB in terms of old generators.
rules := [
  gb[1] = gens[4],
  gb[2] = gens[3],
  gb[3] = (n+k+2)*gens[1] + k*(2*n+3)*(n-k)*gens[3]
]:

# Verifies them in Maple.
eval(rules, {seq(gb[i] = bPs[i], i=1..3), seq(gens[i] = Ps[i], i=1..4)}):
map(evalb@expand, %);

old_module := "annotated_recs_d":
for i to 3 do
  [seq(coeff(op(2, rules[i]), gens[j]), j=1..4)];
  maple2ssr:-pp_in_terms_of(%, [n, k], s, old_module, Ps_names, bPs_names[i]);
  maple2ssr:-write(%, append = evalb(1 < i), file = "ann_s.v")
end do:

extra2 := [n, k], [Sn, Sk], s:
extraA := extra2, mode = tdeg(Sn, Sk):
for i to 3 do
  maple2ssr:-pp_skew_poly(bPs[i], extraA);
  maple2ssr:-write(%, skip_line = evalb(i < 3), file = "ann_s.v")
end do:


### Sequence z_n = \sum_{k=1}^n 1/k^3.

extra1 := [n], [Sn], z:
extraP := extra1, mode = tdeg(Sn):

inhom := (n+1)^3 * (_f(n+1) - _f(n)) - 1:
collect(eval(inhom, n = n + 1) - inhom, _f, factor):
ann_z := to_skew_poly(%, _f, [n], [Sn]):
maple2ssr:-pp_skew_poly(ann_z, extraP):
maple2ssr:-write(%, append = false, skip_line = false, file = "ann_z.v");


### Sequence u_{n,k} = s_{n,k} + z_n.

# Find GB for s in terms of ann_z and Delta_k.  This will show that
# ann_s is a subideal of ann_z (viewed w.r.t. n and k).  Therefore,
# ann_u is taken to be equal to ann_s.

Groebner:-NormalForm(ann_z, bPs, MOrd);
map(Groebner:-NormalForm, bPs, [ann_z, Sk - 1], MOrd);

divide_by_ann_z_and_Delta_k := proc(L, $)
  local R := eval(L, Sk = 1);
  local Q2 := normal((L - R) / (Sk - 1));
  local Q1 := normal(R / ann_z);
  collect([Q1, Q2], {Sn, Sk}, distributed, factor);
end proc:
rules := [seq(divide_by_ann_z_and_Delta_k(bPs[i]), i = 1..3)]:

#old_module := "annotated_recs_d":
for i to 3 do
  #rules[i];
  # What follows breaks as pp_in_terms_of only accepts
  # rational-function coefficients.
#  maple2ssr:-pp_in_terms_of(%, [n, k], s, old_module, ["ann_z", "Delta_k"], bPs_names[i]);
#  maple2ssr:-write(%, append = evalb(1 < i), skip_line = evalb(i < 3), file = "ann_s.v");
# The following is to help write the proof manually.
  gb[i] = rules[i][1] * Ann_z + rules[i][2] * Delta_k
end do;

extra2 := [n, k], [Sn, Sk], w:

maple2ssr:-pp_skew_poly(rules[1][2], extra2, mode = "gb1_cf_of_Delta_k"):
maple2ssr:-write(%, append = false, file = "ann_u.v");

maple2ssr:-pp_skew_poly(rules[2][2], extra2, mode = "gb2_cf_of_Delta_k"):
maple2ssr:-write(%, file = "ann_u.v");

maple2ssr:-pp_skew_poly(rules[3][1], extra2, mode = "gb3_cf_of_Ann_z"):
maple2ssr:-write(%, file = "ann_u.v");

maple2ssr:-pp_skew_poly(rules[3][2], extra2, mode = "gb3_cf_of_Delta_k"):
maple2ssr:-write(%, file = "ann_u.v");

# Express that u_nk has the same annihilator as s_nk.
old_module := "annotated_recs_s":

maple2ssr:-pp_in_terms_of([1], [n, k], u, old_module, ["Sn2_lcomb"], "Sn2"):
maple2ssr:-write(%, file = "ann_u.v");

maple2ssr:-pp_in_terms_of([1], [n, k], u, old_module, ["SnSk_lcomb"], "SnSk"):
maple2ssr:-write(%, file = "ann_u.v");

maple2ssr:-pp_in_terms_of([1], [n, k], u, old_module, ["Sk2_lcomb"], "Sk2"):
maple2ssr:-write(%, file = "ann_u.v");


### Sequence v_{n,k} = u_{n,k} * c_{n,k}.

sys_for_u_nk := map(Ore_algebra:-applyopr, bPs, _f(n, k), Alg):
sys_for_c_nk := Mgfun:-dfinite_expr_to_sys(c_nk, _f(n::shift, k::shift)):
sys_for_v_nk := Mgfun:-`sys*sys`(sys_for_u_nk, sys_for_c_nk):

extra2 := [n, k], [Sn, Sk], v:
extraA := extra2, mode = tdeg(Sn, Sk):
map(to_skew_poly, sys_for_v_nk, _f, [n, k], [Sn, Sk]):
# Sort this Maple set for reproducibility.
sort(
  convert(%, list),
  (a, b) -> Groebner:-TestOrder(a, b, tdeg(Sn, Sk))
):
defs := map(maple2ssr:-pp_skew_poly, %, extraA):
for i to 3 do
  maple2ssr:-write(defs[i], append = evalb(1 < i), file = "ann_v.v")
end do:


### Sequence b_n = \sum_{k=0}^n v_{n,k}.

infolevel[Mgfun]:=5:
CT := Mgfun:-creative_telescoping(LFSol(sys_for_v_nk), n::shift, k::shift)[1]:
ct := eval(CT, {_F = unapply(v_nk, n), _f = unapply(v_nk, n, k)}):
# This summation is over natural boundaries, as shown by inhom = 0 below.
# Just, the computation is more involved than just calling limit.
inhom := eval(ct[2], k = n+5) - eval(ct[2], k = 0):
su_Sum := map(
  proc(s, $)
    local f := op(1, s);
    local lb := op([2, 2, 1], s);
    local ub := op([2, 2, 2], s);
    local pv := op(indets(ub));
    local sv := op([2, 1], s);
    local i;
    s = Sum(f, sv = lb..pv) + add(normal(expand(eval(f, sv = pv + i))), i = 1..ub - pv)
  end proc,
  select(has, indets(inhom, function), Sum)
):
su_binomial := map(
  proc(b, $)
    b = normal(expand(b))
  end proc,
  select(hastype, indets(inhom, function), specfunc(anything, binomial))
):
eval(inhom, su_Sum union su_binomial):
eval(%, binomial = proc(a, b, $) `if`(type(b, negint), 0, 'binomial'(a, b)) end proc):
inhom := collect(%, {Sum, binomial}, normal);

extra1 := [n], [Sn], v:
extra2 := [n, k], [Sn, Sk], v:
extraA := extra2, mode = tdeg(Sn, Sk):
extraP := extra1, mode = "P":
extraQ := extra2, mode = "Q":

maple2ssr:-pp_skew_poly(to_skew_poly(CT[1], _F, [n], [Sn]), extraP):
maple2ssr:-write(%, file = "ann_v.v");

maple2ssr:-pp_skew_poly(to_skew_poly(CT[2], _f, [n, k], [Sn, Sk]), extraQ):
maple2ssr:-write(%, skip_line = false, file = "ann_v.v");

old_module := "annotated_recs_v":
maple2ssr:-pp_in_terms_of([1], [n], b, old_module, ["P_flat"], "Sn4"):
maple2ssr:-write(%, append = false, skip_line = false, file = "ann_b.v");