######################################################################
## Hilbert10.txt Save this file as  Hilbert10.txt  to use it,        #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `Hilbert10.txt`: <Enter>                     #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Robert Dogherty-Bliss and Doron Zeilberger,            #
##Rutgers University ,                                               #
## Send bugs to  DoronZeil at gmail dot com                          # 
######################################################################


print(`First Written: Sept. 2022: tested for Maple 2020`):
print(`Version  : Sept. 2022  `):
print():
print(`This is Hilbert10.txt, A Maple package`):
print(`accompanying Robert Dogherty-Bliss nd Doron Zeilberger's  article: `):
print(` Experimenting with Matiyasevich's Seminal Proof of Hilbret's 10th problem `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Hilbert10.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` EMatrix`):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Hilbert10.txt: A Maple package for  experimenting with Matiyasevich's proof of Hilbert 10th`):

print(`The MAIN procedures are:`):
print(` Bmat, Dio, BackwardsProof, SeqFromRec, SolveD, SolveDd`):
print(`CheckYuriBad, findBadYuriDiag, BackwardsInequalities`):

elif nargs=1 and args[1]=Bmat then
print(`Bmat(b): The matrix associated with the C-finite recurrence with coefficients [b[1], b[2], ..., nops(b)]`):
print(`Try:`):
print(`Bmat([2, -1]);`):
print(`Bmat([a, b, 1]);`):


elif nargs=1 and args[1]=BackwardsProof then
print(`BackwardsProof(b): Returns the simplified difference`):
print(`     P(x[1], x[2], ..., x[nops(b)]) - P(x[0], x[1], x[2], ..., x[nops(b) - 1])`):
print(`where P is the polynomial Dio(b, x) and x[0] is the value obtained by`):
print(`undoing the recurrence with coefficients [b[1], b[2], ..., b[nops(b)]]`):
print(`on the values x[1], x[2], ..., x[nops(b)].`):
print(`In particular, if this is 0, then the diophantine equation is possibly`):
print(`uniquely solved by consecutive terms of the C-finite sequence`):
print(`[[0, 0, ..., 0, 1], [b[1], ..., b[nops(b)]]]`):
print(`Try:`):
print(`BackwardsProof([3, -1])`):
print(`BackwardsProof([3, 2, 1])`):
print(`BackwardsProof([a, b, 1])`):

elif nargs=1 and args[1]=Dio then
print(`Dio(b, x): The polynomial P(x[1],x[2],...x[nops(b)]) such that the solution`):
print(`a(n) to the recurrence with coefficients [b[1], b[2], ..., b[nops(b)]] and`):
print(`initial conditions [0, 0, ..., 0, 1] satisfies`):
print(`P(a(n), a(n + 1), ..., a(n + nops(b) - 1)) = det(B)^n,`):
print(`where B is the matrix associated with the recurrence.`):
print(`Try:`):
print(`Dio([3, -1], x)`):
print(`Dio([3, 2, 1], x)`):


elif nargs=1 and args[1]=DioV then
print(`DioV(b,k): Verbose form of Dio(b,k). Try:`):
print(`DioV(b,2,x);`):

elif nargs=1 and args=EMatrix then

print(`EMatrix(b, e, n): Produce Bmat(b)^n in terms of the first fundamental sequence e(n),`):
print(`i.e., the solution to the recurrence with coefficients b[1], b[2], ..., b[m]`):
print(`with initial conditions 0, 0, ..., 0, 1.`):
print(`Try:`):
print(`EMatrix([b, -1], e, n);`):



elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(S,N): Inputs S=[INI,ope]`):
print(`where INI is the list of initial conditions, a ope a list of`):
print(`size L, say, and a recurrence operator ope, codes a list of`):
print(`size L, finds the first N0 terms of the sequence satisfying`):
print(`the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).`):
print(`For example, for the first 20 Fibonacci numbers,`):
print(`try: SeqFromRec([[1,1],[1,1]],20);`):

elif nargs=1 and args[1]=SolveD then
print(`SolveD (p, vars, K)`):
print(` Input: Polynomial p, list of variables var, and a positive integer K`):
print(` Output: All positive integers K vars[1] <= vars[2] <= ... <= vars[m] <= K`):
print(`         such that p(vars[1], ..., vars[m]) = 0.`):
print(`Try: `):
print(`SolveD(x^2+y^2-z^2,[x,y,z],40);`):


elif nargs=1 and args[1]=SolveDd then
print(`SolveDd(p, vars, K)`):
print(` Input: Polynomial p, list of variables var, and a positive integer K`):
print(` Output: All positive integers K vars[1] < vars[2] < ... <vars[m] <= K`):
print(`and relatively prime to each other `):
print(`         such that p(vars[1], ..., vars[m]) = 0.`):
print(`Try: `):
print(`SolveDd(x^2+y^2-z^2,[x,y,z],40);`):


elif nargs=1 and args[1]=CheckYuriBad then
print(`CheckYuriBad(rec, K)`):
print(`Check if the solutions bounded by K to Dio(rec, x) - 1 are generated by more than`):
print(`one set of initial conditions (efficiently, without generating all`):
print(`possible solutions up to K.)`):
print(`Try:`):
print(`CheckYuriBad([1, 1, 1], 100);`):
print(`CheckYuriBad([4, 4, 1], 100);`):

elif nargs=1 and args[1]=findBadYuriDiag then
print(`findBadYuriDiag(K, M)`):
print(`Find the a <= K such that the solutions to Dio([a, a, 1], x) are probably generated`):
print(`by more than one positive initial condition, looking at solutions bounded by M.`):
print(`Try:`):
print(`findBadYuriDiag(10, 200);`):
print(`findBadYuriDiag(100, 200);`):

elif nargs=1 and args[1]=BackwardsInequalities then
print(`BackwardsInequalities(b, x)`):
print(`A set {H, C} of equalities and inequalities such that, if H implied C,`):
print(`then the going down formula for Dio(b, x) would lead to a single solution.`):
print(`Try:`):
print(`BackwardsInequalities([b, -1], x);`):
print(`BackwardsInequalities([1, 1, 1], x);`):


else
 print(`There is no such thing as`, args):

fi:


end:



with(linalg):


# Bmat(b,k): The associated matrix B for the C-finite recurrence [b[1], b[2], ..., b[k - 1], (-1)^k].
Bmat:=proc(b,k) local i:
    Matrix([b, seq([0 $ (i - 1), 1, 0 $ (nops(b) - i)], i=1..nops(b)-1)]):
end:

#DioV(b,k,x): Verbose form of Dio(b,k). Try:
#DioV(b,2,x);
DioV:=proc(b,k,x) local gu,a,n,i:
gu:=Dio(b,k,x):

print(`A Diophantine Equation satisfied by any`, k, `consecutive values of the solution of the recurrence`):
print(``):
print(a(n)=add(b[i]*a(n-i),i=1..k-1)+(-1)^(k-1)*a(n-k)):
print(``):
print(seq(a[i]=0,i=0..k-2), a[k-1]=1 ):
print(``):
print(`Prop: Let `, P(seq(x[i],i=0..k-1)), `be the following polynomial `):
print(``):
print(gu):
print(``):
print(`and in Maple notation: `):
print(``):
lprint(gu):
print(``):
print(P(seq(a(n+i),i=0..k-1))=0):
print(``):
end:


#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:

IncSeqs := proc(m, K)
    local S, K1, S1, L:

    if m = 0 then
        return {[]}:
    fi:

    S := {}:

    for K1 from 0 to K do
        S1 := IncSeqs(m - 1, K1):
        S := S union {seq([op(L), K1], L in S1)}:
    od:

    S:
end:

# Input: Polynomial p, list of variables var, and a positive integer K
# Output: All positive integers K vars[1] <= vars[2] <= ... <= vars[m] <= K
#         such that p(vars[1], ..., vars[m]) = 0.
SolveD:=proc(p, vars, K)
    local S, sl, L, i:

    S := IncSeqs(nops(vars), K):
    sl := {}:

    for L in S do
        if subs({seq(vars[i] = L[i], i=1..nops(vars))}, p) = 0 then
            sl := {op(sl), L}:
        fi:
    od:

    sl:
end:

# Dio(b, x): The polynomial P(x[1],x[2],...x[nops(b)]) such that the solution
# a(n) to the recurrence with coefficients [b[1], b[2], ..., b[nops(b)]] and
# initial conditions [0, 0, ..., 0, 1] satisfies
#       P(a(n), a(n + 1), ..., a(n + nops(b) - 1)) = det(B)^n,
# where B is the matrix associated with the recurrence.
Dio := proc(b, x)
    local e, M, expr, i:
    M := EMatrix(b, e, n):
    expr := det(M):

    subs({seq(e(n + i)=x[i + 1], i=0..nops(b))}, expr):
end:

# Return a matrix M which transforms the basis {e1(n), e1(n + 1), ..., e1(n + m
# - 1)} into the basis {e1(n), e2(n), ..., em(n)}, where the ei are the
# fundamental solutions to the recurrence with coefficients rec[1], rec[2],
# ..., rec[m], i.e., e1(n) is [[0, 0, ..., 0, 1], [rec[1], rec[2], ...,
# rec[m]]], e2(n) is [[0, 0, ..., 0, 1, 0], ...], and so on.
fundamentalLinCombs := proc(rec)
    local k, m, i:
    m := nops(rec):
    Matrix([seq([seq(-rec[k - i], i=0..k-1), 1, 0 $ m - k - 1], k=0..m-1)]):
end:

# Produce the power of Bmat(b) in terms of the first fundamental sequence,
# i.e., the solution to the recurrence with coefficients b[1], b[2], ..., b[m]
# with initial conditions 0, 0, ..., 0, 1.
EMatrix := proc(b, e, n)
    local transforms, M, eis, i, j, r:
    transforms := fundamentalLinCombs(b):

    # Transform e[i](n + j) into its e[1] form.
    M := []:
    for i from 1 to nops(b) do
        eis := [seq(add(transforms[i][r + 1] * e(n + r + j), r=0..nops(b)-1), j=0..nops(b)-1)]:

        # Replace everything beyond e(n + k - 1) by unfolding the recurrence.
        # XXX: Not efficient bounds.
        for r from 2 * nops(b) to nops(b) by -1 do
            eis := subs(e(n + r) = add(b[i] * e(n + r - i), i=1..nops(b)), eis):
        od:

        M := [op(M), eis]:
    od:

    M := LinearAlgebra[Transpose](Matrix(M)):
    Matrix([seq(convert(M[-i], list), i=1..nops(b))]):
end:

BackwardsProof := proc(b)
    local x, f, B, m, transform, Binv, k, i, d:
    f := Dio(b, x):
    B := Bmat(b):

    Binv := B^(-1):
    m := nops(b):
    transform := {seq(x[m - k] = x[m - k - 1], k=0..m-2)}:
    transform := transform union {x[1] = add(Binv[-1][i] * x[m - i + 1], i=1..m)}:

    d := simplify(subs(transform, f) - f):
end:





StIncSeqs:=proc(m, K)
    local S, K1, S1, L,i:

    if m = 0 then
        return {[]}:
    fi:

if m=1 then
 RETURN({seq([i],i=0..K)}):
fi:

    S := {}:

    for K1 from 0 to K do
        S1 := StIncSeqs(m - 1, K1-1):
   
 for L in S1 do

        S := S union {[op(L), K1]}:
od:

    od:
    S:
end:


# Input: Polynomial p, list of variables var, and a positive integer K
# Output: All positive integers K vars[1] <= vars[2] <= ... <= vars[m] <= K
#         such that p(vars[1], ..., vars[m]) = 0.
SolveDd:=proc(p, vars, K)
    local S, sl, L, i:

    S := StIncSeqs(nops(vars), K):
    sl := {}:

    for L in S do
     if igcd(op(L))=1 then
        if subs({seq(vars[i] = L[i], i=1..nops(vars))}, p) = 0 then
            sl := {op(sl), L}:
        fi:
    fi:
    od:

    sl:
end:


# Check whether the terms toCheck are contained in the first n terms of the
# C-finite sequence C.
contained := (C, L, n) -> member(L, chunks(SeqFromRec(C, n), nops(L))):

findFundamentalSolutions := proc(Ls, rec)
    local ics, L:
    ics := []:

    for L in Ls do
        if not ormap(ic -> contained([ic, rec], L, 100), ics) then
            ics := [op(ics), L]:
        fi:
    od:

    ics:
end:

chunks := proc(L, n)
    local i:
    [seq(L[i..i+n-1], i=1..nops(L)-n)]:
end:

# Example code

(*

    P := Dio([1, 1, 1, -1], x):
    Ls := SolveD(P - 1, [x[1], x[2], x[3], x[4]]):
    findFundamentalSolutions(Ls, [1, 1, 1, -1]):

*)

findYuriLike := proc(K)
    local a, b, c, P, Ls, solns, res, x:

    res := []:

    for a from 1 to K do
    for b from 1 to K do
        P := Dio([a, b, 1], x):
        Ls := SolveD(P - 1, [x[1], x[2], x[3]], 100):
        solns := findFundamentalSolutions(Ls, [a, b, 1]):
        if nops(solns) = 1 then
            res := [op(res), [a, b, 1]]:
        fi:
    od:
    od:

    res:
end:

# Find the equations Dio([a, a, 1], x) that are generated by *more* than 1
# fundamental solution (probably) from a = 1 to K, using the solutions up to M.
findBadYuriDiag := proc(K, M)
    local a, c, P, Ls, solns, res, x:

    res := []:

    for a from 1 to K do
        print(a):
        if CheckYuriBad([a, a, 1], M) then
            res := [op(res), a]:
        fi:
    od:

    res:
end:

# Check if the solutions to Dio(rec, x) - 1 are generated by more than one set
# of initial conditions. (Efficiently, without generating all possible
# solutions up to K.)
CheckYuriBad := proc(rec, K)
    option remember:
    local p, x, L, i, k, solns, res:

    p := Dio(rec, x) - 1:
    solns := {}:
    res := {}:

    L := [0 $ nops(rec)]:

    while L <> [K $ nops(rec)] do
        if subs({seq(x[i] = L[i], i=1..nops(L))}, p) = 0 then
            solns := {op(solns), L}:
            if nops(findFundamentalSolutions(solns, rec)) > 1 then
                # We know it's bad.
                return true:
            fi:
        fi:

        for i from nops(L) to 1 by -1 while L[i] = K do od:
        # i is the index of the first thing we can increment.
        # Note that i will never be zero because of the condition on the while
        # loop.
        L[i] := L[i] + 1:
        for k from i + 1 to nops(L) do
            L[k] := L[i]:
        od:
    od:

    # We think it's good.
    return false:
end:

yuriFundamentalSolutions := proc(a, M)
    local P, Ls, x:
    P := Dio([a, a, 1], x) - 1:
    Ls := SolveD(P, [x[1], x[2], x[3]], M):
    findFundamentalSolutions(Ls, [a, a, 1]):
end:

# Return the inequalities that the backwards transformation should satisfy.
BackwardsInequalities := proc(b, x)
    local f, B, m, transform, Binv, k, i, d, ineqs:
    f := Dio(b, x):
    B := Bmat(b):

    Binv := B^(-1):
    m := nops(b):
    transform := {seq(x[m - k] = x[m - k - 1], k=0..m-2)}:
    transform := transform union {x[1] = add(Binv[-1][i] * x[m - i + 1], i=1..m)}:

    ineqs := {0 <= x[1], seq(x[i] <= x[i + 1], i=1..m-1)}:
    ineqs union {f = 1}, subs(transform, ineqs)
end: