cp's OEIS Frontend

A264737 Primes which divide some term of A000085 (numbers of involutions).

%I A264737 #23 Jul 28 2020 08:48:02
%S A264737 2,5,13,19,23,29,31,43,53,59,61,67,73,79,83,89,97,103,131,137,151,157,
%T A264737 163,173,179,181,191,197,199,211,229,233,239,241,281,293,307,317,347,
%U A264737 359,367,373,379,389,397,409,419,421,431,433,443,449,457,461,463,479,487,491,499
%N A264737 Primes which divide some term of A000085 (numbers of involutions).
%C A264737 Essentially the same as A245177. - _R. J. Mathar_, Nov 25 2015
%H A264737 Robert Israel, <a href="/A264737/b264737.txt">Table of n, a(n) for n = 1..4000</a>
%F A264737 Any individual prime p is easily tested for membership in this set by iterating the recurrence for A000085 mod p, T(n) = T(n-1) + (n-1)T(n-2) modulo p, until either finding a value divisible by p or entering a cycle.
%e A264737 23 divides A000085(11) = 35696 = 2^4 * 23 * 97, so it appears in this set. The sequence A000085 mod 3 cycles: 1,1,2,1,1,2,..., so the prime factor 3 does not appear in this set.
%p A264737 filter:= proc(p) local a,b,c,n,R;
%p A264737   if not isprime(p) then return false fi;
%p A264737   a:= 1; b:= 1;
%p A264737   R[1,1,1]:= 1;
%p A264737   for n from 2 do
%p A264737     c:= a + (n-1)*b mod p;
%p A264737     if c = 0 then return true fi;
%p A264737     b:= a; a:= c;
%p A264737     if R[a,b,(n mod p)] = 1 then return false fi;
%p A264737     R[a,b,(n mod p)]:= 1;
%p A264737   od:
%p A264737 end proc:
%p A264737 select(filter, [2,seq(i,i=3..1000,2)]); # _Robert Israel_, Nov 22 2015
%t A264737 A85 = DifferenceRoot[Function[{y, n}, {(-n - 1) y[n] - y[n + 1] + y[n + 2] == 0, y[1] == 1, y[2] == 2}]];
%t A264737 selQ[p_] := AnyTrue[Range[p - 1], Divisible[A85[#], p]&]; selQ[2] = True;
%t A264737 Reap[For[p = 2, p < 1000, p = NextPrime[p], If[selQ[p], Print[p]; Sow[p] ]]][[2, 1]] (* _Jean-François Alcover_, Jul 28 2020 *)
%Y A264737 Cf. A000085. Essentially a duplicate of A245177.
%K A264737 nonn,easy
%O A264737 1,1
%A A264737 _David Eppstein_, Nov 22 2015