Given a list of primes, how can I plot π(x)?

Question

Given an arbitrarily-long list of primes (or any sequence, I suppose), I'd like to calculate and plot the pairs

(x, count of numbers in my list less than or equal to x)

Bonus points for skipping the 'input a list of primes step' and using How to produce a list of prime numbers in LaTeX :-) But, you can assume that input is an ordered, comma-separated list (i.e. a clist) of prime integers. Of course, any input easier than that is also an option.

MWE

\documentclass{article}
\usepackage{pgfplots}
\pgfplotsset{compat=1.12}

\begin{document}
\begin{tikzpicture}
\begin{axis}[
  title={The prime counting function: $\pi(x)$},
  xlabel=$x$,
  ylabel=$\pi(x)$,
  ]
\addplot[only marks]
table {
   x   primes
   0   0 % 0 is not prime; we have (0 primes ≤ 0)
   1   0 % 1 is not prime; we have (0 primes ≤ 1)
   2   1 % 2 is prime; we now have (1 prime  ≤ 2)
   3   2 % 3 is prime; we now have (2 primes ≤ 3)
   4   2 % 4 is not prime; we have (2 primes ≤ 4)
   5   3 % etc.
   6   3
   7   4
   8   4
   9   4
  10   4
};
\end{axis}
\end{tikzpicture}
\end{document}

Instead of providing π(x) manually like this, I'd like to just provide a list of primes in a text file or a clist. As always, brownie points for generality :-)

Note that, per fun, I'm looking for a pure TeX solution – but all solutions (besides the trivial \includegraphics) are welcome :)

Answer 1

You can use pgfplotstable to create a table that contains one column that stores the output of isprime for every number, and one that counts the number of primes encountered so far by summing up the isprime values:

\documentclass{article}
\usepackage{pgfplots, pgfplotstable}
\pgfplotsset{compat=1.12}

\begin{document}
\pgfplotstablenew[
    create on use/x/.style={
        create col/expr={
            \pgfplotstablerow
        }
    },
    create on use/isprime/.style={
        create col/assign/.code={% Can't use "expr" here because "isint" doesn't work reliably with the FPU engine at the moment
            \pgfmathparse{isprime(\thisrow{x})}%
            \pgfkeyslet{/pgfplots/table/create col/next content}\pgfmathresult%
        }
    },
    create on use/primecount/.style={
        create col/expr={
            \pgfmathaccuma + \thisrow{isprime}
        }
    },
    columns={x, isprime, primecount}
]{25}\loadedtable

\begin{tikzpicture}
\begin{axis}[
  title={The prime counting function: $\pi(x)$},
  xlabel=$x$,
  ylabel=$\pi(x)$,
  ]
\addplot[only marks] table [x=x, y=primecount] {\loadedtable};
\end{axis}
\end{tikzpicture}
\end{document}

---

And here's a way of plotting the counting function for an ordered list, using a PGF math array and a counter pointing to the next element in the list:

\documentclass{article}
\usepackage{pgfplots, pgfplotstable}
\pgfplotsset{compat=1.12}

\begin{document}

\def\mylist{{6, 7, 14, 22, 31, 32, 38, 46, 52, 60, 65, 70, 80, 81, 86, 90, 95, 100, 108, 117, 119, 126, 135, 140, 148, 158, 165, 172, 176, 179}}
\newcounter{listindex}

\pgfplotstablenew[
    create on use/x/.style={
        create col/expr={
            \pgfplotstablerow
        }
    },
    create on use/isinlist/.style={
        create col/assign/.code={
            \pgfmathtruncatemacro\thisx{\thisrow{x}}
            \pgfmathtruncatemacro\nextlistitem{\mylist[\value{listindex}]}
            \ifnum\thisx=\nextlistitem
                \def\result{1}
                \stepcounter{listindex}
            \else
                \def\result{0}
            \fi
            \pgfkeyslet{/pgfplots/table/create col/next content}{\result}%
        }
    },
    create on use/count/.style={
        create col/expr={
            \pgfmathaccuma + \thisrow{isinlist}
        }
    },
    columns={x, isinlist, count}
]{50}\loadedtable

\begin{tikzpicture}
\begin{axis}[
  title={The prime counting function: $\pi(x)$},
  xlabel=$x$,
  ylabel=$\pi(x)$,
  ]
\addplot[only marks] table [x=x, y=count] {\loadedtable};
\end{axis}
\end{tikzpicture}
\end{document}

Answer 2

A sagetex solution:

\documentclass{standalone}
\usepackage{sagetex}
\usepackage{pgfplots}
\pgfplotsset{compat=1.12}
\begin{document}
\begin{sagesilent}
output = r""
output += r"\begin{tikzpicture}"
output += r"\begin{axis}["
output += r"title={The prime counting function: $\pi(x)$},"
output += r"xlabel=$x$,"
output += r"ylabel=$\pi(x)$,"
output += r"]"
output += r"\addplot[only marks] coordinates {"
for i in range(0,20):
    output += r"(%s, %s)"%(i,prime_pi(i))
output += r"};"
output += r"\end{axis}"
output += r"\end{tikzpicture}"  
\end{sagesilent}
\sagestr{output}
\end{document}

This gives the following output in Sagemath Cloud:

All you have to do is change 20 to whatever you want and compile the code to get the revised graph. Note that in Python the last number, 20, doesn't execute. This relies on having Sage on your computer or, easier still, using SagemathCloud (internet connection required). Sage handles the calculations; the y-values are determined by the function prime_pi() which is documented here. No files to read, Sage creates the code on-the-fly.

Answer 3

"Pure-TeX" solution can look like:

\def\primes{2,3,5,7,11,13,17,19,23,29,31,37}

\newcount\tmpnum  \newcount\a
\def\ppi#1{\tmpnum=0 \def\ppiA{#1}\expandafter\ppiB\primes,,\relax}
\def\ppiB#1,{\ifx,#1,\message{use more primes}\def\ppiOUT{0}\else
   \ifnum\ppiA<#1\relax \edef\ppiOUT{\the\tmpnum}\ppiC
   \else \advance\tmpnum by1 \fi
   \expandafter\ppiB\fi
}
\def\ppiC#1,\relax{\fi\fi}

\a=0
\loop
    \ppi\a
    \hbox{\hbox to2em{\hss\the\a:}\hbox to2em{\hss\ppiOUT}}
    \ifnum\a<36 \advance\a by1
    \repeat

\bye

And visualisation:

\newdimen\ystep \ystep=2mm
\vbox{\hrule\hbox{\vrule height28mm depth3mm \kern2mm
\loop
    \ppi\a
    \raise\ppiOUT\ystep\hbox{$\bullet$}\kern.1mm
    \ifnum\a<36 \advance\a by1
    \repeat
\kern2mm\vrule}\hrule}

Note that the graphics is "pure-TeX" too. No PostScript, no Tikz, no \pdfspecial.

Answer 4

Here is a direct approach. The usage is explained in the code comments. I have had an issue with pgfplots, I don't know why passing red, or draw=red, or color=red to addplot has an impact on the line thickness of the plot.

update: downgraded to compat=1.12 as I realized 1.13 is very recent. Also, it seems \addplot+[options] is what I should use. But the result is a bit ... artistic. (last image)

The code is for any (non-negative, ordered) integer sequence given as a comma separated list. Nothing here for computing primes. The sample generates for demonstration a random sequence using \pdfuniformdeviate.

\documentclass[tikz,ignorerest=false, border=12pt]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.12}% 1.13 seems to be very recent

% generated a random strictly increasing sequence of 30 integers
% for the purpose of testing

\pdfsetrandomseed 1234

\makeatletter
\newcommand*\seqA {}%
\edef\@tempa {\pdfuniformdeviate10}%
\edef\seqA {\@tempa}%
% for very very long lists, there are faster ways.
% but let's not bother here.
\count@ 29
\loop
    \edef\@tempa {\the\numexpr\@tempa+\@ne+\pdfuniformdeviate10}%
    \edef\seqA {\seqA, \@tempa}%
\advance\count@\m@ne
\ifnum\count@>\z@
\repeat

\typeout{\string\seqA\space prepared with meaning: \meaning\seqA}
% with pdfrandomseed=1234
% \seqA prepared with meaning macro:->6, 7, 14, 22, 31, 32, 38, 46, 52, 60, 65, 70, 80, 81, 86, 90, 95, 100, 108, 117, 119, 126, 135, 140, 148, 158, 165, 172, 176, 179

% \CumulCnts expandably constructs pairs (x, \pi_S(x)), 0≤ x ≤ xmax,
% for sequence S, given as comma separated increasing list

% It admits optional argument, default xmax=100 to limit x.

% usage: \CumulCnts[optional max x]{\A} will expand to the list of pairs
% inside an \edef or a \csname...\endcsname.

% \edef\cumlA {\CumulCnts[optional max x]{\A}}
% and then use \cumulA

\newcommand*\CumulCnts {}
\def\CumulCnts #1{\expandafter\CumulCnts@i\romannumeral`\^^@#1,\relax,}%
\def\CumulCnts@i #1{\ifx [#1\expandafter\CumulCnts@opt\else
                            \expandafter\CumulCnts@noopt\fi #1}%
\def\CumulCnts@opt [#1,\relax,#2]#3%
   {\expandafter\CumulCnts@ii 
    \the\numexpr #2\expandafter;\romannumeral`\^^@#3,\relax,}%
\def\CumulCnts@noopt {\CumulCnts@ii 100;}%
\def\CumulCnts@ii    {\CumulCnts@iii 0;0;}%
\def\CumulCnts@iii #1;#2;#3;#4#5,{%
    \if\relax #4\expandafter\CumulCnts@finish\fi
    \ifnum #3<#4#5
        \expandafter\@firstoftwo
    \else
        \expandafter\@secondoftwo
    \fi
    {\CumulCnts@c #1;#2;#3;}
    {\CumulCnts@iv #1;#2;#4#5;{#3}}%
}%
\def\CumulCnts@finish \ifnum #1\fi #2#3{#2\relax,}

\def\CumulCnts@iv #1;#2;#3;{%
    \ifnum #1=#3 \expandafter\CumulCnts@v\fi
    (#1, #2)
    \expandafter\CumulCnts@iv\the\numexpr #1+\@ne;#2;#3;%
}
\def\CumulCnts@v #1;#2;#3;{\expandafter\CumulCnts@vi\the\numexpr #2+\@ne;#3;}%
\def\CumulCnts@vi #1;#2;#3{%
    (#2, #1)
    \expandafter\CumulCnts@iii \the\numexpr#2+\@ne;#1;#3;}%

\def\CumulCnts@c #1;#2;#3;{%
    \ifnum #1>#3 \expandafter\CumulCnts@d\fi
    (#1, #2)
    \expandafter\CumulCnts@c\the\numexpr #1+\@ne;#2;#3;%
}% 
\def\CumulCnts@d #1;#2;#3;#4\relax,{}%

\makeatletter  



\begin{document}

\edef\cumulA {\CumulCnts{\seqA}}                       
\typeout {\string\cumulA\space prepared with meaning: \meaning\cumulA}

\begin{tikzpicture}
\begin{axis}[
  title={The counting function: $\pi_S(x)$},
  xlabel=$x$,
  ylabel=$\pi_S(x)$,
]
   \addplot coordinates {\cumulA};
\end{axis}
\end{tikzpicture}

% extending to x≤200
\edef\cumulA {\CumulCnts[200]{\seqA}}                       
\typeout {\string\cumulA\space prepared with meaning: \meaning\cumulA}

\begin{tikzpicture}
\begin{axis}[
  title={The counting function: $\pi_S(x)$},
  xlabel=$x$,
  ylabel=$\pi_S(x)$,
]
   \addplot[red] coordinates {\cumulA};
\end{axis}
\end{tikzpicture}
\end{document}

Adding here for the passers-by what gets written to the log (lines cut by TeX to 79 characters)

First our toy random sequence:

\seqA prepared with meaning: macro:->6, 7, 14, 22, 31, 32, 38, 46, 52, 60, 65, 
70, 80, 81, 86, 90, 95, 100, 108, 117, 119, 126, 135, 140, 148, 158, 165, 172, 
176, 179

Then the result of doing: \edef\cumulA {\CumulCnts[200]{\seqA}} (there is a space token after the final coordinate pair; if people protest the code can be modified to not have it...).

\cumulA prepared with meaning: macro:->(0, 0) (1, 0) (2, 0) (3, 0) (4, 0) (5, 0
) (6, 1) (7, 2) (8, 2) (9, 2) (10, 2) (11, 2) (12, 2) (13, 2) (14, 3) (15, 3) (
16, 3) (17, 3) (18, 3) (19, 3) (20, 3) (21, 3) (22, 4) (23, 4) (24, 4) (25, 4) 
(26, 4) (27, 4) (28, 4) (29, 4) (30, 4) (31, 5) (32, 6) (33, 6) (34, 6) (35, 6)
 (36, 6) (37, 6) (38, 7) (39, 7) (40, 7) (41, 7) (42, 7) (43, 7) (44, 7) (45, 7
) (46, 8) (47, 8) (48, 8) (49, 8) (50, 8) (51, 8) (52, 9) (53, 9) (54, 9) (55, 
9) (56, 9) (57, 9) (58, 9) (59, 9) (60, 10) (61, 10) (62, 10) (63, 10) (64, 10)
 (65, 11) (66, 11) (67, 11) (68, 11) (69, 11) (70, 12) (71, 12) (72, 12) (73, 1
2) (74, 12) (75, 12) (76, 12) (77, 12) (78, 12) (79, 12) (80, 13) (81, 14) (82,
 14) (83, 14) (84, 14) (85, 14) (86, 15) (87, 15) (88, 15) (89, 15) (90, 16) (9
1, 16) (92, 16) (93, 16) (94, 16) (95, 17) (96, 17) (97, 17) (98, 17) (99, 17) 
(100, 18) (101, 18) (102, 18) (103, 18) (104, 18) (105, 18) (106, 18) (107, 18)
 (108, 19) (109, 19) (110, 19) (111, 19) (112, 19) (113, 19) (114, 19) (115, 19
) (116, 19) (117, 20) (118, 20) (119, 21) (120, 21) (121, 21) (122, 21) (123, 2
1) (124, 21) (125, 21) (126, 22) (127, 22) (128, 22) (129, 22) (130, 22) (131, 
22) (132, 22) (133, 22) (134, 22) (135, 23) (136, 23) (137, 23) (138, 23) (139,
 23) (140, 24) (141, 24) (142, 24) (143, 24) (144, 24) (145, 24) (146, 24) (147
, 24) (148, 25) (149, 25) (150, 25) (151, 25) (152, 25) (153, 25) (154, 25) (15
5, 25) (156, 25) (157, 25) (158, 26) (159, 26) (160, 26) (161, 26) (162, 26) (1
63, 26) (164, 26) (165, 27) (166, 27) (167, 27) (168, 27) (169, 27) (170, 27) (
171, 27) (172, 28) (173, 28) (174, 28) (175, 28) (176, 29) (177, 29) (178, 29) 
(179, 30) (180, 30) (181, 30) (182, 30) (183, 30) (184, 30) (185, 30) (186, 30)
 (187, 30) (188, 30) (189, 30) (190, 30) (191, 30) (192, 30) (193, 30) (194, 30
) (195, 30) (196, 30) (197, 30) (198, 30) (199, 30) (200, 30) 

---

With :

\begin{tikzpicture}
\begin{axis}[
  title={The counting function: $\pi_S(x)$},
  xlabel=$x$,
  ylabel=$\pi_S(x)$,
]
   \addplot+[color=red] coordinates {\cumulA};
\end{axis}
\end{tikzpicture}

Answer 5

The R/knitrsolution:

\documentclass[12pt,a5paper,twocolumn]{article}
\usepackage{tikz}
\usepackage{amsmath}
\begin{document}
<<plot,echo=F,dev='tikz'>>=
library(numbers)
library(data.table)
x <- data.table(A=c(0:50))
x$B <- as.numeric(isPrime(x$A))  
x[ , C := cumsum(B)]
plot(x$A,x$C,pch=21, 
main="The prime counting function: $\\boldmath{\\pi(x)}$", 
xlab="$x$", ylab="$\\pi(x)$", lwd=4, 
col=rainbow(start=.7, 51), bg=rainbow(200), 
cex=2, cex.lab=2, cex.main=1.5, cex.axis=1.5)
plot(x$A,x$C, 
main="The prime counting function: $\\boldmath{\\pi(x)}$", 
xlab="$x$", ylab="$\\pi(x)$", lwd=4, type="l",
col="red", cex=2, cex.lab=2, cex.main=1.5, cex.axis=1.5)
@
\end{document}

Answer 6

I hope that you will count Metapost as fun...

prologues := 3; outputtemplate := "%j%c.eps";
% see D.E.Knuth, The Metafont Book, p.173
numeric p[]; boolean n_is_prime; p[1]=2; k:=1;  N:=100;
for n=3 step 2 until infinity:
  n_is_prime := true;
  for j=2 upto k:
    if n mod p[j]=0: n_is_prime := false; fi
    exitif n/p[j] < p[j];
  endfor
  if n_is_prime: p[incr k] := n; exitif k=N; fi
endfor fi

beginfig(1);
path P; P = (0,0)
for x=1 upto N:
  hide(pi := 0;
  for i=1 upto infinity:
    exitif p[i]>x;
    pi := pi + 1;
  endfor) -- (x,pi)
endfor;

numeric u, v;
u = 5;
v = 13;

path xx, yy;
xx = origin -- right scaled xpart point infinity of P scaled u;
yy = origin -- up    scaled ypart point infinity of P scaled v;

for x=10 step 10 until N: 
  draw yy shifted (x*u,0) withcolor .8 white; 
  label(decimal x, (x*u,-8));
endfor

for y=5 step 5 until ypart point infinity of P:
  draw xx shifted (0,y*v) withcolor .8 white;
  label(decimal y, (-8,y*v));
endfor

drawarrow xx;
drawarrow yy;

draw P xscaled u yscaled v withcolor .7[blue,white];
for i=0 upto length P: 
   drawdot point i of P xscaled u yscaled v 
           withpen pencircle scaled 3 
           withcolor .67 blue; 
endfor

endfig;

end.

This uses the same algorithm as my answer to the OP linked question. It works up to N=564 because the 565th prime is the first one greater than 4096, which is infinity in Metapost. If you want more, use the -numbersystem=double option, and redefine infinity to be some suitably large value.