MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: How to do quickest

  • To: mathgroup at smc.vnet.net
  • Subject: [mg123093] Re: How to do quickest
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Wed, 23 Nov 2011 07:05:53 -0500 (EST)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • References: <201111210929.EAA14830@smc.vnet.net>
  • Reply-to: drmajorbob at yahoo.com

The Block construction causes Round[Sqrt[x^3]] to be recomputed each time  
y is mentioned, or so it would seem, and timing bears that out:

Collect[poly =
   Eliminate[{4*m^2 + 6*m*n + n^2 ==
       x, (19*m^2 + 9*m*n + n^2)*Sqrt[m^2 + n^2] == y}, {n}] /.
    Equal -> Subtract, m]

3645 m^12 - 2916 m^10 x + x^6 - 2 x^3 y^2 + y^4 +
  m^6 (270 x^3 - 270 y^2)

Block[{y = Round[Sqrt[x^3]]},
   Reap[Do[x^3 - y^2 != 0 &&
       Or[OddQ[x^6 - 2*x^3*y^2 + y^4],
        Mod[x^6 - 2*x^3*y^2 + y^4, 4] == 0] && !
        IrreduciblePolynomialQ[poly] && Sow@x, {x, 2,
       10^4}]][[2]]] // Timing

{10.98, {{1942, 2878, 3862, 6100, 8380}}}

First@Last@
    Reap[Do[y = Round[Sqrt[x^3]];
      x^3 - y^2 != 0 &&
       Or[OddQ[x^6 - 2*x^3*y^2 + y^4],
        Mod[x^6 - 2*x^3*y^2 + y^4, 4] == 0] && !
        IrreduciblePolynomialQ[poly] && Sow@x, {x, 2, 10^4}]] // Timing

{8.36221, {1942, 2878, 3862, 6100, 8380}}

This is slightly faster (surprisingly?):

First@Last@
    Reap[Do[y = Round[Sqrt[x^3]]; ! IrreduciblePolynomialQ[poly] &&
       x^3 - y^2 != 0 &&
       Or[OddQ[x^6 - 2*x^3*y^2 + y^4],
        Mod[x^6 - 2*x^3*y^2 + y^4, 4] == 0] && Sow@x, {x, 2,
       10^4}]] // Timing

{8.26819, {1942, 2878, 3862, 6100, 8380}}

And that suggests the Eisenstein criterion is unhelpful after all, so:

First@Last@
    Reap[Do[y = Round[Sqrt[x^3]]; ! IrreduciblePolynomialQ[poly] &&
       x^3 - y^2 != 0 && Sow@x, {x, 2, 10^6}]] // Timing

{978.795, {1942, 2878, 3862, 6100, 8380, 11512, 15448, 18694, 31228,
   93844, 111382, 117118, 129910, 143950, 186145, 210025, 375376,
   445528, 468472, 575800, 844596, 950026}}

Bobby

On Tue, 22 Nov 2011 12:40:48 -0600, Andrzej Kozlowski <akoz at mimuw.edu.pl>  
wrote:

> I think the correct use of Eisenstein's criterion is this:
>
>
> Block[{y = Round[Sqrt[x^3]]},
>   Reap[Table[
>      If[x^3 - y^2 != 0 &&
>        Not[Mod[x^6 - 2*x^3*y^2 + y^4, 2] == 0 &&
>          Mod[x^6 - 2*x^3*y^2 + y^4, 4] != 0] && !
>         IrreduciblePolynomialQ[poly], Sow[{x, y}]], {x, 2,
>       1000000}]][[2]]] // Timing
>
> This is because 3645 is not divisible by 2, but  2916 and 270 both are.  
> So if the term x^6-2 x^3 y^2+y^4 is divisible by 2 but not divisible by  
> 4, the polynomial is irreducible and there is no need to test it  
> further. Only when this isn't the case, we need to use  
> IrreduciblePolynomialQ.
>
> I am not sure if adding this will speed up the code. I only have one  
> copy of Mathematica and I need to use it so I can't afford the time to  
> run these programs now, when I am a hurry.
>
> Andrzej
>
>
>
> On 22 Nov 2011, at 19:30, Andrzej Kozlowski wrote:
>
>> My memory of Eisenstein's criterion was wrong (also, I was too much in  
>> a hurry to look it up).  Rather than correcting it I got rid of it  
>> altogether since I think Mathematica probably uses it anyway. I then  
>> get:
>>
>> In[1]:= Collect[
>> poly = Eliminate[{4*m^2 + 6*m*n + n^2 ==
>>      x, (19*m^2 + 9*m*n + n^2)*Sqrt[m^2 + n^2] == y}, {n}] /.
>>   Equal -> Subtract, m]
>>
>> Out[1]= 3645 m^12-2916 m^10 x+m^6 (270 x^3-270 y^2)+x^6-2 x^3 y^2+y^4
>>
>> In[2]:= Block[{y = Round[Sqrt[x^3]]},
>>  Reap[Table[
>>     If[x^3 - y^2 != 0 && Not[IrreduciblePolynomialQ[poly]],
>>      Sow[{x, y}]], {x, 2, 1000000}]][[2]]] // Timing
>>
>>
>> Out[2]=  
>> {1089.54,({1942,85580}	{2878,154396}	{3862,240004}	{6100,476425}	{8380,767125}	{11512,1235168}	{15448,1920032}	{18694,2555956}	{31228,5518439}	{93844,28748141}	{111382,37172564}	{117118,40080716}	{129910,46823500}	{143950,54615700}	{186145,80311375}	{210025,96251275}	{375376,229985128}	{445528,297380512}	{468472,320645728}	{575800,436925600}	{844596,776199807}	{950026,925983476}
>>
>> )}
>>
>> This gets all the numbers but is much slower (I guess it will be better  
>> to add the Eisenstein criterion after all, but of course, in correct  
>> form).
>>
>> Andrzej Kozlowski
>>
>>
>>
>>
>> On 22 Nov 2011, at 18:02, Artur wrote:
>>
>>> Dear Andrzej,
>>> Your procedure omiited some points
>>> 6100
>>> 8380 Best wishes
>>> Artur
>>>
>>> W dniu 2011-11-22 13:23, Andrzej Kozlowski pisze:
>>>> On 22 Nov 2011, at 10:07, Andrzej Kozlowski wrote:
>>>>
>>>>
>>>>> On 22 Nov 2011, at 10:06, Andrzej Kozlowski wrote:
>>>>>
>>>>>
>>>>>> On 21 Nov 2011, at 10:29, Artur wrote:
>>>>>>
>>>>>>
>>>>>>> Dear Mathematica Gurus,
>>>>>>> How to do quickest following procedure (which is very slowly):
>>>>>>>
>>>>>>> qq = {}; Do[y = Round[Sqrt[x^3]];
>>>>>>> If[(x^3 - y^2) != 0,
>>>>>>> kk = m /. Solve[{4 m^2 + 6 m n + n^2 ==
>>>>>>>     x, (19 m^2 + 9 m n + n^2) Sqrt[m^2 + n^2] == y}, {m, n}][[1]];
>>>>>>> ll = CoefficientList[MinimalPolynomial[kk][[1]], #1];
>>>>>>> lll = Length[ll];
>>>>>>> If[lll < 12, Print[{x/(x^3 - y^2)^2, kk, x, y, x^3 - y^2}];
>>>>>>> If[Length[ll] == 3, Print[{kk, x, y}]]]], {x, 2, 1000000}];
>>>>>>> qq
>>>>>>>
>>>>>>>
>>>>>>> (*Best wishes Artur*)
>>>>>>>
>>>>>>>
>>>>>> I think it would be better to send not only the code but also the  
>>>>>> mathematical problem, as there may be a way to do it in a different  
>>>>>> way. Unless I am misunderstanding something, what you are trying to  
>>>>>> do is the same as this:
>>>>>>
>>>>>> In[31]:= Block[{y = Round[Sqrt[x^3]]},
>>>>>> Reap[Table[
>>>>>>  If[(x^3 - y^2) != 0 && Not[IrreduciblePolynomialQ[poly]],
>>>>>>   Sow[{x, y}]], {x, 2, 1000000}]][[2]]] // Timing
>>>>>>
>>>>>> Out[31]= {721.327,{}}
>>>>>>
>>>>>> This ought to be a lot faster than your code, but I have not tried  
>>>>>> to run yours to the end. Also, it is possible that using the  
>>>>>> Eisenstein Test explicitly may be somewhat faster:
>>>>>>
>>>>>> Block[{y = Round[Sqrt[x^3]]},
>>>>>> Reap[Table[
>>>>>> If[x^3 - y^2 != 0 && Mod[x^6 - 2*x^3*y^2 + y^4, 4] == 0 &&
>>>>>>           ! IrreduciblePolynomialQ[poly], Sow[{x, y}]], {x, 2,
>>>>>>  1000000}]][[2]]]
>>>>>>
>>>>>> {}
>>>>>>
>>>>>> but I forgot to use Timing and don't want to wait again,  
>>>>>> particularly that the answer is the empty set.
>>>>>>
>>>>>> Andrzej Kozlowski
>>>>>>
>>>>> I forgot to include the definition of poly:
>>>>>
>>>>> Collect[poly = Eliminate[{4*m^2 + 6*m*n + n^2 == x,
>>>>>    (19*m^2 + 9*m*n + n^2)*Sqrt[m^2 + n^2] == y}, {n}] /. Equal ->  
>>>>> Subtract, m]
>>>>>
>>>>> 3645*m^12 - 2916*m^10*x + m^6*(270*x^3 - 270*y^2) + x^6 -
>>>>> 2*x^3*y^2 + y^4
>>>>>
>>>>> Andrzej Kozlowski
>>>>>
>>>>
>>>> Strange but I run this code with a fresh kernel and got the following  
>>>> answers:
>>>>
>>>> In[1]:=  
>>>> Collect[poly=Eliminate[{4*m^2+6*m*n+n^2==x,(19*m^2+9*m*n+n^2)*Sqrt[m^2+n^2]==y},{n}]/.Equal->Subtract,m]
>>>> Out[1]= 3645 m^12-2916 m^10 x+m^6 (270 x^3-270 y^2)+x^6-2 x^3 y^2+y^4
>>>>
>>>> In[2]:=  
>>>> Block[{y=Round[Sqrt[x^3]]},Reap[Table[If[x^3-y^2!=0&&Mod[x^6-2*x^3*y^2+y^4,4]==0&&!IrreduciblePolynomialQ[poly],Sow[{x,y}]],{x,2,1000000}]][[2]]]//Timing
>>>>
>>>> Out[2]=  
>>>> {766.05,({1942,85580}	{2878,154396}	{3862,240004}	{11512,1235168}	{15448,1920032}	{18694,2555956}	{111382,37172564}	{117118,40080716}	{129910,46823500}	{143950,54615700}	{186145,80311375}	{210025,96251275}	{375376,229985128}	{445528,297380512}	{468472,320645728}	{575800,436925600}	{950026,925983476}
>>>>
>>>> )}
>>>>
>>>>
>>>> I tested the first one and it does seem to be a solution to your  
>>>> problem.
>>>>
>>>> {x, y} = {950026, 925983476};
>>>>
>>>> y == Round[Sqrt[x^3]]
>>>>
>>>> True
>>>>
>>>> x^3 - y^2 != 0
>>>>
>>>> True
>>>>
>>>> kk =
>>>> m /. Solve[{4 m^2 + 6 m n + n^2 ==
>>>>     x, (19 m^2 + 9 m n + n^2) Sqrt[m^2 + n^2] == y}, {m, n}][[1]]
>>>>
>>>> Out[12]= -Sqrt[-(198/5)-(44 I Sqrt[11])/5]
>>>>
>>>> ll = CoefficientList[MinimalPolynomial[kk][[1]], #1];
>>>>
>>>> Length[ll]
>>>>
>>>> 5
>>>>
>>>> I don't know why I got no answers the first time round, perhaps one  
>>>> of the variables had values assigned.
>>>>
>>>> Andrzej
>>>>
>>>>
>>>>
>>>>
>>
>


-- 
DrMajorBob at yahoo.com



  • Prev by Date: Re: Problem displaying user-created .cdf files
  • Next by Date: Re: How to do quickest
  • Previous by thread: Re: How to do quickest
  • Next by thread: Re: How to do quickest