(* LinearSystemSolver. by Manuel Kauers *)

LinearSystemSolver`Private`Version = "0.5 2007-11-01";

BeginPackage["LinearSystemSolver`", {"DiscreteMath`Combinatorica`", "Algebra`PolynomialPowerMod`", "NumberTheory`NumberTheoryFunctions`", "Algebra`Horner`"}]

LinSolve::usage = "";
LinSolveQ::usage = "";
LinSolveUniv::usage = "";

Begin["`Private`"]

PrintLevel = 0;

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 *)

(**
 * 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, n]]*2)/2 ];

(**
 * 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},

    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];

    While[ r1 =!= 0,

      t = First[Timing[
      dr1 = Exponent[r1, x]; dt1 = Exponent[t1, x];
      If[ dr1 < df && drat > dr1 + dt1,
        {n, d} = {r1, t1}; drat = dr1 + dt1;
      ]; ]];
      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;
    ];
    
    {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[n/d];
  ];

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

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

    DebugPrint[0, "algebraic solver."];

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

    deg = Exponent[minpoly, alpha];
    p = 2^28; point = {};
    M = Cancel[#*PolynomialLCM@@(Denominator/@#)]& /@ Together[mat];
    M = Map[Horner[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=PreviousPrime[p])];
        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_] := Module[ {p, P, A, M, sol, lastSol, v, i, k},

    DebugPrint[0, "Q."];

    A = (#/Max[1,GCD@@#])& /@ ((#*LCM@@(Denominator/@#))& /@ mat); 

    p = PreviousPrime[Developer`$MaxMachineInteger];
    M = PolynomialMod[A, p];
    lastSol = {};

    {t0,t1,t2,t3} = {0,0,0,0};

    While[True, 

       DebugPrint[1, p];

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

       Which[ (* initialization and handling of degenerate cases *)
         sol === {}, Return[sol], (* no solution *)
         lastSol === {} || Length[sol] < Length[lastSol], (* init. or all prev. primes unlucky *)
           lastSol = sol; P = p; p = PreviousPrime[p]; M = PolynomialMod[A, p]; Continue[],
         Length[sol] > Length[lastSol], (* this prime unlucky *)
           p = PreviousPrime[p]; M = PolynomialMod[A, p]; Continue[]
       ];

       (* combination *)
       t1 += First[Timing[
       lastSol = Apply[ChineseRemainder[{##}, {P, p}]&, Transpose[{lastSol, sol}, {3, 1, 2}], {2}];
       ]];
       P *= p;

       (* reconstruct & check *)
       sol = {}; p = PreviousPrime[p]; M = PolynomialMod[A, p];
       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}/Second]; Return[sol] ]; (* done. *)
    ];
];

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

     p = PreviousPrime[Developer`$MaxMachineInteger];
     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, "HenselZ."];
     pk = p = PreviousPrime[Developer`$MaxMachineInteger]; mp[x_] := PolynomialMod[x, p];
     q = PreviousPrime[p]; 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_, 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, "univ."]];

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

    While[True,

       p = PreviousPrime[p];
       DebugPrint[If[mod === 0, 1, 2], {mod, p}];

       sol = If[mod === 0, 
               LinSolveUniv[mymat, p, x], 
               NullSpace[mymat /. x -> p, Modulus -> mod]
             ];

       Which[ (* degenerate cases and initialization *)
         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],
           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[ReconstructQ[Coefficient[#, x, j], M]*x^j, {j, 0, Max[1, Exponent[#, x]]}]& /@ v
           ,
             (* ReconstructRatfunVector[v, mod, x, M, Floor[Length[roots]/2]] *)
             v = ReconstructRatfun[#, M, x, mod]& /@ v;
             Cancel[v * PolynomialLCM@@Append[Denominator /@ v, Modulus -> mod], Modulus -> mod]
           ]];
         v = Last[sol];
         If[ MatchQ[v, {$Failed...}], sol = {}; Break[] ];         
         p = PreviousPrime[p];
         q = If[ mod > 0, mod, PreviousPrime[p] ];
         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[]