corecel/math/FerrariSolver.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file corecel/math/FerrariSolver.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Constants.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/cont/Array.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "corecel/math/NumericLimits.hh"
0016 #include "corecel/math/PolyEvaluator.hh"
0017 #include "corecel/math/SoftEqual.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Find positive, real roots for quartic functions using Ferarri.
0024  *
0025  * The quartic solver uses the Ferrari method \citep{polyanin-ferrari-2007,
0026  * https://doi.org/10.1201/9781420010510}, using the underlying quadratic and
0027  * cubic solutions from
0028  * \citet{press-nr-1992, https://dl.acm.org/doi/10.5555/1403886} .
0029  *
0030  * The quartic equation
0031  * \f[
0032    a x^4 + b x^3 + c x^2 + d x + e = 0
0033  * \f]
0034  * has four solutions mathematically, but we only require solutions which are
0035  * both real and positive. This equation is also subject to multiple cases of
0036  * catastrophic precision-limitation-based error both fundamentally and as a
0037  * consequence of the particular algorithm chosen. This solver implements the
0038  * Ferrari method, which is well-established and simple, but more
0039  * prone to numerical error than contemporary methods to be explored such as
0040  * Algorithm 1010 \citep{orellana-alg1010-2020,
0041  * https://doi.org/10.1145/3386241}.
0042  *
0043  * The input argument to an instance of this class is an array \c abcde that
0044  * corresponds to {a, b, c, d, e}. An overload using a four-element array
0045  * \c abcd solves the degenerate case where \f$ e = 0 \f$.
0046  *
0047  * The result is an array of 4 real numbers, where each is either a positive
0048  * valid intersection or the sentinel result \c infinity.
0049  */
0050 class FerrariSolver
0051 {
0052   public:
0053     //!@{
0054     //! \name Type aliases
0055     using Real4 = Array<real_type, 4>;
0056     using Real5 = Array<real_type, 5>;
0057     using result_type = Real4;
0058     //!@}
0059
0060   public:
0061     // Construct with given tolerance
0062     inline CELER_FUNCTION FerrariSolver(real_type tolerance);
0063
0064     //! Construct with default tolerance equal to Orange `Tolerance`.
0065     inline CELER_FUNCTION FerrariSolver() : FerrariSolver{default_tol_} {}
0066
0067     // Solver for fully general case
0068     inline CELER_FUNCTION result_type operator()(Real5 const& abcde) const;
0069
0070     // Solver for surface case
0071     inline CELER_FUNCTION result_type operator()(Real4 const& abcd) const;
0072
0073   private:
0074     //// TYPES ////
0075
0076     using Real2 = Array<real_type, 2>;
0077     using Real3 = Array<real_type, 3>;
0078
0079     //// STATIC DATA ////
0080
0081     //! Default tolerance for quadric solve, taken from Orange `Tolerance`.
0082     static constexpr real_type default_tol_
0083         = (std::is_same_v<real_type, double> ? 1e-5 : 5e-2f);
0084
0085     //! No positive real solution (aka "no intersection")
0086     static constexpr real_type no_solution_
0087         = NumericLimits<real_type>::infinity();
0088
0089     //// DATA ////
0090
0091     // Soft zero for biquadratic and degenerate cubic detection
0092     SoftZero<real_type> const soft_zero_;
0093
0094     //// HELPER FUNCTIONS ////
0095     // Try to place real at given index in list, return next free index
0096     inline CELER_FUNCTION int place_root(
0097         result_type& roots, real_type new_root, int free_index) const;
0098
0099     // Find roots of special reduced quartic which is biquadratic
0100     inline CELER_FUNCTION result_type calc_biquadratic_roots(
0101         real_type qb, real_type p, real_type r) const;
0102
0103     // Find all roots of normalized cubic (unsorted, dominant first)
0104     inline CELER_FUNCTION Real3 real_roots_normalized_cubic(
0105         real_type b, real_type c, real_type d) const;
0106
0107     // Find real quadratic roots
0108     inline CELER_FUNCTION Real2 real_roots_normalized_quadratic(
0109         real_type b, real_type c) const;
0110 };
0111
0112 //---------------------------------------------------------------------------//
0113 // INLINE DEFINITIONS
0114 //---------------------------------------------------------------------------//
0115 /*!
0116  * Construct a solver instance with a specified tolerance for degenerate cases,
0117  * such as the particle starting on the surface.
0118  */
0119 CELER_FUNCTION FerrariSolver::FerrariSolver(real_type tolerance)
0120     : soft_zero_{tolerance}
0121 {
0122 }
0123
0124 //---------------------------------------------------------------------------//
0125 /*!
0126  * Find all positive roots of the polynomial with given a, b, c, d, e.
0127  *
0128  * Replaces negative or complex roots with sentry value no_solution_.
0129  *
0130  * Allows user to specify if the operation takes place on the surface, in which
0131  * case it will solve as a cubic to avoid the root at 0.
0132  */
0133 CELER_FUNCTION auto FerrariSolver::operator()(Real5 const& abcde) const
0134     -> result_type
0135 {
0136     using namespace celeritas::literals;
0137
0138     CELER_EXPECT(abcde[0] != 0);
0139     auto [a, b, c, d, e] = abcde;
0140
0141     // Normalize coefficients
0142     real_type ba = b / a, ca = c / a, da = d / a, ea = e / a;
0143
0144     constexpr real_type half{0.5};
0145     real_type qb = 0.25_r * ba;
0146
0147     // Incomplete quartic
0148     real_type p = PolyEvaluator{-half * ca, 0, 3}(qb);
0149     real_type q = PolyEvaluator{half * da, -ca, 0, 4}(qb);
0150     real_type r = PolyEvaluator{-ea, da, -ca, 0, 3}(qb);
0151
0152     // Edge case: equation is biquadratic
0153     if (soft_zero_(q))
0154     {
0155         return calc_biquadratic_roots(qb, p, r);
0156     }
0157
0158     // One real root of subsidiary cubic
0159     Real3 z = FerrariSolver::real_roots_normalized_cubic(
0160         p, r, p * r - half * ipow<2>(q));
0161     real_type z0 = (z[1] == no_solution_) ? z[0] : max(z[0], max(z[1], z[2]));
0162
0163     real_type s2 = 2 * p + 2 * z0;
0164     if (s2 >= 0)
0165     {
0166         real_type s = std::sqrt(s2);
0167         real_type t;
0168         if (soft_zero_(s))
0169         {
0170             t = ipow<2>(z0) + r;
0171         }
0172         else
0173         {
0174             t = -q / s;
0175         }
0176         auto const [r0, r1] = real_roots_normalized_quadratic(s * half, z0 + t);
0177         auto const [r2, r3]
0178             = real_roots_normalized_quadratic(-s * half, z0 - t);
0179
0180         result_type roots(
0181             no_solution_, no_solution_, no_solution_, no_solution_);
0182         int idx = 0;
0183         idx = place_root(roots, r0 - qb, idx);
0184         idx = place_root(roots, r1 - qb, idx);
0185         idx = place_root(roots, r2 - qb, idx);
0186         place_root(roots, r3 - qb, idx);
0187
0188         return roots;
0189     }
0190     else
0191     {
0192         return result_type(
0193             no_solution_, no_solution_, no_solution_, no_solution_);
0194     }
0195 }
0196
0197 //---------------------------------------------------------------------------//
0198 /*!
0199  * Solve a quartic polynomial where coefficient \c e is known to be 0.
0200  *
0201  * Solves as cubic equation, and does not return the known root of 0.
0202  */
0203 CELER_FUNCTION auto FerrariSolver::operator()(Real4 const& abcd) const
0204     -> result_type
0205 {
0206     auto [a, b, c, d] = abcd;
0207     // Normalize coefficients
0208     real_type ba = b / a, ca = c / a, da = d / a;
0209     Real3 cubic_roots = real_roots_normalized_cubic(ba, ca, da);
0210     auto [z0, z1, z2] = cubic_roots;
0211
0212     result_type roots(no_solution_, no_solution_, no_solution_, no_solution_);
0213
0214     int idx = 0;
0215     idx = place_root(roots, z0, idx);
0216     idx = place_root(roots, z1, idx);
0217     place_root(roots, z2, idx);
0218     return roots;
0219 }
0220
0221 //---------------------------------------------------------------------------//
0222 /*!
0223  * Attempt to put a value into the given list at given index, returning where
0224  * to place the next item.
0225  *
0226  * If the given value is no_solution_ or is not positive, does not place
0227  * the root, and returns the same index for the next one.
0228  */
0229 CELER_FUNCTION int FerrariSolver::place_root(
0230     result_type& roots, real_type new_root, int free_index) const
0231 {
0232     if (!(new_root == no_solution_ || new_root <= 0))
0233     {
0234         roots[free_index] = new_root;
0235         free_index += 1;
0236     }
0237     return free_index;
0238 }
0239
0240 //---------------------------------------------------------------------------//
0241 /*!
0242  * Solve special case of Ferrari where reduced quartic is also biquadratic.
0243  *
0244  * In this special case, the normal solution won't work, and must instead be
0245  * solved as a quadratic equation: The square roots of each quadratic solution
0246  * then go on to form potential quartic solutions, for up to four roots.
0247  */
0248 CELER_FUNCTION auto FerrariSolver::calc_biquadratic_roots(
0249     real_type qb, real_type p, real_type r) const -> result_type
0250 {
0251     auto ir = real_roots_normalized_quadratic(-p, -r);
0252     result_type roots(no_solution_, no_solution_, no_solution_, no_solution_);
0253     int idx = 0;
0254     if (ir[1] != no_solution_ && ir[1] > 0)
0255     {
0256         real_type sqrt_ir1 = std::sqrt(ir[1]);
0257         real_type from_pos1 = sqrt_ir1 - qb;
0258         idx = place_root(roots, from_pos1, idx);
0259         if (from_pos1 > 0)
0260         {
0261             idx = place_root(roots, -sqrt_ir1 - qb, idx);
0262         }
0263     }
0264     if (ir[0] != no_solution_ && ir[0] > 0)
0265     {
0266         real_type sqrt_ir0 = std::sqrt(ir[0]);
0267         real_type from_pos0 = sqrt_ir0 - qb;
0268         idx = place_root(roots, from_pos0, idx);
0269         if (from_pos0 > 0)
0270         {
0271             place_root(roots, -sqrt_ir0 - qb, idx);
0272         }
0273     }
0274     return roots;
0275 }
0276
0277 //---------------------------------------------------------------------------//
0278 /*!
0279  * Solve for the real roots of a cubic function.
0280  *
0281  * Specifically, the cubic function
0282  * \f[
0283    a x^3 + b x^2 + c x + d
0284  * \f]
0285  * where a is assumed to already be 1, and is not provided to the
0286  * function.
0287  *
0288  * \return The real roots of the given cubic equation, with the dominant at
0289  * index 0.
0290  */
0291 CELER_FUNCTION auto FerrariSolver::real_roots_normalized_cubic(
0292     real_type b, real_type c, real_type d) const -> Real3
0293 {
0294     using namespace celeritas::literals;
0295
0296     constexpr real_type half = 0.5_r;
0297     constexpr real_type third = 1_r / 3;
0298     real_type third_b = b * third;
0299
0300     // Intermediate values
0301     real_type q = ipow<2>(third_b) - third * c;
0302     real_type r = half * PolyEvaluator{d, -c, 0, 2}(third_b);
0303
0304     real_type q3 = ipow<3>(q);
0305     real_type r2 = ipow<2>(r);
0306
0307     real_type discrim = r2 - q3;
0308
0309     if (soft_zero_(q) && soft_zero_(r) && soft_zero_(discrim))
0310     {
0311         return Real3(-std::cbrt(d), no_solution_, no_solution_);
0312     }
0313     else if (discrim <= 0)  // All roots real, calculate with trigomonetry
0314     {
0315         real_type theta = std::acos(r / std::sqrt(q3));
0316         real_type n2_root_q = -2_r * std::sqrt(q);
0317         real_type twth_pi = constants::pi * 2_r * third;
0318         real_type third_theta = theta * third;
0319
0320         real_type z0 = n2_root_q * std::cos(third_theta) - third_b;
0321         real_type z1 = n2_root_q * std::cos(third_theta + twth_pi) - third_b;
0322         real_type z2 = n2_root_q * std::cos(third_theta - twth_pi) - third_b;
0323
0324         return Real3(z0, z1, z2);
0325     }
0326     else  // One real and two complex roots, solve for real root with Cardano
0327     {
0328         real_type nr_a = -signum(r)
0329                          * std::cbrt(std::fabs(r) + std::sqrt(discrim));
0330         real_type nr_b = nr_a == 0 ? 0 : q / nr_a;
0331         real_type z0 = nr_a + nr_b - third_b;
0332         return Real3(z0, no_solution_, no_solution_);
0333     }
0334 }
0335
0336 //---------------------------------------------------------------------------//
0337 /*!
0338  * Solve for the real roots of a quadratic function.
0339  *
0340  * Specifically, the quadratic function
0341  * \f[
0342    a x^2 + (hb*2) x + c
0343  * \f]
0344  * where a is assumed to already be 1 and not provided.
0345  *
0346  * \return A pair of roots. If roots are imaginary, returns 2x
0347  * no_solution_.
0348  */
0349 CELER_FUNCTION auto FerrariSolver::real_roots_normalized_quadratic(
0350     real_type hb, real_type c) const -> Real2
0351 {
0352     real_type qb2 = ipow<2>(hb);
0353     if (soft_zero_(qb2 - c))
0354     {
0355         // One critical root
0356         return Real2(-hb, no_solution_);
0357     }
0358     else if (qb2 > c)
0359     {
0360         // Two real roots
0361         real_type ht = std::sqrt(qb2 - c);
0362         return Real2(-hb - ht, -hb + ht);
0363     }
0364     else
0365     {
0366         return Real2(no_solution_, no_solution_);
0367     }
0368 }
0369
0370 //---------------------------------------------------------------------------//
0371 }  // namespace celeritas