Finding the Number of Pythagorean Triples below a bound

[titus_...@yahoo.com]

Hello all,

Given the equation a^2+b^2 = c^2 with (a,b,c) positive integers,
0<a<=b. (For argument's sake, a=b will be allowed.) Let S(10^m) be the
number of solutions with hypotenuse c < bound 10^m. (They need not be
primitive.)

Problem: Find S(10^m).

What would be an efficient Mathematica code for this? The explicit
(a,b,c)'s are not needed, just S(10^m). I also need "m" as high as
possible, maybe m=6,7, with a reasonable run-time, say, less than an
hour. (A kind person gave me a code but told me it may be practical
only up to m=3.)

I stumbled upon some code some time ago, but it only gives the explicit
(a,b,c)'s with c<10^m, not the count:

Block[{a,b}, For[a=1,a<(10^m/Sqrt[2]),a++, For[b=a,b<10^m,b++,
If[Sqrt[a^2+b^2]==Round[Sqrt[a^2+b^2]]&&Sqrt[a^2+b^2]<10^m,
Print[{a,b,Sqrt[a^2+b^2]}]]]]]

for some positive integer "m". Can the above program be modified to
give S(10^m)? Or what program would be fast, for m=6,7?

Any help would be appreciated. Thanks.

-Titus

[Peter Pein]

titus_...@yahoo.com schrieb:

Hi Titus,

using a huge amount of RAM, on can obtain quite short times:

See http://mathworld.wolfram.com/PythagoreanTriple.html, formula (11)

In[1]:= Reduce[{0 < z2^2 - z1^2 && 0 < 2*z1*z2 && z1^2 + z2^2 < 10^m && z1 > 0}, {z1, z2}]
Out[1]= 10^m > 0 && 0 < z1 < Sqrt[10^m]/Sqrt[2] && z1 < z2 < Sqrt[10^m - z1^2]

In[2]:= pythagoCount[m_] :=
Total[Floor[(10^m - 1/2)/Sow[Flatten[Table[(Pick[z^2 + #1^2, (GCD[z, #1] == 1 & ) /@ #1] & )[
Range[z + 1, Floor[Sqrt[10^m - z^2]], 2]], {z, Floor[Sqrt[10^m/2]]}], 1]]]]

In[3]:= {thetime, {allcount, primcount}} =
Timing[Transpose[(Through[{First, Length[hyp = #1[[2,1,1]]] & }[Reap[pythagoCount[#1]]]] & ) /@
Range[7]]]
Out[3]= {11.141*Second,
{{1, 50, 878, 12467, 161431, 1980636, 23471468}, {1, 16, 158, 1593, 15919, 159139, 1591579}}}

In[4]:= N[1/(2*Pi) - primcount/10^Range[7]]
Out[4]= {0.05915494309189534, -0.0008450569081046577, 0.001154943091895344, -0.00014505690810465155,
-0.00003505690810465256, 0.00001594309189534293, -2.9569081046454393*^-6}

For primcount see http://www.research.att.com/~njas/sequences/A101931 and for allcount see http://www.research.att.com/~njas/sequences/A101929.

HTH,
Peter

[titus_...@yahoo.com]

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

[Andrzej Kozlowski]

On 29 Jul 2006, at 07:00, titus_...@yahoo.com wrote:

> Hello all,
>
> Given the equation a^2+b^2 = c^2 with (a,b,c) positive integers,
> 0<a<=b. (For argument's sake, a=b will be allowed.) Let S(10^m) be

> the

> number of solutions with hypotenuse c < bound 10^m. (They need not be
> primitive.)
>
> Problem: Find S(10^m).
>
> What would be an efficient Mathematica code for this? The explicit
> (a,b,c)'s are not needed, just S(10^m). I also need "m" as high as
> possible, maybe m=6,7, with a reasonable run-time, say, less than an
> hour. (A kind person gave me a code but told me it may be practical
> only up to m=3.)
>
> I stumbled upon some code some time ago, but it only gives the
> explicit
> (a,b,c)'s with c<10^m, not the count:
>
> Block[{a,b}, For[a=1,a<(10^m/Sqrt[2]),a++, For[b=a,b<10^m,b++,
> If[Sqrt[a^2+b^2]==Round[Sqrt[a^2+b^2]]&&Sqrt[a^2+b^2]<10^m,
> Print[{a,b,Sqrt[a^2+b^2]}]]]]]
>
> for some positive integer "m". Can the above program be modified to
> give S(10^m)? Or what program would be fast, for m=6,7?
>
> Any help would be appreciated. Thanks.
>
> -Titus
>

Let me start by describing two possible way of attempting to solve
your problem with the help of Mathematica. Afterwords I make a few
remarks about the feasibility of solving your problem in this way for
large values of m.
I will describe two methods: a slow one and a much faster one. The
reason for giving two methods is that in this way we can be more
confident that the answers we are getting are correct ones. Even the
slow method will be vastly faster than the one included in your post.

The first (slow) method makes use of the function
SumOfSquaresRepresentations , which is defined in the package
NumberTheory`NumberTheoryFunctions`. First we load the package:

<< NumberTheory`NumberTheoryFunctions`

?SumOfSquaresRepresentations

SumOfSquaresRepresentations[d, n] gives the list of all \
representations of the positive integer n as a sum of d squares of \
integers.

For example:

SumOfSquaresRepresentations[2, 50^2]

{{-50, 0}, {-48, -14}, {-48, 14},
{-40, -30}, {-40, 30},
{-30, -40}, {-30, 40},
{-14, -48}, {-14, 48},
{0, -50}, {0, 50}, {14, -48},
{14, 48}, {30, -40}, {30, 40},
{40, -30}, {40, 30}, {48, -14},
{48, 14}, {50, 0}}

As you see, this is not quite right, as it includes non-positive
integers as well as the same representation more then once. With
Mathematica it is easy to convert this to the form you want:

pythagoreans1[m_] := Module[{g}, g[n_] := Cases
[SumOfSquaresRepresentations[2,
n], {x_?Positive, y_?Positive} /; x > y]; Flatten[DeleteCases
[g /@
Table[i^2, {i, 1, m}], {}], 1]]

Actually this is not quite the form you wanted because pythagoreans
[m] finds all pythagorean pairs (x,y) with x^2+y^2<=m^2. However, I
find this form most convenient so I will continue using it here. It
can easily be converted to the form you wanted. So we have:

pythagoreans1[50]

{{4, 3}, {8, 6}, {12, 5}, {12, 9}, {15, 8}, {16, 12},
{20, 15}, {24, 7}, {24, 10}, {21, 20}, {24, 18},
{30, 16}, {28, 21}, {35, 12}, {36, 15}, {32, 24},
{40, 9}, {36, 27}, {40, 30}, {48, 14}}

If we do not wish to see the output but only compute the number of
representations we can use the function NumberOFPythagoreans1 defined
as follows:

NumberOfPythagoreans1[m_] := Length[pythagoreans1[m]]

So that:

NumberOfPythagoreans1[50]

20

The method is inefficient and to make it much faster we shall use
some elementary number theory. We know that (e.g., see Hardy and
Wright, "An Introduction to The Theory Of Numbers) that all the
Pythagorean triples (a,b,c) with coprime a and b have the form:
{2*n*m, m^2 - n^2, m^2 + n^2} where n and m are coprime positive
integers and m>n. So we first implement a function
CoprimePythagoreans as follows:

CoprimePythagoreans[p_] := DeleteCases[
With[{q = Sqrt[p]}, Flatten[
Table[If[GCD[n, m] == 1, {2*n*m, m^2 - n^2,
m^2 + n^2}], {n, 1, q}, {m, n + 1,
Sqrt[q^2 - n^2], 2}], 1]], Null]

this computes coprime puthagorean triples in which the "hypotenuse"
is smaller than p, e.g.

CoprimePythagoreans[50]

{{4, 3, 5}, {8, 15, 17}, {12, 35, 37}, {12, 5, 13},
{20, 21, 29}, {24, 7, 25}, {40, 9, 41}}

We can now compute the total number of all pythagorean triples with
hypotenuse less than 50 without actually having to list them:

NumberOfPythagoreans2[k_] := Total[Floor[k/#] & /@ CoprimePythagoreans
[k][[All, 3]]]

NumberOfPythagoreans2[50]

20

Happily this agrees with the previous answer. However, the difference
in performance is pretty huge:

NumberOfPythagoreans1[100]//Timing

{0.427973 Second,52}

v.s.

NumberOfPythagoreans2[100]//Timing

{0.005594 Second,52}

In spite of its impressive speed I doubt that NumberOfPythagoreans2
can cope with your problem for values of m as large as you want (I
have not tried testing it). The problem is not actually speed but
memory. The number of triples that have to be stored is a huge one,
and probably will stretch or exceed the capabilities of even a
powerful computer. This is, of course, something you can test
yourself. The alternative approach would be to try to work out the
number of pythagorean triples without actually storing them as a list
but this is, of course, an entirely different type of problem.

Andrzej Kozlowski

[Andrzej Kozlowski]

When I wrote "an entirely different type of problem" I had, of
course, in mind, an entirely mathematical solution. But almost
immediately after posting my message I realized that there is a very
simple computer solution based on the idea above, that avoids the
storage difficulty altogether. Here it is:

NumberOfPythagoreans[p_] := With[{q = Sqrt[p]}, Sum[If[GCD[n, m] ==
1, Floor[p/(n^2 + m^2)], 0], {n, 1, q},
{m, n + 1, Sqrt[q^2 - n^2], 2}]]

Now,

NumberOfPythagoreans[1000]//Timing

{0.010012 Second,881}

and

NumberOfPythagoreans[10^6]//Timing

{4.70295 Second,1980642}

Which seems to me to solve your problem completely.

Andrzej Kozlowski

[Andrzej Kozlowski]

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

[Carl K. Woll]

titus_...@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

As specified, your problem has polynomials which yield an S value of
Infinity. Besides the trivial example b-a, there are also polynomials
which gives Pell's equation, e.g., b^2-2a^2:

b^2 - 2a^2 == 1

has an infinite number of integer solutions with 0 < a <= b.

Is this intentional? If your problem has other restrictions which force
the number of solutions to be finite you should tell us, as algorithms
for the finite case will probably be very different from algorithms for
the infinite case.

Carl Woll
Wolfram Research

[titus_...@yahoo.com]

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

[Andrzej Kozlowski]

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

[Andrzej Kozlowski]

> 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

[titus_...@yahoo.com]

Hello Andrzej,

The code I used was this:

c1=Sqrt[N[a^2+b^2-k]];

countTriples[m_]:=Block[{a,b,i=0},For[a=1,a<(10^m/Sqrt[2]),a++,
For[b=a,b<10^m,b++,If[c1==Round[c1]&&c1<10^m,i++]]];i]

though m=3 takes about 50 sec (in an old comp) while m=4 will prolly be

around 5000 sec (83 mins).

-Titus

[Andrzej Kozlowski]

>> test[n_] :=

>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>> Integers]
>>

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

[Andrzej Kozlowski]

>>> test[n_] :=

>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>> Integers]
>>>

> Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
> Integers]
>
>

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

[Andrzej Kozlowski]

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

>>>> test[n_] :=

>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>> Integers]
>>>>

>> Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
>> Integers]
>>
>>

> Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
> Integers]
>

[Andrzej Kozlowski]

This is very odd indeed. For example, I get:

c1=Sqrt[N[a^2+b^2-k]];

countTriples[m_]:=Block[{a,b,i=0},For[a=1,a<(10^m/Sqrt[2]),a++,
For[b=a,b<10^m,b++,If[c1==Round[c1]&&c1<10^m,i++]]];i]

In[6]:=
k=3;

countTriples[2]

Invalid comparison
with 0. I attempted.

13

The correct answer is:

f[10^2,-3]

28

Similarly for other values. Your code usually doesn't work on my
machine and when it does it produces wrong answers.

Andrzej Kozlowski

[Andrzej Kozlowski]

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

>>>>> test[n_] :=

>>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>>> Integers]
>>>>>

>>> Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
>>> Integers]
>>>
>>>

>> Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
>> Integers]
>>

