######################################################### # # # Coded by Indranil Ghosh (indranilg49@gmail.com) # # # ######################################################### #Python 2.7.11, OEIS sequence: A265350 from sympy import factorint, prime, primefactors from operator import mul import collections def a007623(n, p=2): return n if n<p else a007623(int(n/p), p+1)*10 + n%p def a001222(n): return 0 if n==1 else a001222(n/primefactors(n)[0]) + 1 def a056169(n): f=factorint(n) return 0 if n==1 else sum([1 for i in f if f[i]==1]) def a275812(n): return a001222(n) - a056169(n) def a275735(n): y=collections.Counter(map(int, list(str(a007623(n)).replace("0", "")))).most_common() return 1 if n==0 else reduce(mul, [prime(y[i][0])**y[i][1] for i in xrange(len(y))]) def a(n): return a275812(a275735(n)) print [n for n in xrange(501) if a(n)>0]