|
%I #23 Jul 28 2020 08:48:02
%S 2,5,13,19,23,29,31,43,53,59,61,67,73,79,83,89,97,103,131,137,151,157,
%T 163,173,179,181,191,197,199,211,229,233,239,241,281,293,307,317,347,
%U 359,367,373,379,389,397,409,419,421,431,433,443,449,457,461,463,479,487,491,499
%N Primes which divide some term of A000085 (numbers of involutions).
%C Essentially the same as A245177. - _R. J. Mathar_, Nov 25 2015
%H Robert Israel, <a href="/A264737/b264737.txt">Table of n, a(n) for n = 1..4000</a>
%F 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 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 filter:= proc(p) local a,b,c,n,R;
%p if not isprime(p) then return false fi;
%p a:= 1; b:= 1;
%p R[1,1,1]:= 1;
%p for n from 2 do
%p c:= a + (n-1)*b mod p;
%p if c = 0 then return true fi;
%p b:= a; a:= c;
%p if R[a,b,(n mod p)] = 1 then return false fi;
%p R[a,b,(n mod p)]:= 1;
%p od:
%p end proc:
%p select(filter, [2,seq(i,i=3..1000,2)]); # _Robert Israel_, Nov 22 2015
%t A85 = DifferenceRoot[Function[{y, n}, {(-n - 1) y[n] - y[n + 1] + y[n + 2] == 0, y[1] == 1, y[2] == 2}]];
%t selQ[p_] := AnyTrue[Range[p - 1], Divisible[A85[#], p]&]; selQ[2] = True;
%t 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 Cf. A000085. Essentially a duplicate of A245177.
%K nonn,easy
%O 1,1
%A _David Eppstein_, Nov 22 2015
|