[Andrzej Kozlowski]

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
>
>
>

>>>>>>> test[n_] :=

>>>>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>>>>> Integers]
>>>>>>>

>>>>> Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
>>>>> Integers]
>>>>>
>>>>>

>>>> Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
>>>> Integers]
>>>>

[Andrzej Kozlowski]

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

>>>>>>>> test[n_] :=

>>>>>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>>>>>> Integers]
>>>>>>>>

>>>>>> Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
>>>>>> Integers]
>>>>>>
>>>>>>

>>>>> Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
>>>>> Integers]
>>>>>

[Andrzej Kozlowski]

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

>>>>>>>>> test[n_] :=

>>>>>>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt
>>>>>>>>> [n],
>>>>>>>>> Integers]
>>>>>>>>>

>>>>>>> Prime[Random[Integer, {2, 20}]]] ­ -1 && Element[Sqrt[n],
>>>>>>> Integers]
>>>>>>>
>>>>>>>

>>>>>> Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
>>>>>> Integers]
>>>>>>

[Andrzej Kozlowski]

{1.0032 Second,282}

Andrzej Kozlowski

>>>>>> test[n_] :=

>>>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>>>> Integers]
>>>>>>

>>>> 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

>>> Prime[Random[Integer, {2, 100}]]] ­ -1 && Element[Sqrt[n],
>>> Integers]
>>>

