OFFSET
2,1
COMMENTS
For the smallest odd prime not generated, see A075227. For information about how often the numerator of these sums is prime, see A075188 and A075189. The Mathematica program also prints the subset that yields the largest prime. For n <=20, the largest prime occurs in a sum of n-2, n-1, or n reciprocals.
LINKS
Chai Wah Wu, Table of n, a(n) for n = 2..441 (terms 2..100 from Martin Fuller)
Martin Fuller, PARI program
EXAMPLE
a(3) =11 because 11 is largest prime numerator in the three sums that yield primes: 1+1/2 = 3/2, 1/2+1/3 = 5/6 and 1+1/2+1/3 = 11/6.
MATHEMATICA
Needs["DiscreteMath`Combinatorica`"]; maxN=20; For[t={}; lst={}; mx=0; i=0; n=2, n<=maxN, n++, While[i<2^n-1, i++; s=NthSubset[i, Range[n]]; k=Numerator[Plus@@(1/s)]; If[PrimeQ[k], If[k>mx, t=s]; mx=Max[mx, k]]]; Print[n, " ", t]; AppendTo[lst, mx]]; lst
Table[Max[Select[Numerator[Total/@Subsets[1/Range[n], {2, 2^n}]], PrimeQ]], {n, 2, 30}] (* The program will take a long time to run. *) (* Harvey P. Dale, Jan 08 2019 *)
PROG
(Haskell)
import Data.Ratio (numerator)
a075226 n = a075226_list !! (n-1)
a075226_list = f 2 [recip 1] where
f x hs = (maximum $ filter ((== 1) . a010051') (map numerator hs')) :
f (x + 1) hs' where hs' = hs ++ map (+ recip x) hs
-- Reinhard Zumkeller, May 28 2013
(PARI) See Fuller link.
(Python)
from math import gcd, lcm
from itertools import combinations
from sympy import isprime
def A075226(n):
m = lcm(*range(1, n+1))
c, mlist = 0, tuple(m//i for i in range(1, n+1))
for l in range(n, -1, -1):
if sum(mlist[:l]) < c:
break
for p in combinations(mlist, l):
s = sum(p)
s //= gcd(s, m)
if s > c and isprime(s):
c = s
return c # Chai Wah Wu, Feb 14 2022
CROSSREFS
KEYWORD
nice,nonn
AUTHOR
T. D. Noe, Sep 08 2002
EXTENSIONS
More terms from Martin Fuller, Jan 19 2008
STATUS
approved