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: [mg68442] 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:44 -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> <B935EC43-291F-4C01-B646-534789B6E908@mimuw.edu.pl> <F7F7DB69-7894-466A-99BC-065835567D5C@mimuw.edu.pl>
  • Sender: owner-wri-mathgroup at wolfram.com

Daniel Lichtblau has already informed me that running the code below  
on a 64 bit computer will not help at all, so there is no point  
trying :-(
However, he also suggested a way around the Compile problem with  
integers larger than machine integers :-)
I am trying it out right now.

Andrzej Kozlowski


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

> I need to add some corrections. There were some some small mistakes  
> in the code I posted earlier, which caused problems with running  
> the compiled code for negative values of k, and which also probably  
> accounted for the slightly incorrect answers which the  
> "approximate" code returned. I attributed the inaccuracy to the use  
> of numerical precision in the test for a number being a perfect  
> square but it seems that the cause was (probably) elsewhere. I  
> don't want to try to make this message too long so I won't bother  
> explaining what I think the mistakes were; but I will give what I  
> think is the correct code. I have decided to separate the code for  
> negative and positive values of k. The code for negative k works  
> also for positive k's but is slightly slower, due to the extra test  
> that needs to be performed. For this reason, and for the sake of  
> greater clarity I have decided to separate the two codes.
>
> The code for positive k:
>
> countTriplesP = 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]]
>
> The code for negative k:
>
> countTriplesN = Compile[{{m, _Integer}, {k, _Integer}}, Module[
>   {i = 1, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2 + k, 4] != 3 && c^2 + k ³ 0, 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]]
>
> Now we get:
>
> countTriplesP[1,6]+countTriplesN[1,-6]//Timing
>
>
> {0.000221 Second,3}
>
>
> countTriplesP[3,2]+countTriplesN[3,-2]//Timing
>
>
> {0.95186 Second,282}
>
>
> countTriplesP[4,2]+countTriplesN[4,-2]//Timing
>
>
> {95.2177 Second,2762}
>
>
>
>
> Note that these values are consistently less by one form the values  
> obtained by Titus and also by me using my earlier "exact" code, but  
> actually I believe that this disparity was due to mistakes in the  
> earlier code. In any case, if we replace the "numerical" test for a  
> perfect square by the much slower "exact" test, the answers will be  
> the same, so the difference of 1 is certainly not due to the use of  
> numerical precision. Anyway, everything works fast and looks  
> perfectly satisfactory but then there is a snag. I decided to run  
> the code for m=5 and k=2, started it and went out for several  
> hours. When I came back I was rather disappointed to see that it  
> was still running and then I saw the message:
>
> countTriplesP[5,2]//Timing
>
> CompiledFunction::"cfn" Numerical error  encountered at instruction  
> 10; proceeding with uncompiled evaluation.
>
> I assume the cause of this is that on 32 bit computers Compile  
> cannot deal with integers larger than 2^32 but we have:
>
>
> c=10^5;
>
> Log[2,c^2]//N
>
> 33.2193
>
> I can't at the moment see any way around this problem except to run  
> uncompiled code (far too long) or try a 64 bit computer.  
> Unfortunately I do not have one and I don't think Titus has one, so  
> if any body has a 64-bit computer with a 64 bit version of  
> Mathematica installed, and a little spare time, I think both of us  
> would like to know how the above code performs in this setting. If  
> it works it will provide a nice bit of advertising for 64 bit  
> computing.
>
> Andrzej Kozlowski
>
>
>
>
> On 4 Aug 2006, at 12:48, Andrzej Kozlowski wrote:
>
>> 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: Finding the Number of Pythagorean Triples below a bound
  • Next by Date: Re: Precision of arguments to FunctionInterpolation
  • Previous by thread: Re: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Next by thread: Re: Re: Re: Re: Finding the Number of Pythagorean Triples below a bound