123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434 |
- // 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.1.2019.10.17
- #pragma once
- #include <Mathematics/Logger.h>
- #include <Mathematics/ArbitraryPrecision.h>
- // The interval [e0,e1] must satisfy e0 <= e1. Expose this define to trap
- // invalid construction where e0 > e1.
- #define GTE_THROW_ON_INVALID_APINTERVAL
- namespace WwiseGTE
- {
- // The APType must be an arbitrary-precision type.
- template <typename APType>
- class APInterval
- {
- public:
- // Construction. This is the only way to create an interval. All such
- // intervals are immutable once created. The constructor
- // APInterval(APType) is used to create the degenerate interval [e,e].
- APInterval()
- :
- mEndpoints{ static_cast<APType>(0), static_cast<APType>(0) }
- {
- static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
- }
- APInterval(APInterval const& other)
- :
- mEndpoints(other.mEndpoints)
- {
- static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
- }
- explicit APInterval(APType e)
- :
- mEndpoints{ e, e }
- {
- static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
- }
- APInterval(APType e0, APType e1)
- :
- mEndpoints{ e0, e1 }
- {
- static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
- #if defined(GTE_THROW_ON_INVALID_APINTERVAL)
- LogAssert(mEndpoints[0] <= mEndpoints[1], "Invalid interval.");
- #endif
- }
- APInterval(std::array<APType, 2> const& endpoint)
- :
- mEndpoints(endpoint)
- {
- static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
- #if defined(GTE_THROW_ON_INVALID_APINTERVAL)
- LogAssert(mEndpoints[0] <= mEndpoints[1], "Invalid interval.");
- #endif
- }
- APInterval& operator=(APInterval const& other)
- {
- static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
- mEndpoints = other.mEndpoints;
- return *this;
- }
- // Member access. It is only possible to read the endpoints. You
- // cannot modify the endpoints outside the arithmetic operations.
- inline APType operator[](size_t i) const
- {
- return mEndpoints[i];
- }
- inline std::array<APType, 2> GetEndpoints() const
- {
- return mEndpoints;
- }
- // Arithmetic operations to compute intervals at the leaf nodes of
- // an expression tree. Such nodes correspond to the raw floating-point
- // variables of the expression. The non-class operators defined after
- // the class definition are used to compute intervals at the interior
- // nodes of the expression tree.
- inline static APInterval Add(APType u, APType v)
- {
- APInterval w;
- w.mEndpoints[0] = u + v;
- w.mEndpoints[1] = w.mEndpoints[0];
- return w;
- }
- inline static APInterval Sub(APType u, APType v)
- {
- APInterval w;
- w.mEndpoints[0] = u - v;
- w.mEndpoints[1] = w.mEndpoints[0];
- return w;
- }
- inline static APInterval Mul(APType u, APType v)
- {
- APInterval w;
- w.mEndpoints[0] = u * v;
- w.mEndpoints[1] = w.mEndpoints[0];
- return w;
- }
- template <typename Dummy = APType>
- inline static
- typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
- Div(APType u, APType v)
- {
- APType const zero = static_cast<APType>(0);
- if (v != zero)
- {
- APInterval w;
- w.mEndpoints[0] = u / v;
- w.mEndpoints[1] = w.mEndpoints[0];
- return w;
- }
- else
- {
- // Division by zero does not lead to a determinate interval.
- // Just return the entire set of real numbers.
- return Reals();
- }
- }
- private:
- std::array<APType, 2> mEndpoints;
- public:
- // FOR INTERNAL USE ONLY. These are used by the non-class operators
- // defined after the class definition.
- inline static APInterval Add(APType u0, APType u1, APType v0, APType v1)
- {
- APInterval w;
- w.mEndpoints[0] = u0 + v0;
- w.mEndpoints[1] = u1 + v1;
- return w;
- }
- inline static APInterval Sub(APType u0, APType u1, APType v0, APType v1)
- {
- APInterval w;
- w.mEndpoints[0] = u0 - v1;
- w.mEndpoints[1] = u1 - v0;
- return w;
- }
- inline static APInterval Mul(APType u0, APType u1, APType v0, APType v1)
- {
- APInterval w;
- w.mEndpoints[0] = u0 * v0;
- w.mEndpoints[1] = u1 * v1;
- return w;
- }
- inline static APInterval Mul2(APType u0, APType u1, APType v0, APType v1)
- {
- APType u0mv1 = u0 * v1;
- APType u1mv0 = u1 * v0;
- APType u0mv0 = u0 * v0;
- APType u1mv1 = u1 * v1;
- return APInterval<APType>(std::min(u0mv1, u1mv0), std::max(u0mv0, u1mv1));
- }
- template <typename Dummy = APType>
- inline static
- typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
- Div(APType u0, APType u1, APType v0, APType v1)
- {
- APInterval w;
- w.mEndpoints[0] = u0 / v1;
- w.mEndpoints[1] = u1 / v0;
- return w;
- }
- template <typename Dummy = APType>
- inline static
- typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
- Reciprocal(APType v0, APType v1)
- {
- APType const one = static_cast<APType>(1);
- APInterval w;
- w.mEndpoints[0] = one / v1;
- w.mEndpoints[1] = one / v0;
- return w;
- }
- template <typename Dummy = APType>
- inline static
- typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
- ReciprocalDown(APType v)
- {
- APType recpv = static_cast<APType>(1) / v;
- APType posinf(0);
- posinf.SetSign(+2);
- return APInterval<APType>(recpv, posinf);
- }
- template <typename Dummy = APType>
- inline static
- typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
- ReciprocalUp(APType v)
- {
- APType recpv = static_cast<APType>(1) / v;
- APType neginf(0);
- neginf.SetSign(-2);
- return APInterval<APType>(neginf, recpv);
- }
- inline static APInterval Reals()
- {
- APType posinf(0), neginf(0);
- posinf.SetSign(+2);
- neginf.SetSign(-2);
- return APInterval(neginf, posinf);
- }
- };
- // Unary operations. Negation of [e0,e1] produces [-e1,-e0]. This
- // operation needs to be supported in the sense of negating a
- // "number" in an arithmetic expression.
- template <typename APType>
- APInterval<APType> operator+(APInterval<APType> const& u)
- {
- return u;
- }
- template <typename APType>
- APInterval<APType> operator-(APInterval<APType> const& u)
- {
- return APInterval<APType>(-u[1], -u[0]);
- }
- // Addition operations.
- template <typename APType>
- APInterval<APType> operator+(APType u, APInterval<APType> const& v)
- {
- return APInterval<APType>::Add(u, u, v[0], v[1]);
- }
- template <typename APType>
- APInterval<APType> operator+(APInterval<APType> const& u, APType v)
- {
- return APInterval<APType>::Add(u[0], u[1], v, v);
- }
- template <typename APType>
- APInterval<APType> operator+(APInterval<APType> const& u, APInterval<APType> const& v)
- {
- return APInterval<APType>::Add(u[0], u[1], v[0], v[1]);
- }
- // Subtraction operations.
- template <typename APType>
- APInterval<APType> operator-(APType u, APInterval<APType> const& v)
- {
- return APInterval<APType>::Sub(u, u, v[0], v[1]);
- }
- template <typename APType>
- APInterval<APType> operator-(APInterval<APType> const& u, APType v)
- {
- return APInterval<APType>::Sub(u[0], u[1], v, v);
- }
- template <typename APType>
- APInterval<APType> operator-(APInterval<APType> const& u, APInterval<APType> const& v)
- {
- return APInterval<APType>::Sub(u[0], u[1], v[0], v[1]);
- }
- // Multiplication operations.
- template <typename APType>
- APInterval<APType> operator*(APType u, APInterval<APType> const& v)
- {
- APType const zero = static_cast<APType>(0);
- if (u >= zero)
- {
- return APInterval<APType>::Mul(u, u, v[0], v[1]);
- }
- else
- {
- return APInterval<APType>::Mul(u, u, v[1], v[0]);
- }
- }
- template <typename APType>
- APInterval<APType> operator*(APInterval<APType> const& u, APType v)
- {
- APType const zero = static_cast<APType>(0);
- if (v >= zero)
- {
- return APInterval<APType>::Mul(u[0], u[1], v, v);
- }
- else
- {
- return APInterval<APType>::Mul(u[1], u[0], v, v);
- }
- }
- template <typename APType>
- APInterval<APType> operator*(APInterval<APType> const& u, APInterval<APType> const& v)
- {
- APType const zero = static_cast<APType>(0);
- if (u[0] >= zero)
- {
- if (v[0] >= zero)
- {
- return APInterval<APType>::Mul(u[0], u[1], v[0], v[1]);
- }
- else if (v[1] <= zero)
- {
- return APInterval<APType>::Mul(u[1], u[0], v[0], v[1]);
- }
- else // v[0] < 0 < v[1]
- {
- return APInterval<APType>::Mul(u[1], u[1], v[0], v[1]);
- }
- }
- else if (u[1] <= zero)
- {
- if (v[0] >= zero)
- {
- return APInterval<APType>::Mul(u[0], u[1], v[1], v[0]);
- }
- else if (v[1] <= zero)
- {
- return APInterval<APType>::Mul(u[1], u[0], v[1], v[0]);
- }
- else // v[0] < 0 < v[1]
- {
- return APInterval<APType>::Mul(u[0], u[0], v[1], v[0]);
- }
- }
- else // u[0] < 0 < u[1]
- {
- if (v[0] >= zero)
- {
- return APInterval<APType>::Mul(u[0], u[1], v[1], v[1]);
- }
- else if (v[1] <= zero)
- {
- return APInterval<APType>::Mul(u[1], u[0], v[0], v[0]);
- }
- else // v[0] < 0 < v[1]
- {
- return APInterval<APType>::Mul2(u[0], u[1], v[0], v[1]);
- }
- }
- }
- // Division operations. If the divisor interval is [v0,v1] with
- // v0 < 0 < v1, then the returned interval is (-infinity,+infinity)
- // instead of Union((-infinity,1/v0),(1/v1,+infinity)). An application
- // should try to avoid this case by branching based on [v0,0] and [0,v1].
- template <typename APType>
- APInterval<APType> operator/(APType u, APInterval<APType> const& v)
- {
- APType const zero = static_cast<APType>(0);
- if (v[0] > zero || v[1] < zero)
- {
- return u * APInterval<APType>::Reciprocal(v[0], v[1]);
- }
- else
- {
- if (v[0] == zero)
- {
- return u * APInterval<APType>::ReciprocalDown(v[1]);
- }
- else if (v[1] == zero)
- {
- return u * APInterval<APType>::ReciprocalUp(v[0]);
- }
- else // v[0] < 0 < v[1]
- {
- return APInterval<APType>::Reals();
- }
- }
- }
- template <typename APType>
- APInterval<APType> operator/(APInterval<APType> const& u, APType v)
- {
- APType const zero = static_cast<APType>(0);
- if (v > zero)
- {
- return APInterval<APType>::Div(u[0], u[1], v, v);
- }
- else if (v < zero)
- {
- return APInterval<APType>::Div(u[1], u[0], v, v);
- }
- else // v = 0
- {
- return APInterval<APType>::Reals();
- }
- }
- template <typename APType>
- APInterval<APType> operator/(APInterval<APType> const& u, APInterval<APType> const& v)
- {
- APType const zero = static_cast<APType>(0);
- if (v[0] > zero || v[1] < zero)
- {
- return u * APInterval<APType>::Reciprocal(v[0], v[1]);
- }
- else
- {
- if (v[0] == zero)
- {
- return u * APInterval<APType>::ReciprocalDown(v[1]);
- }
- else if (v[1] == zero)
- {
- return u * APInterval<APType>::ReciprocalUp(v[0]);
- }
- else // v[0] < 0 < v[1]
- {
- return APInterval<APType>::Reals();
- }
- }
- }
- }
|