[da...@wolfram.com]

> 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

Having been involved in some of the private email I have a bit of an
advantage.

I'd say up to 10^6 or maybe 10^7 the way to go is by looking at the
factorization of c2=c^2+k. This way you only need iterate over c. Also as
Andrzej Kozlowski pointed out, one can remove from consideration all c2
for which, after removing factors of 2, the remaining odd is 3 modulo 4. I
use a prescan to enforce that if factors of 3,7, or 11 are present then so
are their squares (note: I did not check whether this helps for speed but
I suspect it does for large n). The full FactorInteger check enforces that
any factor equal to 3 mod 4 appears to even degree. After that I punt to
NumberTheory`NumberTheoryFunctions` to use
OrderedSumOfSquaresRepresentations. Finally, I found it faster to use this
on c2 with powers of 2 removed. One can show this gives the correct
result.

countTriples[m_, k_] := Module[
{c2, c2odd, total = 0, fax, add, g},
Do[
c2 = c^2 + k;
If[c2 < 5, Continue[]];
c2odd = c2;
While[EvenQ[c2odd], c2odd /= 2];
If[Mod[c2odd, 4] == 3, Continue[]];
g = GCD[c2odd, 231];
If[g ≠ 1 && g^2 ≠ GCD[c2odd, 53361], Continue[]];
fax = Mod[FactorInteger[c2odd], 4];
If[Apply[Or, Map[#[[1]] == 3 && OddQ[#[[2]]] &, fax]], Continue[]];
add = OrderedSumOfSquaresRepresentations[2, c2odd];
total += Length[add];
, {c, 1, m - 1}];
total]

Quick test for agreement with your results:

In[96]:=countTriples[10^3,-5]+countTriples[10^3,5]
Out[96]=426

Here are some timings.

In[54]:=Timing[countTriples[10^4,3]]
Out[54]={0.591 Second,1154}

In[93]:=Timing[countTriples[10^5,3]]
Out[93]={8.562 Second,12115}

In[95]:=Timing[countTriples[10^5,2]]
Out[95]={24.054 Second,9912}

In[94]:=Timing[countTriples[10^6,3]]
Out[94]={131.91 Second,121054}

I would think it will handle 10^7 reasonably well (maybe an hour or two).
But that depends on notgetting bogged in too many relatively expensive
factorizations.

As I noted in email, I would not be surprised if there are further
improvements, maybe by avoiding the OrderedSumOfSquaresRepresentations
usage.

Daniel Lichtblau
Wolfram Research

[Andrzej Kozlowski]

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.
>>>>

>>>>>>>>>>> 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]
>>>>>>>>>>>

>>>>>>>>> Prime[Random[Integer, {2, 20}]]] -1 && Element[Sqrt[n],
>>>>>>>>> Integers]
>>>>>>>>>
>>>>>>>>>

>>>>>>>> Prime[Random[Integer, {2, 100}]]] -1 && Element[Sqrt[n],
>>>>>>>> Integers]
>>>>>>>>

