5

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]
martin
  • 8,678
  • 4
  • 23
  • 70
  • 1
    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. – Nasser Dec 08 '13 at 13:09
  • OK thanks - yes, I know it is a bit convoluted - Could probably take out the Not & the FreeQ!! – martin Dec 08 '13 at 13:13
  • 1
    Yes you could, because I guess it just means Length[PrimeOmega]==3 :P – C. E. Dec 08 '13 at 13:13
  • Yes! Thanks! :) – martin Dec 08 '13 at 13:14
  • 1
    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 – Nasser Dec 08 '13 at 13:25
  • or use Do as in Reap@Do[If[PrimeOmega[n] == 3, Sow[n]], {n, 1, 100}] which is fast also. Not as fancy as your code though but Do is really fast, and I think it qualifies sort of as being functional programming, but may be not pure functional, just a little bit functional. – Nasser Dec 08 '13 at 13:27
  • @Nasser - Great thanks - will try now – martin Dec 08 '13 at 13:34
  • 1
    and if you just want the count, then something like: k = 0; Do[If[PrimeOmega[n] == 3, k++], {n, 1, 1000}]; – Nasser Dec 08 '13 at 13:47
  • @martin It's a bad idea accepting answers after one hour because it discourages others from looking for better solutions than those provided so far. It is better to wait one or two days. – Artes Dec 08 '13 at 15:18
  • @ Artes, I will take that on board for next time. I know this is your area of expertise, and I would like to have seen your contribution - & I realise I shouldn't have accepted an answer so readily, but on a positive note, I have been really pleased by the quality of some of the solutions, an should be very interested if you have anything to add to the post :) – martin Dec 09 '13 at 14:03

3 Answers3

7

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}

ybeltukov
  • 43,673
  • 5
  • 108
  • 212
  • why the decimal point? in count3[100000.] ? – Nasser Dec 08 '13 at 13:59
  • @Nasser It is a bit faster because division integer by integer produce fractions. – ybeltukov Dec 08 '13 at 14:01
  • 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? – martin Dec 08 '13 at 14:15
  • 1
    @martin See the general solution in my update – ybeltukov Dec 08 '13 at 14:29
  • @ 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! :) – martin Dec 09 '13 at 14:00
3

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}
Coolwater
  • 20,257
  • 3
  • 35
  • 64
  • 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. – Artes Dec 08 '13 at 14:16
  • 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? – martin Dec 08 '13 at 14:16
  • 1
    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])]}] – Coolwater Dec 08 '13 at 14:19
  • Great - thanks :) – martin Dec 08 '13 at 14:21
  • 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];

     – Coolwater Dec 08 '13 at 14:22
  • @ Coolwater - many thanks for the update - it would be nice to include a reference to the almost-prime formulas in your answer. I will include it in this comment for anyone interested: Almost Primes – martin Dec 09 '13 at 14:07
2

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

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • cnt appears very baroque; were you practicing alternative methods? You could write: cnt[k_, n_] := PrimeOmega@Range@n ~Count~ k – Mr.Wizard Dec 19 '13 at 07:07
  • 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... – ubpdqn Dec 20 '13 at 00:41
  • 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. – Mr.Wizard Dec 20 '13 at 02:38