(* LinearSystemSolver. by Manuel Kauers *) LinearSystemSolver`Private`Version = "0.15 2009-12-17"; (* for MMA6 from version 0.6 on *) (* for MMA7 from version 0.9 on *) SetAttributes[HornerForm, Listable]; BeginPackage["LinearSystemSolver`", {"Combinatorica`"}] 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`"] 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_] := 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[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_] := Module[ {u, myvars, data, den, terms, lt, res, deg, p}, 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, 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])]; data = Transpose[Table[PolynomialMod[data[[All, i]], points[[i]]], {i, 1, Length[points]}]]; (* reconstruct one by one *) data = Table[ReconstructQ0[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, 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[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, 2 (* 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[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 = 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[]