cp's OEIS Frontend

A266815 Primes whose sum of reciprocal of digits is a prime.

%I A266815 #15 May 16 2019 08:53:22
%S A266815 11,2441,3313,3331,4241,4421,12163,12613,13313,13331,16231,16363,
%T A266815 16633,21163,21613,26113,31663,32233,32323,32611,33113,33223,33311,
%U A266815 48281,48821,61231,61363,62131,62311,63211,63361,88241,112121,114643,116443,122263,123323,126223
%N A266815 Primes whose sum of reciprocal of digits is a prime.
%C A266815 Subset of A034708.
%H A266815 Paolo P. Lava, <a href="/A266815/b266815.txt">Table of n, a(n) for n = 1..1000</a>
%e A266815 11: 1/1 + 1/1 = 2;
%e A266815 2441: 1/2 + 1/4 + 1/4 + 1/1 = 2;
%e A266815 3313: 1/3 + 1/3 + 1/1 + 1/3 = 2; etc.
%p A266815 P:=proc(q) local a,k,n,ok;
%p A266815 for n from 1 to q do if isprime(n) then ok:=1; a:=0; for k from 0 to ilog10(n) do
%p A266815 if trunc(n/10^k) mod 10>0 then a:=a+1/(trunc(n/10^k) mod 10) else ok:=0; break; fi; od;
%p A266815 if ok=1 and type(a,integer) then if isprime(a) then print(n) fi; fi; fi; od; end: P(10^9);
%p A266815 # Alternative:
%p A266815 N:= 8: # to get all terms of up to N digits
%p A266815 S1:= proc(t,k,N)
%p A266815     option remember;
%p A266815     if t = 0 then {[]}
%p A266815     elif k = 0 then {}
%p A266815     else
%p A266815       `union`(seq(map(p -> [op(p),k$m], procname(t - m*2520/k, k-1,N-m)),
%p A266815          m = 0 .. min(N, floor(t*k/2520))))
%p A266815     fi
%p A266815 end proc:
%p A266815 targets:= 2520*select(isprime,[$2..N]):
%p A266815 Dlists:= select(p -> convert(p,`+`) mod 3 <> 0, `union`(seq(S1(t,9,N),t=targets))):
%p A266815 g:= proc(L) local i,m;
%p A266815    m:= nops(L);
%p A266815    op(select(isprime, map(t -> add(t[i]*10^(i-1),i=1..m), combinat:-permute(L))));
%p A266815 end proc:
%p A266815 sort(convert(map(g, Dlists),list)); # _Robert Israel_, Feb 12 2016
%t A266815 Select[Prime@ Range@ 12000, If[MemberQ[#, 0], False, PrimeQ@ Total[1/#]] &@ IntegerDigits@ # &] (* _Michael De Vlieger_, Feb 12 2016 *)
%o A266815 (PARI) isok(n) = if (isprime(n), my(d = digits(n)); vecmin(d) && (denominator(s=sum(k=1, #d, 1/d[k])) == 1) && isprime(s)) \\ _Michel Marcus_, Feb 12 2016
%Y A266815 Cf. A000040, A034708.
%K A266815 nonn,base
%O A266815 1,1
%A A266815 _Paolo P. Lava_, Feb 12 2016