[go: up one dir, main page]

login
A173656
Primes p such that p^2 divides P(p), where P = Perrin sequence A001608.
4
OFFSET
1,1
COMMENTS
It is not known if this sequence is infinite.
The squares are in A013998.
No other terms below 10^10. - Max Alekseyev, Aug 27 2023
LINKS
G. L. Honaker, Jr. and Chris Caldwell, Prime Curios! 521
Wikipedia, Perrin number
EXAMPLE
521 is in the sequence since its square 271441 is the factor of A001608(521).
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 */
CROSSREFS
Sequence in context: A167734 A122715 A153180 * A015291 A028484 A057699
KEYWORD
bref,hard,more,nonn
AUTHOR
STATUS
approved