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: [mg68426] 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:02 -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>
  • Sender: owner-wri-mathgroup at wolfram.com

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