
MATHEMATICA

lst = {}; a = 3; b = 0; c = 2; Do[P = b + a; If[PrimeQ[n] && Divisible[P, n^2], AppendTo[lst, n]]; a = b; b = c; c = P, {n, 3, 2*10^5}]; lst
lst = {}; PowerMod2[mat_, n_, m_] := Mod[Fold[Mod[If[#2 == 1, #1.#1.mat, #1.#1], m] &, mat, Rest@IntegerDigits[n, 2]], m]; LinearRecurrence2[coeffs_, init_, n_, m_] := Mod[First@PowerMod2[Append[RotateRight /@ Most@IdentityMatrix@Length[coeffs], coeffs], n, m].init, m] /; n >= Length[coeffs]; Do[n = Power[p, 2]; If[PrimeQ[p] && LinearRecurrence2[{1, 1, 0}, {3, 0, 2}, n, n] == 0, AppendTo[lst, p]], {p, 1, 2*10^5, 2}]; lst


PROG

(PARI)
/* Assuming b13998 containing second column of b013998.txt */
A013998 = readvec(b13998);
for (k=1, #A013998, if (issquare(A013998[k])==1, print(k, " ", A013998[k])));
/* Hugo Pfoertner, Sep 01 2017 */
