Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

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

Search the Archive

Re: Re: Re: Finding the Number of Pythagorean Triples below a bound

  • To: mathgroup at smc.vnet.net
  • Subject: [mg68431] Re: [mg68382] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Sat, 5 Aug 2006 03:46:14 -0400 (EDT)
  • References: <eaeqa3$53v$1@smc.vnet.net><200607300848.EAA25171@smc.vnet.net> <eakfgm$rl6$1@smc.vnet.net> <200608020923.FAA28520@smc.vnet.net> <79C36C70-E091-4A82-8EC5-0EDC743D081D@mimuw.edu.pl> <200608031007.GAA15743@smc.vnet.net> <76C59C34-5A55-4318-B0AF-3A572E71B421@mimuw.edu.pl> <C56DFF0B-6F0B-4E64-B72A-D805B3C6C063@mimuw.edu.pl> <7AA8804E-73F2-4D7B-BCAB-633CE33CD1D0@mimuw.edu.pl> <F7520FE0-16DB-4F8F-A39E-DF64749B85F0@mimuw.edu.pl>
  • Sender: owner-wri-mathgroup at wolfram.com

I have good news: the code I just posted can be compiled and then it  
becomes really fast ;-)

countT = Compile[{{m, _Integer}, {k, _Integer}}, Module[
   {i = 0, c2, diff, sdiff},
   Do [
    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
     Do [
       diff = c^2 + k - b^2;
          sdiff = Sqrt[N[diff]];
          If [sdiff >= b && sdiff == Round[sdiff], i++];
       , {b, 1, c2}]];
, {c, 1, 10^m - 1}];
i]]


countT[3,2]+countT[3,-2]//Timing


{1.0032 Second,282}

I will try the case m=4 and m=5 and send you the results. I am not  
promising to do this very soon, just in case ;-)

Andrzej Kozlowski



On 4 Aug 2006, at 12:09, Andrzej Kozlowski wrote:

> Here is the fastest code I have seen so far that works (it is based  
> on the one that Daniel sent you). I have corrected and enhanced   
> and enhanced it. It gives approximate answers for reasons that I  
> explained in earlier postings (use of machine arithmetic to test  
> for perfect integers). Of course it is easy to replace the code by  
> exact code by replacing the numerical test for a perfect square by  
> an exact one.
>
>
> countTriples[m_,k_] := Module[
>   {i=0, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]];
>     Do [
>       diff = c^2+k-b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff>=b&&sdiff==Round[sdiff],i++];
>       , {b,1,c2}]];
> ,{c,1,10^m-1}];
> i]
>
>
> countTriples[3,2]+countTriples[3,-2]//Timing
>
>
> {12.3746 Second,282}
>
> The correct answer is 283.
>
> This code should easily deal with the case m=4 (I have not yet  
> tried it) and I think even m=5 should now be within reach.
>
> Andrzej Kozlowski
>
>
>
> On 4 Aug 2006, at 11:27, Andrzej Kozlowski wrote:
>
>> The "improvement" below which I sent a little earlier is wrong  
>> (even though it returned correct answers). Obviously the point is  
>> that c^2+k  can't be of the form 4n + 3, but there is no reason  
>> why a^2+b^2-k can't be of that form. Since my code does not  
>> explicitly select c it can't make use of this additional  
>> improvement. A different code, which uses explicit choices of  
>> (say) a and c and tests for b being a perfect square could exploit  
>> this fact and perhaps gain extra speed. It should not be difficult  
>> to write such a code along the lines I have been using.
>>
>> Andrzej
>>
>>
>> On 4 Aug 2006, at 10:47, Andrzej Kozlowski wrote:
>>
>>>
>>> On 3 Aug 2006, at 16:42, Andrzej Kozlowski wrote:
>>>
>>>>
>>>> On 3 Aug 2006, at 12:07, Andrzej Kozlowski wrote:
>>>>
>>>>> On 2 Aug 2006, at 20:01, Andrzej Kozlowski wrote:
>>>>>
>>>>>>
>>>>>> On 2 Aug 2006, at 11:23, titus_piezas at yahoo.com wrote:
>>>>>>
>>>>>>> Hello all,
>>>>>>>
>>>>>>> My thanks to Peter and Andrzej, as well as those who privately
>>>>>>> emailed
>>>>>>> me.
>>>>>>>
>>>>>>> To recall, the problem was counting the number of solutions to a
>>>>>>> bivariate polynomial equal to a square,
>>>>>>>
>>>>>>> Poly(a,b) = c^2
>>>>>>>
>>>>>>> One form that interested me was the Pythagorean-like equation:
>>>>>>>
>>>>>>> a^2 + b^2 = c^2 + k
>>>>>>>
>>>>>>> for {a,b} a positive integer, 0<a<=b, and k any small  
>>>>>>> integer. I was
>>>>>>> wondering about the density of solutions to this since I knew  
>>>>>>> in the
>>>>>>> special case of k=0, let S(N) be the number of primitive  
>>>>>>> solutions
>>>>>>> with
>>>>>>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>>>>>>>
>>>>>>> For k a squarefree integer, it is convenient that any  
>>>>>>> solution is
>>>>>>> also
>>>>>>> primitive. I used a simple code that allowed me to find S 
>>>>>>> (10^m) with
>>>>>>> m=1,2,3 for small values of k (for m=4 took my code more than  
>>>>>>> 30 mins
>>>>>>> so I aborted it). The data is given below:
>>>>>>>
>>>>>>> Note: Values are total S(N) for *both* k & -k:
>>>>>>>
>>>>>>> k = 2
>>>>>>> S(N) = 4, 30, 283
>>>>>>>
>>>>>>> k = 3
>>>>>>> S(N) = 3, 41, 410
>>>>>>>
>>>>>>> k = 5
>>>>>>> S(N) = 3, 43, 426
>>>>>>>
>>>>>>> k = 6
>>>>>>> S(N) = 3, 36, 351
>>>>>>>
>>>>>>> Question: Does S(N)/N for these also converge? For example,  
>>>>>>> for the
>>>>>>> particular case of k = -6, we have
>>>>>>>
>>>>>>> S(N) = 2, 20, 202
>>>>>>>
>>>>>>> which looks suspiciously like the ratio might be converging.
>>>>>>>
>>>>>>> Anybody know of a code for this that can find m=4,5,6 in a  
>>>>>>> reasonable
>>>>>>> amount of time?
>>>>>>>
>>>>>>>
>>>>>>> Yours,
>>>>>>>
>>>>>>> Titus
>>>>>>>
>>>>>>
>>>>>>
>>>>>> Here is a piece code which utilises the ideas I have described in
>>>>>> my previous posts:
>>>>>>
>>>>>> ls = Prime /@ Range[3, 10];
>>>>>>
>>>>>> test[n_] :=
>>>>>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>>>> Integers]
>>>>>>
>>>>>> f[P_, k_] := Sum[If[(w =
>>>>>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>>>>>           Floor[Sqrt[P^2 - a^2]]}]
>>>>>>
>>>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>>>>
>>>>>> We can easily confirm the results of your computations,.e.g.  
>>>>>> for k=2.
>>>>>>
>>>>>>
>>>>>> Table[g[i,2],{i,3}]
>>>>>>
>>>>>>
>>>>>> {4,30,283}
>>>>>>
>>>>>>
>>>>>> Since you have not revealed the "simple code" you have used it is
>>>>>> hard to tell if the above one is any better. It is however,
>>>>>> certainly capable of solving the problem for m=4:
>>>>>>
>>>>>>
>>>>>> g[4,2]//Timing
>>>>>>
>>>>>>
>>>>>> {4779.39 Second,2763}
>>>>>>
>>>>>> The time it took on my 1 Gigahertz PowerBook was over 70 minutes,
>>>>>> which is longer than you thought "reasonable", so I am still not
>>>>>> sure if this is any improvement on what you already have. The  
>>>>>> time
>>>>>> complexity of this algorithm seems somewhat larger than  
>>>>>> exponential
>>>>>> so I would expect that it will take about 6 hours to deal with  
>>>>>> n=5
>>>>>> on my computer, and perhaps 2 weeks to deal with n=6.
>>>>>>
>>>>>> Andrzej Kozlowski
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>> I mistakenly copied and pasted a wrong (earlier) definition of f.
>>>>> Here is the correct one:
>>>>>
>>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>>>      1 && (w = a^2 +
>>>>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>>>>>
>>>>> The definition of g is as before:
>>>>>
>>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>>>
>>>>> Andrzej Kozlowski
>>>>>
>>>>
>>>>
>>>> Below is a faster version of the above code. (It owes a  
>>>> significant improvement to Daniel Lichtblau, which I borrowed  
>>>> from him without his knowledge ;-))
>>>>
>>>> test[n_] :=
>>>>     JacobiSymbol[n,
>>>>   Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
>>>> Integers]
>>>>
>>>>
>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>>     1 && (w = a^2 +
>>>>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>>>> +k]]}]]
>>>>
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>>
>>>> The improvement is in the upper bound on a in the sum. Since a  
>>>> is the smaller of the two squares whose sum is equal to P^2+k it  
>>>> can't be larger than  Floor[Sqrt[(P^2+k)/2]].
>>>>
>>>> Note that you can improve the performance by loosing some  
>>>> accuracy if you use a cruder test for a perfect square:
>>>>
>>>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>>>>    ]
>>>>
>>>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>>     1 && (w = a^2 +
>>>>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>>>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>>>> +k]]}]]
>>>>
>>>>
>>>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
>>>>
>>>>
>>>> Let's compare the two cases.
>>>>
>>>> In[7]:=
>>>> g[3,2]//Timing
>>>>
>>>> Out[7]=
>>>> {89.554 Second,283}
>>>>
>>>> In[8]:=
>>>> g1[3,2]//Timing
>>>>
>>>> Out[8]=
>>>> {37.376 Second,283}
>>>>
>>>> So we see that we get the same answer and the second approach is  
>>>> considerably faster. However:
>>>>
>>>> In[9]:=
>>>> g[1,6]//Timing
>>>>
>>>> Out[9]=
>>>> {0.008863 Second,3}
>>>>
>>>> In[10]:=
>>>> g1[1,6]//Timing
>>>>
>>>> Out[10]=
>>>> {0.005429 Second,5}
>>>>
>>>>
>>>> The correct answer is 3 (returned by the first method). The  
>>>> faster method found two false solutions. This should not matter  
>>>> if you are interested only in approximate answers (as you seem  
>>>> to be) but it is worth keeping in mind.
>>>>
>>>> Andrzej Kozlowski
>>>
>>>
>>> I have noticed one more obvious improvement. We can replace test by:
>>>
>>> test[n_] :=
>>>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>>>   Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
>>> Integers]
>>>
>>> and test1 by
>>>
>>> test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round[w]
>>>    ]
>>>
>>> We are simply making use of the easily to prove fact that an  
>>> integer of the form 4 k + 1 can never be the sum of two squares.  
>>> There is a noticeable improvement in the performance of g:
>>>
>>>
>>> g[3,2]//Timing
>>>
>>> {58.0786 Second,283}
>>>
>>> However, the performance of g1 seems to actually decline slightly:
>>>
>>>
>>> g1[3,2]//Timing
>>>
>>> {40.8776 Second,283}
>>>
>>>
>>> However, there are fewer cases of "false solutions":
>>>
>>> In[22]:=
>>> g1[1,6]//Timing
>>>
>>> Out[22]=
>>> {0.006694 Second,4}
>>>
>>>
>>> I am still not sure what is the most efficient use of  
>>> JacobiSymbol in this kind of problems. In the first version of my  
>>> code I used a test involving the first few odd primes, afterwards  
>>> I switched to using just one random prime in taken form a certain  
>>> range. Evaluation of JacobiSympol[m,p] != -1 is certainly much  
>>> faster than that of Element[Sqrt[m],Integers], but it is not  
>>> clear what is the most efficient number of primes to use, which  
>>> are the best primes and whether it is better to choose them at  
>>> random or just use a fixed selection.
>>> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster,  
>>> but it will sometimes produce "false squares".
>>>
>>> Andrzej Kozlowski
>>>
>>>
>>
>


  • Prev by Date: Re: Re: Re: Re: "No more memory available" -- a recurring problem
  • Next by Date: Re: Re: Re: Re: "No more memory available" -- a recurring problem
  • Previous by thread: Re: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Next by thread: Re: Re: Re: Finding the Number of Pythagorean Triples below a bound