[Andrzej Kozlowski]

On 5 Aug 2006, at 22:32, da...@wolfram.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
>

> If[g 狎 1 && g^2 狎 GCD[c2odd, 53361], Continue[]];

Congratulations, this is rather spectacular, I think.

I have to confess, I was very sceptical about the prospects of this
method; now I think it was because I forgot how much good mathematics
has gone into making FactorInteger as fast as it is!
Actually, it is vastly faster not only than my compiled code but also
than the C code that can be found here:
http://pat7.com/jp/221k.c
Which only goes to show, what all of us here have always known
anyway ;-) namely, that to solve mathematical problems it is better
to use mathematical methods and Mathematica than general purpose
programming languages like C.

Andrzej Kozlowski

[Andrzej Kozlowski]

I got

test1[5,2]//Timing

{5136.24 Second,8929}

countTriplesP[4,2]

983

8929+983

Andrzej Kozlowski

>>>>>>>>>> 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]
>>>>>>>>>>

>>>>>>>> Prime[Random[Integer, {2, 20}]]] -1 && Element[Sqrt[n],
>>>>>>>> Integers]
>>>>>>>>
>>>>>>>>

>>>>>>> Prime[Random[Integer, {2, 100}]]] -1 && Element[Sqrt[n],
>>>>>>> Integers]
>>>>>>>

