Post here to introduce yourself and say hello.
Dozens Demigod
icarus
Dozens Demigod
Joined: Apr 11 2006, 12:29 PM
Ah, well in that case we can take the index n of A246547(n) and there, we have such a function.

We can write a more efficient program than the Mathematica in that sequence. This one gives 10^7 terms nearly instantly, but the programs there search all numbers, testing them for PrimePowerQ and Not-PrimeQ, requiring a couple minutes to find the same terms. We can use logarithm approach. Since we are not interested in primes, we consider the range p = prime(1 .. π(√n)), where π(x) is the prime counting function. Then for each prime p, write the array p^k for 2 ≤ k ≤ log_p n, sorting the aggregated results. Thereby, we derive all terms of A246547 less than or equal to n.

Code: Select all

``````With[{n = 15000},
Union@ Flatten@
Table[Array[p^# &, Floor@ Log[p, n] - 1, 2], {p,
Prime@ Range@ PrimePi@ Sqrt@ n}] ]``````
This is a good fast program so I am going to contribute it there! It generated 80070 terms (A246547 for all n ≤ 10^12) in less than a second!

Now having a dataset that large would then allow us to merely count the number of terms less than or equal to an arbitrary number. Suppose we choose 9. Then the index of 9 in {4, 8, 9, 16, 25, 27, ...} is 3. In Mathematica we would use LengthWhile[s, # <= 9 &] and that would yield 3. If we tried 12, we could see that there still are 3 terms in that sequence less than 12, so f(12) = 3.

Written in uncial (dozenal), here are the 102; terms less than 100,000; :

Code: Select all

``4, 8, 9, 14, 21, 23, 28, 41, 54, 69, a1, a5, a8, 121, 183, 194, 201, 247, 261, 368, 381, 441, 509, 5a1, 681, 714, 92b, 961, b81, 10a1, 1228, 1323, 1331, 1341, 1481, 1761, 1985, 2021, 21a1, 2454, 2721, 2a15, 2b01, 3101, 3741, 3969, 3b77, 3ba1, 4701, 48a8, 5541, 5aa1, 6181, 6761, 6a61, 705b, 7481, 8581, 9061, 9401, 9594, 9887, 9b21, aa41, b221, b483, 10a21, 11241, 12145, 12321, 13461, 14181, 14641, 152a7, 153a1, 16661, 16b61, 16b68, 19141, 19681, 1a561, 1ab01, 21921, 24941, 25391, 259a1, 26421, 27501, 29081, 29741, 2a209, 30561, 31b14, 32281, 33a75, 34041, 35a61, 36601, 384a1, 39265, 39841, 3a017, 3a421, 40401, 41821, 46661, 47b81, 48841, 4a1a1, 500bb, 534a1, 55881, 58101, 59821, 5a5a1, 60141, 62701, 63501, 63a28, 65b41, 68621, 6b161, 70a81, 721a5, 736a1, 77261, 79081, 7924b, 80981, 85721, 86623, 866a1, 8b601, 90601, 93641, 956a1, 98801, 9aa2b, a0a41, a2ba1, a4081, a6261, ab431, b0941, b5301, b7621``
Of these, two dozen two are regular to the dozen; these are the powers m of 2 and 3 less than 100,000 such that m ≥ 2.

Code: Select all

``4, 8, 9, 14, 23, 28, 54, 69, a8, 183, 194, 368, 509, 714, 1228, 1323, 2454, 3969, 48a8, 9594, b483, 16b68, 2a209, 31b14, 63a28, 86623``
We can see that the non-regular terms end in -1, -5, -7, or -b as they are coprime to the dozen (the highest common factor of the term and the dozen is 1).

Alternatively, we can take the integer part of log_p n for any prime p ≤ √n then add all the terms and get f(n)

{16,10,6,5,4,3,3,3,2,2,2,2,2,2,2,2,2,2,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1} = 146

Thus the function is easy and relatively efficient:

Code: Select all

``f[n_] := Total@ Table[Floor@ Log[p, n] - 1, {p, Prime@ Range@ PrimePi@ Sqrt@ n}]``
This computes f(10^12) = 80070 in under a second.