Skip to content

Instantly share code, notes, and snippets.

/gist:478bd74a511575deba83 Secret
Created Feb 12, 2016

Embed
What would you like to do?
Clear[idAlgSum]
Options[idAlgSum] = {StopBeforeGB -> False};
idAlgSum[z_, degTry_Integer?Positive, coefSubstitutions_List,
OptionsPattern[]] :=
Module[{r, \[Delta]z, p, \[Delta]p, q, \[Delta]q, t, u, conds, vars,
sols, solsP, solsQ},
r = MinimalPolynomial[z];
\[Delta]z = Length[CoefficientList[r[t], t]] - 1;
\[Delta]p = degTry;
\[Delta]q = Ceiling[\[Delta]z/\[Delta]p];
Print["Degrees: ", {\[Delta]z, \[Delta]p, \[Delta]q}];
p[t_] :=
Subscript[a, 0] + Sum[Subscript[a, k] t^k, {k, \[Delta]p}];
q[t_] :=
Subscript[b, 0] + Sum[Subscript[b, k] t^k, {k, \[Delta]q}];
vars = Table[Subscript[a, k], {k, 0, \[Delta]p}]~Join~
Table[Subscript[b, k], {k, 0, \[Delta]q}];
Print["Variables: ", vars];
conds = CoefficientList[Resultant[p[t], q[u - t], t] - r[u], u];
conds = conds /. coefSubstitutions;
Print["Conditions (Length ", Length@conds, "): ", conds // Column];
If[OptionValue[StopBeforeGB], Return[conds]];
conds = GroebnerBasis[conds, vars];
Print["Basis (Length ", Length[conds], "): ", conds // Column];
sols = Solve[Thread[conds == 0], vars];
Print["Solutions: ", sols // TableForm];
sols = Cases[#, Except[Null]] &@Table[
Module[{coefValues = vars /. Join[sol, coefSubstitutions],
whichIntegers},
whichIntegers =
NumericQ[#] && Abs[# - Round[#]] < 1.*^-8 + 1.*^-8 Abs[#] & /@
coefValues;
(*Print[coefValues];Print[whichIntegers];*)
If[And @@ whichIntegers
, sol
, Print["Not an integer solution: ", coefValues];]
]
, {sol, sols}];
Print["Solutions (reduced): ", sols // TableForm];
sols = Join[#, coefSubstitutions] & /@ sols;
Cases[#, Except[Null]] &@Flatten[#, 2] &@Table[
Module[{
solP = t /. Solve[p[t] == 0 /. sol, t],
solQ = t /. Solve[q[t] == 0 /. sol, t]
},
Table[
If[N[Abs[(rootP + rootQ)/z - 1], 400] < 1.*^-200,
{rootP, rootQ}]
, {rootP, solP}, {rootQ, solQ}]
], {sol, sols}]
]
idAlgSum[Sqrt[2] + Sqrt[3], 2, {Subscript[b, 2] -> 1,
Subscript[a, 2] -> 1, Subscript[b, 1] -> 0}]
idAlgSum[Sqrt[2] + Sqrt[3] + Sqrt[5], 2, {Subscript[b, 1] -> 0,
Subscript[a, 2] -> 1}]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.