That partial difference equation is

dS(x,y) = (P(x/y,y-1) - P(y-1, sqrt(y-1)))(P(y, sqrt(y)) - P(y-1, sqrt(y-1)))

where

P(x,y) = [x] - S(x,y) - 1

with the constraint that if y>sqrt(x), then P(x,y) = P(x,sqrt(x)), which is necessary to handle the fact that [x/p] - 1 is the count of composites that have a prime p as a factor for a positive integer x, as long as p<=sqrt(x).

And finally, S(x,y) is the sum of dS from 2 to y or sqrt(x) if y>sqrt(x), as it has the same constraint, and S(x,1) = 0, by definition.

It's easy enough to move from the partial difference equation dS(x,y) to a partial differential equation:

ΔS(x,y) = (P(x/y,y-Δy) - P(y-Δy, sqrt(y-Δy)))(P(y, sqrt(y)) - P(y-Δy, sqrt(y-Δy)))

and divide by Δy to get

ΔS(x,y)/Δy = (P(x/y,y-Δy) - P(y-Δy, sqrt(y-Δy)))(P(y, sqrt(y)) - P(y-Δy, sqrt(y-Δy)))/Δy

and let Δy approach 0 to get

S'

_{y}(x,y) = (P(x/y,y) - P(y, sqrt(y))) P'(y, sqrt(y))

and now I also have

P(x,y) = x - S(x,y) - 1

and differentiating with respect to y, gives

P'

_{y}(x,y) = -S'

_{y}

so I have

P'

_{y}(x,y) = -(P(x/y,y) - P(y, sqrt(y))) P'(y, sqrt(y))

and a continuous function.

Now you can numerically integrate that and I'm not really up on numerical integration so I made a really basic program and did so, and found that I got values closer to pi(x) than li(x) but a bit further than R(x), but I did just a rough go at it.

To date, I've yet to find anyone to help me with the numerical integration!

If that partial differential integrates to a close value to pi(x) then that could be rather important news.

James