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: [mg68415] Re: [mg68382] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Fri, 4 Aug 2006 03:59:29 -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>
  • Sender: owner-wri-mathgroup at wolfram.com

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


  • Prev by Date: Re: mathematica 4.0 + NET/Link doesn't work?
  • Next by Date: Re: Programming with .NET / hand over of variables
  • Previous by thread: 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