5
$\begingroup$

Is there any way I can speed up this prime factor counting function? (I am looking for all numbers in a range with 3 prime factors (counted with multiplicity).)

Omega3[n_Integer] := \[Not] FreeQ[PrimeOmega[n], _?(# == 3 &)]
Omega3Count[n_] := Count[Range@n, _?Omega3]
|improve this question
$\endgroup$
  • 1
    $\begingroup$ you can start by skipping the prime numbers themselves. maybe that will speed it up a tiny bit. Your code is too advanced for me, having hard time knowing wHat this do Not[FreeQ[PrimeOmega[n], _?(# == 3 &)]] ? btw, bad idea to use UpperCase first letter for your function names. They look like build-in commands. $\endgroup$ – Nasser Dec 8 '13 at 13:09
  • $\begingroup$ OK thanks - yes, I know it is a bit convoluted - Could probably take out the Not & the FreeQ!! $\endgroup$ – martin Dec 8 '13 at 13:13
  • 1
    $\begingroup$ Yes you could, because I guess it just means Length[PrimeOmega]==3 :P $\endgroup$ – C. E. Dec 8 '13 at 13:13
  • 1
    $\begingroup$ Why not just use a simple table? Table[If[PrimeOmega[n] == 3, n, Sequence @@ {}], {n, 1, 100}] gives {8, 12, 18, 20, 27, 28, 30, 42, 44, 45, 50, 52, 63, 66, 68, 70, 75, 76, 78, 92, 98, 99} and it seems faster than what you have using a quick test. May be you can double check $\endgroup$ – Nasser Dec 8 '13 at 13:25
  • 1
    $\begingroup$ and if you just want the count, then something like: k = 0; Do[If[PrimeOmega[n] == 3, k++], {n, 1, 1000}]; $\endgroup$ – Nasser Dec 8 '13 at 13:47
3
$\begingroup$

I found same answer as ybeltukov, but a little improvment using cubic root (i see now the difference is actually significant (130 times faster than omega3count)):

co2[k_]:=Sum[1,{n,PrimePi[Power[k, (3)^-1]]},
{m,n,PrimePi[k/Prime[n]^2]},{l,m,PrimePi[k/(Prime[n]Prime[m])]}]

Result:

Timing[Omega3Count[310123]]
{14.383000000000001`,78591}

Versus

Timing[co2[310123]]
{0.10900000000000176`,78591}
|improve this answer
$\endgroup$
  • $\begingroup$ It's not 130 times faster, perhaps 3 times at most. One thing affecting timings is the order of computation, Prime and PrimePi use caching and sieving, for more information see e.g. What is so special about Prime?, +1. $\endgroup$ – Artes Dec 8 '13 at 14:16
  • $\begingroup$ I am having a little trouble comprehending the syntax of this - what values to m, n and k take? Apologies if I am being a bit slow here ... How would I implement this as Omega4Count, for example? $\endgroup$ – martin Dec 8 '13 at 14:16
  • 1
    $\begingroup$ like this: omega4[k_] := Sum[1, {n, PrimePi[Power[k, (4)^-1]]}, {m, n, PrimePi[k/Prime[n]^3]}, {l, m, PrimePi[k/(Prime[n]^2 Prime[m])]}, {j, l, PrimePi[k/(Prime[n] Prime[m] Prime[l])]}] $\endgroup$ – Coolwater Dec 8 '13 at 14:19
  • $\begingroup$ Great - thanks :) $\endgroup$ – martin Dec 8 '13 at 14:21
  • $\begingroup$ Also, If you call the function (e.g omega4) with many different integers, it will speed up very much, if you replace PrimePi with PrimePiM that is defined by PrimePiM[n_]:=PrimePiM[n]=PrimePi[n]; $\endgroup$ – Coolwater Dec 8 '13 at 14:22
6
$\begingroup$

There is a nice combination of Prime and PrimePi:

count3[n_] := Sum[1, {i, PrimePi[n]}, {j, i, PrimePi[n/Prime[i]]}, 
     {k, j, PrimePi[n/Prime[i]/Prime[j]]}];

count3[100000.] // AbsoluteTiming

{0.157486, 25556}

It is ~30 times faster:

Omega3Count[100000] // AbsoluteTiming

{4.445524, 25556}

Update

A general solution (with Coolwater's improvement)

count[k_, n_] := Sum[1, ##] & @@ 
     Transpose[{#, Prepend[Most[#], 1], PrimePi@Prepend[Prime[First[#]]^(1 - k)
           Rest@FoldList[Times, n, Prime@First[#]/Prime@Most[#]], n^(1/k)]}] 
             &@Table[Unique[], {k}];

count[3, 100000] // AbsoluteTiming

{0.048649, 25556}

|improve this answer
$\endgroup$
  • $\begingroup$ why the decimal point? in count3[100000.] ? $\endgroup$ – Nasser Dec 8 '13 at 13:59
  • $\begingroup$ @Nasser It is a bit faster because division integer by integer produce fractions. $\endgroup$ – ybeltukov Dec 8 '13 at 14:01
  • $\begingroup$ I am having a little trouble comprehending the syntax of this - what values to i, j and k take? Apologies if I am being a bit slow here ... How would I implement this as Omega4Count, for example? $\endgroup$ – martin Dec 8 '13 at 14:15
  • 1
    $\begingroup$ @martin See the general solution in my update $\endgroup$ – ybeltukov Dec 8 '13 at 14:29
  • $\begingroup$ @ ybeltukov, thank you very much for your update. I don't think it would be very nice of me to re-award the correct answer status now, but this is a really fantastic answer to the post. I wish I could give more + points for this answer. I wasn't aware of almost-prime formulas before now - so it has really opened my eyes to something new. Thanks again! :) $\endgroup$ – martin Dec 9 '13 at 14:00
2
$\begingroup$

EDIT Using Sow and Reap for general function. Mush less efficient than ybeltukov:

cnt[k_, n_] := 
 Last@Reap[Sow[1, PrimeOmega@#] & /@ Range[n], k, Total@#2 &]

Timing:

cnt[3, 100000] // AbsoluteTiming

yields:

{2.263500, {25556}}

Reassuringly same result...

ORIGINALANSWER

You could use Pick:

f[u_] := Pick[Range[u], PrimeOmega /@ Range[u], 3]

Timing[f[100]] yields:

{0., {8, 12, 18, 20, 27, 28, 30, 42, 44, 45, 50, 52, 63, 66, 68, 70,
75, 76, 78, 92, 98, 99}}

The timing for 10000: 0.187500

|improve this answer
$\endgroup$
  • $\begingroup$ cnt appears very baroque; were you practicing alternative methods? You could write: cnt[k_, n_] := PrimeOmega@Range@n ~Count~ k $\endgroup$ – Mr.Wizard♦ Dec 19 '13 at 7:07
  • $\begingroup$ Thank you for feedback. Agree my approach was not deep thinking or efficient...a combination of my own limitations and being time poor but still wanting to participate... $\endgroup$ – ubpdqn Dec 20 '13 at 0:41
  • $\begingroup$ Thanks for accepting my feedback. What you wrote is actually quite advanced and has strong potential in other applications. It can conserve memory compared to Count (especially if Scan is used in place of Map) and it is easily adaptable to efficient individual counts of several different elements. Here however it seems to be overkill when Count is available. $\endgroup$ – Mr.Wizard♦ Dec 20 '13 at 2:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.