%I #28 Aug 12 2024 04:33:57
%S 1,3,7,13,31,39,91,93,127,217,307,381,403,921,961,1093,1209,1651,1723,
%T 2149,2801,2821,3279,3541,3937,3991,4953,5113,5169,7651,8011,8191,
%U 8403,9517,10303,10623,11811,11973,12061,12493,15339,17293,19531,19607,22399
%N Numbers that have all their divisors in A002191 (possible values for sigma(n), A000203).
%C Since 2 does not belong to A002191, all terms are odd.
%C All primes p that are in A023195 (Prime numbers that are the sum of the divisors of some n), are also in this sequence; and the prime factors of all terms can only belong to A023195.
%C Up to 10^7, only one term is a prime power: 961=31^2 (being a square, see A038688, A228061 and A243810).
%H Amiram Eldar, <a href="/A243765/b243765.txt">Table of n, a(n) for n = 1..2000</a>
%e The divisors of 3 are 1 and 3 that both belong to A002191, 1 as sigma(1) and 3 as sigma(2).
%e The divisors of 39 are 1, 3, 13 and 39 all of which belong to A002191, 13 as sigma(9) 39 as sigma(18).
%p N:= 10^6: # to get all terms up to N
%p A002191:= select(`<=`,{seq(numtheory[sigma](i),i=1..N)},N):
%p A243765:= select(t -> numtheory[divisors](t) subset A002191, A002191); # _Robert Israel_, Jun 16 2014
%o (PARI) list(lim) = select(n->n<=lim, Set(vector(lim\=1, n, sigma(n))));
%o isok(n, lists) = {fordiv (n, d, if (!vecsearch(lists, d), return(0))); return(1);}
%o lista(nn) = {lists = list(nn); for(n=1, nn, if (isok(n, lists), print1(n, ", ")););}
%Y Cf. A000203, A002191, A023195.
%Y Cf. A045572 (analog sequence with the sum of proper divisors instead).
%K nonn
%O 1,2
%A _Michel Marcus_, Jun 10 2014