cp's OEIS Frontend

A380858 a(n) is the number of primes p <= n such that p^(p + n) == p (mod p + n).

%I A380858 #18 Mar 13 2025 08:56:37
%S A380858 0,0,2,1,2,1,1,2,1,3,2,3,1,3,2,3,2,4,1,3,1,3,1,6,0,6,1,4,2,7,1,3,0,6,
%T A380858 3,6,1,5,2,5,2,8,1,5,1,5,1,8,0,6,2,5,1,9,0,8,1,5,3,12,1,8,1,7,2,11,1,
%U A380858 8,2,8,2,10,1,6,0,9,1,12,1,7,1,5,1,13,0,9,3,6,1,15
%N A380858 a(n) is the number of primes p <= n such that p^(p + n) == p (mod p + n).
%H A380858 Robert Israel, <a href="/A380858/b380858.txt">Table of n, a(n) for n = 1..10000</a>
%e A380858 a(3) = 2 because 2^(2+3) = 32 mod (2+3) is equal to 2 and 3^(3+3) = 729 mod (3+3) is equal to 3;
%e A380858 a(4) = 1 because 2^(2+4) = 64 mod (2+4) is equal to 4, but not is equal to 2, and 3^(3+4) = 2187 mod (3+4) is equal to 3.
%p A380858 P:= NULL: R:= NULL:
%p A380858 for n from 1 to 100 do
%p A380858   if isprime(n) then P:= P,n fi;
%p A380858   R:= R, nops(select(p -> p &^ (p+n) mod (p+n) = p, [P]));
%p A380858 od:
%p A380858 R; # _Robert Israel_, Mar 12 2025
%t A380858 a[n_] := Count[Range[n], _?(PrimeQ[#] && PowerMod[#, # + n, # + n] == # &)]; Array[a, 100] (* _Amiram Eldar_, Feb 06 2025 *)
%o A380858 (Magma) [#[p: p in PrimesUpTo(n) | p^(p+n) mod (p+n) eq p]: n in [1..90]];
%o A380858 (PARI) a(n) = my(nb=0); forprime(p=2, n, if (Mod(p, p+n)^(p+n) == p, nb++)); nb; \\ _Michel Marcus_, Feb 06 2025
%Y A380858 Cf. A000040, A371883.
%K A380858 nonn,look
%O A380858 1,3
%A A380858 _Juri-Stepan Gerasimov_, Feb 06 2025