[Andrzej Kozlowski]

On 6 Aug 2006, at 10:01, Titus Piezas wrote:

> Hello Dan and Andrzej,

>
> >I'd say up to 10^6 or maybe 10^7 the way to go is by looking at the
> >factorization of c2=c^2+k. This way you only need iterate over c.
>

> Yes, turns out the problem was equivalent to counting the number of
> representations of N as
>
> a^2+b^2 = N
>
> where 0<a<=b (I should have phrased it that way, but that is the
> benefit of high-sight) which I assume can be given by the function
> OrderedSumOfSquaresRepresentations[2,N]? Hence the code you have
> given me, in essence, is, say for k = 5,
>
> Sum[OrderedSumOfSquaresRepresentations[2,c^2+5],{c,1,10^m-1}]
>
> roughly speaking? (With finer details of removing from
> consideration those c^2+5 that are of from 4m+3, etc.)

Well, no, it is rather:

Sum[Length[OrderedSumOfSquaresRepresentations[2,c^2+5]],{c,1,10^m-1}]

In fact I forgot that there was OrderedSumOfSquaresRepresentations as
well as SumOfSquaresRepresentations, which is the main reason why I
did not think it would be very useful. Another reason is that, I
expected, computing Length of many lists and then adding these
Lengths up should be slower than just counting the total number of
elements.

