(PARI) a(p) = {my(v=List(vector(p-1, i, i)), m=p, h=0); for(i=1, p-1, h+=1/i; if(numerator(h)%p==0, for(j=i*p, i*p+p-1, listput(v, j)))); while(m<#v, h+=sum(i=v[m-1]+1, v[m]-1, 1/i); for(i=m, m+p-1, h+=1/v[i]; if(numerator(h)%p==0, for(j=v[i]*p, v[i]*p+p-1, listput(v, j)))); m+=p); v[#v]; } lista(nn) = {my(w=vector(nn)); forprime(p=2, nn, w[p]=a(p)); forprime(p=2, sqrtint(nn), for(e=2, logint(nn, p), w[p^e]=p^(e-1)*(w[p]+1)-1)); for(n=6, nn, f=factor(n); w[n]=vecmax(vector(#f[, 1], i, w[f[i, 1]^f[i, 2]]))); w; } Note: if a(p) is calculated for more than one hour, it is possible that a(p) = 0 and the calculation may continue indefinitely. Quickly calculate whether A002805(k) is not divided by prime p: is(k, p) = {my(m=k\p); Mod(numerator(sum(i=1, m\p, 1/i)/p + sum(j=1+m-m%p, m, 1/j)), p) == 0; } n = 2: [1] n = 3: [1, 2, 6, 7, 8, 21, 22, 23, 66, 67, 68] n = 4: [1, 2, 3] n = 5: [1, 2, 3, 4, 20, 21, 22, 23, 24, 100, 101, 102, 103, 104, 120, 121, 122, 123, 124] n = 6: [1, 2, 6, 7, 8, 21, 22, 23, 66, 67, 68] n = 7: [1, 2, 3, 4, 5, 6, 42, 43, 44, 45, 46, 47, 48, 294, 295, 296, 297, 298, 299, 300, 336, 337, 338, 339, 340, 341, 342, 2065, 2066, 2067, 2068, 2069, 2070, 2071, 2093, 2094, 2095, 2096, 2097, 2098, 2099, 2359, 2360, 2361, 2362, 2363, 2364, 2365, 2387, 2388, 2389, 2390, 2391, 2392, 2393, 14672, 14673, 14674, 14675, 14676, 14677, 14678, 16730, 16731, 16732, 16733, 16734, 16735, 16736, 102725, 102726, 102727, 102728, 102729, 102730, 102731, 117117, 117118, 117119, 117120, 117121, 117122, 117123, 117145, 117146, 117147, 117148, 117149, 117150, 117151, 719096, 719097, 719098, 719099, 719100, 719101, 719102] n = 8: [1, 2, 3, 4, 5, 6, 7] n = 9: [1, 2, 3, 4, 5, 6, 7, 8, 18, 19, 20, 21, 22, 23, 24, 25, 26, 63, 64, 65, 66, 67, 68, 69, 70, 71, 198, 199, 200, 201, 202, 203, 204, 205, 206] n = 10: [1, 2, 3, 4, 20, 21, 22, 23, 24, 100, 101, 102, 103, 104, 120, 121, 122, 123, 124]