######################################################################
## 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  : Oct. 2023  `):
print():
print(`This is Hilbert10.txt, A Maple package`):
print(`accompanying Robert Dogherty-Bliss, Charles Kenney, and  Doron Zeilberger's  article: `):
print(`Using Ideas in  Matiyasevich's Seminal Proof of the Undecidability of the General `):
print(`Diophantine Equation  in order to concoct lots and lots of Decidable Ones`):



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(`-------------------------`):

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():

 print(`-------------------------`):

print(`For general help, and a list of the STORY functions,`):
print(` type "ezraSt();". For specific help type "ezra(procedure_name);" `):
print():

 print(`-------------------------`):



ezraSt:=proc()
if args=NULL then

print(`The STORY  procedures are`):
print(` DioThe, DioThe3, verboseProof`):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` Box, Coe, ConjSol, EMatrix, MinSolEven, MinSolOdd `):

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(` allSolns, Bmat, Dio, BackwardsProof, FindPromising3, SeqFromRec,  SolveDd, SolveDi, SolveDp, IsGood, IsGoodNorm `):
print(`CheckYuriBad, findBadYuriDiag, BackwardsInequalities`):


elif nargs=1 and args[1]=allSolns then
print(`allSolns(L): all solutions to the Diophantine equation generated by L. Try`):
print(`allSolns([2,3,1]);`):

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);`):

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]=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]=Box then
print(`Box(L): inputs a list of non-negative integers L, of length k, say, outputs the set [0,L[1]]x ... [0,L[k]]. Try:`):
print(`Box([2,2,2]);`):


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]=Coe then
print(`Coe(f,var): inputs a polynomial f in the list of variables var and outputs the set of coefficients. Try:`):
print(`Coe(a*x+b*y+c*x*y,[x,y]);`):


elif nargs=1 and args[1]=ConjSol then
print(`ConjSol(R,K): Inputs a list R of integers  of length k, say, and a positive integer K, finds the set of inital conditons `):
print(`such that the solutions of`):
print(`Dio(([0$(k-1),1],R],[seq(x[i],i=1..k)])=1, are unions of consecutive k-segments of SeqFromRec(INI,20) for that set. `):
print(`If IsGood(R,K) then it is a singleton. If it is not a union of such tuples then it returns false. Try:`):
print(`ConjSol([1,1,1,-1],100);`):

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]=DioThe then
print(`DioThe(R,x,K): inputs a list of integers R, of length k, say, representing an order-k recurrence, outputs a paper about a diophantine equation satisfied by ceratin recursive sequenes`):
print(`satisfying this recurrence (including [0$(k-1),1]). It also conjectures that these are the only ones. Try:`):
print(`DioThe([3,-1],x,50);`):
print(`DioThe([2,3,1],x,50);`):


elif nargs=1 and args[1]=DioThe3 then
print(`#DioThe3(a1,a2,x,r,K): inputs two positive integers a1, a2, and an integer r,  outputs a paper about all solutions of the diphantine equation`):
print(`in x[1],x[2],x[3], all non-negative integers  of the diophantine equation. K is a positive integer for the semi-rigorous investigation`):
print(`Dio([a1,a2,1],x)=r. Try:`):
print(`DioThe3(2,3,x,1,10);`):



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]=FindPromising3 then
print(`FindPromising3(K1,K2,r): inputs a small positive integer K1, and a fairly large positive integer K2, and a positive integer r, outputs all pairs 1<=a,b<=K1 such that`):
print(`MinSolOdd(A,r,K2)=MinSolOdd(A,r,2*K2), Try:`):
print(`FindPromising3(3,50,1);`):

elif nargs=1 and args[1]=IsGood then
print(`IsGood(R,K): Inputs a list R of integers  of length k, say, and a positive integer K, say, explores whether the only solutions of`):
print(`Dio(([0$(k-1),1],R],[seq(x[i],i=1..k)])=1, are consecutive k-segments of SeqFromRec([0$(k-1),R],20). Try:`):
print(`IsGood([3,-1],100);`):

