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

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)$,
]
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 , I'm looking for a pure TeX solution – but all solutions (besides the trivial \includegraphics) are welcome :)

• /pgfplots/hist/cumulative? – Symbol 1 Jan 10 '16 at 3:07
• @Symbol1 Ooh, that looks promising. I can't quite it to work though: \begin{axis}[hist=cumulative] \addplot[only marks] table [y index=0] { 2 3 5 7 11 … 1811 }; \end{axis} does not work. – Sean Allred Jan 10 '16 at 3:16
• Would the downvoter care to leave a helpful comment on how this post can be improved? – Sean Allred Jan 11 '16 at 3:24
• Use a specialized program for ploting, like gnuplot or perhaps the plotting facilities of asymptote if LaTeX look-and-feel is important? – vonbrand Jan 13 '16 at 1:24
• @vonbrand Note that this is for fun :-) I've never used gnuplot or asymptote, but it would be neat to compare it to the other answers. – Sean Allred Jan 13 '16 at 1:30

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}

\begin{tikzpicture}
\begin{axis}[
title={The prime counting function: $\pi(x)$},
xlabel=$x$,
ylabel=$\pi(x)$,
]
\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}

\begin{tikzpicture}
\begin{axis}[
title={The prime counting function: $\pi(x)$},
xlabel=$x$,
ylabel=$\pi(x)$,
]
\end{axis}
\end{tikzpicture}
\end{document}


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.

• Out of curiosity, what version of Python does Sage rely on? If it's 2.x, I imagine xrange should be used instead. – Sean Allred Jan 13 '16 at 1:32
• The documentation says 2.6x. I just tried using xrange and there was no difference in number of points plotted. – DJP Jan 13 '16 at 2:02
• There wouldn't be – xrange is just more memory-efficient. It uses a generator instead of creating the whole list on the stack. – Sean Allred Jan 13 '16 at 2:03

"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
\expandafter\ppiB\fi
}
\def\ppiC#1,\relax{\fi\fi}

\a=0
\loop
\ppi\a
\hbox{\hbox to2em{\hss\the\a:}\hbox to2em{\hss\ppiOUT}}
\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
\repeat
\kern2mm\vrule}\hrule}


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

• When I saw your answer I thought you would be printing it with raw PostScript :) – Manuel Jan 11 '16 at 0:42
• @Manuel Don't be silly – that wouldn't be device-independent :-) Awesome stuff, wipet! – Sean Allred Jan 23 '16 at 19:51

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}%
\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)$,
]
\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)$,
]
\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)$,
]
\end{axis}
\end{tikzpicture}


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}


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.

• Metapost is fun :-) +1 – Sean Allred Jan 11 '16 at 16:09