Piecewise algebraic B-spline surfaces (ABS surfaces) are capable of modeling globally smooth shapes of arbitrary topology. These can be potentially applied in geometric modeling, scientific visualization, computer animation and mathematical illustration. However, real-time ray casting the surface is still an obstacle for interactive applications, due to the large amount of numerical root findings of nonlinear polynomial systems that are required. In this paper, we present a GPU-based real-time ray casting method for ABS surfaces. To explore the powerful parallel computing capacity of contemporary GPUs, we adopt iterative numerical root-finding algorithms, e.g., the Newton–Raphson and regula falsi algorithms, rather than recursive ones. To facilitate convergence of the Newton–Raphson or regula falsi algorithm, their initial guesses are determined through rasterization of the isotopic isosurface, and the isosurface is generated based on regular criteria for surface domain subdivision. Meanwhile, polar surfaces are adopted to identify single roots or to isolate different roots, i.e., ray and surface intersections. As an important geometric feature, the silhouette curve is elaborately computed to floating-point accuracy, which can be applied in further anti-aliasing processes. The experimental results show that the proposed method can render thousands of piecewise algebraic surface patches of degrees 6–9 in real time.