######################################################### # # # Coded by Indranil Ghosh (indranilg49@gmail.com) # # # ######################################################### #Python 2.7.11, OEIS sequence: A276076 from sympy import factorial as f, factorint, prime, primepi from operator import mul def a007623(n, p=2): return n if n<p else a007623(int(n/p), p+1)*10 + n%p def a003961(n): f=factorint(n) return 1 if n==1 else reduce(mul, [prime(primepi(i) + 1)**f[i] for i in f]) def a276009(n): x=str(a007623(n)) y="".join([str(int(i) - 1) if int(i)>0 else '0' for i in x])[::-1] return sum([int(y[i])*f(i + 1) for i in xrange(len(y))]) def a257684(n): x=str(a007623(n))[:-1] y="".join([str(int(i) - 1) if int(i)>0 else '0' for i in x])[::-1] return 0 if n==1 else sum([int(y[i])*f(i + 1) for i in xrange(len(y))]) def a275732(n): x=str(a007623(n))[::-1] return 1 if n==0 or x.count("1")==0 else reduce(mul, [prime(i + 1) for i in xrange(len(x)) if x[i]=='1']) def a275733(n): return 1 if n==0 else a275732(n)*a003961(a275733(a257684(n))) def a(n): return 1 if n==0 else a275733(n)*a(a276009(n)) print [a(n) for n in xrange(201)]