
PROG

(PARI)
N=10^10;
default(primelimit, N);
M = [0, 1, 0, 0, 0; 0, 0, 1, 0, 0; 0, 0, 0, 1, 0; 0, 0, 0, 0, 1; 1, 0, 1, 1, 1];
a(n)=lift( trace( Mod(M, n)^n ) );
ta(n)=lift( trace( Mod(M, n) ) );
{ for (n=2, N,
if ( isprime(n), next() );
if ( a(n)==ta(n), print1(n, ", "); );
); }
/* Matt McIrvin, after Joerg Arndt's program for A013998, May 23, 2013 */
