From your question, you want to find every cube that a ray passes through, and enumerate them.
The direction he's looking in will be signified by $(d_x,d_y,d_z)\newcommand{\sgn}{\operatorname{sgn}}$, note that I assume that the direction is a vector starting at the center of the coordinate system, and not from the player.
We can simply enumerate coordinates on the ray with some interval, and discard duplicates, we just have to make sure that the interval is small enough such that we miss no blocks.
If we want this interval to be small enough such that it always ends up in the current or an adjacent block, we want $e$ to be the minumum distance to one of the lines in the grid between the boxes. (If the point is on a line, we discard that one) So let $e$ be that distance then we simply need to normalize the vector to have that size.
Finding $e$ is tricky, we need to find the distance from the point to many lines in a grid. For the following I assume that $q$ is the side length between the boxes, and that $x,y,z$ is the coordinate of the point.
We will perform an algorithm 3 times, one for each dimension. Each time we perform it, we are given a number, totaling at $3$ numbers. $e$ is the smallest of those $3$ numbers. The first time we will have $h=x,p=y$, the second time $h=x,p=z$ and the third time will have $h=y,p=z$.
The algorithm is as follows:
Let $\lfloor x\rfloor_q$ be $x$ rounded down with resolution $q$, and let $[x]_q$ be defined as the smallest of the following two numbers: $(x-\lfloor x\rfloor_q)^2$ and $(q-x+\lfloor x\rfloor_q)^2$
If $[h]_q=[p]_q=0$ then the number returned by the algorithm is $q$.
Otherwise the number returned by the algorithm is $\sqrt{[h]_q+[p]_q}$
Let $s=\dfrac1e\sqrt{(d_x)^2+(d_y)^2+(d_z)^2}$, and define:
$$
\begin{align}
\delta_x&=\frac{d_x}s\\
\delta_y&=\frac{d_y}s\\
\delta_z&=\frac{d_z}s
\end{align}
$$
Now simply given a point, calculate the next point with the given formula:
$$
\begin{align}
\text{next $x$}&=x+\delta_x\\
\text{next $y$}&=y+\delta_y\\
\text{next $z$}&=z+\delta_z
\end{align}
$$
Where $x,y,z$ are the coordinates of the previous point. Each of these points will be inside a box on the ray, and you wont miss any because we normalized the direction. Note that we might get the same box several times in a row, so if that's a problem you have to check for duplicates. The starting point is of course the players position.
Best Answer
By $\,\, \times \,\, $ I denote cross-product and by $(\,\, \cdot \,\, )$ I denote dot product. The length of a vector is $|\vec{u}| = \sqrt{(\vec{u} \cdot \vec{u})}$.
Look at the picture. $S$ is the source of the light ray and the direction of the ray is denoted by $\vec{v}$. The three points are $A, \, B, \, C$. Say $O$ is the origin of your coordinate system. Now, the unit normal vector to the plane determined by points $A, \, B, \, C$ can be calculated for example as $$\vec{n} = \frac{\vec{AB} \times \vec{AC}}{\big|\vec{AB} \times \vec{AC}\big|}$$ Vector $\vec{OK}$ in the right side of the picture is the orthogonal projection of vector $\vec{v} = \vec{OM}$ onto $\vec{n}$ and is thus defined as $$\vec{OK} = (\vec{v} \cdot \vec{n})\, \vec{n}$$ so twice the vector $\vec{OK}$ is $$\vec{OL} = 2 \,\vec{OK} = 2\, (\vec{v} \cdot \vec{n}) \, \vec{n}$$ Thus, the reflected vector is $$\vec{w} = \vec{LM} = \vec{OM} - \vec{OL} = \vec{v} - 2\, (\vec{v} \cdot \vec{n}) \, \vec{n}$$
Now the intersection point $R$ between the light ray and the plane $A, B, C$. The vector equation for the plane $A, B, C$ can be written as $$\vec{n} \cdot \big(\vec{OX} - \vec{OA}\big) = 0$$ where $X$ is any point on the plane and $\vec{OX}$ is its position (coordinate) vector, i.e. the vector pointing from the origin $O$ to the point $X$. The equation for the light ray is $$\vec{OX} = \vec{OS} + t \, \vec{v}$$ for any point $X$ on the light ray ($t\geq 0$ in order to have a ray). Thus the intersection point $R$ should satisfy both equations, i.e. \begin{align} &\vec{n} \cdot \big(\vec{OR} - \vec{OA}\big) = 0\\ &\vec{OR} = \vec{OS} + t \, \vec{v} \end{align} so substitute the second equation in the first and \begin{align} 0 &= \vec{n} \cdot \big(\vec{OS} + t \, \vec{v} - \vec{OA}\big) = \vec{n} \cdot \big(\vec{OS} - \vec{OA} + t \, \vec{v} \big) = \vec{n} \cdot \big(\vec{AS} + t \, \vec{v} \big) = \big(\vec{n} \cdot \vec{AS}\big) + t \big(\vec{n} \cdot \vec{v} \big) \end{align} and when you solve for $t$ you get $$t = \frac{- \big(\vec{n} \cdot \vec{AS}\big)}{\big(\vec{n} \cdot \vec{v} \big)} = \frac{\big(\vec{n} \cdot \vec{SA}\big)}{\big(\vec{n} \cdot \vec{v} \big)}$$ notice $\vec{SA} = - \vec{AS}$. Therefore, $$\vec{OR} = \vec{OS} + \frac{\big(\vec{n} \cdot \vec{SA}\big)}{\big(\vec{n} \cdot \vec{v} \big)} \, \vec{v}$$
Thus, the reflected ray is given by the following point $R$ and vector $\vec{w}$ \begin{align} \vec{OR} &= \vec{OS} + \frac{\big(\vec{n} \cdot \vec{SA}\big)}{\big(\vec{n} \cdot \vec{v} \big)} \, \vec{v}\\ \vec{w} &= \vec{v} - 2\, (\vec{v} \cdot \vec{n}) \, \vec{n}\\ \vec{n} &= \frac{\vec{AB} \times \vec{AC}}{\big|\vec{AB} \times \vec{AC}\big|} \end{align} You can permute the points $A, B, C$ any way you want, you will get the same result.