>
> >Here are some timings.
> >
> >In[54]:=Timing[countTriples[10^4,3]]
> >Out[54]={0.591 Second,1154}
> >
> >In[93]:=Timing[countTriples[10^5,3]]
> >Out[93]={8.562 Second,12115}
> >

> >In[94]:=Timing[countTriples[10^6,3]]
> >Out[94]={131.91 Second,121054}

> There is something odd with the first value. James' count (up to
> m=5) vs this gives the sequences,
>
> {1219, 12115}
> {1154, 12115, 121054}
>
> Why the disparity? If the ratio S(N)/N would converge, it should
> start as 0.121... (Note however, that James' table has bound 10^m
> *not* 10^m-1.}

My code also gives:

countTriplesP[4,3]

1219

Changing 10^m-1 to 10^m should give a value at least as large, and in
fact gives the same answer.

which suggests that something still needs to be looked carefully at
in Daniel's code. I can't do it now but will think about it when I can.

Andrzej

[Andrzej Kozlowski]

I have looked carefully at Daniel's code and I found I needn't have
bothered ;-
On my PowerBook Daniel's code gives:

countTriples[10^4,3]//Timing

{0.94218 Second,1219}

which shows that my computer is about half the speed of his but also
that somehow he probably pasted the wrong output into his posting?
Everything about the code seems perfect.

In fact, I am sure we could probably make it even faster if we looked
inside the code for OrderedSumOfSquaresRepresentations and changed
it so that it would count the number of representations rather than
make a list of them. But I am not sure that the extra speed we would
gain would be worth the extra effort.

Andrzej Kozlowski

[da...@wolfram.com]

>
> On 6 Aug 2006, at 10:01, Titus Piezas wrote:
>
>> Hello Dan and Andrzej,
>>
>> >I'd say up to 10^6 or maybe 10^7 the way to go is by looking at the
>> >factorization of c2=c^2+k. This way you only need iterate over c.
>>
>> Yes, turns out the problem was equivalent to counting the number of
>> representations of N as
>>
>> a^2+b^2 = N
>>
>> where 0<a<=b (I should have phrased it that way, but that is the
>> benefit of high-sight) which I assume can be given by the function
>> OrderedSumOfSquaresRepresentations[2,N]? Hence the code you have
>> given me, in essence, is, say for k = 5,
>>
>> Sum[OrderedSumOfSquaresRepresentations[2,c^2+5],{c,1,10^m-1}]
>>
>> roughly speaking? (With finer details of removing from
>> consideration those c^2+5 that are of from 4m+3, etc.)
>
> Well, no, it is rather:
>
> Sum[Length[OrderedSumOfSquaresRepresentations[2,c^2+5]],{c,1,10^m-1}]
>
> In fact I forgot that there was OrderedSumOfSquaresRepresentations as
> well as SumOfSquaresRepresentations, which is the main reason why I
> did not think it would be very useful. Another reason is that, I
> expected, computing Length of many lists and then adding these
> Lengths up should be slower than just counting the total number of
> elements.

Length[] is quite fast (constant time) to compute, so that aspect at least
should not be an issue.

>> >Here are some timings.
>> >
>> >In[54]:=Timing[countTriples[10^4,3]]
>> >Out[54]={0.591 Second,1154}
>> >
>> >In[93]:=Timing[countTriples[10^5,3]]
>> >Out[93]={8.562 Second,12115}
>> >
>> >In[94]:=Timing[countTriples[10^6,3]]
>> >Out[94]={131.91 Second,121054}
>> There is something odd with the first value. James' count (up to
>> m=5) vs this gives the sequences,
>>
>> {1219, 12115}
>> {1154, 12115, 121054}
>>
>> Why the disparity? If the ratio S(N)/N would converge, it should
>> start as 0.121... (Note however, that James' table has bound 10^m
>> *not* 10^m-1.}
>
> My code also gives:
>
>
> countTriplesP[4,3]
>
> 1219
>
> Changing 10^m-1 to 10^m should give a value at least as large, and in
> fact gives the same answer.
>
>
> which suggests that something still needs to be looked carefully at
> in Daniel's code. I can't do it now but will think about it when I can.
>
> Andrzej
>

