By accumulating/integrating the result from a variant of the Fourier transform of the von Mangoldt function I believe it should be possible:

```
Clear[f]
scale = 1000000;
f = Range[scale];
f[[1]] = N@MangoldtLambda[1];
Monitor[Do[
f[[i]] = N@MangoldtLambda[i] + f[[i - 1]], {i, 2, scale}], i]
xres = .002;
xlist = Exp[Range[0, Log[scale], xres]];
tmax = 60;
tres = .015;
s = 1/2;
Monitor[errList1 =
Table[((xlist^(-s + I t).(f[[Floor[xlist]]] - xlist)))*(s +
I t), {t, Range[0, 60, tres]}];, t]
Print["Variant of the Fourier transform of the von Mangoldt function"]
g1 = ListLinePlot[Accumulate[Re[errList1]/Length[xlist]],
DataRange -> {0, 60}, PlotRange -> {-.3, 300}, Axes -> True,
Filling -> Axis];
Show[g1, ImageSize -> Large]
```

I have put a picture of the result here:

http://mobiusfunction.files.wordpress.com/2012/07/zeta-zero-counting-function.png

The scale is wrong though and I don't know how to do it by integration, but the
steps at x-values of imaginary parts of Riemann zeta zeros should all be roughly
the same.

Edit 12.8.2012:
The integration of the formula in the program above is not too hard to do:

```
Clear[xlist, s, t, g];
Integrate[((xlist^(-s + I t)*(g - xlist)))*(-s + I t), t]
Print["g=f[[Floor[xlist]]]"]
Clear[f]
scale = 1000000;
f = Range[scale];
f[[1]] = N@MangoldtLambda[1];
Monitor[Do[
f[[i]] = N@MangoldtLambda[i] + f[[i - 1]], {i, 2, scale}], i]
xres = .002;
xlist = Exp[Range[0, Log[scale], xres]];
xlist[[1]] = 1.001;
tmax = 60;
tres = .015;
s = 1/2;
Monitor[errList1 =
Table[Total[((f[[Floor[xlist]]] - xlist) xlist^(-s +
I t) (I + (I s + t) Log[xlist]))/Log[xlist]^2], {t,
Range[0, 60, tres]}];, t]
Print["Variant of the Fourier transform of the von Mangoldt function"]
g1 = ListLinePlot[Re[errList1]/Length[xlist] - 1,
DataRange -> {0, 60}, PlotRange -> {-0, 3.7}, Axes -> True];
Show[g1, ImageSize -> Full]
```

I have put picture of the result here:

http://mobiusfunction.files.wordpress.com/2012/08/zeta-zero-counting-august-2012.png

I don't know if it is entirely correct and the scaling is wrong.
The minus one seems to make the counting start at the first zeta zero.

Edit 5.2.2014:

Using this as a starting point:

$$\Lambda(n)=\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}.$$

Where $\Lambda(n)$ is the von Mangoldt function.

Edit 4.2.2014:

The Dirichlet series:

$$f(t) = (1-\sum\limits_{n=1}^{n=k} \frac{1}{1.2\log(k)} \frac{1}{n} \zeta(1/2+i \cdot t)\sum\limits_{d|n} \frac{\mu(d)}{d^{(1/2+i \cdot t-1)}})^{12}$$

leads to a sequence of numbers that when accumulated gives the plot:

Mathematica code for the plot above:

```
scale = 400;
Print["Counting to 60"]
Monitor[g1 =
ListLinePlot[
Accumulate[(1 -
Table[Re[
Zeta[1/2 - I*k]*
Total[Table[
Total[MoebiusMu[Divisors[n]]/Divisors[n]^(1/2 - I*k - 1)]/
n/Log[scale]/1.2, {n, 1, scale}]]], {k, 0 + 1/1000, 60,
N[1/6]}])^12], DataRange -> {0, 60},
PlotRange -> {-0.15, 15}], Floor[k]]
```

Edit 13.2.2014:

```
(*program start*)
scale = 300;
Print["Counting to 60"]
Monitor[g1 =
ListLinePlot[
0.69*Accumulate[
Table[Exp[-Re[
Zeta[1/2 - I*k]*
Total[Table[
Total[MoebiusMu[Divisors[n]]/
Divisors[n]^(1/2 - I*k - 1)]/n, {n, 1, scale}]]]], {k,
0 + 1/1000, 60, N[1/6]}]], DataRange -> {0, 60},
PlotRange -> {-0.15, 15}];, Floor[k]]
Show[g1, ListPlot[Table[{N[Im[ZetaZero[n]]], n}, {n, 1, 13}],
PlotStyle -> Black, Filling -> Axis]]
(*program end*)
```

The plot below from the program above is something like:

$$f(k)=\exp \left(-\Re\left(\zeta \left(\frac{1}{2}-i k\right) \sum _{n=1}^{\text{scale}} \frac{\sum\limits_{d|n}\frac{\mu (d(n))}{d(n)^{-i k+\frac{1}{2}-1}}}{n}\right)\right)$$