// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2020
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2019.08.13

#pragma once

#include <Mathematics/Math.h>

// Let f(t,A) = sin(t*A)/sin(A).  The slerp of quaternions q0 and q1 is
//   slerp(t,q0,q1) = f(1-t,A)*q0 + f(t,A)*q1.
// Let y = 1-cos(A); we allow A in [0,pi], so y in [0,1].  As a function of y,
// a series representation for f(t,y) is
//   f(t,y) = sum_{i=0}^{infinity} c_{i}(t) y^{i}
// where c_0(t) = t, c_{i}(t) = c_{i-1}(t)*(i^2 - t^2)/(i*(2*i+1)) for i >= 1.
// The c_{i}(t) are polynomials in t of degree 2*i+1.  The document
// https://www.geometrictools/com/Documentation/FastAndAccurateSlerp.pdf
// derives an approximation
//   g(t,y) = sum_{i=0}^{n-1} c_{i}(t) y^{i} + (1+u_n) c_{n}(t) y^n
// which has degree 2*n+1 in t and degree n in y.
//
// Given q0 and q1 such that cos(A) = dot(q0,q1) in [0,1], in which case
// A in [0,pi/2], let qh = (q0+q1)/|q0 + q1| = slerp(1/2,q0,q1).  Note that
// |q0 + q1| = 2*cos(A/2) because
//   sin(A/2)/sin(A) = sin(A/2)/(2*sin(A/2)*cos(A/2)) = 1/(2*cos(A/2))
// The angle between q0 and qh is the same as the angle between qh and q1,
// namely, A/2 in [0,pi/4].  It may be shown that
//                    +--
//   slerp(t,q0,q1) = | slerp(2*t,q0,qh),     0 <= t <= 1/2
//                    | slerp(2*t-1,qh,q1), 1/2 <= t <= 1
//                    +--
// The slerp functions on the right-hand side involve angles in [0,pi/4], so
// the approximation is more accurate for those evaluations than using the
// original.
//
// TODO: The constants in GetEstimate are those of the published paper.
// Modify these to match the aforementioned GT document.

namespace WwiseGTE
{
    template <typename Real>
    class ChebyshevRatio
    {
    public:
        // Compute f(t,A) = sin(t*A)/sin(A), where t in [0,1], A in [0,pi/2],
        // cosA = cos(A), f0 = f(1-t,A), and f1 = f(t,A).
        static void Get(Real t, Real cosA, Real& f0, Real& f1)
        {
            if (cosA < (Real)1)
            {
                // The angle A is in (0,pi/2].
                Real A = std::acos(cosA);
                Real invSinA = (Real)1 / std::sin(A);
                f0 = std::sin(((Real)1 - t) * A) * invSinA;
                f1 = std::sin(t * A) * invSinA;
            }
            else
            {
                // The angle theta is 0.
                f0 = (Real)1 - t;
                f1 = (Real)t;
            }
        }

