Fast Way to Find Pythagorean Decompositions

Question

I am looking for a computationally-efficient way to solve s Diophantine equation related to Pythagorean triples.

In particular, the problem is the following:

Given a number $c$, what are the numbers $b_1,b_2,...,b_n$ such that both $c^2-b_k$ and $c^2+b_k$ are square numbers.

For example, given the number $c=65$, we have the values

$$[4056, 3696, 3000, 2016]$$

for the sequence $b_n$, since $65^2-4056$, $65^2+4056$, $65^2-3696$, $65^2+3696$, etc. are all perfect integer squares.

The naive approach that my program is using is relatively simple:

1. For all integer $a$ such that $1\leq a<c$, test if $2c^2-a^2$ is a perfect square.
2. If it is, then $c^2-a^2$ is in $b_n$.

Define a function $f(k)$ for integer $k$ such that $f(k)$ is the number of elements in the corresponding sequence $b_n$. For example, $f(65)=4$, as shown.

I decided to investigate the local maximums of the $f$, finding that they were successively

$$5,25,65,325,1105,5525,27625,32045,....$$

To my surprise, this sequence appears in the OEIS at http://oeis.org/A088959. The title of this integer sequence is "Lowest numbers which are d-Pythagorean decomposable, i.e. square is expressible as sum of two distinct squares in more ways than for any smaller number."

I noticed that the values in this entry have been computed up to a high value (the 70th term is computed as $49360435020710806319207602548125$), and I am confused as to how they did this. The naive approach I laid out would take billions of years on this number.

Thus, I am asking for a faster way to generate the sequence $b_n$ for some given integer $c$.

Answer 1

the primitive solutions to $$ u^2 + v^2 = 2 w^2, $$ that is $$ \gcd(u,v,w) = 1, $$ are given by $$ u = r^2 - 2 r s - s^2, \; \; v = r^2 + 2 r s - s^2, \; \; w = r^2 + s^2, $$ where $$ \gcd(r,s) = 1 $$ but $r,s$ are not both odd.

There is more that can be said... if you are testing $c$ that is divisible by any prime such as $3,7, 11, 19, 23,$ i.e. $p \equiv 3 \pmod 4,$ then you are wasting time, as the other numbers are divisible by the same prime.

Solving $r^2 + s^2 = w$ for prime $w \equiv 1 \pmod 4$ is fairly understandable. For composite $w,$ see Hardy Muskat Williams.