Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Re: Finding the Number of Pythagorean Triples below a bound

• To: mathgroup at smc.vnet.net
• Subject: [mg68290] Re: [mg68242] Re: Finding the Number of Pythagorean Triples below a bound
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Mon, 31 Jul 2006 04:33:42 -0400 (EDT)
• References: <eaeqa3$53v$1@smc.vnet.net> <200607300848.EAA25171@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

On 30 Jul 2006, at 10:48, titus_piezas at yahoo.com wrote:

> Hello all,
>
> To be specific, I am looking for some code applicable to the more
> general bivariate polynomial,
>
> Poly(a,b) = c^2
>
> where {a,b} are positive integers, 0 < a <=b, and of which the
> Pythagorean triples are just a special case. The problem is to find
> S(10^m) which is the number of solutions with c < 10^m for as  
> "high" as
> m=5,6 with a reasonable run-time of, say, an hour or less.
>
> Anybody knows of such code?
>
> -Titus
>


I doubt that one can construct efficient code that will deal with the  
general case. There have been several examples of related problems  
posted to this list in the past, for example look at the threads:

http://forums.wolfram.com/mathgroup/archive/2005/Oct/msg00594.html
http://forums.wolfram.com/mathgroup/archive/2005/Dec/msg00262.html

One idea that is useful when testing if a number is a perfect square  
is to perform the test module a group of suitable chosen primes  
rather than just taking the square root.  Let me illustrate this by  
an example. Let's start with a list of primes; I will take just the  
first 20 odd ones.

ls = Prime[Range[3, 20]];

We now define three tests for an integer to be a perfect square:

test1[n_] := Element[Sqrt[n], Integers]

test2[n_] :=
      Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],  
Integers]

test3[n_] := Not[MemberQ[JacobiSymbol[n, ls], -1]]

Of course test3 will sometimes produce "false" perfect squares. Note  
also the order of the subtests which make up test2 - this is  
important. test2 will first perform the JacobiSymbol subtest; only if  
this is passed will the second subtest be carried out. Now we compare  
the performance of each test:

testList = Range[10^6];


a=Select[testList,test1];//Timing

{155.21 Second,Null}


b=Select[testList,test2];//Timing


{58.8689 Second,Null}


c=Select[testList,test3];//Timing


{56.188 Second,Null}


We see that test2 is nearly three times faster than test1, while  
test3 is only somewhat faster than test2. Now let's compare the  
lenghts of lists we obtained:


Map[Length,{a,b,c}]

{1000,1000,1001}

We see that we got one "pseudo-square" in list c. The pseudo-square is:


Complement[c,a][[1]]

860574

If one is interested only in an approximate number of perfect squares  
test3 might be just as good as test2, but on this showing it does not  
seem to offer much gain in performance. This could be different in  
other situations. Also, using a smaller collection of test primes  
will probably, in the above situation, improve performance of test2.  
All such considerations seem to be dependent on  the particular  
problem at hand and have to be tested empirically, which is one  
reason why it is difficult to suggest any general code for this sort  
of problems.

Andrzej Kozlowski

• References:
  • Re: Finding the Number of Pythagorean Triples below a bound
    • From: titus_piezas@yahoo.com