Spaces:
Runtime error
Runtime error
| template <typename T> | |
| DEVICE | |
| inline bool solve_quadratic(T a, T b, T c, T *t0, T *t1) { | |
| // From https://github.com/mmp/pbrt-v3/blob/master/src/core/pbrt.h#L419 | |
| T discrim = square(b) - 4 * a * c; | |
| if (discrim < 0) { | |
| return false; | |
| } | |
| T root_discrim = sqrt(discrim); | |
| T q; | |
| if (b < 0) { | |
| q = -0.5f * (b - root_discrim); | |
| } else { | |
| q = -0.5f * (b + root_discrim); | |
| } | |
| *t0 = q / a; | |
| *t1 = c / q; | |
| if (*t0 > *t1) { | |
| swap_(*t0, *t1); | |
| } | |
| return true; | |
| } | |
| template <typename T> | |
| DEVICE | |
| inline int solve_cubic(T a, T b, T c, T d, T t[3]) { | |
| if (fabs(a) < 1e-6f) { | |
| if (solve_quadratic(b, c, d, &t[0], &t[1])) { | |
| return 2; | |
| } else { | |
| return 0; | |
| } | |
| } | |
| // normalize cubic equation | |
| b /= a; | |
| c /= a; | |
| d /= a; | |
| T Q = (b * b - 3 * c) / 9.f; | |
| T R = (2 * b * b * b - 9 * b * c + 27 * d) / 54.f; | |
| if (R * R < Q * Q * Q) { | |
| // 3 real roots | |
| T theta = acos(R / sqrt(Q * Q * Q)); | |
| t[0] = -2.f * sqrt(Q) * cos(theta / 3.f) - b / 3.f; | |
| t[1] = -2.f * sqrt(Q) * cos((theta + 2.f * T(M_PI)) / 3.f) - b / 3.f; | |
| t[2] = -2.f * sqrt(Q) * cos((theta - 2.f * T(M_PI)) / 3.f) - b / 3.f; | |
| return 3; | |
| } else { | |
| T A = R > 0 ? -pow(R + sqrt(R * R - Q * Q * Q), T(1./3.)): | |
| pow(-R + sqrt(R * R - Q * Q * Q), T(1./3.)); | |
| T B = fabs(A) > 1e-6f ? Q / A : T(0); | |
| t[0] = (A + B) - b / T(3); | |
| return 1; | |
| } | |
| } | |