Some hand analysis first.
Suppose we have numbers A<B<C<D, there are obviously no two equal numbers. We first pick arbitrary number N. Then let A=N-1, B=N+1, so A*B = N² - 1. Now we can use C*A = (2N-1)²-1 = 4N²-4N = 4N(N-1). This sets C=4N Then C*B = 4N*(N+1) = (2N+1)² - 1. Finally, we require D. If we take D=(4N²-1)*C we get DC = (4N²-1)*(16N²) = 64N⁴-16N² = (8N²-1)² - 1, DB = (4N²-1)*4N*(N+1) = 16N⁴+16N³-4N²-4N = (4N²+2N-1)² - 1 and DA = (4N²-1)*4N*(N-1) = 16N⁴-16N³-4N²+4N = (4N²-2N-1)² - 1.
This is sufficient to prove there are infinitely many 4 element sets.
These are NOT all of the 4 element options.
For example, there are also 1,8,15,528 or 2, 12, 24, 2380.
I believe you can start with any number (at all), then pick any number that generates almost square with the first one (there are infinitely many), then pick any third number that generates almost square with the first two (I guess still infinitely many options "but way fewer than before" - you know what I mean...), then pick the unique fourth number. As the fourth number is unique, there is obviously no fifth. But I don't know how to prove this, this is mainly how my test code seems to spit out results (ignoring limited search range. 1e6 in example; I tried this with up to 1e8 but it looks the same)
Code (matlab, most likely works in octave too but I haven't tried it there). findSQS does the actual looping and calling itself, while printing final possibilities. Initial thingy outside the function is to ensure we can start findSQS with arbitrary number, such as 1 or 2.
sqm = (1:1e6).^2 - 1;
for initCand = 1:1000
initList = sqm / initCand;
initList(initList ~= round(initList)) = [];
initList(initList <= initCand) = [];
[elems] = findSQS(initCand, initList, initList);
end
function [nextElems, nextList] = findSQS(elems, sqm, remList)
if (isempty(remList)) % no more numbers to add, print what we have.
nextElems = elems
nextList = [];
return
end
for i = 1 : length(remList) % pick any number still on the list, starting with first one.
NE = remList(i);
nextElems = [elems, NE]; % add element
nextList = sqm / NE; % all almost squares with the current one, after we
nextList(nextList ~= round(nextList)) = []; % delete non-integers
nextList(nextList < NE) = []; % and remove smaller numbers (remove duplicate sets).
nextList = intersect(remList, nextList); % Find intersection of that list with numbers not ruled out by previous elements.
[nextElems, nextList] = findSQS(nextElems, sqm, nextList); % And call itself again.
end
end
