Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Update

I would now like to produce a plot that explores the difference between the squarefree count and the expected mean (without using SquareFreeQ) - something along the lines of:

a = PrimeOmega[Range[1000]];
b = PrimeNu[Range[1000]];
ListLinePlot[{#/Zeta[2] - Accumulate[Flatten[Inner[If[#1 === #2, 1, 0] 
&, a, b, List]]]}, Filling -> Axis]

or

a = PrimeOmega[Range[100]];
b = PrimeNu[Range[100]];
Plot[{x/Zeta[2] - Accumulate[Flatten[Inner[If[#1 === #2, 1, 0] &, a, b, List]]]}, 
{x, 0, 100}, Filling -> Axis]

In the hope of achieving something similar to:

Plot[{{x/Log[x]} - PrimePi[x]}, {x, 0, 100}, Filling -> Axis]

Or alternatively, slightly alter Artes' cf plot (see below).

Instead of:

cf = {#, Count[Range[#, # + 500], _?SquareFreeQ]} & /@ Range[0, 100000, 500]

Using something like:

a = PrimeOmega[Range[100000]];
b = PrimeNu[Range[100000]];
cf = {#, Count[Range[#, # + 500], _?Accumulate[Flatten[Inner[If[#1 === #2, 1, 0] 
&, a, b, List]]]]} & /@ Range[0, 100000, 500]

... Much clumsier, I realise, but I would like to avoid using SquareFreeQ if I can.

Old

I would like mathematica to output instances where PrimeOmega[x]=PrimeNu[x] up to a given range. (Clearly in this example, the output will be the set of square free numbers.)

a = PrimeOmega[Range[1000000]];
b = PrimeNu[Range[1000000]];
Count[Transpose[{a, b}], {1, 1}]
a = PrimeOmega[Range[1000000]];
b = PrimeNu[Range[1000000]];
Count[Transpose[{a, b}], {2, 2}] .....

works, but is (a) very slow, and (b) very clumsy.

I was wondering if there was a more succinct way of putting it?

In addition, I would like to ultimately create a plot for what will essentially be a counting function for PrimeOmega[x]&&PrimeNu[x]=={1, 1},{2, 2}... etc. Is there any way of doing this?

share|improve this question
add comment

5 Answers

up vote 5 down vote accepted

Instead of using both PrimeOmega and PrimeNu I'd rather use only SquareFreeQ.

Let's compare appropriate timings:

First @ AbsoluteTiming[ a = PrimeOmega[Range[300000]];
                b = PrimeNu[Range[300000]]; 
                Inner[If[#1 === #2, True, False] &, a, b, List];] 

First @ AbsoluteTiming[SquareFreeQ /@ Range[300000];]
19.748000
 1.521000 

and of course:

Inner[ If[ #1 === #2, True, False] &, a, b, List] == (SquareFreeQ /@ Range[300000])
True

Edit

If we are to find numbers which satisfy PrimeOmega[x] == PrimeNu[x]] in a given range we can use Select, e.g.

Select[ Range[10^6, 10^6 + 11], SquareFreeQ]
{1000001, 1000002, 1000003, 1000005, 1000006, 1000007, 1000009, 
 1000010, 1000011, 1000013, 1000014, 1000015, 1000018, 1000019}

to count them we use:

Count[ Range[10^6, 10^6 + 20], _?SquareFreeQ]
14

of course we might use Count[Range[10^6, 10^6 + 20], _?(PrimeOmega[#] == PrimeNu[#] &)] instead but the latter is slower.
Let's define appropriate counting function:

cf = {#, Count[Range[#, # + 500], _?SquareFreeQ]} & /@ Range[0, 100000, 500];

Namely we count square free numbers in every range: {0, 500}, {500, 1000},...,{99500, 100000}.

Let's plot the counting function cf:

With[{ mcf = Mean @ cf[[All, 2]]}, 
     ListPlot[ cf, AxesOrigin -> {0, 280}, PlotRange -> {280, 320}, AspectRatio -> 1/5, 
               PlotMarkers -> Automatic, Filling -> mcf, 
               Epilog -> {Darker @ Green, Line[{{0, mcf}, {100000, mcf}}]}]]

enter image description here

share|improve this answer
    
That's great, but I want to extend it to non squarefree also. –  martin Oct 9 '13 at 18:01
    
To ensure that mapping SquareFreeQ is better, try to find timings of generating lists a and b. It takes: AbsoluteTiming[a = PrimeOmega[Range[300000]]; b = PrimeNu[Range[300000]];] yields 19.615000 while AbsoluteTiming[Inner[...]] only 0.503000. –  Artes Oct 9 '13 at 18:02
    
@martin What do you mean by non-squarefree??? I've just provided what you've been looking for. The main problem is generating lists of PrimeOmega and PrimeNu, instead you can just play with SquareFreeQ. Did you miss anything? –  Artes Oct 9 '13 at 18:09
    
@martin Could you explain what you mean by extension to non-squarefree numbers, and clarify what kind of plot you are looking for? –  Artes Oct 9 '13 at 23:10
    
Sorry, yes - was looking for this kind of thing: a = PrimeOmega[Range[1000]];b = PrimeNu[Range[1000]]; ListPlot[{Accumulate[ Flatten[Inner[If[#1 === #2, 1, 0] &, a, b, List]]], Accumulate[Flatten[Inner[If[#1 === #2 + 1, 1, 0] &, a, b, List]]], Accumulate[Flatten[Inner[If[#1 === #2 + 2, 1, 0] &, a, b, List]]]}] –  martin Oct 9 '13 at 23:17
show 9 more comments

Is this what you want:

Reap[Sow[#, PrimeOmega[#] == PrimeNu[#]] & /@ Range[300000];, True][[2,1]] // Length

OR using just SquareFreeQ which is obviously faster

Reap[Sow[#, SquareFreeQ[#]] & /@ Range[300000];, True][[2, 1]] // Length

Here is another way based on Inner

a = PrimeOmega @ Range[300000]; b = PrimeNu @ Range[300000];
Inner[If[#1 === #2, True, ## &[]] &, a, b, List] // Length
share|improve this answer
    
Outputs list, rather than counting number of times PrimeOmega[n] is equal to PrimeNu[n] within a given range. Output should be single number. - Sorry, obviously I am not being very clear! –  martin Oct 9 '13 at 18:03
    
@martin. See my edit. –  RunnyKine Oct 9 '13 at 18:05
    
Yes! Perfect! Many thanks! –  martin Oct 9 '13 at 18:24
add comment

In general, to count how many entries in your two lists are the same, you can use

a = PrimeOmega[Range[1000000]];
b = PrimeNu[Range[1000000]];

Length@a - (Unitize[a - b] // Total)

607926

The last line, the counting line, will be fast.

However for the underlying problem, I will rephrase what @Artes has pointed out already: PrimeOmega[x] == PrimeNu[x] if and only if x is square-free. So SquareFree will be a faster way to solve the whole problem, as Artes has shown.


Response to updated question

Perhaps this is what you're after:

a = PrimeOmega[Range[100000]];
b = PrimeNu[Range[100000]];
sums = Accumulate[1 - Unitize[a - b]];

Plot[x/Zeta[2] - sums[[Min[1 + Floor[x], Length@sums]]], {x, 0, 100000},
 Filling -> Axis]

Mathematica graphics

share|improve this answer
    
OK, yes, got that thanks. –  martin Oct 10 '13 at 17:31
add comment

Update

OK - First part of update solved:

a = PrimeOmega[Range[100000]];
b = PrimeNu[Range[100000]];
c = Range[100000];
ListLinePlot[{c/Zeta[2] - Accumulate[Flatten[Inner[If[#1 === #2, 1, 0] 
&, a, b, List]]]}, PlotStyle -> Red]

enter image description here

... Just unsure how to modify Artes' count plot.

Old

Sorry Artes, yes - was looking for this kind of thing:

a = PrimeOmega[Range[100]];
b = PrimeNu[Range[100]];
ListPlot[{Accumulate[Flatten[Inner[If[#1 === #2, 1, 0] &, a, b, List]]],    
Accumulate[Flatten[Inner[If[#1 === #2 + 1, 1, 0] &, a, b, List]]],
Accumulate[Flatten[Inner[If[#1 === #2 + 2, 1, 0] &, a, b, List]]]}]

enter image description here

share|improve this answer
    
You can edit your questions or answers to improve them. But I think you've included this plot in another answer, i.e. here, posting duplicate answers or questions is not (in general) appropriate. –  Artes Oct 9 '13 at 23:38
    
Nearly the same - only slightly different: –  martin Oct 9 '13 at 23:56
    
1st plot is square free, second has one square + a few other things, third, ...& so on –  martin Oct 10 '13 at 0:03
    
I meant you posted another question where it was more clearly stated what were your expectations. From the above question it is not clear what is the essence. –  Artes Oct 10 '13 at 0:05
    
Other thread was asking about plot comparing semiprime & prime counting functions; this one was about squarefree, numbers with square prime factors, numbers with cubed prime factors ,etc. counting functions –  martin Oct 10 '13 at 0:09
show 1 more comment

Late, but a slightly different approach avoidingSquareFreeQ, PrimeOmega, and PrimeNu. See Sloane's A143658, the number of squarefree integers not exceeding $2^n$. Fast, since the sum is only to the square root of the upper limit; however, it only has samples at the powers of 2.

SloanesA143658[n_] := 
   Module[{t = 2^n}, 
          ParallelSum[MoebiusMu[k] Floor[t/k^2], {k, 1, Floor[Sqrt[t]]}]]

Plot the count as dots with the theoretical line, $2^n/\zeta(2)$, in red.

With[{nmax = 20}, 
   ListLogPlot[Table[{n, SloanesA143658[n]}, {n, 1, nmax}], 
               Frame -> True, 
               FrameLabel -> {"Exponent  n", "Number of SquareFree"}, 
               PlotLabel -> "SquareFree Numbers <= 2^n", BaseStyle -> {FontSize -> 14}, 
               Epilog -> {Red, Line[Table[{n, Log[2^n/Zeta[2]]}, {n, 1, nmax}]]}
   ]]

Plot the relative error of the theoretical value versus the count.

With[{nmax = 20}, 
   ListLogPlot[Table[{n, Abs[2^n/Zeta[2] - SloanesA143658[n]]/2^n}, {n, 1, nmax}], 
               Frame -> True, Joined -> True, 
               FrameLabel -> {"Exponent  n", "Relative Error"}, 
               PlotLabel -> "SquareFree Numbers <= 2^n", BaseStyle -> {FontSize -> 14}
   ]]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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