OFFSET
1,1
COMMENTS
The general case is (b,b+2*k) to give |sopfr(b)-sopfr(b+2*k)|=2*k=p+2*k-p for prime pair (p,p+2*k). For various k compare all pairs (p,p+2*k)<10^m to determine if the percentage of matching (b,b+2*k)is a maximum for some value of k. For k=1 and m=6, there are 15 of 8169 twin primes that had matching (b,b+2), or only 0.1836%, about 1 out of 545.
LINKS
Chai Wah Wu, Table of n, a(n) for n = 1..222
EXAMPLE
sopfr(845)=5+13+13=31 and sopfr(847)=7+11+11=29 and (29,31) are twin primes.
MATHEMATICA
sopfr[n_]:=Total[Flatten[Table[#[[1]], #[[2]]]&/@FactorInteger[n]]]; Select[ Range[ 673*10^4], CompositeQ[#]&&CompositeQ[#+2]&&AllTrue[ {sopfr[ #], sopfr[#+2]}, PrimeQ]&&Abs[sopfr[#]-sopfr[#+2]]==2&] (* Requires Mathematica version 10 or later *) (* Harvey P. Dale, Apr 16 2018 *)
PROG
(PARI) sopfr(n)=my(f=factor(n)); sum(i=1, #f~, f[i, 1]*f[i, 2]);
is(m, n)=my(p, q); !isprime(m) && !isprime(n) && isprime(p=sopfr(m)) && isprime(q=sopfr(n)) && abs(p-q)==2;
for(n=1, 10^6, if(is(n, n+2), print1(n", ")))
\\ Charles R Greathouse IV, Sep 08 2014
(Python)
from sympy import isprime, factorint
A247048_list = []
for n in range(1, 10**6):
....if not(isprime(n) or isprime(n+2)):
........m = sum(p*e for p, e in factorint(n).items())
........if isprime(m):
............m2 = sum(p*e for p, e in factorint(n+2).items())
............if ((m2 == m+2) or (m == m2+2)) and isprime(m2):
................A247048_list.append(n)
# Chai Wah Wu, Sep 25 2014
CROSSREFS
KEYWORD
nonn
AUTHOR
J. M. Bergot, Sep 10 2014
EXTENSIONS
More terms from Charles R Greathouse IV, Sep 10 2014
STATUS
approved