Abstract: The discrete ordinates (S_N) method proposed by Carlson constitutes an angular discretization for neutron transport equation. Owing to its superior computational accuracy and efficiency compared to alternative numerical approaches, the S_N method has been extensively adopted in neutron transport applications, particularly in nuclear reactor design and radiation shielding optimization. Although the S_N method demonstrates satisfactory performance and yields reliable results in most practical scenarios, it inherently suffers from ray effects, a well-documented numerical artifact that fundamentally limits its solution accuracy in the problems involving localized neutron source and wea scattering media. These angular discretization artifacts induce non-physical distortions in computed flux distributions, specifically generating artificial flux overestimations along preferential angular directions while suppressing flux magnitudes in geometrically unaligned regions. To address this defect, numerous of ray effect mitigation methods were developed. The virtual source method is formulated through the strategic introduction of pseudo-source terms into the S_N equations, establishing equivalence with spherical harmonics (P_N) formulations to eliminate ray effects. While this approach achieves angular flux correction by enforcing moment-matching conditions between S_N and P_N operators, its computational formalism becomes increasingly intricate due to higher-order spherical harmonic expansions, and suffers from convergence reliability concerns when handling multigroup problems with strongly anisotropic scattering coupling. The first collision source method, widely implemented in S_N codes to mitigate ray effects, operates by employing alternative transport methodologies, such as ray tracing, Monte Carlo simulations, or P_N expansions, to compute the uncollided flux component most severely impacted by ray effects, thereby reducing numerical distortions. However, this approach necessitates hybrid code architectures that integrate supplementary transport solvers with S_N frameworks, substantially increasing computational complexity and programming challenges, while inheriting inherent limitations from the auxiliary methods: ray tracing struggles with reflective boundary conditions, Monte Carlo suffers from expensive computational costs, and P_N method exhibit numerical oscillations in angularly anisotropic problems. In summary, ray effect severely constrains the computational accuracy of the S_N method and have not yet been effectively resolved. In this paper, it is noticed that the S_N method fundamentally constitutes a deterministic numerical integration scheme in angular dimension, whose characteristic ray effects arise from the inherent precision limitations of low-order angular quadrature techniques, particularly evident in non-smooth angular flux distributions with pronounced directional dependencies. Motivated by the capacity of Monte Carlo method to achieve precise integration of non-smooth functions, a novel random S_N method based on quasi-Monte Carlo quadrature sets was proposed. Preliminary numerical results demonstrate that with equivalent numbers of discrete angular directions, the proposed method exhibits superior accuracy to conventional quadrature sets in strongly angularly anisotropic problems, but underperforms in weakly anisotropic scenarios. A collision-coupled rS_N-S_N methodology that synergistically integrates the advantages of stochastic and deterministic quadrature sets was further developed, with numerical verification confirming its capability to achieve enhanced computational precision at reduced computational time and resource expenditure compared to traditional implementations.