For culling faraway particles: as lightxbulb said, the keyword to search for is "acceleration structure". There are a variety of options: BVH, octree, kd-tree, uniform grid, or some combination of these. Which one will be most suitable is going to depend on the details of your situation—how large/dense is the simulation, how much you care about performance vs memory vs implementation complexity, and so on.
They all work by subdividing space into volumes in some fashion and keeping track of which particles are in each volume. Then, when you're marching you only need to keep track of which volume you're in, and that will tell you which particles can potentially contribute to the density at that point. For instance if it's a uniform grid, you only need to look at the current cell and possibly some of its immediate neighbors if they're within the maximum particle radius of your current point.
Getting a smoothed mesh with surface normals etc. is usually done by defining a density field for the particle system. This is a function on 3D space that consists of summing up a falloff function for each particle over a defined radius:
$$
D(\vec r) = \sum_{i \in \text{particles}} f \left(\frac{|\vec r - \vec r_i|}{R_i} \right)
$$
where $f(u)$ is some function that falls off smoothly from $f(0) = 1$ to $f(1) = 0$ (for example, a simple cubic falloff: $f(u) = 1 - 3u^2 + 2u^3$), and $\vec r_i, R_i$ are the position and radius of the $i$th particle. Then you set some threshold density value and let the surface be defined by all the points where $D(\vec r) = \text{threshold}$.
For getting the surface normal: the surface normal is just the (normalized) gradient vector of the density field. This can be calculated by finite differences between nearby points in the field:
$$
\begin{aligned}
\nabla D(\vec r) &= \left( \frac{\partial D}{\partial x}, \frac{\partial D}{\partial y}, \frac{\partial D}{\partial z} \right) \\
&\approx \frac{1}{\epsilon} \left( D(\vec r + \epsilon \vec x) - D(\vec r), D(\vec r + \epsilon \vec y) - D(\vec r), D(\vec r + \epsilon \vec z) - D(\vec r) \right)
\end{aligned}
$$
or analytically by adding up gradient contributions from each particle.
$$
\begin{aligned}
\nabla D(\vec r) &= \sum_i \nabla \left( f \left(\frac{|\vec r - \vec r_i|}{R_i} \right) \right) \\
&= \sum_i \text{normalize}(\vec r - \vec r_i) \frac{1}{R_i} \frac{df}{du}
\end{aligned}
$$
For getting the surface in the form of a mesh: the keyword to search for is "isosurface extraction". Marching cubes is one example of such an algorithm, as lightxbulb mentioned. It evaluates the density field at evenly spaced grid points, uses those values to estimate where the intersection of the surface with the cubes' edges will be, and constructs from that a tiny mesh within each cube. Stitching all those meshes together yields a representation of the full surface. Naturally, the accuracy of it depends on the size of the grid used.
There is also a variant that uses tetahedra instead of cubes, and another similar algorithm called "dual contouring". Here's a nice article showing how these work interactively (in 2D). There's also a great deal of literature on how to improve accuracy, performance, etc. with these kinds of methods.
