Order of the linear difference equation that solves the Diophantine equation x^2 + (x+n)^2 = y^2.
3, 3, 3, 3, 3, 3, 7, 3, 3, 3, 3, 3, 3, 7, 3, 3, 7, 3, 3, 3, 7, 3, 7, 3, 3, 3, 3, 7, 3, 3, 7, 3, 3, 7, 7, 3, 3, 3, 3, 3, 7, 7, 3, 3, 3, 7, 7, 3, 11, 3, 7, 3, 3, 3, 3, 7, 3, 3, 3, 3, 3, 7, 7, 3, 3, 3, 3, 7, 7, 7, 7, 3, 7, 3, 3, 3, 7, 3, 7, 3, 3, 7, 3, 7, 7, 3, 3
1
For instance, the equation x^2 + (x+1)^2 = y^2 has solutions that obey a third-order linear difference equation x(n+1) = 7*x(n) - 7*x(n-1) + x(n-2) with initial terms (0, 3, 20). It appears that all primes of the form 8k+-1 satisfy a seventh-order linear difference equation; all other primes satisfy a third-order linear difference equation. See sequence S000643 for the difference equation and the initial terms for n up to 1000.
T. D. Noe, Plot of 1000 terms
T. D. Noe, Table of 1000 terms
(Mma) PerfectSquareQ[n_] := JacobiSymbol[n, 13] =!= -1 && JacobiSymbol[n, 19] =!= -1 && JacobiSymbol[n, 17] =!= -1 && JacobiSymbol[n, 23] =!= -1 && IntegerQ[Sqrt[n]]; FindCoef[seq_List, i_, n_] := Module[{m, a, s, noSoln = {}}, a = Take[seq, {i, i + 2 n - 1}]; m = Reverse /@ Partition[Most[a], n, 1]; If[Det[m] == 0, noSoln, s = LinearSolve[m, Take[a, -n]]; If[And @@ IntegerQ /@ s, s, noSoln]]]; Table[td = {}; n = 0; c = {}; While[If[PerfectSquareQ[n^2 + (n + d)^2], AppendTo[td, n]; If[Length[td] > 2 && OddQ[Length[td]], c = FindCoef[td, 1, (Length[td] - 1)/2]]]; c == {}, n++]; Length[c], {d, 100}]
Cf. A129837, A201916, A201917, S000640, S000643.
nonn,hard,nice
T. D. Noe, May 21 2015