elif nargs=1 and args[1]=IsGoodNorm then
print(`IsGoodNorm(R,x,N): inputs a list of integers R representing a recurrence, and a suggested norm, decides whether it seems to be a good norm for proving that all the solutions`):
print(`of the diophantine equation ispired by it is indeed given by it`):
print(`If it is positive, or noegative but very close to zero it is a good (non-rigorous!) indication that it is a good norm. Try:`):
print(`IsGoodNorm([1,1,1],x,x[1]^2+x[2]^2+x[3]^2);`):



elif nargs=1 and args[1]=MinSolEven then
print(`MinSolEven(A,r,K): inputs a list of positive integers A of length k-1, say, where k is even, and an integer K, finds all solutions (x[1],x[2],x[3],...,x[k]) `):
print(`where 0<=x[1],x[2],x[3], ...x[k]<=K that satisfy the `):
print(`the Diophantine equation Dio([op(A),1],x)=r, and in addition the strict inequality`):
print(`a1*x[k-1]+a[2]*x[k-2]+...+a[k-1]*x[1]-x[k]<0. Try:`):
print(`MinSolEven([1,1,1],1,5);`):

elif nargs=1 and args[1]=MinSolOdd then
print(`MinSolOdd(A,r,K): inputs a list of positive integers A of length k-1, say, where k is odd, and an integer K, finds all solutions (x[1],x[2],x[3],...,x[k]) `):
print(`where 0<=x[1],x[2],x[3], ...x[k]<=K that satisfy the `):
print(`the Diophantine equation Dio([op(A),1],x)=r, and in addition the strict inequality`):
print(`x[k]-a1*x[k-1]-a[2]*x[k-2]-...-a[k-1]*x[1]<0. Try:`):
print(`MinSolOdd([1,1],1,10);`):


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]=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]=SolveDi then
print(`SolveDi (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(`SolveDi(x^2+y^2-z^2,[x,y,z],40);`):




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


elif nargs=1 and args[1]=verboseProof then
print(`verboseProof(rec): Prints a proof that all solutions 0 <= x <= y <= z`):
print(`to the equation P(x, y, z) = 1, where P is obtained from Dio(rec, x),`):
print(`come from applying rec to finitely many initial solutions.`):
print(`Try:`):
print(`verboseProof([1, 1, 1])`):
print(`verboseProof([2, 1, 1])`):
print(`verboseProof([10, 1, 1])`):

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);`):


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.
SolveDi:=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 := SolveDi(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 := SolveDi(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 := SolveDi(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:

Coe1:=proc(f, x)
    local i;
    {seq(expand(coeff(f, x, i)), i = ldegree(f, x) .. degree(f, x))};
end:

Coe:=proc(f, x)
    local gu,  x1, gu1;
    if nops(x) = 1 then RETURN(Coe1(f, x[1])); end if;
    x1 := [op(2 .. nops(x), x)];
    gu := Coe(f, x1);
    {seq(op(Coe1(gu1, x[1])), gu1 in gu)};
end:



#IsGood(R,K): Inputs a list R of integers  of length k, say, and a positive integer K, say, explores whether the only solutions of
#Dio(([0$(k-1),1],R],[seq(x[i],i=1..k)])=1, are consecutive k-segments of SeqFromRec([0$(k-1),R],20). Try:
#IsGood([3,-1],100);
IsGood:=proc(R,K) local P,x,L,c,k,lu,i,n1,gu,i1:
k:=nops(R):
L:=SeqFromRec([[0$(k-1),1],R],20):
P:=Dio(R,x):

c:=expand({seq(subs({seq(x[i]=L[n1+i-1],i=1..k)},P),n1=1..nops(L)-k)}):
if nops(c)<>1 then
print(c):
  RETURN(FAIL):
fi:
c:=c[1]:

lu:=SolveDi(P-c,[seq(x[i],i=1..k)],K):

for i from 1 while L[i+k-1]<=K do od:
L:=[op(1..i+k-2,L)]:
gu:={seq([seq(L[n1+i1],i1=0..k-1)],n1=1..nops(L)-k+1)}:
evalb(lu=gu):
end:




#ConjSol1(R,K): Inputs a list R of integers  of length k, say, and a positive integer K, finds the set of inital conditons 
#such that the solutions of
#Dio(([0$(k-1),1],R],[seq(x[i],i=1..k)])=1, are unions of consecutive k-segments of SeqFromRec(INI,20) for that set. 
#If IsGood(R,K) then it is a singleton. If it is not a union of such tuples then it returns false. Try:
#ConjSol1([1,1,1,-1],100);
ConjSol1:=proc(R,K) local P,x,L,c,k,lu,i,n1,gu,i1,GU,INI:
k:=nops(R):
L:=SeqFromRec([[0$(k-1),1],R],40):
P:=Dio(R,x):

c:=expand({seq(subs({seq(x[i]=L[n1+i-1],i=1..k)},P),n1=1..nops(L)-k)}):
if nops(c)<>1 then
print(c):
  RETURN(FAIL):
fi:

c:=c[1]:

lu:=SolveDi(P-c,[seq(x[i],i=1..k)],K):

GU:={[0$(k-1),1]}:

for i from 1 while L[i+k-1]<=K do od:
L:=[op(1..i+k-2,L)]:
gu:={seq([seq(L[n1+i1],i1=0..k-1)],n1=1..nops(L)-k+1)}:
lu:=lu minus gu:

if lu={} then
 RETURN(GU):
fi:

while lu<>{} do
 INI:=lu[1]:
 L:=SeqFromRec([INI,R],40):

for i from 1 while L[i+k-1]<=K do od:
L:=[op(1..i+k-2,L)]:
gu:={seq([seq(L[n1+i1],i1=0..k-1)],n1=1..nops(L)-k+1)}:

if gu minus lu<>{} then
 RETURN(false):
else
lu:=lu minus gu:
GU:=GU union {INI}:
fi:
od:
GU:
end:




#ConjSol(R,K): Inputs a list R of integers  of length k, say, and a positive integer K, finds the set of inital conditons 
#such that the solutions of
#Dio(([0$(k-1),1],R],[seq(x[i],i=1..k)])=1, are unions of consecutive k-segments of SeqFromRec(INI,20) for that set. 
#If IsGood(R,K) then it is a singleton. If it is not a union of such tuples then it returns false. Try:
#ConjSol([1,1,1,-1],100);
ConjSol:=proc(R,K) local P,gu,gu1,mu,x,i,k,n1,L:
k:=nops(R):
P:=Dio(R,x)-1:
 gu:=ConjSol1(R,K):
 if gu=FAIL or gu=false then
  RETURN(gu):
 fi:

 gu1:=ConjSol1(R,2*K):

if gu<>gu1 then
 RETURN(false):
else
 gu1:=ConjSol1(R,3*K):
if gu<>gu1 then
 RETURN(false):
else

for mu in gu do
L:=SeqFromRec([mu,R],100):
if {seq(subs({seq(x[i]=L[n1+i-1],i=1..k)},P),n1=1..nops(L)-k)}<>{0} then
 RETURN(FAIL):
fi:
od:

RETURN(gu):
fi:
fi:
end:





#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):


L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:


#DioThe(R,x,K): inputs a list of integers R, of length k, say, representing an order-k recurrence, outputs a paper about a diophantine equation satisfied by ceratin recursive sequenes
#satisfying this recurrence (including [0$(k-1),1]). It also conjectures that these are the only ones. Try:
#DioThe([-3,1],x,50);
#DioThe([2,3,1],x,50);
DioThe:=proc(R,x,K) local k,gu,mu,f,a,t0,j,i1,P,i,j1:
t0:=time():
k:=nops(R):
gu:=ConjSol(R,K):
if gu=FAIL or gu=false then
 RETURN(gu):
fi:

P:=Dio(R,x):

print(`On Solutions of the Diophantine equation`, P=1 ):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

if nops(gu)=1 then
print(`The following set of increasing positive integers`, [seq(x[i],i=1..k)], `  are definitely solutions of the diophantine equation`):
else
print(`The following`, nops(gu),  `sets of increasing positive integers`, [seq(x[i],i=1..k)], `  are definitely solutions of the diophantine equation`):
fi:

print(``):
print(P=1):
print(``):

for i from 1 to nops(gu) do
 mu:=gu[i]:
 f[i]:=CtoR([mu,R],x):
print(``):
print(`-----------------------------------------------------`):
print(``):

print(`Define the sequence `, a[i][j], ` for j from 1 to infinity as the solution of the recurrence `):
print(``):
print(a[i][j]=add(R[i1]*a[i][j-i1],i1=1..k)):
print(``):
print(`Subject to the initial conditions `):
print(``):
print(seq(a[i][j1-1]=mu[j1],j1=1..nops(mu))):
print(``):
print(`Or equivalently in terms of the generating function`):
print(``):
print(f[i]=Sum(a[i][j]*x^j,j=0..infinity)):
print(``):
print(`Then the following increasing `, k, ` -tuples `):
print(``):
print([seq(x[i1],i1=1..k)]=[seq(a[i][j+i1-1],i1=1..k)]):
print(``):
print(`for j from 0 to infinity `):

print(``):
print(`are solutions of the diophantine equation`, P=1):
od:

print(`Furthermore, we conjecture that these are all the increasing positive solutions`):

print(`-----------------------------------------`):
print(``):
print(`This ends thie article that took`, time()-t0, `to create. `):

end:



#IsGoodNorm(R,x,N): inputs a list of integers R representing a recurrence, and a suggested norm, decides whether it is a good norm for proving that all the solutions
#of the diophantine equation ispired by it is indeed given by it. 
#It if is positive or noegative but very close to zero it is a good (non-rigorous!) indication that it is a good norm
#Try:
#IsGoodNorm([1,1,1],x,x[1]^2+x[2]^2+x[3]^2);
IsGoodNorm:=proc(R,x,N) local k,P,B,i,j,N1:
Digits:=100:
k:=nops(R):
P:=Dio(R,x):
B:=inverse(Bmat(R)):
if expand(subs({seq(x[k-i+1]=add(B[i,j]*x[k-j+1],j=1..k),i=1..k)},P)-P)<>0 then
 RETURN(FAIL):
fi:

N1:=expand(N-subs({seq(x[k-i+1]=add(B[i,j]*x[k-j+1],j=1..k),i=1..k)},N)):


#gu:=Optimization[Minimize](N1,{P=0,seq(x[k-i+1]>=x[k-i],i=1..k-1),x[1]>=0}):
Optimization[Minimize](N1,{P=0,seq(x[k-i+1]>=x[k-i],i=1..k-1),x[1]>=0},initialpoint={seq(x[i]=i,i=1..k)})[1]:

#if gu[1]>=-evalf(10^(-Digits/4)) then
# true:
#else
#false:
#fi:

end:

recToBackwards := proc(L, x)
    local d, k:
    d := nops(L):
    simplify((x[d] - add(L[-k - 1] * x[k], k=1..d-1)) / L[d]):
end:

# find the absolute minimum of f(t, s) on the region Q(t, s) <= 0 within the
# unit cube. we assume that Q(t, s) is linear in t and s, and has negative
# coefficients on both. thus the region is bounded by Q(t, s) = 0, t = 1, and s
# = 1.
infolevel[findAbsoluteMin] := 1:
findAbsoluteMin := proc(f, Q, t, s)
    local expr, a, b, region_1, region_2, region_3, min_1, min_2, min_3, crits,
           crit_values, crit, val, min_crit:
    # examine the boundary Q(t, s) = 0.
    expr := solve(Q, s):
    a := solve(expr = 1, t):
    b := solve(expr = 0, t):
    region_1 := subs(s = solve(Q, s), f):
    min_1 := simplify(minimize(region_1, t=a..b)):

    # Q(t, s) should have negative coefficients on s and t, so Q(t, s) <= 0 is
    # bounded by the edges s = 1 and t = 1.
    region_2 := subs(s = 1, f):
    min_2 := simplify(minimize(region_2, t=0..1)):

    region_3 := subs(t = 1, f):
    min_3 := simplify(minimize(region_3, s=0..1)):

    # check for critical points.
    crits := solve({diff(f, s), diff(f, t)}, {s, t}):

    crit_values := []:
    for crit in crits do
        val := evalf(subs(crit, Q)):
        if val <= 0 then
            userinfo(1, findAbsoluteMin, `making a floating point guess:`, val <= 0):
            crit_values := [op(crit_values), simplify(subs(crit, f))]:
        fi:
    od:
    min_crit := min(crit_values):

    min(min_1, min_2, min_3, min_crit):
end:

infolevel[verboseProof] := 1:
verboseProof:=proc(rec)
    local P, deg, f, B, Q, pos_min, order_min, pos_bound, order_bound, s, t, x, y, z:
    if nops(rec) <> 3 then
        userinfo(1, verboseProof, `only third order recurrences are allowed`):
        return:
    elif min(rec) <= 0 then
        userinfo(1, verboseProof, `only nonnegative coefficients are allowed`):
        return:
    elif rec[3] <> 1 then
        userinfo(1, verboseProof, `last coefficient must be 1`):
        return:
    fi:

    P := Dio(rec, x):
    P := subs({x[1] = x, x[2] = y, x[3] = z}, P):
    deg := degree(P, z):
    f := P / z^deg:
    f := expand(subs({x = t * z, y = s * z}, f)):

    B := recToBackwards(rec, [x, y, z]):
    Q := expand(subs({x = t * z, y = s * z}, B / z)):

    pos_min := findAbsoluteMin(f, Q, t, s):
    order_min := findAbsoluteMin(f, -Q + t, t, s):

    if evalf(pos_min) < 0 then
        userinfo(1, verboseProof, `positive proof didn't work!`):
        return:
    elif evalf(order_min) < 0 then
        userinfo(1, verboseProof, `order proof didn't work!`):
        return:
    fi:

    print(`THEOREM. The nonnegative, increasing solutions to the diophantine equation`):
    print(P = 1):
    print(`are generated by applying the recurrence`):
    print(rec):
    print(`to finitely many initial solutions.`):

    print(`PROOF. Let P be the polynomial P on the left-hand side.`):
    print(`Note that it is invariant under the recurrence:`):
    print(`P - P(shift) is`, expand(P - subs({x = y, y = z, z = rec[1] * z + rec[2] * y + rec[3] * x}, P))):

    print(`The backwards shift formula to get the previous term from the triple (x, y, z) is`):
    print(B):
    print(`We will show that this backwards shift gives a smaller increasing solution for sufficiently large z.`):

    print(`Divide both sides of our equation by`, z^deg, `and make the change of variables {t = x / z, s = y / z}`):
    print(`This gives`):
    print(f = 1 / z^deg):
    print(`where (t, s) is in the unit square.`):

    print(`Let (x, y, z) be a generic solution. Then the following inequalities are routine calculus exercises.`):

    print(`the inequality`, B > 0, `also known as`, Q > 0, `holds for`):
    pos_bound := 1 / pos_min^(1/deg):
    print(z >= pos_bound):
    print(`more explicitly for`):
    print(z >= evalf(pos_bound)):

    print(`the inequality`, B < x, `also known as`, Q < t, `holds for`):
    order_bound := 1 / order_min^(1/deg):
    print(z >= order_bound):
    print(`more explicitly for`):
    print(z >= evalf(order_bound)):

    print(`We only need to look for solutions with`, z < max(evalf(order_bound), evalf(pos_bound))):
    print(`and there are finitely many of these.`):
    print(`Q.E.D.`):
end:





#DioThe3(a1,a2,x,r,K): inputs two positive integers a1, a2  and an intgeger r, outputs a paper about all solutions of the diphantine equation
#in x[1],x[2],x[3], all non-negative integers  of the diophantine equation. K is a positive integer for the semi-rigorous investigation
#Dio([a1,a2,1],x)=1. Try:
#DioThe3(2,3,x,1,10);
DioThe3:=proc(a1,a2,x,r,K) local gu,mu,f,a,t0,j,i1,P,i,j1:
t0:=time():


P:=Dio([a1,a2,1],x):
gu:=MinSolOdd([a1,a2],r,K):
print(`On Solutions of the Diophantine equation`, P=r ):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):



if gu={} then
print(`There are no solutions in non-negative triples [x[1],x[2],x[3]] of the Diophantine equation `):

print(``):
print(P=r):
print(``):
print(`---------------------------------------`):
print(``):
print(`Proof: It is readily seen that if [x[1],x[2],x[3]] is a solution of the diophantine equation `):
print(``):
print(P=r):
print(``):
print(`Then so is  the triple`):
print(``):
print([x[3]-a1*x[2]-a2*x[1],x[1],x[2]]):
print(``):
print(`So the mapping `):
print(``):
print(`[x[1],x[2],x[3]] -> `, [x[3]-a1*x[2]-a2*x[1],x[1],x[2]]):
print(``):
print(`maps solution triples to other solution triples`):
print(``):
print(`If you start with any solution triple with non-negative entries sooner or later you would have the scenario that for the FIRST time, the first component, x[1], is negative`):
print(``):
print(`It is readily seen that there do not exist such  non-negative integers [x[1],x[2],x[3]]  that satisfy the above diophantine equation and in addition satisfy the inequality`):
print(``):
print(x[3]-a1*x[2]-a2*x[1]<0):
print(``):
print(`hence there are no solutions at all.`):
print(``):
print(` QED. `):
print(``):
print(`-----------------------------------------`):
print(``):
print(`This ends thie article that took`, time()-t0, `to create. `):

RETURN():
fi:
 if nops(gu)=1 then
print(`The following set of triples of non-negative integers`, [seq(x[i],i=1..3)], `  are the ONLLY solutions (with non-negative integers) of the diophantine equation`):
else
print(`The following`, nops(gu),  `sets of triples of non-negative integers`, [seq(x[i],i=1..3)], `  are the ONLY solutions (with non-negative integers) of the diophantine equation`):
fi:

print(``):
print(P=r):
print(``):

for i from 1 to nops(gu) do
 mu:=gu[i]:
 f[i]:=CtoR([mu,[a1,a2,1]],x):
print(``):
print(`-----------------------------------------------------`):
print(``):

print(`Define the sequence `, a[i][j], ` for j from 0 to infinity as the solution of the recurrence `):
print(``):
print(a[i][j]=a1*a[i][j-1]+a2*a[i][j-2]+ a[i][j-3]):
print(``):
print(`Subject to the initial conditions `):
print(``):
print(seq(a[i][j1-1]=mu[j1],j1=1..nops(mu))):
print(``):
print(`Or equivalently in terms of the generating function`):
print(``):
print(f[i]=Sum(a[i][j]*x^j,j=0..infinity)):
print(``):
print(`Then the following  triples of non-negative solutions `):
print(``):
print([seq(x[i1],i1=1..3)]=[seq(a[i][j+i1-1],i1=1..3)]):
print(``):
print(`for j from 0 to infinity `):

print(``):
print(`are solutions of the diophantine equation`, P=r):
od:

print(``):
print(`---------------------------------------`):
print(``):
print(`Proof: It is readily seen that if [x[1],x[2],x[3]] is a solution of the diophantine equation `):
print(``):
print(P=r):
print(``):
print(`Then so is  the triple`):
print(``):
print([x[3]-a1*x[2]-a2*x[1],x[1],x[2]]):
print(``):
print(`So the mapping `):
print(``):
print(`[x[1],x[2],x[3]] -> `, [x[3]-a1*x[2]-a2*x[1],x[1],x[2]]):
print(``):
print(`maps solution triples to other solution triples`):
print(``):
print(`If you start with any solution triple with non-negative entries sooner or later you would have the scenario that for the FIRST time, the first component, x[1], is negative`):
print(``):
print(`It is readily seen that the only such triples of non-negative integers [x[1],x[2],x[3]]  that satisfy the above diophantine equation and in addition satisfy the inequality`):
print(``):
print(x[3]-a1*x[2]-a2*x[1]<0):
print(``):
print(`are the following triples`):
print(``):
print(gu):
print(``):
print(`Going forward gives the above proposed solutions `):
print(``):
print(` QED. `):
print(``):
print(`-----------------------------------------`):
print(``):
print(`This ends thie article that took`, time()-t0, `to create. `):

end:

 









#Box(L): inputs a list of non-negative integers L, of length k, say, outputs the set [0,L[1]]x ... [0,L[k]]. Try:
#Box([2,2,2]);
Box:=proc(L) local k,i,gu,L1,gu1:
option remember:
k:=nops(L):
if k=0 then
 RETURN({[]}):
fi:

L1:=[op(1..k-1,L)]:
gu:=Box(L1):

{seq(seq([op(gu1),i],i=0..L[k]),gu1 in gu)}:

end:





#MinSolOdd(A,r,K): inputs a list of positive integers A of length k-1, say, where k is odd, and an integer K, finds all solutions (x[1],x[2],x[3],...,x[k]) 
#where 0<=x[1],x[2],x[3], ...x[k]<=K that satisfy the 
#the Diophantine equation Dio([op(A),1],x)=r, and in addition the strict inequality
#x[k]-a1*x[k-1]-a[2]*x[k-2]-...-a[k-1]*x[1]<0. Try:
#MinSolOdd([1,1],1,5);
MinSolOdd:=proc(A,r,K) local k,x,P,i1,gu,lu,i,lu1:
k:=nops(A)+1:
if k mod 2<>1 then
print(`The list`, A , `must be of even length `):
 RETURN(FAIL):
fi:

P:=Dio([op(A),1],x)-r:



if expand(subs({x[1]=x[k]-add(A[i]*x[k-i],i=1..k-1),seq(x[i]=x[i-1],i=2..k)},P)-P)<>0 then
 print(`Something bad happened`):
 RETURN(FAIL):
fi:



lu:=Box([K$k]):


gu:={}:

for lu1 in lu do
     if subs({seq(x[i1]=lu1[i1],i1=1..k)},P)=0 and lu1[k]-add(A[i1]*lu1[k-i1],i1=1..k-1)<0 then
       gu:=gu union {lu1}:
      fi:
od:
gu:
end:



#MinSolEven(A,r,K): inputs a list of positive integers A of length k-1, say, where k is even, and an integer K, finds all solutions (x[1],x[2],x[3],...,x[k]) 
#where 0<=x[1],x[2],x[3], ...x[k]<=K that satisfy the 
#the Diophantine equation Dio([op(A),1],x)=r, and in addition the strict inequality
#a1*x[k-1]+a[2]*x[k-2]+...+a[k-1]*x[1]-x[k]<0. Try:
#MinSolEven([1,1,1],1,5);
MinSolEven:=proc(A,r,K) local k,x,P,i1,gu,lu,i,lu1:
k:=nops(A)+1:
if k mod 2<>0 then
print(`The list`, A , `must be of odd length `):
 RETURN(FAIL):
fi:

P:=Dio([op(A),-1],x)-r:



if expand(subs({x[1]=add(A[i]*x[k-i],i=1..k-1)-x[k],seq(x[i]=x[i-1],i=2..k)},P)-P)<>0 then
 print(`Something bad happened`):
 RETURN(FAIL):
fi:



lu:=Box([K$k]):


gu:={}:

for lu1 in lu do
     if subs({seq(x[i1]=lu1[i1],i1=1..k)},P)=0 and add(A[i1]*lu1[k-i1],i1=1..k-1)-lu1[k]<0 then
       gu:=gu union {lu1}:
      fi:
od:
gu:
end:



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

    S := Box([K$nops(vars)]) :
    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:


#FindPromising3(K1,K2,r): inputs a small positive integer K1, and a fairly large positive integer K2, and a positive integer r, outputs all pairs 1<=a,b<=K1 such that
#MinSolOdd([a,b],r,K2)=MinSolOdd([a,b],r,2*K2), Try:
#FindPromising3(3,50,1);
FindPromising3:=proc(K1,K2,r) local gu,a,b:
gu:={}:

for a from 1 to K1 do
 for b from 1 to K1 do
  if MinSolOdd([a,b],r,K2)=MinSolOdd([a,b],r,2*K2) then
    gu:=gu union {[a,b]}:
  fi:
 od:
od:
gu:
end:

infolevel[allSolns] := 1:
allSolns:=proc(rec)
    local P, deg, f, B, Q, pos_min, order_min, pos_bound, order_bound, s, t, x,
    y, z, max_bound, solns, z0, y0, x0, reduced_solns, soln, a, b, c:
    if nops(rec) <> 3 then
        userinfo(1, allSolns, `only third order recurrences are allowed`):
        return FAIL:
    elif min(rec) <= 0 then
        userinfo(1, allSolns, `only nonnegative coefficients are allowed`):
        return FAIL:
    elif rec[3] <> 1 then
        userinfo(1, allSolns, `last coefficient must be 1`):
        return FAIL:
    fi:

    P := Dio(rec, x):
    P := subs({x[1] = x, x[2] = y, x[3] = z}, P):
    deg := degree(P, z):
    f := P / z^deg:
    f := expand(subs({x = t * z, y = s * z}, f)):

    B := recToBackwards(rec, [x, y, z]):
    Q := expand(subs({x = t * z, y = s * z}, B / z)):

    pos_min := findAbsoluteMin(f, Q, t, s):

    order_min := findAbsoluteMin(f, -Q + t, t, s):

    if evalf(pos_min) < 0 then
        userinfo(1, allSolns, `positive proof didn't work!`):
        return FAIL:
    elif evalf(order_min) < 0 then
        userinfo(1, allSolns, `order proof didn't work!`):
        return FAIL:
    fi:

    pos_bound := 1 / pos_min^(1/deg):
    order_bound := 1 / order_min^(1/deg):
    max_bound := ceil(max(order_bound, pos_bound)):

    solns := {}:
    userinfo(1, allSolns, `looking for monotonically increasing solutions up to`, max_bound):
    for z0 from 0 to max_bound do
    for y0 from 0 to z0 do
    for x0 from 0 to y0 do
        if simplify(subs({x = x0, y = y0, z = z0}, P) - 1) = 0 then
            solns := solns union {[x0, y0, z0]}:
        fi:
    od:
    od:
    od:

    reduced_solns := {}:

    while nops(solns) > 0 do
        soln := solns[1]:
        reduced_solns := reduced_solns union {soln}:
        a, b, c := soln[1], soln[2], soln[3]:
        while c < max_bound do
            solns := solns minus {[a, b, c]}:
            a, b, c := b, c, rec[1] * c + rec[2] * b + rec[3] * a:
        od:
    od:

    reduced_solns:
end: