Re: Re: Re: Re: Finding the Number of Pythagorean Triples below a bound
- To: mathgroup at smc.vnet.net
- Subject: [mg68468] Re: [mg68442] Re: [mg68382] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Sun, 6 Aug 2006 02:56:58 -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> <200608050746.DAA14183@smc.vnet.net> <71350A02-470E-4AFE-87B5-DCFACFC6483B@mimuw.edu.pl>
- Sender: owner-wri-mathgroup at wolfram.com
A small correction: 5136.24 seconds is about 1 hour and 25 minutes;
hence the total time taken was about 1 hour 26.5 minute, not 1 hour
45 minutes as I wrote earlier.
Andrzej Kozlowski
On 5 Aug 2006, at 15:40, Andrzej Kozlowski wrote:
> I have succeeded in computing the number of solutions of a^2
> +b^2==c^2+2 where c<10^5. This was done by means of the follwoign
> code:
>
> test1 = 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, 10.^(m - 1), 10.^m - 1}];
> i]]
>
> I got
>
>
> test1[5,2]//Timing
>
>
> {5136.24 Second,8929}
>
> To this one should add the outcome of
>
> countTriplesP[4,2]
>
> 983
>
> where countTriplesP is defined in my earlier post. So the actual
> number of solutions is
>
>
> 8929+983
>
> 9912
> This means that compiled Mathematica code can quite successfully
> deal with the case m=5 of Titus's problem. This computation took 1
> hour and 43 minutes on a 1 Gigahertz PowerBook G4, to which we have
> to add about 1.5 minute for solving the 10^4 case. That gives the
> total time of 1.45 minutes, which seems to me pretty impressive,
> and shows, I think, that Mathematica is quite usable for hard
> numerical computations, provided one one is lucky enough to have
> ones code compile properly. I really hope some work at WRI is going
> into making Compile work better (in more cases and more predictably).
>
>
> Andrzej Kozlowski
>
> On 5 Aug 2006, at 09:46, Andrzej Kozlowski wrote:
>
>> 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
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
- 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: 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