I think I see the problem. Looking at the In/Out numbers I notice the 10^4
case was run much earlier than the others. I am guessing it was run when I
had not fully debugged the code and was getting results that were too
small (pecifically, I was missing the cases with factors of form 4*k+3
that appear to even powers). I may simply have failed to rerun that one
before cutting-and-pasting into my note.

If you run the code I sent and get the correct result that will confirm
the matter. I'll try it later as right now I am not on a machine with a
reliable Mathematica installation.

Daniel

[Titus Piezas]

Hello Dan and Andrzej,

>I'd say up to 10^6 or maybe 10^7 the way to go is by looking at the
>factorization of c2=c^2+k. This way you only need iterate over c.

Yes, turns out the problem was equivalent to counting the number of representations of N as

a^2+b^2 = N

where 0<a<=b (I should have phrased it that way, but that is the benefit of high-sight) which I assume can be given by the function OrderedSumOfSquaresRepresentations[2,N]? Hence the code you have given me, in essence, is, say for k = 5,

Sum[OrderedSumOfSquaresRepresentations[2,c^2+5],{c,1,10^m-1}]

roughly speaking? (With finer details of removing from consideration those c^2+5 that are of from 4m+3, etc.)

>Here are some timings.
>
>In[54]:=Timing[countTriples[10^4,3]]
>Out[54]={0.591 Second,1154}
>
>In[93]:=Timing[countTriples[10^5,3]]
>Out[93]={8.562 Second,12115}
>
>In[94]:=Timing[countTriples[10^6,3]]
>Out[94]={131.91 Second,121054}

There is something odd with the first value. James' count (up to m=5) vs this gives the sequences,

{1219, 12115}
{1154, 12115, 121054}

Why the disparity? If the ratio S(N)/N would converge, it should start as 0.121... (Note however, that James' table has bound 10^m *not* 10^m-1.}

-Tito

[Andrzej Kozlowski]

On 7 Aug 2006, at 07:40, Andrzej Kozlowski wrote:

>
> On 6 Aug 2006, at 13:33, Andrzej Kozlowski wrote:
>
>>

>> On 6 Aug 2006, at 10:01, Titus Piezas wrote:
>>

>>> Hello Dan and Andrzej,
>>>
>>>> I'd say up to 10^6 or maybe 10^7 the way to go is by looking at the
>>>> factorization of c2=c^2+k. This way you only need iterate over c.
>>>
>>> Yes, turns out the problem was equivalent to counting the number
>>> of representations of N as
>>>
>>> a^2+b^2 = N
>>>
>>> where 0<a<=b (I should have phrased it that way, but that is the
>>> benefit of high-sight) which I assume can be given by the function
>>> OrderedSumOfSquaresRepresentations[2,N]? Hence the code you have
>>> given me, in essence, is, say for k = 5,
>>>
>>> Sum[OrderedSumOfSquaresRepresentations[2,c^2+5],{c,1,10^m-1}]
>>>
>>> roughly speaking? (With finer details of removing from
>>> consideration those c^2+5 that are of from 4m+3, etc.)
>>

>> Well, no, it is rather:
>>
>> Sum[Length[OrderedSumOfSquaresRepresentations[2,c^2+5]],{c,1,10^m-1}]
>>
>> In fact I forgot that there was OrderedSumOfSquaresRepresentations
>> as well as SumOfSquaresRepresentations, which is the main reason
>> why I did not think it would be very useful. Another reason is
>> that, I expected, computing Length of many lists and then adding
>> these Lengths up should be slower than just counting the total
>> number of elements.
>>
>>>

>>>> Here are some timings.
>>>>
>>>> In[54]:=Timing[countTriples[10^4,3]]
>>>> Out[54]={0.591 Second,1154}
>>>>
>>>> In[93]:=Timing[countTriples[10^5,3]]
>>>> Out[93]={8.562 Second,12115}
>>>>
>>>> In[94]:=Timing[countTriples[10^6,3]]
>>>> Out[94]={131.91 Second,121054}
>>> There is something odd with the first value. James' count (up to
>>> m=5) vs this gives the sequences,
>>>
>>> {1219, 12115}
>>> {1154, 12115, 121054}
>>>
>>> Why the disparity? If the ratio S(N)/N would converge, it should
>>> start as 0.121... (Note however, that James' table has bound 10^m
>>> *not* 10^m-1.}
>>

