5
$\begingroup$

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?

$\endgroup$
0

5 Answers 5

5
$\begingroup$

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

$\endgroup$
14
  • $\begingroup$ That's great, but I want to extend it to non squarefree also. $\endgroup$
    – martin
    Commented Oct 9, 2013 at 18:01
  • $\begingroup$ 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. $\endgroup$
    – Artes
    Commented Oct 9, 2013 at 18:02
  • $\begingroup$ @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? $\endgroup$
    – Artes
    Commented Oct 9, 2013 at 18:09
  • $\begingroup$ @martin Could you explain what you mean by extension to non-squarefree numbers, and clarify what kind of plot you are looking for? $\endgroup$
    – Artes
    Commented Oct 9, 2013 at 23:10
  • $\begingroup$ 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]]]}] $\endgroup$
    – martin
    Commented Oct 9, 2013 at 23:17
3
$\begingroup$

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
$\endgroup$
3
  • $\begingroup$ 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! $\endgroup$
    – martin
    Commented Oct 9, 2013 at 18:03
  • $\begingroup$ @martin. See my edit. $\endgroup$
    – RunnyKine
    Commented Oct 9, 2013 at 18:05
  • $\begingroup$ Yes! Perfect! Many thanks! $\endgroup$
    – martin
    Commented Oct 9, 2013 at 18:24
3
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ OK, yes, got that thanks. $\endgroup$
    – martin
    Commented Oct 10, 2013 at 17:31
1
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Artes
    Commented Oct 9, 2013 at 23:38
  • $\begingroup$ Nearly the same - only slightly different: $\endgroup$
    – martin
    Commented Oct 9, 2013 at 23:56
  • $\begingroup$ 1st plot is square free, second has one square + a few other things, third, ...& so on $\endgroup$
    – martin
    Commented Oct 10, 2013 at 0:03
  • $\begingroup$ 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. $\endgroup$
    – Artes
    Commented Oct 10, 2013 at 0:05
  • $\begingroup$ 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 $\endgroup$
    – martin
    Commented Oct 10, 2013 at 0:09
1
$\begingroup$

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}
   ]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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