Celebrating Ten Trinadays Today!

Post here to introduce yourself and say hello.
icarus
Dozens Demigod
icarus
Dozens Demigod
Joined: Apr 11 2006, 12:29 PM

Jul 6 2018, 07:22 PM #25

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.
Quote
Like
Share