>> My code also gives:
>>
>>
>> countTriplesP[4,3]
>>
>> 1219
>>
>> Changing 10^m-1 to 10^m should give a value at least as large, and
>> in fact gives the same answer.
>>
>>
>> which suggests that something still needs to be looked carefully at
>> in Daniel's code. I can't do it now but will think about it when I
>> can.
>>
>> Andrzej
>
>

> I have looked carefully at Daniel's code and I found I needn't have
> bothered ;-
> On my PowerBook Daniel's code gives:
>
>
> countTriples[10^4,3]//Timing
>
> {0.94218 Second,1219}
>
> which shows that my computer is about half the speed of his but also
> that somehow he probably pasted the wrong output into his posting?
> Everything about the code seems perfect.
>
> In fact, I am sure we could probably make it even faster if we looked
> inside the code for OrderedSumOfSquaresRepresentations and changed
> it so that it would count the number of representations rather than
> make a list of them. But I am not sure that the extra speed we would
> gain would be worth the extra effort.
>
> Andrzej Kozlowski
>

As Daniel pointed out, my last comment above was based on forgetting
how fast Length is in Mathematica - now "I am sure that we
probably" ;-) could not make the code significantly faster by
tinkling with OrderedSumOfSquaresRepresentations.

However, there are two additional comments, that I would like to
make, which I think are of some interest. First, I think I know how
Daniel got his wrong answer

In[54]:=Timing[countTriples[10^4,3]]
Out[54]={0.591 Second,1154}

If you forget to load the package

<< NumberTheory`NumberTheoryFunctions`

And run Daniel's code above you will get the answer 1154. This is
what must have happened.

Secondly, I think I have found another speed improvement to Daniel's
code. This is my new version:

countTriples1[m_, k_] := Module[

{c2, c2odd, total = 0, fax, add, g},
Do[
c2 = c^2 + k;
If[c2 < 5, Continue[]];
c2odd = c2;
While[EvenQ[c2odd], c2odd /= 2];
If[Mod[c2odd, 4] == 3, Continue[]];
g = GCD[c2odd, 231];

If[g ­ 1 && g^2 ­ GCD[c2odd, 53361], Continue[]];

fax = Mod[FactorInteger[c2odd], 4];
If[Apply[Or, Map[#[[1]] == 3 && OddQ[#[[2]]] &, fax]], Continue
[]];

add = If[PrimeQ[c2odd], 1, Length
[OrderedSumOfSquaresRepresentations[2, \
c2odd]]];
total += add;

, {c, 1, m - 1}];
total]

With Daniel's code I get on my PowerBook:

In[4]:=
countTriples[10^4,3]//Timing

Out[4]=
{2.23016 Second,1219}

With the new code I get:

In[5]:=
countTriples1[10^4,3]//Timing

Out[5]=
{0.728627 Second,1219}

Which is more than 3 times faster. The improvement is again based on
a bit of mathematics and a some powerful Mathematica technology. The
mathematics is that it is known and not hard to prove that a prime
number the form 4m+1 there has exactly one representation as a sum
of two squares. Of course to use this fact we need PrimeQ to be as
fast as it is; there is no way you could use this fact without having
a function like that at your disposal.

It seems the code is now pretty close to the optimum.

Andrzej Kozlowski