        // Compute estimates for f(t,y) = sin(t*A)/sin(A), where t in [0,1],
        // A in [0,pi/2], y = 1 - cos(A), f0 is the estimate for f(1-t,y), and
        // f1 is the estimate for f(t,y).  The approximating function is a
        // polynomial of two variables. The template parameter N must be in
        // {1..16}.  The degree in t is 2*N+1 and the degree in Y is N.
        template <int N>
        static void GetEstimate(Real t, Real y, Real & f0, Real & f1)
        {
            static_assert(1 <= N && N <= 16, "Invalid degree.");

            // The ASM output shows that the constants/ in these arrays are
            // loaded to XMM registers as literal values, and only those
            // constants required for the specified degree D are loaded.
            // That is, the compiler does a good job of optimizing the code.

            Real const onePlusMu[16] =
            {
                (Real)1.62943436108234530,
                (Real)1.73965850021313961,
                (Real)1.79701067629566813,
                (Real)1.83291820510335812,
                (Real)1.85772477879039977,
                (Real)1.87596835698904785,
                (Real)1.88998444919711206,
                (Real)1.90110745351730037,
                (Real)1.91015881189952352,
                (Real)1.91767344933047190,
                (Real)1.92401541194159076,
                (Real)1.92944142668012797,
                (Real)1.93413793373091059,
                (Real)1.93824371262559758,
                (Real)1.94186426368404708,
                (Real)1.94508125972497303
            };

            Real const a[16] =
            {
                (N != 1 ? (Real)1 : onePlusMu[0]) / ((Real)1 * (Real)3),
                (N != 2 ? (Real)1 : onePlusMu[1]) / ((Real)2 * (Real)5),
                (N != 3 ? (Real)1 : onePlusMu[2]) / ((Real)3 * (Real)7),
                (N != 4 ? (Real)1 : onePlusMu[3]) / ((Real)4 * (Real)9),
                (N != 5 ? (Real)1 : onePlusMu[4]) / ((Real)5 * (Real)11),
                (N != 6 ? (Real)1 : onePlusMu[5]) / ((Real)6 * (Real)13),
                (N != 7 ? (Real)1 : onePlusMu[6]) / ((Real)7 * (Real)15),
                (N != 8 ? (Real)1 : onePlusMu[7]) / ((Real)8 * (Real)17),
                (N != 9 ? (Real)1 : onePlusMu[8]) / ((Real)9 * (Real)19),
                (N != 10 ? (Real)1 : onePlusMu[9]) / ((Real)10 * (Real)21),
                (N != 11 ? (Real)1 : onePlusMu[10]) / ((Real)11 * (Real)23),
                (N != 12 ? (Real)1 : onePlusMu[11]) / ((Real)12 * (Real)25),
                (N != 13 ? (Real)1 : onePlusMu[12]) / ((Real)13 * (Real)27),
                (N != 14 ? (Real)1 : onePlusMu[13]) / ((Real)14 * (Real)29),
                (N != 15 ? (Real)1 : onePlusMu[14]) / ((Real)15 * (Real)31),
                (N != 16 ? (Real)1 : onePlusMu[15]) / ((Real)16 * (Real)33)
            };

            Real const b[16] =
            {
                (N != 1 ? (Real)1 : onePlusMu[0]) * (Real)1 / (Real)3,
                (N != 2 ? (Real)1 : onePlusMu[1]) * (Real)2 / (Real)5,
                (N != 3 ? (Real)1 : onePlusMu[2]) * (Real)3 / (Real)7,
                (N != 4 ? (Real)1 : onePlusMu[3]) * (Real)4 / (Real)9,
                (N != 5 ? (Real)1 : onePlusMu[4]) * (Real)5 / (Real)11,
                (N != 6 ? (Real)1 : onePlusMu[5]) * (Real)6 / (Real)13,
                (N != 7 ? (Real)1 : onePlusMu[6]) * (Real)7 / (Real)15,
                (N != 8 ? (Real)1 : onePlusMu[7]) * (Real)8 / (Real)17,
                (N != 9 ? (Real)1 : onePlusMu[8]) * (Real)9 / (Real)19,
                (N != 10 ? (Real)1 : onePlusMu[9]) * (Real)10 / (Real)21,
                (N != 11 ? (Real)1 : onePlusMu[10]) * (Real)11 / (Real)23,
                (N != 12 ? (Real)1 : onePlusMu[11]) * (Real)12 / (Real)25,
                (N != 13 ? (Real)1 : onePlusMu[12]) * (Real)13 / (Real)27,
                (N != 14 ? (Real)1 : onePlusMu[13]) * (Real)14 / (Real)29,
                (N != 15 ? (Real)1 : onePlusMu[14]) * (Real)15 / (Real)31,
                (N != 16 ? (Real)1 : onePlusMu[15]) * (Real)16 / (Real)33
            };

            Real term0 = (Real)1 - t, term1 = t;
            Real sqr0 = term0 * term0, sqr1 = term1 * term1;
            f0 = term0;
            f1 = term1;
            for (int i = 0; i < N; ++i)
            {
                term0 *= (b[i] - a[i] * sqr0) * y;
                term1 *= (b[i] - a[i] * sqr1) * y;
                f0 += term0;
                f1 += term1;
            }
        }
    };
}