There is a much faster way that doesn't even require testing. You need a bit of elementary number theory to find that way, but here goes:
If n² + (n+1)²
is a perfect square, that means there is an m
such that
m² = n² + (n+1)² = 2n² + 2n + 1
<=> 2m² = 4n² + 4n + 1 + 1
<=> 2m² = (2n+1)² + 1
<=> (2n+1)² - 2m² = -1
Equations of that type are easily solved, starting from the "smallest" (positive) solution
1² - 2*1² = -1
of
x² - 2y² = -1
corresponding to the number 1 + √2
, you obtain all further solutions by multiplying that with a power of the primitive solution of
a² - 2b² = 1
which is (1 + √2)² = 3 + 2*√2
.
Writing that in matrix form, you obtain all solutions of x² - 2y² = -1
as
|x_k| |3 4|^k |1|
|y_k| = |2 3| * |1|
and all x_k
are necessarily odd, thus can be written as 2*n + 1
.
The first few solutions (x,y)
are
(1,1), (7,5), (41,29), (239,169)
corresponding to (n,m)
(0,1), (3,5), (20,29), (119,169)
You can get the next (n,m)
solution pair via
(n_(k+1), m_(k+1)) = (3*n_k + 2*m_k + 1, 4*n_k + 3*m_k + 2)
starting from (n_0, m_0) = (0,1)
.
Quick Haskell code since I don't speak MatLab:
Prelude> let next (n,m) = (3*n + 2*m + 1, 4*n + 3*m + 2) in take 20 $ iterate next (0,1)
[(0,1),(3,5),(20,29),(119,169),(696,985),(4059,5741),(23660,33461),(137903,195025)
,(803760,1136689),(4684659,6625109),(27304196,38613965),(159140519,225058681)
,(927538920,1311738121),(5406093003,7645370045),(31509019100,44560482149)
,(183648021599,259717522849),(1070379110496,1513744654945),(6238626641379,8822750406821)
,(36361380737780,51422757785981),(211929657785303,299713796309065)]
Prelude> map (\(n,m) -> (n^2 + (n+1)^2 - m^2)) it
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
Edit by EitanT:
Here's the MATLAB code to calculate the first N
numbers:
res = zeros(1, N);
nm = [0, 1];
for k = 1:N
nm = nm * [3 4; 2 3] + [1, 2];
res(k) = nm(1);
end
The resulting array res
should hold the values of n
that satisfy the condition of the perfect square.