## Prime factor counting function

5

1

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]


1you 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 – 2013-12-08T13:09:39.080

OK thanks - yes, I know it is a bit convoluted - Could probably take out the Not & the FreeQ!! – martin – 2013-12-08T13:13:02.773

1Yes you could, because I guess it just means Length[PrimeOmega]==3 :P – C. E. – 2013-12-08T13:13:53.803

Yes! Thanks! :) – martin – 2013-12-08T13:14:27.687

1Why 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 – 2013-12-08T13:25:15.050

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 – 2013-12-08T13:27:54.720

@Nasser - Great thanks - will try now – martin – 2013-12-08T13:34:43.187

1and if you just want the count, then something like: k = 0; Do[If[PrimeOmega[n] == 3, k++], {n, 1, 1000}]; – Nasser – 2013-12-08T13:47:12.660

@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 – 2013-12-08T15:18:29.117

@ 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 – 2013-12-09T14:03:17.080

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}


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 – 2013-12-08T14:16:19.573

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 – 2013-12-08T14:16:30.097

1like 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 – 2013-12-08T14:19:37.067

Great - thanks :) – martin – 2013-12-08T14:21:09.387

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 – 2013-12-08T14:22:01.607

@ 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 – 2013-12-09T14:07:53.903

6

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}

why the decimal point? in count3[100000.] ? – Nasser – 2013-12-08T13:59:12.023

@Nasser It is a bit faster because division integer by integer produce fractions. – ybeltukov – 2013-12-08T14:01:59.950

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 – 2013-12-08T14:15:39.607

1@martin See the general solution in my update – ybeltukov – 2013-12-08T14:29:00.063

@ 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 – 2013-12-09T14:00:25.690

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...

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

cnt appears very baroque; were you practicing alternative methods? You could write: cnt[k_, n_] := PrimeOmega@Range@n ~Count~ k – Mr.Wizard – 2013-12-19T07:07:28.183

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 – 2013-12-20T00:41:30.880

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 – 2013-12-20T02:38:41.477