
(* LinearSystemSolver. by Manuel Kauers *)

LinearSystemSolver`Private`Version = "0.23 2012-06-06";
(* for MMA6 from version 0.6 on *)
(* for MMA7 from version 0.9 on *)
(* for MMA8 from version 0.19 on *)

SetAttributes[HornerForm, Listable];

If[$VersionNumber >= 6,
BeginPackage["LinearSystemSolver`"]
,
BeginPackage["LinearSystemSolver`", {"DiscreteMath`Combinatorica`", "NumberTheory`NumberTheoryFunctions`"}];

(* The command NextPrime takes only one parameter in Mathematica Version < 6. *)
Unprotect[NumberTheory`NumberTheoryFunctions`NextPrime];
NumberTheory`NumberTheoryFunctions`NextPrime[x_?NumericQ, n_Integer?Positive] := 
  Nest[NumberTheory`NumberTheoryFunctions`NextPrime, x, n];
NumberTheory`NumberTheoryFunctions`NextPrime[x_?NumericQ, 0] := NextPrime[x];
NumberTheory`NumberTheoryFunctions`NextPrime[x_?NumericQ, n_Integer?Negative] := 
  Nest[NumberTheory`NumberTheoryFunctions`PreviousPrime, x, -n];
Protect[NumberTheory`NumberTheoryFunctions`NextPrime];
];

LinSolveAlg::usage = "Compute ker(A) for a given A in Q(alpha)^(n times m) with alpha algebraic over Q. " <>
 "Ex: LinSolveAlg[{{1,alpha,3},{alpha,2,3}}, alpha, alpha^2-5]";
LinSolveQ::usage = "Compute ker(A) for a given A in Q^(n times m). Ex: LinSolveQ[{{1,2,3},{4,5,6}}]";
LinSolveUniv::usage = 
 "Compute ker(A) for a given A in Q[x]^(n times m). Ex: LinSolveUniv[{{1,x,3},{x,4,x}}, modulus, x]";
LinSolveHensel::usage = 
 "Solve an inhomogeneous system A.x=B for A in Z^(n times n) and B in Z^(n times m). " <>
 "Ex: LinSolveHensel[A, B]";
Lift::usage =
 "Given a list of homomorphic images {phi[1],phi[2],...}, and a list of moduli {r[1],r[2],...}, compute a " <>
 "rational function rat such that rat mod r[i] == phi[i] for all i. The phi[i] may be themselves rational " <>
 "functions. In this case, the function is applied to corresponding coefficients. " <>
 "Ex: Lift[{...},{...},{y,z},x,modulus] for rational function reconstruction " <>
 "or Lift[{...},{...},{x,y,z}] for rational number reconstruction."
  
Begin["`Private`"]

If[$VersionNumber < 6, 
  Quiet[a_, ___] := a;
];

PrintLevel = -1;

DebugPrint[u_, v__] := If[u <= PrintLevel, Print[v]];

(** ------------------------------------------------------------------- **)
(* Arithmetic *)

List2Poly[coeffs_List, var_] := coeffs . var^(Range[Length[coeffs]]-1);
Poly2List[p_, var_, deg_] := Module[ {l}, 
   l = CoefficientList[p, var]; 
   If[ Length[l] <= deg + 1, 
       Return[Join[l, Table[0, {deg + 1 - Length[l]}]]],
       Return[Take[l, deg + 1]]
   ]];

(**
 * polynomial division
 **)
QuoRem[p_, q_, x_, m_] :=
  Module[ {p0, q0, quo, c, i, j},

    p0 = CoefficientList[p, x]; q0 = CoefficientList[q, x]; quo = {};

    If[ Length[p0] < Length[q0], Return[ {0, p} ] ];

    While[ Length[p0] >= Length[q0],

      c = PolynomialMod[Last[p0]/Last[q0], m];
      PrependTo[quo, c];

      j = Length[p0] - Length[q0] + 1;
      Do[ p0[[j++]] = p0[[j]] - c q0[[i]], {i, 1, Length[q0]} ]; 
      p0 = Most[p0];
    ];

    Return[ {List2Poly[quo, x], List2Poly[p0, x]} ]
  ];

