<< RISC`LinearSystemSolver`
<< "BimonicRecurrences.m";
temp = ""; (* directory for storing temporary results *)
twist = ""; (* path to twist knot recurrences *)
tw = -3; (* Cases 0 and 1 are trivial, start with -1, total range is from -14 to 15. *)
While[tw > -15,
Print["=========================================================="];
Print["Starting with tw = ", tw];
Print["=========================================================="];
file = "tw" <> ToString[tw];
rec = Collect[Get[twist <> "rec.twist.knot." <> ToString[tw] <> ".m"] /. f[a_] :> f[a]/(q^a - 1), f[_], Together];
inord = Max[Cases[rec, f[n + a_.] :> a, Infinity]];
Print["Input order is ", inord];
ord = inord + 1; deg = -1/2; sol = 0;
While[sol === 0,
  deg = Floor[3/2 * deg] + 1;
  If[deg > 209, ord++; deg = 0];
  sol = BimonicRecurrence[rec, ord, deg, Return -> "nsol"];
];
While[sol > 0,
  deg--;
  sol = BimonicRecurrence[rec, ord, deg, Return -> "nsol"];
];
Print["ord = ", ord, ", deg = ", deg, ", sol = ", sol];
qn = 0;
While[sol === 0 && qn < 200,
  qn++;
  sol = BimonicRecurrence[rec, ord, deg, Power -> qn, Return -> "nsol"];
];
If[sol === 0, qn = Automatic; deg++];
sol = BimonicRecurrence[rec, ord, deg, Power -> qn, Return -> "nsol"];
Print["Have configuration: ord = ", ord, ", deg = ", deg, ", qn = ", qn, ", sol = ", sol];
If[sol === 0, Print["Error: 0 solutions."]; Abort[]];
beg = 0; end = 0; min = Infinity; minpos = {};
While[beg <= ord,
  sol1 = BimonicRecurrence[rec, ord, deg, Power -> qn, Begin -> beg, End -> end, Return -> "nsol"];
  If[sol1 === sol,
    dim = First[BimonicRecurrence[rec, ord, deg, Power -> qn, Begin -> beg, End -> end, Return -> "dims"]];
    If[dim < min, min = dim; minpos = {beg, end}];
    beg++; end = beg;
  ,
    end++;
    If[end > ord, beg = ord + 1];
  ];
];
If[minpos === {}, Print["Error with minpos."]; Abort[]];
{beg, end} = minpos;
Print["Minimized range: beg = ", beg, ", end = ", end];
vecs = vecs1 = {Table[0, {ord - inord - 1}]}; zeros = {}; vals = {sol};
todegs[vec_] := Table[{0, deg}, {ord - inord - 1}] + (({1, -1}#)& /@ Partition[Join[vec, Reverse[vec]], 2]);
(*While[Max[vals] > 0,
  vecs1 = vecs;
  If[Length[vecs] === Max[vals] === 1, Break[]];
  vecs = Union[Flatten[Function[v, (v + #)& /@ IdentityMatrix[Length[v]]] /@ vecs, 1]];
  vecs = Select[vecs, Max[#] <= Ceiling[deg / 2] &];
  vecs = Select[vecs, Function[v, Max[Min[v - #]& /@ zeros] < 0]];
  Print["Testing ", Length[vecs], " vecs."];
  vals = Function[vec,
    BimonicRecurrence[rec, ord, todegs[vec], Power -> qn, Begin -> beg, End -> end, Return -> "nsol"]] /@ vecs;
  zeros = Join[zeros, vecs[[Flatten[Position[vals, 0]]]]];
  vecs = vecs[[Flatten[Position[vals, Max[vals]]]]];
  Print["Starting again with ", Length[vecs], " vecs."];
];*)
vec = First[vecs1];
sol = BimonicRecurrence[rec, ord, todegs[vec], Power -> qn, Begin -> beg, End -> end, Return -> "nsol"];
Print["Have ", sol, " solutions."];
pos = 1;
While[pos <= ord - inord - 1,
(*
  vec[[pos]]++;
  sol = BimonicRecurrence[rec, ord, todegs[vec], Power -> qn, Begin -> beg, End -> end, Return -> "nsol"];
  If[sol === 0, vec[[pos]]--; pos++];
*)
  (* Additional symmetry for twist knots. *)
  vec[[pos]]++; If[pos < Length[vec], vec[[pos+1]]++];
  sol = BimonicRecurrence[rec, ord, todegs[vec], Power -> qn, Begin -> beg, End -> end, Return -> "nsol"];
  If[sol === 0, vec[[pos]]--; If[pos < Length[vec], vec[[pos+1]]--]; pos = pos + 2];
];
deg = todegs[vec];
sol = BimonicRecurrence[rec, ord, deg, Power -> qn, Begin -> beg, End -> end, Return -> "nsol"];
Print["Have degs = ", deg, " with ", sol, " solutions."];
If[sol === 0, Print["Error: 0 solutions."]; Abort[]];
r = 500; test = "";
While[StringQ[test],
  Print["Computing evaluation points up to ", r];
  BimonicRecurrence[rec, ord, deg, Power -> qn, Begin -> beg, End -> end,
    Evaluate -> Range[Max[r - 500 + 1, 3], r], FileName -> temp <> file];
  test = Catch[InterpolateAndReconstruct[temp <> file, {NextPrime[2^31, -2]}, Print -> True]];
  r += 500;
];
qdeg = Max[Get[temp <> file <> "_qdegs.m"]] + 20;
Print["qdeg = ", qdeg];
pp = 3; test = "";
While[StringQ[test],
  Print["Compute data for NextPrime[2^31, -", pp, "]"];
  BimonicRecurrence[rec, ord, deg, Power -> qn, Begin -> beg, End -> end,
    Modulus -> NextPrime[2^31, -pp], Evaluate -> Range[3, qdeg], FileName -> temp <> file];
  InterpolateAndReconstruct[temp <> file, {NextPrime[2^31, -pp]}];
  test = Catch[RationalCoefficients[temp <> file]];
  pp++;
];
reconst = Get[temp <> file <> "_reconst.m"];
Print["reconst: ", {Dimensions[#], ByteCount[#], N[Max[Abs[Cases[#, _Integer, Infinity]]]]} &[reconst]];
ansatz = BimonicRecurrence[rec, ord, deg, Power -> qn, Return -> "ansatz"];
lincomb = ansatz /. C[i_] :> reconst[[i]];
Put[lincomb, twist <> "lincomb.twist." <> ToString[tw] <> ".m"];
bmrec = Collect[lincomb /. f[a_] :> (rec /. n -> a) /. q^a_ :> q^Expand[a], f[_]];
cf = (Coefficient[ansatz, f[n]] /. C[1] -> 1) * Coefficient[rec, f[n]];
cf = Together[cf / q^(Exponent[cf, q^n, Min] * n)];
cf = Together[cf / (PolynomialGCD @@ CoefficientList[cf, q^n])];
bmrec = Collect[bmrec, f[_], Together[# / cf]&];
bmrec = Collect[bmrec, f[_], 
  Function[c, Factor[c] /. a_Plus :> Collect[a, q^n, With[{d = Exponent[#, q, Min]}, q^d*Expand[#/q^d]] &]]];
Put[bmrec, twist <> "bmrec.twist." <> ToString[tw] <> ".m"];
Run["rm " <> temp <> file <> "_*"];
tw = -tw + If[tw < 0, 1, 0];
]; (* end While *)
Quit[];
