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: [mg68432] 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:16 -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> <B45555BF-3093-4487-A3EA-0981783CA75A@mimuw.edu.pl>
  • Sender: owner-wri-mathgroup at wolfram.com

My latest code has just managed:

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


{93.9638 Second,2762}


That is less than two minutes on a 1 gigahertz computer. The correct  
answer is actually 2763 (by my earlier computation using the exact  
test) so we have lost one solution but gained more than 50 fold  
improvement in performance!
The case m=5 is now certainly feasible, although I am not sure if I  
wish my not very powerful PowerBook to be occupied for so long, as I  
need to use Mathematica for other tasks. Perhaps I can now leave this  
to others.

Andrzej Kozlowski

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

> 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: Finding the Number of Pythagorean Triples below a bound
  • 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