MyPolynomialQuotient[p_, q_, x_, m_] := First[QuoRem[p, q, x, m]];
MyPolynomialRemainder[p_, q_, x_, m_] := Last[QuoRem[p, q, x, m]];

(** ------------------------------------------------------------------- **)
(* Reconstruction *)

Lift[rats_, points_, vars_, x_, p_Integer] :=
  Module[ {u, myvars, data, den, terms, lt, res, deg},
    If[ Length[rats] =!= Length[points], Throw["illegal argument"] ];
    myvars = Append[vars, u];
    data = PolynomialSupport[#, myvars, Modulus -> p]& /@ 
          (Numerator[#] + u Denominator[#]&) /@ Together[rats, Modulus -> p]; 
    If[ Length[Union[Length /@ data]] =!= 1, Throw["unexpected zeros."] ]; 
    terms = (#/(#/.Thread[myvars->1]))&[data[[1]]];
    (* normalize *)
    lt = Position[terms, Last[Sort[terms]]][[1, 1]];
    data = Transpose[#/#[[lt]]& /@ (data /. Thread[myvars -> 1])];
    (* find denominator and clear it *)
    den = Sum[Prime[500+i]data[[i]], {i, 1, Length[terms]}];
    den = InterpolatingPolynomial[Transpose[{points, den}], x, Modulus -> p];
    den = ReconstructRatfun[den, points, x, p];
    deg = 1 + Exponent[Numerator[den], x]; den = Denominator[den];
    data = Transpose[Table[data[[All, i]]*(den /. x -> points[[i]]), {i, 1, Length[points]}]];
    (* interpolation *)
    (#1/#2)&@@CoefficientList[Table[Expand[InterpolatingPolynomial[Transpose[Take[#, deg]& /@ {points, data[[i]]}], x, Modulus -> p], Modulus -> p], {i, 1, Length[terms]}].terms, u]
  ];
Lift[rats_, points_, vars_, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {u, myvars, data, den, terms, lt, res, deg, p, nondefective},
    recon = "recon" /. {opts} /. {"recon" -> ReconstructQ0};
    If[ Length[rats] =!= Length[points], Throw["illegal argument"] ];
    myvars = Append[vars, u];
    data = Table[PolynomialSupport[(Numerator[#]+u Denominator[#]&)[Together[rats[[i]], Modulus -> points[[i]]]], myvars, Modulus -> points[[i]]], {i, 1, Length[points]}];
    If[ Length[Union[Length /@ data]] =!= 1, 
       nondefective = Position[Length /@ data, Max[Length /@ data]];
       Print["discarding ", Length[points] - Length[nondefective], " images."];
       Return[Lift[Extract[rats, nondefective], Extract[points, nondefective], vars, opts]];
    ]; 
    terms = (#/(#/.Thread[myvars->1]))&[data[[1]]];
    (* normalize *)
    lt = Position[terms, Last[Sort[terms]]][[1, 1]];
    data = Transpose[#/#[[lt]]& /@ (data /. Thread[myvars -> 1])];
    data = Transpose[Table[PolynomialMod[data[[All, i]], points[[i]]], {i, 1, Length[points]}]];
    (* reconstruct one by one *)
    data = Table[recon[ChineseRemainder[data[[i]], points], Times@@points], {i, 1, Length[terms]}];
    (* data *= LCM@@(Denominator /@ data); *)
    (#1/#2)&@@CoefficientList[data.terms, u]
  ];

PolynomialSupport[poly_, vars_, opts___] := (* multivariate polynomial expression to list of monomials *)
  Module[ {x, rest, terms},
    If[ poly === 0, Return[ {} ] ];
    x = First[vars]; rest = Rest[vars];
    terms = CoefficientList[poly, x, opts];
    If[ Length[rest] > 0, terms = PolynomialSupport[#, rest, opts]& /@ terms];
    DeleteCases[Flatten[terms*x^Range[0, Exponent[poly, x, opts]]], 0]
  ];

(**
 * Reconstruction of rational numbers
 **)
ReconstructQ[n_, pq_] := (#[[1]]/#[[2]])&[First[LatticeReduce[{{n, 1}, {pq, 0}}]]];

ReconstructQ0[n_, pq_] := If[ n===0, 0, (((#[[2,2]]/#[[1,2,2]])&)[Internal`HGCD[pq, Mod[n,pq]]]*2)/2 ];

(* use this if you expect u : v to be close to  #digits in numerator : #digits in denominator. *)
ReconstructQbiased[u_, v_][n_, q_] := If[ n==0, 0,
  Module[ {r, s, rr, ss, quo},
    {r, s, rr, ss} = {q, 0, n, 1};
    While[ v Log[N[1 + Abs[rr]]] > u Log[N[1 + Abs[ss]]],
      quo = Quotient[r, rr];
      {r, s, rr, ss} = {rr, ss, r - quo rr, s - quo ss}
    ];
    rr/ss
  ]];

(* this gives the smallest of all *)
ReconstructQ2[n_, q_] := If[ n==0, 0,
  Module[ {r, s, rr, ss, quo, u, v},
    {u, v} = {n, 1}; 
    {r, s, rr, ss} = {q, 0, Mod[n, q], 1};
    While[ rr != 0, 
      quo = Quotient[r, rr];
      {r, s, rr, ss} = {rr, ss, r - quo rr, s - quo ss};
      If[ s != 0 && Abs[r*s] < Abs[u*v], {u, v} = {r, s} ]
    ];
    u/v
  ]];

(**
 * Reconstruction of rational functions
 **)
Do[T[i]=0,{i,1,6}];
ReconstructRatfun[g_, mod_, x_, p_] := 
  Module[ {t, r},
    {t, r} = Timing[MyReconstructRatfun[g, mod, x, p]];
    T[1] += t;
    Return[r];    
  ];

MyReconstructRatfun[g_, mod_, x_, p_] :=
  Module[ {t, f, r0, r1, t0, t1, n, d, q, c, df, drat, dr1, dt1, unique},

    c = PolynomialMod[#, p]&;
    If[ c[g] === 0, Return[ 0 ] ];

    f = If[ ListQ[mod], Product[x - mod[[i]], {i, 1, Length[mod]}], mod ];
    df = Exponent[f, x];
    {r0, r1, t0, t1} = {f, g, 0, 1}; {n, d} = {r1, t1};
    drat = Exponent[n, x] + Exponent[d, x];
    unique = True;

    While[ r1 =!= 0,

      t = First[Timing[
      dr1 = Exponent[r1, x]; dt1 = Exponent[t1, x];
      If[ dr1 < df && t1 =!= 1 && drat == dr1 + dt1, unique = False];
      If[ dr1 < df && drat > dr1 + dt1,
        {n, d} = {r1, t1}; drat = dr1 + dt1;
        unique = True;
      ]; 
      ]];
      T[5]+=t;

      {t, q} = Timing[MyPolynomialQuotient[r0, r1, x, p]];
      T[2]+=t;
      {t, {r0, r1, t0, t1}} = Timing[{r1, c[r0 - q r1], t1, c[t0 - q t1]}];
      T[4]+=t;
    ];
    If[ !unique, Throw["ambigous reconstruction result"] ];
    
    {t, {n, d}} = Timing[CoefficientList[#, x]& /@ {n, d}];
    T[3]+=t;
    {t, {n, d}} = Timing[Map[c, {n, d}/Last[d], {2}]];
    T[6]+=t;
    {n, d} = List2Poly[#, x]& /@ {n, d};


    Return[Together[n/d, Modulus -> p]];
  ];

(** ------------------------------------------------------------------- **)
(* matrix operations *)

Strassen[A_, B_, p_Integer] := Module[ {n, n2, AA, BB, M},

  n = First[Dimensions[A]];
  Which[
   n <= 127, PolynomialMod[A.B, p],
   OddQ[n], Strassen[PadRight[A, {n+1, n+1}], PadRight[B, {n+1, n+1}], p][[1;;n, 1;;n]],
   True, n2 = n/2;
    AA[1, 1] = A[[1;;n2,1;;n2]];   AA[1, 2] = A[[1;;n2,n2+1;;n]];
    AA[2, 1] = A[[n2+1;;n,1;;n2]]; AA[2, 2] = A[[n2+1;;n,n2+1;;n]];
    BB[1, 1] = B[[1;;n2,1;;n2]];   BB[1, 2] = B[[1;;n2,n2+1;;n]];
    BB[2, 1] = B[[n2+1;;n,1;;n2]]; BB[2, 2] = B[[n2+1;;n,n2+1;;n]];
    M[1] = Strassen[AA[1,1] + AA[2,2], BB[1,1] + BB[2,2], p];
    M[2] = Strassen[AA[2,1] + AA[2,2], BB[1,1], p];
    M[3] = Strassen[AA[1,1], BB[1,2] - BB[2,2], p];
    M[4] = Strassen[AA[2,2], BB[2,1] - BB[1,1], p];
    M[5] = Strassen[AA[1,1] + AA[1,2], BB[2,2], p];
    M[6] = Strassen[AA[2,1] - AA[1,1], BB[1,1] + BB[1,2], p];
    M[7] = Strassen[AA[1,2] - AA[2,2], BB[2,1] + BB[2,2], p];
    Join[
     Apply[Join, Transpose[{M[1] + M[4] - M[5] + M[7], M[3] + M[5]}, {2,1}], {1}],
     Apply[Join, Transpose[{M[2] + M[4], M[1] - M[2] + M[3] + M[6]}, {2,1}], {1}]
    ]
  ]
];

StrassenInverse[A_, p_Integer] := Module[ {n, B, Inv},

  n = First[Dimensions[A]];
  B = Table[RandomInteger[p], {n}, {n}];
  Strassen[B, StrassenInverseRec[Strassen[A, B, p], p], p]
];
StrassenInverseRec[A_, p_Integer] := Module[ {n, n2, P, C, AA},

  n = First[Dimensions[A]];
  Which[
   n <= 2047, Inverse[A, Modulus -> p],
   OddQ[n], StrassenInverseRec[Append[Append[#, 0]& /@ A, Append[Table[0, {n}], 1]], p][[1;;n, 1;;n]],
   True, n2 = n/2;
    AA[1, 1] = A[[1;;n2,1;;n2]];   AA[1, 2] = A[[1;;n2,n2+1;;n]];
    AA[2, 1] = A[[n2+1;;n,1;;n2]]; AA[2, 2] = A[[n2+1;;n,n2+1;;n]];
    P[1] = StrassenInverseRec[AA[1,1], p];
    P[2] = Strassen[AA[2,1], P[1], p];
    P[3] = Strassen[P[1], AA[1,2], p];
    P[4] = Strassen[AA[2,1], P[3], p];
    P[5] = P[4] - AA[2, 2];
    P[6] = StrassenInverseRec[P[5], p];
    C[1, 2] = Strassen[P[3], P[6], p];
    C[2, 1] = Strassen[P[6], P[2], p];
    C[1, 1] = P[1] - Strassen[P[3], C[2, 1], p];
    C[2, 2] = -P[6];
    Join[ Join[C[1, 1], C[1, 2], 2], Join[C[2, 1], C[2, 2], 2]]
  ]
];

(** ------------------------------------------------------------------- **)
(* Dense Solvers *)

(* linear solver for algebraic number fields; coeffs must be in alpha *)
LinSolveAlg[mat_, alpha_Symbol, minpoly_, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {p, deg, point, sol, lastP, lastSol, mod, lift, i, j, M},

    DebugPrint[0, "Q(alpha)."];

    If[ Variables[minpoly] =!= {alpha}, Throw["illegal minpoly specification."] ];
    If[ !MemberQ[{{alpha}, {}}, Variables[mat]], Throw["matrix contains illegal symbols."] ];

    deg = Exponent[minpoly, alpha];
    p = Developer`$MaxMachineInteger; point = {};
    M = Cancel[#*PolynomialLCM@@(Denominator/@#)]& /@ Together[mat];
    M = Map[HornerForm[PolynomialRemainder[#, minpoly, alpha]]&, M, {2}];

    lastP = 0; lastSol = {};

    While[True,

      (* find a modulus where minpoly is square free and factorizes into linear factors *)
      While[ Head[point] =!= Times || Length[point] < deg, 
        point = PolynomialMod[minpoly/Coefficient[minpoly, alpha, deg], (p=NextPrime[p,-1])];
        point = Factor[point, Modulus -> p];
      ];
      DebugPrint[1, "modulus ", p, " "];
      point = alpha - List@@point;
      mod = {#}& /@ Thread[alpha -> point];

      (* solve in homomorphic image *)
      sol = NullSpace[#, Modulus -> p]& /@ (M /. mod); (* (Hom[M, p, #]& /@ mod); (this is slower) *)
      If[ Min[Length /@ sol] == 0, Return[ {} ]]; (* no solution *)

      (* interpolate *)
      sol = Map[InterpolatingPolynomial[Transpose[{point, #}], alpha, Modulus -> p]&, 
                 Transpose[sol, {3, 1, 2}], {2}];

      (* check for degenerate cases; covers also special handling of first iteration *)
      If[ Length[sol] > Length[lastSol], lastSol = sol; lastP = p; Continue[] ];
      If[ Length[sol] < Length[lastSol], Continue[] ];

      (* combine with previous solution *)
      sol = Apply[Sum[ChineseRemainder[{Coefficient[#1, alpha, i], Coefficient[#2, alpha, i]}, 
                                       {lastP, p}]alpha^i, {i, 0, deg - 1}]&, 
                  Transpose[{lastSol, sol}, {3, 1, 2}], {2}];      

      lastP *= p; lastSol = sol;      

      (* lift and check for termination *)
      lift = {};
      Do[
        AppendTo[lift, Sum[ReconstructQ[Coefficient[#, alpha, j], lastP]alpha^j, {j, 0, deg-1}]& /@ sol[[i]]];
        If[ !MatchQ[Expand[PolynomialRemainder[#, minpoly, alpha]]& /@ (M . Last[lift]), {0...} ],
          lift = {}; Break[] (* not a basis *)
        ];
      , {i, 1, Length[sol]}];

      If[ Length[lift] == Length[sol], Return[lift] ];
    ];
  ];

(* linear solver for Q via chinese remaindering *)
LinSolveQ[mat_?MatrixQ] := Module[ {A}, 

    A = (#/Max[1,GCD@@#])& /@ ((#*LCM@@(Denominator/@#))& /@ mat); 
    LinSolveQ[PolynomialMod[A, #]&]
];
LinSolveQ[A_Symbol] := LinSolveQ[A[#]&];
LinSolveQ[A_Function] := Module[ {p, P, M, sol, lastSol, v, i, k},

    DebugPrint[0, "Q."];

    lastSol = {}; p = Developer`$MaxMachineInteger;
    {t0,t1,t2,t3} = {0,0,0,0};

    While[True, 

       p = NextPrime[p, -1];
       DebugPrint[1, p];
       t3 += First[Timing[M = A[p]]];

       If[ Length[lastSol] > 0, (* lift Z_P to Q and check for termination *)
         sol = {}; 
         t2 += First[Timing[
         Do[
           v = ReconstructQ0[#, P]& /@ lastSol[[i]];
           If[ !MatchQ[PolynomialMod[M . PolynomialMod[v, p], p], {0..}], sol = {}; Break[] ];
           AppendTo[sol, v];
         , {i, 1, Length[lastSol]}];       
         ]];
         If[ Length[sol] > 0, DebugPrint[1, {t0,t1,t2,t3}/.Second->1]; Return[sol] ]; (* STOP. *)
       ];

       (* include next hom. image *)
       t0 += First[Timing[sol = NullSpace[M, Modulus -> p]]];

       Which[ (* initialization and handling of degenerate cases *)
         sol === {}, Return[sol], (* no solution. STOP. *)
         lastSol === {} || Length[sol] < Length[lastSol], (* init. or all prev. primes unlucky *)
           lastSol = sol; P = p, 
         Length[sol] > Length[lastSol], (* this prime is unlucky, discard it. *)
           Print["unlucky prime discarded."]; Null,
         True, (* combine with previus solution *)
           t1 += First[Timing[
           lastSol = Apply[ChineseRemainder[{##}, {P, p}]&, Transpose[{lastSol, sol}, {3, 1, 2}], {2}];
           ]];
           P *= p;
       ];
    ];
];

(* linear solver for Q using hensel lifting *)
LinSolveHensel[mat_] :=
  Module[ {rank, regularCols, pi, piInv, A, At, A0, A0t, B, p, i, s, X},

     p = NextPrime[Developer`$MaxMachineInteger,-1];
     A = (#/Max[1,GCD@@#])& /@ ((#*LCM@@(Denominator/@#))& /@ mat);
     At = Transpose[A];
     A0 = DeleteCases[RowReduce[A, Modulus -> p], {0...}];
     A0t = Transpose[A0];
     rank = Length[A0];
     If[ rank === Length[At], Return[ {} ] ]; (* no solution *)
     (* move regular columns to front *)
     regularCols = Position[A0t, {0..., 1, 0...}];
     pi = Join[regularCols, Complement[{#}& /@ Range[Length[A0t]], regularCols]];
     piInv = Table[0, {Length[pi]}]; Do[piInv[[First[pi[[i]]]]]={i}, {i, 1, Length[pi]}];
     A = Transpose[Extract[At, pi]];
     (* find a row subset with maximal rank *)
     Do[
       If[ MatrixRank[Extract[A, s = Transpose[Subsets[Range[Length[A]], {rank}, {i}]]], Modulus -> p] === rank,
         A = Extract[A, s]; Break[]
       ];
     , {i, 1, Binomial[Length[A], rank]}];
     (* call solver *)
     X = LinSolveHensel[Take[#, rank]& /@ A, Drop[#, rank]& /@ -A];
     Return[Transpose[Extract[Join[X, IdentityMatrix[Length[A0t] - rank]], piInv]]];
  ];

(* hensel solver with linear convergence for inhomogeneous systems over Z *)
LinSolveHensel[A_, B_] := (* A must be a regular square matrix *)
  Module[ {pk, p, Ainv, Aq, Bq, mp, mq, X, k, u},

     DebugPrint[0, "Z."];
     pk = p = NextPrime[Developer`$MaxMachineInteger, -1]; mp[x_] := PolynomialMod[x, p];
     q = NextPrime[p,-1]; mq[x_] := PolynomialMod[x, q]; Aq = mq[A]; Bq = mq[B];

     Ainv = Inverse[A, Modulus -> p];
     u = LinearSolve[A, B, Modulus -> p]; k = 0;

     While[True,

        u -= pk mp[Ainv.Internal`ExactQuotient[A.u - B, pk]]; pk *= p; k++; (* lift *)
        X = Map[ReconstructQ0[#, pk]&, u, {2}]; (* reconstruct *)
        If[ MatchQ[mq[Aq.mq[X] - Bq], {{0...}...}], Break[] ]; (* check *)
     ];

     Return[X];
  ];

(* special solver for univariate polynomials as matrix entries *)
LinSolveUniv[mat_, x_] := LinSolveUniv[mat, 0, x];
LinSolveUniv[mat_, mod_, x_] :=
  Module[ {p, q, mymat, myTimes, myPlus, myPower, d,
           sol, lastSol, roots, x0, n, m, u, v, i, j, M, Mmat, deg, heads, tails, ker},

    p = If[ mod == 0, (* prime seed *)
            Developer`$MaxMachineInteger, 
            Floor[Developer`$MaxMachineInteger/2]
        ]; 

    mymat = mat;

    If[mod === 0, DebugPrint[0, "Q(x)."]];

    lastSol = {}; roots = {}; M = 1;

    While[True,

       p = If[ mod == 0, NextPrime[p, -1], p + 1 ];
       DebugPrint[If[mod === 0, 1, 2], {mod, p}];

       sol = If[mod === 0, 
               LinSolveUniv[mymat, p, x], 
               Quiet[Check[NullSpace[mymat /. x -> p, Modulus -> mod], $Failed, NullSpace::nmod], {NullSpace::nmod}]
             ];

       Which[ (* degenerate cases and initialization *)
         sol === $Failed,
           Print["unlucky prime discarded."]; Continue[],
         sol === {},
           Return[{}],
         mod === 0 && Length[roots] == 0,
           Null, (* nothing *)
         Length[sol] < Length[lastSol] || Length[roots] == 0,
           lastSol = sol; roots = {p}; M = If[mod === 0, p, x - p]; Continue[],
         Length[sol] > Length[lastSol],
           Print["unlucky prime discarded."]; Continue[]
       ];

       If[ mod === 0,
         (* chinese remaindering *)
         If[ Length[roots] == 0,
           lastSol = sol
         ,
           {n, m} = Dimensions[sol];
           lastSol = Table[
                {u, v} = #[[i, j]]& /@ {lastSol, sol}; 
                {u, v} = Poly2List[#, x, Max[1, Exponent[u, x], Exponent[v, x]]]& /@ {u, v};
               (ChineseRemainder[#, {M, p}]& /@ Transpose[{u, v}])
                 . x^Range[0, Length[u] - 1]
             , {i, 1, n}, {j, 1, m}];
         ]
       ,
         (* interpolation step *)
         lastSol = Expand[lastSol + (sol - (lastSol /. x -> p))*(M/Times@@(p-roots)), Modulus -> mod];
       ];
       AppendTo[roots, p];
       If[ mod === 0,
         M *= p,
         M = Expand[M*(x-p), Modulus -> mod]
       ];

       (* reconstruction *)
       (* If[ OddQ[Length[roots]], Continue[] ]; *)
       If[ mod != 0 && Mod[Length[roots], 15] != 0, Continue[] ]; (** ************************** **)

       sol = {};
       Do[
         v = lastSol[[i]];
         AppendTo[sol, 
           If[ mod === 0,
             Sum[ReconstructQ0[Coefficient[#, x, j], M]*x^j, {j, 0, Max[0, Exponent[#, x]]}]& /@ v
           ,
             (* ReconstructRatfunVector[v, mod, x, M, Floor[Length[roots]/2]] *)
             v = Catch[ReconstructRatfun[#, M, x, mod]]& /@ v;
             v = v /. _String -> 1;
             Cancel[v * PolynomialLCM@@Append[Denominator /@ v, Modulus -> mod], Modulus -> mod]
           ]];
         v = Last[sol];
         If[ MatchQ[v, {$Failed...}], sol = {}; Break[] ];         
         p = If[ mod == 0, NextPrime[p, -1], p + 1 ];
         q = If[ mod > 0, mod, NextPrime[p, -1] ];
         If[ !MatchQ[
               PolynomialMod[PolynomialMod[mymat /. x -> p, q] . PolynomialMod[v /. x -> p, q], q],
               {0...}],
           sol = {}; Break[];
         ];
       , {i, 1, Length[lastSol]}];

       If[ Length[sol] > 0, Return[sol] ];
    ];    
  ];

End[]

EndPackage[]
