APInterval.h 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434
  1. // David Eberly, Geometric Tools, Redmond WA 98052
  2. // Copyright (c) 1998-2020
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // https://www.boost.org/LICENSE_1_0.txt
  5. // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
  6. // Version: 4.1.2019.10.17
  7. #pragma once
  8. #include <Mathematics/Logger.h>
  9. #include <Mathematics/ArbitraryPrecision.h>
  10. // The interval [e0,e1] must satisfy e0 <= e1. Expose this define to trap
  11. // invalid construction where e0 > e1.
  12. #define GTE_THROW_ON_INVALID_APINTERVAL
  13. namespace WwiseGTE
  14. {
  15. // The APType must be an arbitrary-precision type.
  16. template <typename APType>
  17. class APInterval
  18. {
  19. public:
  20. // Construction. This is the only way to create an interval. All such
  21. // intervals are immutable once created. The constructor
  22. // APInterval(APType) is used to create the degenerate interval [e,e].
  23. APInterval()
  24. :
  25. mEndpoints{ static_cast<APType>(0), static_cast<APType>(0) }
  26. {
  27. static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
  28. }
  29. APInterval(APInterval const& other)
  30. :
  31. mEndpoints(other.mEndpoints)
  32. {
  33. static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
  34. }
  35. explicit APInterval(APType e)
  36. :
  37. mEndpoints{ e, e }
  38. {
  39. static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
  40. }
  41. APInterval(APType e0, APType e1)
  42. :
  43. mEndpoints{ e0, e1 }
  44. {
  45. static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
  46. #if defined(GTE_THROW_ON_INVALID_APINTERVAL)
  47. LogAssert(mEndpoints[0] <= mEndpoints[1], "Invalid interval.");
  48. #endif
  49. }
  50. APInterval(std::array<APType, 2> const& endpoint)
  51. :
  52. mEndpoints(endpoint)
  53. {
  54. static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
  55. #if defined(GTE_THROW_ON_INVALID_APINTERVAL)
  56. LogAssert(mEndpoints[0] <= mEndpoints[1], "Invalid interval.");
  57. #endif
  58. }
  59. APInterval& operator=(APInterval const& other)
  60. {
  61. static_assert(WwiseGTE::is_arbitrary_precision<APType>::value, "Invalid type.");
  62. mEndpoints = other.mEndpoints;
  63. return *this;
  64. }
  65. // Member access. It is only possible to read the endpoints. You
  66. // cannot modify the endpoints outside the arithmetic operations.
  67. inline APType operator[](size_t i) const
  68. {
  69. return mEndpoints[i];
  70. }
  71. inline std::array<APType, 2> GetEndpoints() const
  72. {
  73. return mEndpoints;
  74. }
  75. // Arithmetic operations to compute intervals at the leaf nodes of
  76. // an expression tree. Such nodes correspond to the raw floating-point
  77. // variables of the expression. The non-class operators defined after
  78. // the class definition are used to compute intervals at the interior
  79. // nodes of the expression tree.
  80. inline static APInterval Add(APType u, APType v)
  81. {
  82. APInterval w;
  83. w.mEndpoints[0] = u + v;
  84. w.mEndpoints[1] = w.mEndpoints[0];
  85. return w;
  86. }
  87. inline static APInterval Sub(APType u, APType v)
  88. {
  89. APInterval w;
  90. w.mEndpoints[0] = u - v;
  91. w.mEndpoints[1] = w.mEndpoints[0];
  92. return w;
  93. }
  94. inline static APInterval Mul(APType u, APType v)
  95. {
  96. APInterval w;
  97. w.mEndpoints[0] = u * v;
  98. w.mEndpoints[1] = w.mEndpoints[0];
  99. return w;
  100. }
  101. template <typename Dummy = APType>
  102. inline static
  103. typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
  104. Div(APType u, APType v)
  105. {
  106. APType const zero = static_cast<APType>(0);
  107. if (v != zero)
  108. {
  109. APInterval w;
  110. w.mEndpoints[0] = u / v;
  111. w.mEndpoints[1] = w.mEndpoints[0];
  112. return w;
  113. }
  114. else
  115. {
  116. // Division by zero does not lead to a determinate interval.
  117. // Just return the entire set of real numbers.
  118. return Reals();
  119. }
  120. }
  121. private:
  122. std::array<APType, 2> mEndpoints;
  123. public:
  124. // FOR INTERNAL USE ONLY. These are used by the non-class operators
  125. // defined after the class definition.
  126. inline static APInterval Add(APType u0, APType u1, APType v0, APType v1)
  127. {
  128. APInterval w;
  129. w.mEndpoints[0] = u0 + v0;
  130. w.mEndpoints[1] = u1 + v1;
  131. return w;
  132. }
  133. inline static APInterval Sub(APType u0, APType u1, APType v0, APType v1)
  134. {
  135. APInterval w;
  136. w.mEndpoints[0] = u0 - v1;
  137. w.mEndpoints[1] = u1 - v0;
  138. return w;
  139. }
  140. inline static APInterval Mul(APType u0, APType u1, APType v0, APType v1)
  141. {
  142. APInterval w;
  143. w.mEndpoints[0] = u0 * v0;
  144. w.mEndpoints[1] = u1 * v1;
  145. return w;
  146. }
  147. inline static APInterval Mul2(APType u0, APType u1, APType v0, APType v1)
  148. {
  149. APType u0mv1 = u0 * v1;
  150. APType u1mv0 = u1 * v0;
  151. APType u0mv0 = u0 * v0;
  152. APType u1mv1 = u1 * v1;
  153. return APInterval<APType>(std::min(u0mv1, u1mv0), std::max(u0mv0, u1mv1));
  154. }
  155. template <typename Dummy = APType>
  156. inline static
  157. typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
  158. Div(APType u0, APType u1, APType v0, APType v1)
  159. {
  160. APInterval w;
  161. w.mEndpoints[0] = u0 / v1;
  162. w.mEndpoints[1] = u1 / v0;
  163. return w;
  164. }
  165. template <typename Dummy = APType>
  166. inline static
  167. typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
  168. Reciprocal(APType v0, APType v1)
  169. {
  170. APType const one = static_cast<APType>(1);
  171. APInterval w;
  172. w.mEndpoints[0] = one / v1;
  173. w.mEndpoints[1] = one / v0;
  174. return w;
  175. }
  176. template <typename Dummy = APType>
  177. inline static
  178. typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
  179. ReciprocalDown(APType v)
  180. {
  181. APType recpv = static_cast<APType>(1) / v;
  182. APType posinf(0);
  183. posinf.SetSign(+2);
  184. return APInterval<APType>(recpv, posinf);
  185. }
  186. template <typename Dummy = APType>
  187. inline static
  188. typename std::enable_if<has_division_operator<Dummy>::value, APInterval>::type
  189. ReciprocalUp(APType v)
  190. {
  191. APType recpv = static_cast<APType>(1) / v;
  192. APType neginf(0);
  193. neginf.SetSign(-2);
  194. return APInterval<APType>(neginf, recpv);
  195. }
  196. inline static APInterval Reals()
  197. {
  198. APType posinf(0), neginf(0);
  199. posinf.SetSign(+2);
  200. neginf.SetSign(-2);
  201. return APInterval(neginf, posinf);
  202. }
  203. };
  204. // Unary operations. Negation of [e0,e1] produces [-e1,-e0]. This
  205. // operation needs to be supported in the sense of negating a
  206. // "number" in an arithmetic expression.
  207. template <typename APType>
  208. APInterval<APType> operator+(APInterval<APType> const& u)
  209. {
  210. return u;
  211. }
  212. template <typename APType>
  213. APInterval<APType> operator-(APInterval<APType> const& u)
  214. {
  215. return APInterval<APType>(-u[1], -u[0]);
  216. }
  217. // Addition operations.
  218. template <typename APType>
  219. APInterval<APType> operator+(APType u, APInterval<APType> const& v)
  220. {
  221. return APInterval<APType>::Add(u, u, v[0], v[1]);
  222. }
  223. template <typename APType>
  224. APInterval<APType> operator+(APInterval<APType> const& u, APType v)
  225. {
  226. return APInterval<APType>::Add(u[0], u[1], v, v);
  227. }
  228. template <typename APType>
  229. APInterval<APType> operator+(APInterval<APType> const& u, APInterval<APType> const& v)
  230. {
  231. return APInterval<APType>::Add(u[0], u[1], v[0], v[1]);
  232. }
  233. // Subtraction operations.
  234. template <typename APType>
  235. APInterval<APType> operator-(APType u, APInterval<APType> const& v)
  236. {
  237. return APInterval<APType>::Sub(u, u, v[0], v[1]);
  238. }
  239. template <typename APType>
  240. APInterval<APType> operator-(APInterval<APType> const& u, APType v)
  241. {
  242. return APInterval<APType>::Sub(u[0], u[1], v, v);
  243. }
  244. template <typename APType>
  245. APInterval<APType> operator-(APInterval<APType> const& u, APInterval<APType> const& v)
  246. {
  247. return APInterval<APType>::Sub(u[0], u[1], v[0], v[1]);
  248. }
  249. // Multiplication operations.
  250. template <typename APType>
  251. APInterval<APType> operator*(APType u, APInterval<APType> const& v)
  252. {
  253. APType const zero = static_cast<APType>(0);
  254. if (u >= zero)
  255. {
  256. return APInterval<APType>::Mul(u, u, v[0], v[1]);
  257. }
  258. else
  259. {
  260. return APInterval<APType>::Mul(u, u, v[1], v[0]);
  261. }
  262. }
  263. template <typename APType>
  264. APInterval<APType> operator*(APInterval<APType> const& u, APType v)
  265. {
  266. APType const zero = static_cast<APType>(0);
  267. if (v >= zero)
  268. {
  269. return APInterval<APType>::Mul(u[0], u[1], v, v);
  270. }
  271. else
  272. {
  273. return APInterval<APType>::Mul(u[1], u[0], v, v);
  274. }
  275. }
  276. template <typename APType>
  277. APInterval<APType> operator*(APInterval<APType> const& u, APInterval<APType> const& v)
  278. {
  279. APType const zero = static_cast<APType>(0);
  280. if (u[0] >= zero)
  281. {
  282. if (v[0] >= zero)
  283. {
  284. return APInterval<APType>::Mul(u[0], u[1], v[0], v[1]);
  285. }
  286. else if (v[1] <= zero)
  287. {
  288. return APInterval<APType>::Mul(u[1], u[0], v[0], v[1]);
  289. }
  290. else // v[0] < 0 < v[1]
  291. {
  292. return APInterval<APType>::Mul(u[1], u[1], v[0], v[1]);
  293. }
  294. }
  295. else if (u[1] <= zero)
  296. {
  297. if (v[0] >= zero)
  298. {
  299. return APInterval<APType>::Mul(u[0], u[1], v[1], v[0]);
  300. }
  301. else if (v[1] <= zero)
  302. {
  303. return APInterval<APType>::Mul(u[1], u[0], v[1], v[0]);
  304. }
  305. else // v[0] < 0 < v[1]
  306. {
  307. return APInterval<APType>::Mul(u[0], u[0], v[1], v[0]);
  308. }
  309. }
  310. else // u[0] < 0 < u[1]
  311. {
  312. if (v[0] >= zero)
  313. {
  314. return APInterval<APType>::Mul(u[0], u[1], v[1], v[1]);
  315. }
  316. else if (v[1] <= zero)
  317. {
  318. return APInterval<APType>::Mul(u[1], u[0], v[0], v[0]);
  319. }
  320. else // v[0] < 0 < v[1]
  321. {
  322. return APInterval<APType>::Mul2(u[0], u[1], v[0], v[1]);
  323. }
  324. }
  325. }
  326. // Division operations. If the divisor interval is [v0,v1] with
  327. // v0 < 0 < v1, then the returned interval is (-infinity,+infinity)
  328. // instead of Union((-infinity,1/v0),(1/v1,+infinity)). An application
  329. // should try to avoid this case by branching based on [v0,0] and [0,v1].
  330. template <typename APType>
  331. APInterval<APType> operator/(APType u, APInterval<APType> const& v)
  332. {
  333. APType const zero = static_cast<APType>(0);
  334. if (v[0] > zero || v[1] < zero)
  335. {
  336. return u * APInterval<APType>::Reciprocal(v[0], v[1]);
  337. }
  338. else
  339. {
  340. if (v[0] == zero)
  341. {
  342. return u * APInterval<APType>::ReciprocalDown(v[1]);
  343. }
  344. else if (v[1] == zero)
  345. {
  346. return u * APInterval<APType>::ReciprocalUp(v[0]);
  347. }
  348. else // v[0] < 0 < v[1]
  349. {
  350. return APInterval<APType>::Reals();
  351. }
  352. }
  353. }
  354. template <typename APType>
  355. APInterval<APType> operator/(APInterval<APType> const& u, APType v)
  356. {
  357. APType const zero = static_cast<APType>(0);
  358. if (v > zero)
  359. {
  360. return APInterval<APType>::Div(u[0], u[1], v, v);
  361. }
  362. else if (v < zero)
  363. {
  364. return APInterval<APType>::Div(u[1], u[0], v, v);
  365. }
  366. else // v = 0
  367. {
  368. return APInterval<APType>::Reals();
  369. }
  370. }
  371. template <typename APType>
  372. APInterval<APType> operator/(APInterval<APType> const& u, APInterval<APType> const& v)
  373. {
  374. APType const zero = static_cast<APType>(0);
  375. if (v[0] > zero || v[1] < zero)
  376. {
  377. return u * APInterval<APType>::Reciprocal(v[0], v[1]);
  378. }
  379. else
  380. {
  381. if (v[0] == zero)
  382. {
  383. return u * APInterval<APType>::ReciprocalDown(v[1]);
  384. }
  385. else if (v[1] == zero)
  386. {
  387. return u * APInterval<APType>::ReciprocalUp(v[0]);
  388. }
  389. else // v[0] < 0 < v[1]
  390. {
  391. return APInterval<APType>::Reals();
  392. }
  393. }
  394. }
  395. }