myw = 1 + z^(1/4) + z^(1/2) + z^(3/4) + z + z^(5/4)
fExponents = Exponent[myw, z, List]
fdegree = Length[fExponents];
matrixFlag = False;
While[! matrixFlag && fdegree < 10,
{
pdegree = fdegree;
Quiet[Table[
Subscript[s, n] =., {k, 0, fdegree}, {n, 20 k, 20 k + pdegree}]];
theCoefficients =
Table[Subscript[s, n], {k, 0, fdegree}, {n, 20 k, 20 k + pdegree}];
thePolyTerms =
Table[theCoefficients[[i, j]] z^(j - 1), {i, 1, fdegree + 1}, {j,
1, pdegree + 1}];
thePolys = (Plus @@ # &) /@ thePolyTerms;
theFunctionForm =
Plus @@ Flatten[
Table[thePolys[[i]] w^(i - 1), {i, 1, fdegree + 1}]];
expandedForm = Expand[theFunctionForm /. w -> myw];
elist = Exponent[expandedForm, z, List];
theEquations =
Table[Coefficient[expandedForm, z, elist[[i]]], {i, 1,
Length[elist]}];
theSolution =
Solve[Table[
Coefficient[expandedForm, z, elist[[i]]] == 0, {i, 1,
Length[elist]}], Flatten[theCoefficients], Method -> "Reduce"];
val1 = theSolution // Flatten;
val2 = Association[val1];
clist = Flatten[theCoefficients];
iList = Intersection[Flatten[theCoefficients], Keys[val2]];
theNormalForm =
Normal[CoefficientArrays[#, Flatten[theCoefficients]]] & /@
theEquations;
theCMatrix = #[[2]] & /@ theNormalForm;
theMRank = MatrixRank[theCMatrix];
Print["number of equations", Length[theEquations]];
Print["number of variables: ", Length[Flatten[theCoefficients]]];
Print["Matrix Rank: ", MatrixRank[theCMatrix]];
If[theMRank < Length[Flatten[theCoefficients]],
{
Print[
"Rank less than number of variables. Computing function . . . \
"];
theFreeVals = Complement[clist, iList];
Print["free vars: ", theFreeVals];
For[i = 1, i <= Length[theFreeVals], i++,
Subscript[theFreeVals[[i, 1]], theFreeVals[[i, 2]]] = 1;
];
Print["Function: ",
polyForm[(theFunctionForm /. theSolution)[[1]], w]];
matrixFlag = True;
}
,
fdegree++;
];
}
];