Re: Re: Re: Finding the Number of Pythagorean Triples below a bound
- To: mathgroup at smc.vnet.net
- Subject: [mg68441] 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:42 -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>
- Sender: owner-wri-mathgroup at wolfram.com
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
>>>>>
>>>>>
>>>>
>>>
>>
>
- References:
- Re: Finding the Number of Pythagorean Triples below a bound
- From: titus_piezas@yahoo.com
- Re: Re: Finding the Number of Pythagorean Triples below a bound
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: Finding the Number of Pythagorean Triples below a bound