On Tue, Sep 20, 2011 at 3:34 AM, Bill Gosper <billgosper@gmail.com> wrote:
RootAp1[z_, n_] := Block[{r = RootApproximant[z, n]}, If[Log[Abs[(MinimalPolynomial[r] /. Plus -> Times)@1]] > 9*Precision[z], 1, r]]
[This replaces previous version, which had one of those delicious bugs that disappears when you insert Print statements. Thank you for not playing with it and sending me hate mail.-]
Ssexy[q_]:=Catch[Block[{rts,foos=1,v=Variables[q][[1]],Q},rts=#[[1,2]]&/@NSolve[q==0,WorkingPrecision->Ceiling[Log[Abs[Discriminant[q,v]]]]];If[FreeQ[Q=Solve[0==q],Root],Q, Do[foos*=Select[v*Factor[q,Extension->{RootAp1[rts[[j]]+rts[[k]],3],RootAp1[rts[[j]]*rts[[k]],3]}],Exponent[#,v]==2&];If[Exponent[foos,v]==6,Throw[ToRadicals[Solve[0==foos]]]],{j,5},{k,j+1,6}];foos=Factor[q]; Do[If[foos=!=(foos=Factor[q,Extension->(RootAp1[#,2]&/@CoefficientList[(#-rts[[i]])*(#-rts[[j]])*(#-rts[[k]]),#])]),Throw[Solve[foos==0]]],{i,4},{j,i+1,5},{k,j+1,6}]];Q]]
This is seriously untested, longer, faster on most cases, but several minutes on Ssexy[2925951033851274156588135512485165232256823853056 - 2697817290737324449800236848640264992467435520 #1 + 9932351343021963689693473396732415411860992 #1^2 - 1881654619801628210127689611068977299937 #1^3 + 11450425009897563891465337536606118710 #1^4 + 3847649781964086608961673413540069 #1^5 + 378818692265664781682717625943 #1^6]
The slowness was due to WorkingPrecision failing to track the problem size. Now takes .5 sec. That's right, decupling the precision sped it up 200 fold! But now In[51]:= Timing[Ssexy[1-2 x^2-2 x^3+x^6]]
Out[51]= {3.64753,{{x->1/3 (27/2-3/2 Sqrt[3 (27-8 Sqrt[2])])^(1/3)+(1/2 (9+Sqrt[3 (27-8 Sqrt[2])]))^(1/3)/3^(2/3)},{x->-(1/6) (1+I Sqrt[3]) (27/2-3/2 Sqrt[3 (27-8 Sqrt[2])])^(1/3)-((1-I Sqrt[3]) (1/2 (9+Sqrt[3 (27-8 Sqrt[2])]))^(1/3))/(2 3^(2/3))},{x->-(1/6) (1-I Sqrt[3]) (27/2-3/2 Sqrt[3 (27-8 Sqrt[2])])^(1/3)-((1+I Sqrt[3]) (1/2 (9+Sqrt[3 (27-8 Sqrt[2])]))^(1/3))/(2 3^(2/3))},{x->(1/2 (9+Sqrt[3 (27+8 Sqrt[2])]))^(1/3)/3^(2/3)-2^(5/6)/(3 (9+Sqrt[3 (27+8 Sqrt[2])]))^(1/3)},{x->-(((1+I Sqrt[3]) (1/2 (9+Sqrt[3 (27+8 Sqrt[2])]))^(1/3))/(2 3^(2/3)))+(1-I Sqrt[3])/(2^(1/6) (3 (9+Sqrt[3 (27+8 Sqrt[2])]))^(1/3))},{x->-(((1-I Sqrt[3]) (1/2 (9+Sqrt[3 (27+8 Sqrt[2])]))^(1/3))/(2 3^(2/3)))+(1+I Sqrt[3])/(2^(1/6) (3 (9+Sqrt[3 (27+8 Sqrt[2])]))^(1/3))}}} takes seven times longer than the hairy case! (Apparently the hairy case is at least somewhat "promiscuous".) Interlacing the 6choose2 quadratic(cubic) tests with the 6choose3 cubic(quadratic) tests ought speed up the unlucky cases more than slowing down the lucky ones. Anyway, play with this new Ssexy if you're that sort of person. It would be cool if a function this short solves all solvable (numerically explicit) sextics. --rwg
On Mon, Sep 19, 2011 at 5:41 AM, Tito Piezas <tpiezas@gmail.com> wrote:
Dear Bill,
In addition to that problematic sextic you found, can you also test the solver with a "promiscuous" sextic? It can be solved in 4 ways. Define,
F(x) = 1 - 5 x + 9 x^3 - 2 x^4 - 3 x^5 + x^6 = 0
then F(x) decomposes as quadratics in 3 ways,
1. F(x) = Resultant[x^2 + a x + (a^2 + a - 3), a^3 + 3a^2 - a - 5, a] 2. F(x) = Resultant[x^2 + b x + (b^2 + 2b - 2), b^3 + 3b^2 - b - 1, b] 3. F(x) = Resultant[x^2 - x + c, c^3 + 5c^2 + 5c - 1, c]
plus it factors into cubics over the extension, Sqrt[37]],
4. Factor[F(x), Extension -> Sqrt[37]]
I'm curious which of the four the "ssexy" program will find.
Sincerely,
- Tito
None! As written, Ssexy first tries Mma Solve, which immediately scores with In[10]:= Decompose[1 - 5 x + 9 x^3 - 2 x^4 - 3 x^5 + x^6, x]
Out[10]= {1 + 5 x - 5 x^2 + x^3, -x + x^2}
which I guess is your case 3. Experimentally disabling the call to Solve, it scored nine times looking for three quadratics with cubic extensions, and then found the sqrt 37 split very promptly, which gave massively simpler expressions. Maybe I should switch the order of checking. Also, if the three quadratics case always hits several (9?) times or never, it might be possible to give up after 1/severalth of the search. --rwg [...]