diff --git a/.gitmodules b/.gitmodules index 1eeade25..a403d43d 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,13 +2,15 @@ [submodule "core"] path = lol-core url = ../../lolengine/lol.git - branch = core/master + branch = core +# A very large repository with Windows binaries [submodule "external"] path = external url = ../../lolengine/ext-binaries.git branch = master +# Some libraries we embed in the engine [submodule "imgui"] path = src/3rdparty/imgui url = ../../lolengine/ext-imgui.git @@ -21,14 +23,6 @@ path = src/3rdparty/lua url = ../../lolengine/ext-lua.git branch = lol -[submodule "pegtl"] - path = src/3rdparty/pegtl - url = ../../lolengine/ext-pegtl.git - branch = lol -[submodule "mingw-std-threads"] - path = src/3rdparty/mingw-std-threads - url = ../../lolengine/ext-mingw-std-threads.git - branch = master [submodule "cpp-httplib"] path = src/3rdparty/cpp-httplib url = ../../lolengine/ext-cpp-httplib.git diff --git a/doc/tutorial/11_fractal.cpp b/doc/tutorial/11_fractal.cpp index 84082e00..273e2c9f 100644 --- a/doc/tutorial/11_fractal.cpp +++ b/doc/tutorial/11_fractal.cpp @@ -99,7 +99,7 @@ public: m_aabb.aa = m_position; m_aabb.bb = vec3((vec2)m_window_size, 0); - if (has_threads()) + if (thread::has_threads()) { // Spawn worker threads and wait for their readiness for (int i = 0; i < MAX_THREADS; i++) @@ -111,7 +111,7 @@ public: ~Fractal() { - if (has_threads()) + if (thread::has_threads()) { // Signal worker threads for completion and wait for them to quit for (int i = 0; i < MAX_THREADS; i++) @@ -291,7 +291,7 @@ public: for (int i = 0; i < m_size.y; i += MAX_LINES * 2) { - if (has_threads()) + if (thread::has_threads()) m_jobqueue.push(i); else DoWork(i); @@ -473,7 +473,7 @@ public: if (m_dirty[m_frame]) { - if (has_threads()) + if (thread::has_threads()) { for (int i = 0; i < m_size.y; i += MAX_LINES * 2) m_donequeue.pop(); diff --git a/src/lol/base/features.h b/legacy/features.h similarity index 100% rename from src/lol/base/features.h rename to legacy/features.h diff --git a/lol-core b/lol-core index 29336e3d..f30c17f1 160000 --- a/lol-core +++ b/lol-core @@ -1 +1 @@ -Subproject commit 29336e3d2fcd112d6f66bc5d8569876ccb22b6ec +Subproject commit f30c17f180002170d3d25bcc133b6792175c8bb7 diff --git a/src/3rdparty/Makefile.am b/src/3rdparty/Makefile.am index e1528f34..817401d1 100644 --- a/src/3rdparty/Makefile.am +++ b/src/3rdparty/Makefile.am @@ -13,13 +13,10 @@ liblol_lua_a_CPPFLAGS = $(AM_CPPFLAGS) -DLUA_ANSI $(disable_cflags_lua) include lol-lua.am -EXTRA_DIST += $(imgui_sources) $(mingw_std_threads_sources) $(pegtl_sources) \ - $(cpp_httplib_sources) +EXTRA_DIST += $(imgui_sources) $(cpp_httplib_sources) EXTRA_DIST += lol-lua.vcxproj lol-lua.vcxproj.filters include lol-cpp-httplib.am -include lol-mingw-std-threads.am -include lol-pegtl.am # XXX: this is included by the parent Makefile instead #include lol-imgui.am diff --git a/src/3rdparty/mingw-std-threads b/src/3rdparty/mingw-std-threads deleted file mode 160000 index 86ca7dc8..00000000 --- a/src/3rdparty/mingw-std-threads +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 86ca7dc89659a556dcd1b6f954ed8f2799e23e2f diff --git a/src/3rdparty/pegtl b/src/3rdparty/pegtl deleted file mode 160000 index 5e9e1e8c..00000000 --- a/src/3rdparty/pegtl +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 5e9e1e8c2792293efced6998ae8e28ac46d93c03 diff --git a/src/Makefile.am b/src/Makefile.am index e0e18984..13b83dc6 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -31,15 +31,15 @@ liblol_core_headers = \ lol/lua.h \ \ lol/base/all.h \ - lol/base/avl_tree.h lol/base/features.h lol/base/tuple.h lol/base/types.h \ - lol/base/array.h lol/base/assert.h lol/base/string.h lol/base/map.h \ - lol/base/enum.h lol/base/log.h \ + lol/base/avl_tree.h ../legacy/features.h lol/base/tuple.h lol/base/types.h \ + lol/base/array.h lol/base/assert.h lol/base/map.h lol/base/enum.h \ + lol/base/log.h \ \ lol/math/all.h \ - lol/math/functions.h lol/math/vector.h lol/math/half.h lol/math/real.h \ - lol/math/geometry.h lol/math/interp.h lol/math/rand.h lol/math/arraynd.h \ + lol/math/functions.h lol/math/vector.h lol/math/half.h \ + lol/math/geometry.h lol/math/interp.h lol/math/arraynd.h \ lol/math/constants.h lol/math/matrix.h lol/math/ops.h \ - lol/math/transform.h lol/math/polynomial.h lol/math/bigint.h \ + lol/math/transform.h lol/math/bigint.h \ lol/math/noise/gradient.h lol/math/noise/perlin.h \ lol/math/noise/simplex.h \ \ @@ -53,8 +53,7 @@ liblol_core_headers = \ lol/engine/tickable.h \ \ lol/sys/all.h \ - lol/sys/init.h lol/sys/file.h lol/sys/getopt.h lol/sys/thread.h \ - lol/sys/timer.h \ + lol/sys/init.h lol/sys/file.h \ \ lol/image/all.h \ lol/image/pixel.h lol/image/color.h lol/image/image.h \ @@ -94,7 +93,7 @@ liblol_core_sources = \ base/assert.cpp base/features.cpp base/log.cpp base/string.cpp \ \ math/vector.cpp math/matrix.cpp math/transform.cpp math/half.cpp \ - math/geometry.cpp math/real.cpp \ + math/geometry.cpp \ \ gpu/shader.cpp gpu/indexbuffer.cpp gpu/vertexbuffer.cpp \ gpu/framebuffer.cpp gpu/texture.cpp gpu/renderer.cpp \ @@ -117,7 +116,7 @@ liblol_core_sources = \ mesh/mesh.cpp mesh/mesh.h \ mesh/primitivemesh.cpp mesh/primitivemesh.h \ \ - sys/init.cpp sys/file.cpp sys/hacks.cpp sys/getopt.cpp \ + sys/init.cpp sys/file.cpp sys/hacks.cpp \ \ image/resource.cpp image/resource-private.h \ image/image.cpp image/image-private.h image/kernel.cpp image/pixel.cpp \ diff --git a/src/gpu/shader.cpp b/src/gpu/shader.cpp index 2589b690..11a774e0 100644 --- a/src/gpu/shader.cpp +++ b/src/gpu/shader.cpp @@ -25,7 +25,7 @@ # undef WIN32_LEAN_AND_MEAN #endif -#include "tao/pegtl.hpp" +#include #include "lolgl.h" diff --git a/src/lol-core.vcxproj b/src/lol-core.vcxproj index 74f0fc21..edb857a5 100644 --- a/src/lol-core.vcxproj +++ b/src/lol-core.vcxproj @@ -178,7 +178,6 @@ - @@ -202,7 +201,6 @@ - @@ -256,10 +254,9 @@ - + - @@ -299,9 +296,6 @@ - - - @@ -309,10 +303,7 @@ - - - diff --git a/src/lol-core.vcxproj.filters b/src/lol-core.vcxproj.filters index 3f37d0d7..579a2be2 100644 --- a/src/lol-core.vcxproj.filters +++ b/src/lol-core.vcxproj.filters @@ -210,9 +210,6 @@ math - - math - math @@ -236,9 +233,6 @@ sys - - sys - sys @@ -356,7 +350,7 @@ lol\base - + lol\base @@ -365,9 +359,6 @@ lol\base - - lol\base - lol\base @@ -482,15 +473,6 @@ lol\math - - lol\math - - - lol\math - - - lol\math - lol\math @@ -512,18 +494,9 @@ lol\sys - - lol\sys - lol\sys - - lol\sys - - - lol\sys - mesh diff --git a/src/lol/base/all.h b/src/lol/base/all.h index a395de71..2ae2adde 100644 --- a/src/lol/base/all.h +++ b/src/lol/base/all.h @@ -12,7 +12,8 @@ #pragma once -#include +#include + #include #include #include diff --git a/src/lol/base/string.h b/src/lol/base/string.h deleted file mode 100644 index d32efd0f..00000000 --- a/src/lol/base/string.h +++ /dev/null @@ -1,46 +0,0 @@ -// -// Lol Engine -// -// Copyright © 2010—2018 Sam Hocevar -// © 2013—2015 Benjamin “Touky” Huet -// -// Lol Engine is free software. It comes without any warranty, to -// the extent permitted by applicable law. You can redistribute it -// and/or modify it under the terms of the Do What the Fuck You Want -// to Public License, Version 2, as published by the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The string tools -// ---------------- -// Contains some utilities to work with std::string objects. -// - -#include - -namespace lol -{ - -/* Split a string along a single separator */ -array split(std::string const &s, char sep = '\n'); - -/* Split a string along multiple separators */ -array split(std::string const &s, std::string const &seps); - -/* Check whether a string starts or ends with a given substring */ -bool starts_with(std::string const &s, std::string const &prefix); -bool ends_with(std::string const &s, std::string const &suffix); - -/* Convert a string to lowercase or uppercase */ -std::string tolower(std::string const &s); -std::string toupper(std::string const &s); - -/* Format a string, printf-style */ -std::string format(char const *format, ...) LOL_ATTR_FORMAT(1, 2); -std::string vformat(char const *format, va_list ap); - -} /* namespace lol */ - diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h deleted file mode 100644 index e55f4126..00000000 --- a/src/lol/math/polynomial.h +++ /dev/null @@ -1,430 +0,0 @@ -// -// Lol Engine -// -// Copyright: (c) 2010-2014 Sam Hocevar -// This program is free software; you can redistribute it and/or -// modify it under the terms of the Do What The Fuck You Want To -// Public License, Version 2, as published by Sam Hocevar. See -// http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The polynomial class -// -------------------- -// -// The data structure is a simple dynamic array of scalars, with the -// added guarantee that the leading coefficient is always non-zero. -// - -#include - - -namespace lol -{ - -template -struct LOL_ATTR_NODISCARD polynomial -{ - /* The zero polynomial */ - explicit inline polynomial() {} - - /* A constant polynomial */ - explicit inline polynomial(T const &a) - { - m_coefficients.push(a); - reduce_degree(); - } - - /* Create a polynomial from a list of coefficients */ - explicit polynomial(std::initializer_list const &init) - { - for (auto a : init) - m_coefficients.push(a); - - reduce_degree(); - } - - /* Factory for Chebyshev polynomials */ - static polynomial chebyshev(int n) - { - /* Use T0(x) = 1, T1(x) = x, Tn(x) = 2 x Tn-1(x) - Tn-2(x) */ - std::function coeff = [&](int i, int j) - { - if (i > j || i < 0 || ((j ^ i) & 1)) - return (int64_t)0; - if (j < 2) - return (int64_t)1; - return 2 * coeff(i - 1, j - 1) - coeff(i, j - 2); - }; - - polynomial ret; - for (int k = 0; k <= n; ++k) - ret.m_coefficients.push(T(coeff(k, n))); - return ret; - } - - /* We define the zero polynomial (with no coefficients) as having - * degree -1 on purpose. */ - inline int degree() const - { - return (int)m_coefficients.count() - 1; - } - - /* Set one of the polynomial’s coefficients */ - void set(int n, T const &a) - { - ASSERT(n >= 0); - - if (n > degree() && !a) - return; - - while (n > degree()) - m_coefficients.push(T(0)); - - m_coefficients[n] = a; - reduce_degree(); - } - - /* Evaluate polynomial at a given value. This method can also - * be used to compose polynomials, i.e. using another polynomial - * as the value instead of a scalar. */ - template LOL_ATTR_NODISCARD U eval(U x) const - { - U ret(leading()); - for (int i = degree() - 1; i >= 0; --i) - ret = ret * x + U(m_coefficients[i]); - return ret; - } - - polynomial derive() const - { - /* No need to reduce the degree after deriving. */ - polynomial ret; - for (int i = 1; i <= degree(); ++i) - ret.m_coefficients.push(m_coefficients[i] * T(i)); - return ret; - } - - array roots() const - { - /* For now we can only solve polynomials of degrees 0, 1, 2 or 3. */ - ASSERT(degree() >= 0 && degree() <= 3, - "roots() called on polynomial of degree %d", degree()); - - if (degree() == 0) - { - /* p(x) = a > 0 */ - return array {}; - } - else if (degree() == 1) - { - /* p(x) = ax + b */ - T const &a = m_coefficients[1]; - T const &b = m_coefficients[0]; - return array { -b / a }; - } - else if (degree() == 2) - { - /* p(x) = ax² + bx + c */ - T const &a = m_coefficients[2]; - T const &b = m_coefficients[1]; - T const &c = m_coefficients[0]; - - T const k = b / (a + a); - T const delta = k * k - c / a; - - if (delta < T(0)) - { - return array {}; - } - else if (delta > T(0)) - { - T const sqrt_delta = sqrt(delta); - return array { -k - sqrt_delta, -k + sqrt_delta }; - } - else - { - return array { -k }; - } - } - else if (degree() == 3) - { - static T const pi = acos(T(-1)); - - /* p(x) = ax³ + bx² + cx + d */ - T const &a = m_coefficients[3]; - T const &b = m_coefficients[2]; - T const &c = m_coefficients[1]; - T const &d = m_coefficients[0]; - - /* Find k, m, and n such that p(x - k) / a = x³ + mx + n - * q(x) = p(x - k) / a - * = x³ + (b/a-3k) x² + ((c-2bk)/a+3k²) x + (d-ck+bk²)/a-k³ - * - * So we want k = b/3a and thus: - * q(x) = x³ + (c-bk)/a x + (d-ck+2bk²/3)/a - * - * k = b / 3a - * m = (c - bk) / a - * n = (d - ck + 2bk²/3) / a */ - T const k = b / (T(3) * a); - T const m = (c - b * k) / a; - T const n = ((T(2) / T(3) * b * k - c) * k + d) / a; - - /* Let x = u + v, such that 3uv = -m, then: - * q(u + v) = u³ + v³ + n - * - * We then need to solve the following system: - * u³v³ = -m³/27 - * u³ + v³ = -n - * - * which gives : - * Δ = n² + 4m³/27 - * u³ - v³ = √Δ - * u³ + v³ = -n - * - * u³,v³ = (-n ± √Δ) / 2 - */ - T const delta = n * n + T(4) / T(27) * m * m * m; - - /* Because 3uv = -m and m is not complex - * angle(u³) + angle(v³) must equal 0. - * - * This is why we compute u³ and v³ by norm and angle separately - * instead of using a std::complex class */ - T u_norm, u3_angle; - T v_norm, v3_angle; - - if (delta < 0) - { - v_norm = u_norm = sqrt(m / T(-3)); - - u3_angle = atan2(sqrt(-delta), -n); - v3_angle = -u3_angle; - } - else - { - T const sqrt_delta = sqrt(delta); - - u_norm = cbrt(abs(n - sqrt_delta) / T(2)); - v_norm = cbrt(abs(n + sqrt_delta) / T(2)); - - u3_angle = (n >= sqrt_delta) ? pi : 0; - v3_angle = (n <= -sqrt_delta) ? 0 : -pi; - } - - T solutions[3]; - - for (int i : { 0, 1, 2 }) - { - T u_angle = (u3_angle + T(2 * i) * pi) / T(3); - T v_angle = (v3_angle - T(2 * i) * pi) / T(3); - - solutions[i] = u_norm * cos(u_angle) + v_norm * cos(v_angle); - } - - if (a == d && b == c && b == T(3)*a) // triple solution - { - return array { solutions[0] - k }; - } - - // if root of the derivative is also root of the current polynomial, we have a double root. - for (auto root : (polynomial{ c, T(2) * b, T(3) * a }).roots()) - { - if (eval(root) == T(0)) - return array { (solutions[0] + solutions[2]) / T(2) - k, - solutions[1] - k }; - } - - // we have 3 or 1 root depending on delta sign - if (delta > 0) - { - return array { solutions[0] - k }; - } - else // if (delta < 0) 3 real solutions - { - return array { solutions[0] - k, - solutions[1] - k, - solutions[2] - k }; - } - } - - /* It is an error to reach this point. */ - ASSERT(false); - return array {}; - } - - /* Access individual coefficients. This is read-only and returns a - * copy because we cannot let the user mess with the integrity of - * the structure (i.e. the guarantee that the leading coefficient - * remains non-zero). */ - LOL_ATTR_NODISCARD inline T operator[](ptrdiff_t n) const - { - if (n < 0 || n > degree()) - return T(0); - - return m_coefficients[n]; - } - - /* Return the leading coefficient */ - LOL_ATTR_NODISCARD inline T leading() const - { - return (*this)[degree()]; - } - - /* Boolean operations */ - bool operator !() const - { - return degree() < 0; - } - - operator bool() const - { - return !!*this; - } - - /* Unary plus */ - polynomial operator +() const - { - return *this; - } - - /* Unary minus */ - polynomial operator -() const - { - polynomial ret; - for (auto a : m_coefficients) - ret.m_coefficients.push(-a); - return ret; - } - - /* Add two polynomials */ - polynomial &operator +=(polynomial const &p) - { - int min_degree = lol::min(p.degree(), degree()); - - for (int i = 0; i <= min_degree; ++i) - m_coefficients[i] += p[i]; - - for (int i = min_degree + 1; i <= p.degree(); ++i) - m_coefficients.push(p[i]); - - reduce_degree(); - return *this; - } - - polynomial operator +(polynomial const &p) const - { - return polynomial(*this) += p; - } - - /* Subtract two polynomials */ - polynomial &operator -=(polynomial const &p) - { - return *this += -p; - } - - polynomial operator -(polynomial const &p) const - { - return polynomial(*this) += -p; - } - - /* Multiply polynomial by a scalar */ - polynomial &operator *=(T const &k) - { - for (auto &a : m_coefficients) - a *= k; - reduce_degree(); - return *this; - } - - polynomial operator *(T const &k) const - { - return polynomial(*this) *= k; - } - - /* Divide polynomial by a scalar */ - polynomial &operator /=(T const &k) - { - return *this *= (T(1) / k); - } - - polynomial operator /(T const &k) const - { - return *this * (T(1) / k); - } - - /* Multiply polynomial by another polynomial */ - polynomial &operator *=(polynomial const &p) - { - return *this = *this * p; - } - - polynomial operator *(polynomial const &p) const - { - polynomial ret; - polynomial const &q = *this; - - if (p.degree() >= 0 && q.degree() >= 0) - { - int n = p.degree() + q.degree(); - for (int i = 0; i <= n; ++i) - ret.m_coefficients.push(T(0)); - - for (int i = 0; i <= p.degree(); ++i) - for (int j = 0; j <= q.degree(); ++j) - ret.m_coefficients[i + j] += p[i] * q[j]; - - ret.reduce_degree(); - } - - return ret; - } - - /* Divide a polynomial by another one. There is no /= variant because - * the return value contains both the quotient and the remainder. */ - tuple, polynomial> operator /(polynomial p) const - { - ASSERT(p.degree() >= 0); - - tuple, polynomial> ret; - polynomial "ient = ret.m1; - polynomial &remainder = ret.m2; - - remainder = *this / p.leading(); - p /= p.leading(); - - for (int n = remainder.degree() - p.degree(); n >= 0; --n) - { - quotient.set(n, remainder.leading()); - for (int i = 0; i < p.degree(); ++i) - remainder.m_coefficients[n + i] -= remainder.leading() * p[i]; - (void)remainder.m_coefficients.pop(); - } - - return ret; - } - -private: - /* Enforce the non-zero leading coefficient rule. */ - void reduce_degree() - { - while (m_coefficients.count() && !m_coefficients.last()) - (void)m_coefficients.pop(); - } - - /* The polynomial coefficients */ - array m_coefficients; -}; - -template -polynomial operator *(T const &k, polynomial const &p) -{ - /* We assume k * p is the same as p * k */ - return p * k; -} - -} /* namespace lol */ - diff --git a/src/lol/math/rand.h b/src/lol/math/rand.h deleted file mode 100644 index 6256b90a..00000000 --- a/src/lol/math/rand.h +++ /dev/null @@ -1,126 +0,0 @@ -// -// Lol Engine -// -// Copyright: (c) 2010-2013 Sam Hocevar -// This program is free software; you can redistribute it and/or -// modify it under the terms of the Do What The Fuck You Want To -// Public License, Version 2, as published by Sam Hocevar. See -// http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The Random number generators -// ---------------------------- -// - -#include -#include - -namespace lol -{ - -/* Random number generators */ -template LOL_ATTR_NODISCARD static inline T rand(); -template LOL_ATTR_NODISCARD static inline T rand(T a); -template LOL_ATTR_NODISCARD static inline T rand(T a, T b); - -/* One-value random number generators */ -template LOL_ATTR_NODISCARD static inline T rand(T a) -{ - return a ? rand() % a : T(0); -} - -template<> LOL_ATTR_NODISCARD inline half rand(half a) -{ - float f = (float)std::rand() / (float)RAND_MAX; - return (half)(a * f); -} - -template<> LOL_ATTR_NODISCARD inline float rand(float a) -{ - float f = (float)std::rand() / (float)RAND_MAX; - return a * f; -} - -template<> LOL_ATTR_NODISCARD inline double rand(double a) -{ - double f = (double)std::rand() / (double)RAND_MAX; - return a * f; -} - -template<> LOL_ATTR_NODISCARD inline ldouble rand(ldouble a) -{ - ldouble f = (ldouble)std::rand() / (ldouble)RAND_MAX; - return a * f; -} - -/* Two-value random number generator -- no need for specialisation */ -template LOL_ATTR_NODISCARD static inline T rand(T a, T b) -{ - return a + rand(b - a); -} - -/* Default random number generator */ -template LOL_ATTR_NODISCARD static inline T rand() -{ - switch (sizeof(T)) - { - case 1: - return static_cast(std::rand() & 0x7f); - case 2: - { - uint16_t ret = std::rand(); - if (RAND_MAX < 0x7fff) - ret = (ret << 7) ^ std::rand(); - return static_cast(ret & 0x7fffu); - } - case 4: - { - uint32_t ret = std::rand(); - if (RAND_MAX >= 0xffff) - ret = (ret << 16) ^ std::rand(); - else - { - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - } - return static_cast(ret & 0x7fffffffu); - } - case 8: - { - uint64_t ret = std::rand(); - if (RAND_MAX >= 0xffff) - { - ret = (ret << 16) ^ std::rand(); - ret = (ret << 16) ^ std::rand(); - ret = (ret << 16) ^ std::rand(); - } - else - { - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - ret = (ret << 8) ^ std::rand(); - } - return static_cast(ret & (~(uint64_t)0 >> 1)); - } - default: - ASSERT(false, "rand() doesn’t support types of size %d\n", - (int)sizeof(T)); - return 0; - } -} - -template<> LOL_ATTR_NODISCARD inline half rand() { return rand(1.f); } -template<> LOL_ATTR_NODISCARD inline float rand() { return rand(1.f); } -template<> LOL_ATTR_NODISCARD inline double rand() { return rand(1.0); } -template<> LOL_ATTR_NODISCARD inline ldouble rand() { return rand(1.0); } - -} /* namespace lol */ - diff --git a/src/lol/math/real.h b/src/lol/math/real.h deleted file mode 100644 index 47ecb2de..00000000 --- a/src/lol/math/real.h +++ /dev/null @@ -1,368 +0,0 @@ -// -// Lol Engine -// -// Copyright © 2010—2019 Sam Hocevar -// -// Lol Engine is free software. It comes without any warranty, to -// the extent permitted by applicable law. You can redistribute it -// and/or modify it under the terms of the Do What the Fuck You Want -// to Public License, Version 2, as published by the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The Real class -// -------------- -// - -#include - -#include -#include -#include - -namespace lol -{ - -/* This is OUR namespace. Don't let Windows headers mess with it. */ -#undef min -#undef max - -/* - * The base class for reals. The only real reason for making this a template - * class is so we can have implicit constructors ("real x = 1" works) but - * avoid accidental implicit conversions ("int x = 1; sqrt(x)" will never - * call real::sqrt). - */ -template -class LOL_ATTR_NODISCARD Real -{ -public: - typedef T bigit_t; - typedef int64_t exponent_t; - - Real() = default; - - Real(float f); - Real(double f); - Real(long double f); - Real(int32_t i); - Real(uint32_t i); - Real(int64_t i); - Real(uint64_t i); - - Real(char const *str); - - LOL_ATTR_NODISCARD bool is_zero() const { return m_mantissa.size() == 0; } - LOL_ATTR_NODISCARD bool is_negative() const { return m_sign; } - LOL_ATTR_NODISCARD bool is_nan() const { return m_nan; } - LOL_ATTR_NODISCARD bool is_inf() const { return m_inf; } - - LOL_ATTR_NODISCARD operator float() const; - LOL_ATTR_NODISCARD operator double() const; - LOL_ATTR_NODISCARD operator long double() const; - LOL_ATTR_NODISCARD operator int32_t() const; - LOL_ATTR_NODISCARD operator uint32_t() const; - LOL_ATTR_NODISCARD operator int64_t() const; - LOL_ATTR_NODISCARD operator uint64_t() const; - - Real operator +() const; - Real operator -() const; - Real operator +(Real const &x) const; - Real operator -(Real const &x) const; - Real operator *(Real const &x) const; - Real operator /(Real const &x) const; - Real const &operator +=(Real const &x); - Real const &operator -=(Real const &x); - Real const &operator *=(Real const &x); - Real const &operator /=(Real const &x); - - LOL_ATTR_NODISCARD bool operator ==(Real const &x) const; - LOL_ATTR_NODISCARD bool operator !=(Real const &x) const; - LOL_ATTR_NODISCARD bool operator <(Real const &x) const; - LOL_ATTR_NODISCARD bool operator >(Real const &x) const; - LOL_ATTR_NODISCARD bool operator <=(Real const &x) const; - LOL_ATTR_NODISCARD bool operator >=(Real const &x) const; - - LOL_ATTR_NODISCARD bool operator !() const; - LOL_ATTR_NODISCARD operator bool() const; - - /* Comparison functions */ - template friend Real min(Real const &a, Real const &b); - template friend Real max(Real const &a, Real const &b); - template friend Real clamp(Real const &x, - Real const &a, Real const &b); - - /* Trigonometric functions */ - template friend Real sin(Real const &x); - template friend Real cos(Real const &x); - template friend Real tan(Real const &x); - template friend Real asin(Real const &x); - template friend Real acos(Real const &x); - template friend Real atan(Real const &x); - template friend Real atan2(Real const &y, Real const &x); - - /* Hyperbolic functions */ - template friend Real sinh(Real const &x); - template friend Real cosh(Real const &x); - template friend Real tanh(Real const &x); - - /* Exponential and logarithmic functions */ - template friend Real exp(Real const &x); - template friend Real exp2(Real const &x); - template friend Real erf(Real const &x); - template friend Real log(Real const &x); - template friend Real log2(Real const &x); - template friend Real log10(Real const &x); - - /* Floating-point functions */ - template friend Real frexp(Real const &x, typename Real::exponent_t *exp); - template friend Real ldexp(Real const &x, typename Real::exponent_t exp); - template friend Real modf(Real const &x, Real *iptr); - template friend Real nextafter(Real const &x, Real const &y); - - /* Power functions */ - template friend Real inverse(Real const &x); - template friend Real sqrt(Real const &x); - template friend Real cbrt(Real const &x); - template friend Real pow(Real const &x, Real const &y); - template friend Real gamma(Real const &x); - - /* Rounding, absolute value, remainder etc. */ - template friend Real ceil(Real const &x); - template friend Real copysign(Real const &x, Real const &y); - template friend Real floor(Real const &x); - template friend Real fabs(Real const &x); - template friend Real round(Real const &x); - template friend Real fmod(Real const &x, Real const &y); - - /* Functions inherited from GLSL */ - template friend Real abs(Real const &x); - template friend Real fract(Real const &x); - template friend Real degrees(Real const &x); - template friend Real radians(Real const &x); - - /* Additional functions */ - template friend Real franke(Real const &x, Real const &y); - template friend Real peaks(Real const &x, Real const &y); - - std::string str(int ndigits = 150) const; - std::string xstr() const; - - /* Additional operators using base C++ types */ -#define __LOL_REAL_OP_HELPER_GENERIC(op, type) \ - inline Real operator op(type x) const { return *this op (Real)x; } \ - inline Real const &operator op##=(type x) { return *this = (*this op x); } -#define __LOL_REAL_OP_HELPER_FASTMULDIV(op, type) \ - inline Real operator op(type x) const \ - { \ - Real tmp = *this; return tmp op##= x; \ - } \ - inline Real const &operator op##=(type x) \ - { \ - /* If multiplying or dividing by a power of two, take a shortcut */ \ - if (!is_zero() && x && !(x & (x - 1))) \ - { \ - while (x >>= 1) \ - m_exponent += 1 op 2 - 1; /* 1 if op is *, -1 if op is / */ \ - } \ - else \ - *this = *this op (Real)x; \ - return *this; \ - } -#define __LOL_REAL_OP_HELPER_INT(type) \ - __LOL_REAL_OP_HELPER_GENERIC(+, type) \ - __LOL_REAL_OP_HELPER_GENERIC(-, type) \ - __LOL_REAL_OP_HELPER_FASTMULDIV(*, type) \ - __LOL_REAL_OP_HELPER_FASTMULDIV(/, type) -#define __LOL_REAL_OP_HELPER_FLOAT(type) \ - __LOL_REAL_OP_HELPER_GENERIC(+, type) \ - __LOL_REAL_OP_HELPER_GENERIC(-, type) \ - __LOL_REAL_OP_HELPER_GENERIC(*, type) \ - __LOL_REAL_OP_HELPER_GENERIC(/, type) - - __LOL_REAL_OP_HELPER_INT(int) - __LOL_REAL_OP_HELPER_INT(unsigned int) - __LOL_REAL_OP_HELPER_INT(int64_t) - __LOL_REAL_OP_HELPER_INT(uint64_t) - __LOL_REAL_OP_HELPER_FLOAT(float) - __LOL_REAL_OP_HELPER_FLOAT(double) - __LOL_REAL_OP_HELPER_FLOAT(long double) - - /* Constants */ - static Real const R_0(); - static Real const& R_1(); - static Real const& R_2(); - static Real const& R_3(); - static Real const& R_4(); - static Real const& R_10(); - - static Real const& R_E(); - static Real const& R_LOG2E(); - static Real const& R_LOG10E(); - static Real const& R_LN2(); - static Real const& R_LN10(); - static Real const& R_PI(); - static Real const& R_PI_2(); - static Real const& R_PI_3(); - static Real const& R_PI_4(); - static Real const& R_TAU(); - static Real const& R_1_PI(); - static Real const& R_2_PI(); - static Real const& R_2_SQRTPI(); - static Real const& R_SQRT2(); - static Real const& R_SQRT3(); - static Real const& R_SQRT1_2(); - - static Real const R_INF(); - static Real const R_NAN(); - - static Real const& R_MIN(); - static Real const& R_MAX(); - -private: - std::vector m_mantissa; - exponent_t m_exponent = 0; - bool m_sign = false, m_nan = false, m_inf = false; - -public: - static int DEFAULT_BIGIT_COUNT; - - static inline int bigit_bits() { return 8 * (int)sizeof(bigit_t); } - inline int bigit_count() const { return (int)m_mantissa.size(); } - inline int total_bits() const { return bigit_count() * bigit_bits(); } -}; - -template -std::ostream& operator <<(std::ostream &s, Real const &x); - -/* - * Mandatory forward declarations of template specialisations - */ -template<> real::Real(float f); -template<> real::Real(double f); -template<> real::Real(long double f); -template<> real::Real(int32_t i); -template<> real::Real(uint32_t i); -template<> real::Real(int64_t i); -template<> real::Real(uint64_t i); -template<> real::Real(char const *str); - -template<> LOL_ATTR_NODISCARD real::operator float() const; -template<> LOL_ATTR_NODISCARD real::operator double() const; -template<> LOL_ATTR_NODISCARD real::operator long double() const; -template<> LOL_ATTR_NODISCARD real::operator int32_t() const; -template<> LOL_ATTR_NODISCARD real::operator uint32_t() const; -template<> LOL_ATTR_NODISCARD real::operator int64_t() const; -template<> LOL_ATTR_NODISCARD real::operator uint64_t() const; -template<> real real::operator +() const; -template<> real real::operator -() const; -template<> real real::operator +(real const &x) const; -template<> real real::operator -(real const &x) const; -template<> real real::operator *(real const &x) const; -template<> real real::operator /(real const &x) const; -template<> real const &real::operator +=(real const &x); -template<> real const &real::operator -=(real const &x); -template<> real const &real::operator *=(real const &x); -template<> real const &real::operator /=(real const &x); -template<> LOL_ATTR_NODISCARD bool real::operator ==(real const &x) const; -template<> LOL_ATTR_NODISCARD bool real::operator !=(real const &x) const; -template<> LOL_ATTR_NODISCARD bool real::operator <(real const &x) const; -template<> LOL_ATTR_NODISCARD bool real::operator >(real const &x) const; -template<> LOL_ATTR_NODISCARD bool real::operator <=(real const &x) const; -template<> LOL_ATTR_NODISCARD bool real::operator >=(real const &x) const; -template<> LOL_ATTR_NODISCARD bool real::operator !() const; -template<> LOL_ATTR_NODISCARD real::operator bool() const; - -template Real min(Real const &a, Real const &b); -template Real max(Real const &a, Real const &b); -template Real clamp(Real const &x, - Real const &a, Real const &b); - -template Real sin(Real const &x); -template Real cos(Real const &x); -template Real tan(Real const &x); -template Real asin(Real const &x); -template Real acos(Real const &x); -template Real atan(Real const &x); -template Real atan2(Real const &y, Real const &x); -template Real sinh(Real const &x); -template Real cosh(Real const &x); -template Real tanh(Real const &x); -template Real exp(Real const &x); -template Real exp2(Real const &x); -template Real erf(Real const &x); -template Real log(Real const &x); -template Real log2(Real const &x); -template Real log10(Real const &x); -template Real frexp(Real const &x, typename Real::exponent_t *exp); -template Real ldexp(Real const &x, typename Real::exponent_t exp); -template Real modf(Real const &x, Real *iptr); -template Real nextafter(Real const &x, Real const &y); -template Real inverse(Real const &x); -template Real sqrt(Real const &x); -template Real cbrt(Real const &x); -template Real pow(Real const &x, Real const &y); -template Real gamma(Real const &x); -template Real ceil(Real const &x); -template Real copysign(Real const &x, Real const &y); -template Real floor(Real const &x); -template Real fabs(Real const &x); -template Real round(Real const &x); -template Real fmod(Real const &x, Real const &y); -template Real abs(Real const &x); -template Real fract(Real const &x); -template Real degrees(Real const &x); -template Real radians(Real const &x); -template Real franke(Real const &x, Real const &y); -template Real peaks(Real const &x, Real const &y); - -template<> real min(real const &a, real const &b); -template<> real max(real const &a, real const &b); -template<> real clamp(real const &x, real const &a, real const &b); - -template<> real sin(real const &x); -template<> real cos(real const &x); -template<> real tan(real const &x); -template<> real asin(real const &x); -template<> real acos(real const &x); -template<> real atan(real const &x); -template<> real atan2(real const &y, real const &x); -template<> real sinh(real const &x); -template<> real cosh(real const &x); -template<> real tanh(real const &x); -template<> real exp(real const &x); -template<> real exp2(real const &x); -template<> real erf(real const &x); -template<> real log(real const &x); -template<> real log2(real const &x); -template<> real log10(real const &x); -template<> real frexp(real const &x, real::exponent_t *exp); -template<> real ldexp(real const &x, real::exponent_t exp); -template<> real modf(real const &x, real *iptr); -template<> real nextafter(real const &x, real const &y); -template<> real inverse(real const &x); -template<> real sqrt(real const &x); -template<> real cbrt(real const &x); -template<> real pow(real const &x, real const &y); -template<> real gamma(real const &x); -template<> real ceil(real const &x); -template<> real copysign(real const &x, real const &y); -template<> real floor(real const &x); -template<> real fabs(real const &x); -template<> real round(real const &x); -template<> real fmod(real const &x, real const &y); -template<> real abs(real const &x); -template<> real fract(real const &x); -template<> real degrees(real const &x); -template<> real radians(real const &x); -template<> real franke(real const &x, real const &y); -template<> real peaks(real const &x, real const &y); - -template<> std::string real::str(int ndigits) const; -template<> std::string real::xstr() const; - -} /* namespace lol */ - diff --git a/src/lol/sys/all.h b/src/lol/sys/all.h index 765c54d3..cd76273d 100644 --- a/src/lol/sys/all.h +++ b/src/lol/sys/all.h @@ -12,9 +12,6 @@ #pragma once -#include -#include /* requires thread.h */ -#include #include #include diff --git a/src/lol/sys/getopt.h b/src/lol/sys/getopt.h deleted file mode 100644 index ff7d5430..00000000 --- a/src/lol/sys/getopt.h +++ /dev/null @@ -1,41 +0,0 @@ -// -// Lol Engine -// -// Copyright © 2002—2016 Sam Hocevar -// -// Lol Engine is free software. It comes without any warranty, to -// the extent permitted by applicable law. You can redistribute it -// and/or modify it under the terms of the Do What the Fuck You Want -// to Public License, Version 2, as published by the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The getopt functions -// -------------------- -// - -namespace lol -{ - -class getopt -{ -public: - getopt(int argc, char ** _argv); - getopt(int argc, char * const * _argv); - ~getopt(); - - void add_opt(int short_opt, char const *long_opt, bool has_arg); - int parse(); - - int index; - char *arg; - -private: - std::unique_ptr m_private; -}; - -} - diff --git a/src/lol/sys/thread.h b/src/lol/sys/thread.h deleted file mode 100644 index f52e320d..00000000 --- a/src/lol/sys/thread.h +++ /dev/null @@ -1,222 +0,0 @@ -// -// Lol Engine -// -// Copyright © 2010—2019 Sam Hocevar -// -// Lol Engine is free software. It comes without any warranty, to -// the extent permitted by applicable law. You can redistribute it -// and/or modify it under the terms of the Do What the Fuck You Want -// to Public License, Version 2, as published by the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The Threading classes -// --------------------- -// - -#include - -#include -#include -#include - -/* XXX: workaround for a bug in Visual Studio 2012 and 2013! - * https://connect.microsoft.com/VisualStudio/feedback/details/747145 */ -#if defined(_MSC_VER) && (_MSC_VER < 1900) -# define LOL_VISUAL_STUDIO_BUG_747145_WORKAROUND 1 -# ifdef WIN32_LEAN_AND_MEAN -# include -# else -# define WIN32_LEAN_AND_MEAN 1 -# include -# undef WIN32_LEAN_AND_MEAN -# endif -#endif - -/* XXX: workaround for missing std::thread in mingw */ -#if _GLIBCXX_MUTEX && !_GLIBCXX_HAS_GTHREADS && _WIN32 -# include "mingw.thread.h" -# include "mingw.mutex.h" -# include "mingw.condition_variable.h" -# undef near /* Fuck Microsoft */ -# undef far /* Fuck Microsoft again */ -#endif - -namespace lol -{ - -// This is like std::mutex but we can add debug information to it -class mutex -{ -public: - inline void lock() { m_mutex.lock(); } - - inline bool try_lock() { return m_mutex.try_lock(); } - - inline void unlock() { m_mutex.unlock(); } - -private: - std::mutex m_mutex; -}; - -// A FIFO queue for threads -template -class queue -{ -public: - int size() const { return m_count; } - - // Will block the thread if another has already locked - void push(T value) - { - std::unique_lock uni_lock(m_mutex); - - if (has_threads()) - { - // Wait for the mutex availability or non-fullness - m_full_cond.wait(uni_lock, [&]{ return m_count < CAPACITY; }); - } - - do_push(value); /* Push value */ - - if (has_threads()) - { - // Release lock and notify empty condition var (in that order) - uni_lock.unlock(); - m_empty_cond.notify_one(); - } - } - - // Will not block if another has already locked - bool try_push(T value) - { - std::unique_lock uni_lock(m_mutex, std::defer_lock); - - if (has_threads()) - { - // Try to lock, bail out if we fail - if (!uni_lock.try_lock()) - return false; - } - - if (m_count == CAPACITY) - return false; // unlocks uni_lock - - do_push(value); - - if (has_threads()) - { - // Release lock and notify empty condition var (in that order) - uni_lock.unlock(); - m_empty_cond.notify_one(); - } - - return true; - } - - // Will block the thread if another has already locked - T pop() - { - ASSERT(has_threads(), "Pop should only be used with threads. Use try_pop instead."); - - // Wait for the mutex availability or non-emptiness - std::unique_lock uni_lock(m_mutex); - m_empty_cond.wait(uni_lock, [&]{return m_count > 0; }); - - T ret = do_pop(); - - // Release lock and notify full condition var (in that order) - uni_lock.unlock(); - m_full_cond.notify_one(); - - return ret; - } - - // Will not block if another has already locked - bool try_pop(T &ret) - { - std::unique_lock uni_lock(m_mutex, std::defer_lock); - - if (has_threads()) - { - if (!uni_lock.try_lock()) - return false; - } - - if (m_count == 0) - return false; - - ret = do_pop(); /* Pop value */ - - if (has_threads()) - { - // Release lock and notify full condition var (in that order) - uni_lock.unlock(); - m_full_cond.notify_one(); - } - - return true; - } - - // Inner methods for actual update -private: - void do_push(T &value) - { - m_values[(m_start + m_count) % CAPACITY] = value; - m_count++; - } - - T& do_pop() - { - int idx = m_start; - m_start = (m_start + 1) % CAPACITY; - m_count--; - return m_values[idx]; - } - -private: - static int const CAPACITY = N; - T m_values[CAPACITY]; - int m_start = 0, m_count = 0; - - std::mutex m_mutex; - std::condition_variable m_empty_cond, m_full_cond; -}; - -// Base class for threads -class thread -{ -public: - thread(std::function fn) - : m_function(fn) - { - m_thread = std::thread(trampoline, this); - } - - ~thread() - { -#if LOL_VISUAL_STUDIO_BUG_747145_WORKAROUND - m_thread.detach(); -#else - m_thread.join(); -#endif - } - -private: - static void trampoline(thread *that) - { - that->m_function(that); -#if LOL_VISUAL_STUDIO_BUG_747145_WORKAROUND - ExitThread(0); -#endif - } - - std::thread m_thread; - std::function m_function; -}; - -} /* namespace lol */ - diff --git a/src/math/real.cpp b/src/math/real.cpp deleted file mode 100644 index 639d92d7..00000000 --- a/src/math/real.cpp +++ /dev/null @@ -1,1698 +0,0 @@ -// -// Lol Engine -// -// Copyright © 2010—2019 Sam Hocevar -// -// Lol Engine is free software. It comes without any warranty, to -// the extent permitted by applicable law. You can redistribute it -// and/or modify it under the terms of the Do What the Fuck You Want -// to Public License, Version 2, as published by the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#include - -#include -#include -#include -#include -#include -#include - -namespace lol -{ - -/* - * First handle explicit specialisation of our templates. - */ - -template<> int real::DEFAULT_BIGIT_COUNT = 16; - -/* - * Initialisation order is not important because everything is - * done on demand, but here is the dependency list anyway: - * - fast_log() requires R_1 - * - log() requires R_LN2 - * - inverse() require R_2 - * - exp() requires R_0, R_1, R_LN2 - * - sqrt() requires R_3 - */ - -static real fast_log(real const &x); - -static real load_min(); -static real load_max(); -static real load_pi(); - -/* These getters do not need caching, their return values are small */ -template<> real const real::R_0() { return real(); } -template<> real const real::R_INF() { real ret; ret.m_inf = true; return ret; } -template<> real const real::R_NAN() { real ret; ret.m_nan = true; return ret; } - -#define LOL_CONSTANT_GETTER(name, value) \ - template<> real const& real::name() \ - { \ - static real ret; \ - static int prev_bigit_count = -1; \ - /* If the default bigit count has changed, we must recompute - * the value with the desired precision. */ \ - if (prev_bigit_count != DEFAULT_BIGIT_COUNT) \ - { \ - ret = (value); \ - prev_bigit_count = DEFAULT_BIGIT_COUNT; \ - } \ - return ret; \ - } - -LOL_CONSTANT_GETTER(R_1, real(1.0)); -LOL_CONSTANT_GETTER(R_2, real(2.0)); -LOL_CONSTANT_GETTER(R_3, real(3.0)); -LOL_CONSTANT_GETTER(R_10, real(10.0)); - -LOL_CONSTANT_GETTER(R_MIN, load_min()); -LOL_CONSTANT_GETTER(R_MAX, load_max()); - -LOL_CONSTANT_GETTER(R_LN2, fast_log(R_2())); -LOL_CONSTANT_GETTER(R_LN10, log(R_10())); -LOL_CONSTANT_GETTER(R_LOG2E, inverse(R_LN2())); -LOL_CONSTANT_GETTER(R_LOG10E, inverse(R_LN10())); -LOL_CONSTANT_GETTER(R_E, exp(R_1())); -LOL_CONSTANT_GETTER(R_PI, load_pi()); -LOL_CONSTANT_GETTER(R_PI_2, R_PI() / 2); -LOL_CONSTANT_GETTER(R_PI_3, R_PI() / R_3()); -LOL_CONSTANT_GETTER(R_PI_4, R_PI() / 4); -LOL_CONSTANT_GETTER(R_TAU, R_PI() + R_PI()); -LOL_CONSTANT_GETTER(R_1_PI, inverse(R_PI())); -LOL_CONSTANT_GETTER(R_2_PI, R_1_PI() * 2); -LOL_CONSTANT_GETTER(R_2_SQRTPI, inverse(sqrt(R_PI())) * 2); -LOL_CONSTANT_GETTER(R_SQRT2, sqrt(R_2())); -LOL_CONSTANT_GETTER(R_SQRT3, sqrt(R_3())); -LOL_CONSTANT_GETTER(R_SQRT1_2, R_SQRT2() / 2); - -#undef LOL_CONSTANT_GETTER - -/* - * Now carry on with the rest of the Real class. - */ - -template<> real::Real(int32_t i) { new(this) real((double)i); } -template<> real::Real(uint32_t i) { new(this) real((double)i); } -template<> real::Real(float f) { new(this) real((double)f); } - -template<> real::Real(int64_t i) -{ - new(this) real((uint64_t)lol::abs(i)); - m_sign = i < 0; -} - -template<> real::Real(uint64_t i) -{ - new(this) real(); - if (i) - { - /* Only works with 32-bit bigits for now */ - static_assert(sizeof(bigit_t) == 4, "bigit_t must be 32-bit"); - - int delta = 1; - while ((i >> 63) == 0) - { - i <<= 1; - ++delta; - } - i <<= 1; /* Remove implicit one */ - - m_exponent = 64 - delta; - m_mantissa.resize(DEFAULT_BIGIT_COUNT); - m_mantissa[0] = (bigit_t)(i >> 32); - if (bigit_count() > 1) - m_mantissa[1] = (bigit_t)i; - } -} - -template<> real::Real(double d) -{ - union { double d; uint64_t x; } u = { d }; - - m_sign = bool(u.x >> 63); - - exponent_t exponent = (u.x << 1) >> 53; - - switch (exponent) - { - case 0x00: /* +0 / -0 */ - break; - case 0x7ff: /* Inf/NaN (FIXME: handle NaN!) */ - m_inf = true; - break; - default: - /* Only works with 32-bit bigits for now */ - static_assert(sizeof(bigit_t) == 4, "bigit_t must be 32-bit"); - m_exponent = exponent - ((1 << 10) - 1); - m_mantissa.resize(DEFAULT_BIGIT_COUNT); - m_mantissa[0] = (bigit_t)(u.x >> 20); - if (bigit_count() > 1) - m_mantissa[1] = (bigit_t)(u.x << 12); - break; - } -} - -template<> real::Real(long double f) -{ - /* We don’t know the long double layout, so we get rid of the - * exponent, then load it into a real in two steps. */ - int exponent; - f = frexpl(f, &exponent); - new(this) real(double(f)); - *this += double(f - (long double)*this); - m_exponent += exponent; -} - -template<> real::operator float() const { return (float)(double)*this; } -template<> real::operator int32_t() const { return (int32_t)(double)floor(*this); } -template<> real::operator uint32_t() const { return (uint32_t)(double)floor(*this); } - -template<> real::operator uint64_t() const -{ - uint32_t msb = (uint32_t)ldexp(*this, -32); - uint64_t ret = ((uint64_t)msb << 32) - | (uint32_t)(*this - ldexp((real)msb, 32)); - return ret; -} - -template<> real::operator int64_t() const -{ - /* If number is positive, convert it to uint64_t first. If it is - * negative, switch its sign first. */ - return is_negative() ? -(int64_t)-*this : (int64_t)(uint64_t)*this; -} - -template<> real::operator double() const -{ - union { double d; uint64_t x; } u; - - /* Get sign */ - u.x = (is_negative() ? 1 : 0) << 11; - - /* Compute new exponent (FIXME: handle Inf/NaN) */ - int64_t e = m_exponent + ((1 << 10) - 1); - - if (is_zero()) - u.x <<= 52; - else if (e < 0) /* if exponent underflows, set to zero */ - u.x <<= 52; - else if (e >= 0x7ff) - { - u.x |= 0x7ff; - u.x <<= 52; - } - else - { - u.x |= e; - - /* Store mantissa if necessary */ - u.x <<= 32; - if (bigit_count() > 0) - u.x |= m_mantissa[0]; - u.x <<= 20; - if (bigit_count() > 1) - { - u.x |= m_mantissa[1] >> 12; - /* Rounding */ - u.x += (m_mantissa[1] >> 11) & 1; - } - } - - return u.d; -} - -template<> real::operator long double() const -{ - double hi = double(*this); - double lo = double(*this - hi); - return (long double)(hi) + (long double)(lo); -} - -/* - * Create a real number from an ASCII representation - */ -template<> real::Real(char const *str) -{ - real ret = 0; - exponent_t exponent = 0; - bool hex = false, comma = false, nonzero = false, negative = false, finished = false; - - for (char const *p = str; *p && !finished; p++) - { - switch (*p) - { - case '-': - case '+': - if (p != str) - break; - negative = (*p == '-'); - break; - case '.': - if (comma) - finished = true; - comma = true; - break; - case 'x': - case 'X': - /* This character is only valid for 0x... and 0X... numbers */ - if (p != str + 1 || str[0] != '0') - finished = true; - hex = true; - break; - case 'p': - case 'P': - if (hex) - exponent += atoi(p + 1); - finished = true; - break; - case 'e': - case 'E': - if (!hex) - { - exponent += atoi(p + 1); - finished = true; - break; - } - LOL_ATTR_FALLTHROUGH - case 'a': case 'b': case 'c': case 'd': case 'f': - case 'A': case 'B': case 'C': case 'D': case 'F': - case '0': case '1': case '2': case '3': case '4': - case '5': case '6': case '7': case '8': case '9': - if (nonzero) - { - /* Multiply ret by 10 or 16 depending the base. */ - if (!hex) - { - real x = ret + ret; - ret = x + x + ret; - } - ret.m_exponent += hex ? 4 : 1; - } - if (*p != '0') - { - ret += (*p >= 'a' && *p <= 'f') ? (int)(*p - 'a' + 10) - : (*p >= 'A' && *p <= 'F') ? (int)(*p - 'A' + 10) - : (int)(*p - '0'); - nonzero = true; - } - if (comma) - exponent -= hex ? 4 : 1; - break; - default: - finished = true; - break; - } - } - - if (hex) - ret.m_exponent += exponent; - else if (exponent) - ret *= pow(R_10(), (real)exponent); - - if (negative) - ret = -ret; - - *this = ret; -} - -template<> real real::operator +() const -{ - return *this; -} - -template<> real real::operator -() const -{ - real ret = *this; - ret.m_sign ^= true; - return ret; -} - -template<> real real::operator +(real const &x) const -{ - if (x.is_zero()) - return *this; - - if (is_zero()) - return x; - - /* Ensure both arguments are positive. Otherwise, switch signs, - * or replace + with -. */ - if (is_negative()) - return -(-*this + -x); - - if (x.is_negative()) - return *this - (-x); - - /* Ensure *this has the larger exponent (no need for the mantissa to - * be larger, as in subtraction). Otherwise, switch. */ - if (m_exponent < x.m_exponent) - return x + *this; - - int64_t e1 = m_exponent; - int64_t e2 = x.m_exponent; - - int64_t bigoff = (e1 - e2) / bigit_bits(); - int64_t off = e1 - e2 - bigoff * bigit_bits(); - - /* FIXME: ensure we have the same number of bigits */ - if (bigoff > bigit_count()) - return *this; - - real ret; - ret.m_mantissa.resize(bigit_count()); - ret.m_exponent = m_exponent; - - uint64_t carry = 0; - for (int i = bigit_count(); i--; ) - { - carry += m_mantissa[i]; - if (i - bigoff >= 0) - carry += x.m_mantissa[i - bigoff] >> off; - - if (off && i - bigoff > 0) - carry += (x.m_mantissa[i - bigoff - 1] << (bigit_bits() - off)) & 0xffffffffu; - else if (i - bigoff == 0) - carry += (uint64_t)1 << (bigit_bits() - off); - - ret.m_mantissa[i] = (uint32_t)carry; - carry >>= bigit_bits(); - } - - /* Renormalise in case we overflowed the mantissa */ - if (carry) - { - carry--; - for (int i = 0; i < bigit_count(); ++i) - { - uint32_t tmp = ret.m_mantissa[i]; - ret.m_mantissa[i] = ((uint32_t)carry << (bigit_bits() - 1)) - | (tmp >> 1); - carry = tmp & 1u; - } - ret.m_exponent++; - } - - return ret; -} - -template<> real real::operator -(real const &x) const -{ - if (x.is_zero()) - return *this; - - if (is_zero()) - return -x; - - /* Ensure both arguments are positive. Otherwise, switch signs, - * or replace - with +. */ - if (is_negative()) - return -(-*this + x); - - if (x.is_negative()) - return (*this) + (-x); - - /* Ensure *this is larger than x */ - if (*this < x) - return -(x - *this); - - exponent_t e1 = m_exponent; - exponent_t e2 = x.m_exponent; - - exponent_t bigoff = (e1 - e2) / bigit_bits(); - exponent_t off = e1 - e2 - bigoff * bigit_bits(); - - /* FIXME: ensure we have the same number of bigits */ - if (bigoff > bigit_count()) - return *this; - - real ret; - ret.m_mantissa.resize(bigit_count()); - ret.m_exponent = m_exponent; - - /* int64_t instead of uint64_t to preserve sign through shifts */ - exponent_t carry = 0; - for (int i = 0; i < bigoff; ++i) - { - carry -= x.m_mantissa[bigit_count() - 1 - i]; - /* Emulates a signed shift */ - carry >>= bigit_bits(); - carry |= carry << bigit_bits(); - } - if (bigoff < bigit_count()) - carry -= x.m_mantissa[bigit_count() - 1 - bigoff] & (((exponent_t)1 << off) - 1); - carry /= (exponent_t)1 << off; - - for (int i = bigit_count(); i--; ) - { - carry += m_mantissa[i]; - if (i - bigoff >= 0) - carry -= x.m_mantissa[i - bigoff] >> off; - - if (off && i - bigoff > 0) - carry -= (x.m_mantissa[i - bigoff - 1] << (bigit_bits() - off)) & 0xffffffffu; - else if (i - bigoff == 0) - carry -= (uint64_t)1 << (bigit_bits() - off); - - ret.m_mantissa[i] = (bigit_t)carry; - carry >>= bigit_bits(); - carry |= carry << bigit_bits(); - } - - carry += 1; - - /* Renormalise if we underflowed the mantissa */ - if (carry == 0) - { - /* How much do we need to shift the mantissa? FIXME: this could - * be computed above */ - off = 0; - for (int i = 0; i < bigit_count(); ++i) - { - if (!ret.m_mantissa[i]) - { - off += bigit_bits(); - continue; - } - - /* “~tmp > tmp” checks that the MSB is not set */ - for (bigit_t tmp = ret.m_mantissa[i]; ~tmp > tmp; tmp <<= 1) - off++; - break; - } - if (off == total_bits()) - ret.m_mantissa.resize(0); - else - { - off++; /* Shift once more to get rid of the leading 1 */ - ret.m_exponent -= off; - - bigoff = off / bigit_bits(); - off -= bigoff * bigit_bits(); - - for (int i = 0; i < bigit_count(); ++i) - { - bigit_t tmp = 0; - if (i + bigoff < bigit_count()) - tmp |= ret.m_mantissa[i + bigoff] << off; - if (off && i + bigoff + 1 < bigit_count()) - tmp |= ret.m_mantissa[i + bigoff + 1] >> (bigit_bits() - off); - ret.m_mantissa[i] = tmp; - } - } - } - - return ret; -} - -template<> real real::operator *(real const &x) const -{ - real ret; - - /* The sign is easy to compute */ - ret.m_sign = is_negative() ^ x.is_negative(); - - /* If any operand is zero, return zero. FIXME: 0 * Inf? */ - if (is_zero() || x.is_zero()) - return ret; - - ret.m_mantissa.resize(bigit_count()); - ret.m_exponent = m_exponent + x.m_exponent; - - /* Accumulate low order product; no need to store it, we just - * want the carry value */ - uint64_t carry = 0, hicarry = 0, prev; - for (int i = 0; i < bigit_count(); ++i) - { - for (int j = 0; j < i + 1; j++) - { - prev = carry; - carry += (uint64_t)m_mantissa[bigit_count() - 1 - j] - * (uint64_t)x.m_mantissa[bigit_count() - 1 + j - i]; - if (carry < prev) - hicarry++; - } - carry >>= bigit_bits(); - carry |= hicarry << bigit_bits(); - hicarry >>= bigit_bits(); - } - - /* Multiply the other components */ - for (int i = 0; i < bigit_count(); ++i) - { - for (int j = i + 1; j < bigit_count(); j++) - { - prev = carry; - carry += (uint64_t)m_mantissa[bigit_count() - 1 - j] - * (uint64_t)x.m_mantissa[j - 1 - i]; - if (carry < prev) - hicarry++; - } - prev = carry; - carry += m_mantissa[bigit_count() - 1 - i]; - carry += x.m_mantissa[bigit_count() - 1 - i]; - if (carry < prev) - hicarry++; - ret.m_mantissa[bigit_count() - 1 - i] = carry & ~(bigit_t)0; - carry >>= bigit_bits(); - carry |= hicarry << bigit_bits(); - hicarry >>= bigit_bits(); - } - - /* Renormalise in case we overflowed the mantissa */ - if (carry) - { - carry--; - for (int i = 0; i < bigit_count(); ++i) - { - bigit_t tmp = ret.m_mantissa[i]; - ret.m_mantissa[i] = ((bigit_t)carry << (bigit_bits() - 1)) - | (tmp >> 1); - carry = tmp & 1u; - } - ++ret.m_exponent; - } - - return ret; -} - -template<> real real::operator /(real const &x) const -{ - return *this * inverse(x); -} - -template<> real const &real::operator +=(real const &x) -{ - real tmp = *this; - return *this = tmp + x; -} - -template<> real const &real::operator -=(real const &x) -{ - real tmp = *this; - return *this = tmp - x; -} - -template<> real const &real::operator *=(real const &x) -{ - real tmp = *this; - return *this = tmp * x; -} - -template<> real const &real::operator /=(real const &x) -{ - real tmp = *this; - return *this = tmp / x; -} - -template<> bool real::operator ==(real const &x) const -{ - /* If NaN is involved, return false */ - if (is_nan() || x.is_nan()) - return false; - - /* If both zero, they are equal; if either is zero, they are different */ - if (is_zero() || x.is_zero()) - return is_zero() && x.is_zero(); - - /* FIXME: handle NaN/Inf */ - return m_exponent == x.m_exponent && m_mantissa == x.m_mantissa; -} - -template<> bool real::operator !=(real const &x) const -{ - return !(is_nan() || x.is_nan() || *this == x); -} - -template<> bool real::operator <(real const &x) const -{ - /* If NaN is involved, return false */ - if (is_nan() || x.is_nan()) - return false; - - /* Ensure we are positive */ - if (is_negative()) - return -*this > -x; - - /* If x is zero or negative, we can’t be < x */ - if (x.is_negative() || x.is_zero()) - return false; - - /* If we are zero, we must be < x */ - if (is_zero()) - return true; - - /* Compare exponents */ - if (m_exponent != x.m_exponent) - return m_exponent < x.m_exponent; - - /* Compare all relevant bits */ - for (int i = 0; i < bigit_count(); ++i) - if (m_mantissa[i] != x.m_mantissa[i]) - return m_mantissa[i] < x.m_mantissa[i]; - - return false; -} - -template<> bool real::operator <=(real const &x) const -{ - return !(is_nan() || x.is_nan() || *this > x); -} - -template<> bool real::operator >(real const &x) const -{ - /* If NaN is involved, return false */ - if (is_nan() || x.is_nan()) - return false; - - /* Ensure we are positive */ - if (is_negative()) - return -*this < -x; - - /* If x is zero, we’re > x iff we’re non-zero since we’re positive */ - if (x.is_zero()) - return !is_zero(); - - /* If x is strictly negative, we’re > x */ - if (x.is_negative()) - return true; - - /* If we are zero, we can’t be > x */ - if (is_zero()) - return false; - - /* Compare exponents */ - if (m_exponent != x.m_exponent) - return m_exponent > x.m_exponent; - - /* Compare all relevant bits */ - for (int i = 0; i < bigit_count(); ++i) - if (m_mantissa[i] != x.m_mantissa[i]) - return m_mantissa[i] > x.m_mantissa[i]; - - return false; -} - -template<> bool real::operator >=(real const &x) const -{ - return !(is_nan() || x.is_nan() || *this < x); -} - -template<> bool real::operator !() const -{ - return !(bool)*this; -} - -template<> real::operator bool() const -{ - /* A real is "true" if it is non-zero AND not NaN */ - return !is_zero() && !is_nan(); -} - -template<> real min(real const &a, real const &b) -{ - return (a < b) ? a : b; -} - -template<> real max(real const &a, real const &b) -{ - return (a > b) ? a : b; -} - -template<> real clamp(real const &x, real const &a, real const &b) -{ - return (x < a) ? a : (x > b) ? b : x; -} - -template<> real inverse(real const &x) -{ - real ret; - - /* If zero, return infinite */ - if (x.is_zero()) - return copysign(real::R_INF(), x); - - /* Use the system’s float inversion to approximate 1/x */ - union { float f; uint32_t x; } u = { 1.0f }; - u.x |= x.m_mantissa[0] >> 9; - u.f = 1.0f / u.f; - - ret.m_mantissa.resize(x.bigit_count()); - ret.m_mantissa[0] = u.x << 9; - ret.m_sign = x.m_sign; - ret.m_exponent = -x.m_exponent + (u.x >> 23) - 0x7f; - - /* FIXME: 1+log2(bigit_count) steps of Newton-Raphson seems to be enough for - * convergence, but this hasn’t been checked seriously. */ - for (int i = 1; i <= x.bigit_count(); i *= 2) - ret = ret * (real::R_2() - ret * x); - - return ret; -} - -template<> real sqrt(real const &x) -{ - /* if zero, return x (FIXME: negative zero?) */ - if (x.is_zero()) - return x; - - /* if negative, return NaN */ - if (x.is_negative()) - return real::R_NAN(); - - int tweak = x.m_exponent & 1; - - /* Use the system’s float inversion to approximate 1/sqrt(x). First - * we construct a float in the [1..4[ range that has roughly the same - * mantissa as our real. Its exponent is 0 or 1, depending on the - * parity of x’s exponent. The final exponent is 0, -1 or -2. We use - * the final exponent and final mantissa to pre-fill the result. */ - union { float f; uint32_t x; } u = { 1.0f }; - u.x += tweak << 23; - u.x |= x.m_mantissa[0] >> 9; - u.f = 1.0f / sqrtf(u.f); - - real ret; - ret.m_mantissa.resize(x.bigit_count()); - ret.m_mantissa[0] = u.x << 9; - - ret.m_exponent = -(x.m_exponent - tweak) / 2 + (u.x >> 23) - 0x7f; - - /* FIXME: 1+log2(bigit_count()) steps of Newton-Raphson seems to be enough for - * convergence, but this hasn’t been checked seriously. */ - for (int i = 1; i <= x.bigit_count(); i *= 2) - { - ret = ret * (real::R_3() - ret * ret * x); - --ret.m_exponent; - } - - return ret * x; -} - -template<> real cbrt(real const &x) -{ - /* if zero, return x */ - if (x.is_zero()) - return x; - - int tweak = x.m_exponent % 3; - if (tweak < 0) - tweak += 3; - - /* Use the system’s float inversion to approximate cbrt(x). First - * we construct a float in the [1..8[ range that has roughly the same - * mantissa as our real. Its exponent is 0, 1 or 2, depending on the - * value of x. The final exponent is 0 or 1 (special case). We use - * the final exponent and final mantissa to pre-fill the result. */ - union { float f; uint32_t x; } u = { 1.0f }; - u.x += tweak << 23; - u.x |= x.m_mantissa[0] >> 9; - u.f = powf(u.f, 1.f / 3); - - real ret; - ret.m_mantissa.resize(x.bigit_count()); - ret.m_mantissa[0] = u.x << 9; - ret.m_exponent = (x.m_exponent - tweak) / 3 + (u.x >> 23) - 0x7f; - ret.m_sign = x.m_sign; - - /* FIXME: 1+log2(bigit_count()) steps of Newton-Raphson seems to be enough - * for convergence, but this hasn’t been checked seriously. */ - real third = inverse(real::R_3()); - for (int i = 1; i <= x.bigit_count(); i *= 2) - { - ret = third * (x / (ret * ret) + (ret * 2)); - } - - return ret; -} - -template<> real pow(real const &x, real const &y) -{ - /* Shortcuts for degenerate cases */ - if (!y) - return real::R_1(); - if (!x) - return real::R_0(); - - /* Small integer exponent: use exponentiation by squaring */ - int64_t int_y = (int64_t)y; - if (y == (real)int_y) - { - real ret = real::R_1(); - real x_n = int_y > 0 ? x : inverse(x); - - while (int_y) /* Can be > 0 or < 0 */ - { - if (int_y & 1) - ret *= x_n; - x_n *= x_n; - int_y /= 2; - } - - return ret; - } - - /* If x is positive, nothing special to do. */ - if (x > real::R_0()) - return exp(y * log(x)); - - /* XXX: manpage for pow() says “If x is a finite value less than 0, - * and y is a finite noninteger, a domain error occurs, and a NaN is - * returned”. We check whether y is closer to an even number or to - * an odd number and return something reasonable. */ - real round_y = round(y); - bool is_odd = round_y / 2 == round(round_y / 2); - return is_odd ? exp(y * log(-x)) : -exp(y * log(-x)); -} - -/* A fast factorial implementation for small numbers. An optional - * step argument allows to compute double factorials (i.e. with - * only the odd or the even terms. */ -static real fast_fact(int x, int step = 1) -{ - if (x < step) - return 1; - - if (x == step) - return x; - - unsigned int start = (x + step - 1) % step + 1; - real ret(start); - uint64_t multiplier = 1; - - for (int i = start, exponent = 0;;) - { - if (i >= x) - return ldexp(ret * multiplier, exponent); - - i += step; - - /* Accumulate the power of two part in the exponent */ - unsigned int tmp = i; - while ((tmp & 1) == 0) - { - tmp >>= 1; - exponent++; - } - - /* Accumulate the other factors in the multiplier */ - if (multiplier * tmp / tmp != multiplier) - { - ret *= multiplier; - multiplier = 1; - } - multiplier *= tmp; - } -} - -template<> real gamma(real const &x) -{ - /* We use Spouge’s formula. FIXME: precision is far from acceptable, - * especially with large values. We need to compute this with higher - * precision values in order to attain the desired accuracy. It might - * also be useful to sort the ck values by decreasing absolute value - * and do the addition in this order. */ - int a = (int)ceilf(logf(2) / logf(2 * F_PI) * x.total_bits()); - - real ret = sqrt(real::R_PI() * 2); - real fact_k_1 = real::R_1(); - - for (int k = 1; k < a; k++) - { - real a_k = (real)(a - k); - real ck = pow(a_k, (real)((float)k - 0.5)) * exp(a_k) - / (fact_k_1 * (x + (real)(k - 1))); - ret += ck; - fact_k_1 *= (real)-k; - } - - ret *= pow(x + (real)(a - 1), x - (real::R_1() / 2)); - ret *= exp(-x - (real)(a - 1)); - - return ret; -} - -template<> real fabs(real const &x) -{ - real ret = x; - ret.m_sign = false; - return ret; -} - -template<> real abs(real const &x) -{ - return fabs(x); -} - -template<> real fract(real const &x) -{ - return x - floor(x); -} - -template<> real degrees(real const &x) -{ - /* FIXME: need to recompute this for different mantissa sizes */ - static real mul = real(180) * real::R_1_PI(); - - return x * mul; -} - -template<> real radians(real const &x) -{ - /* FIXME: need to recompute this for different mantissa sizes */ - static real mul = real::R_PI() / real(180); - - return x * mul; -} - -static real fast_log(real const &x) -{ - /* This fast log method is tuned to work on the [1..2] range and - * no effort whatsoever was made to improve convergence outside this - * domain of validity. It can converge pretty fast, provided we use - * the following variable substitutions: - * y = sqrt(x) - * z = (y - 1) / (y + 1) - * - * And the following identities: - * ln(x) = 2 ln(y) - * = 2 ln((1 + z) / (1 - z)) - * = 4 z (1 + z^2 / 3 + z^4 / 5 + z^6 / 7...) - * - * Any additional sqrt() call would halve the convergence time, but - * would also impact the final precision. For now we stick with one - * sqrt() call. */ - real y = sqrt(x); - real z = (y - real::R_1()) / (y + real::R_1()), z2 = z * z, zn = z2; - real sum = real::R_1(); - - for (int i = 3; ; i += 2) - { - real newsum = sum + zn / (real)i; - if (newsum == sum) - break; - sum = newsum; - zn *= z2; - } - - return z * sum * 4; -} - -template<> real log(real const &x) -{ - /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M), - * with the property that M is in [1..2[, so fast_log() applies here. */ - if (x.is_negative() || x.is_zero()) - return real::R_NAN(); - - real tmp(x); - tmp.m_exponent = 0; - return real(x.m_exponent) * real::R_LN2() + fast_log(tmp); -} - -template<> real log2(real const &x) -{ - /* Strategy for log2(x): see log(x). */ - if (x.is_negative() || x.is_zero()) - return real::R_NAN(); - - real tmp(x); - tmp.m_exponent = 0; - return real(x.m_exponent) + fast_log(tmp) * real::R_LOG2E(); -} - -template<> real log10(real const &x) -{ - return log(x) * real::R_LOG10E(); -} - -static real fast_exp_sub(real const &x, real const &y) -{ - /* This fast exp method is tuned to work on the [-1..1] range and - * no effort whatsoever was made to improve convergence outside this - * domain of validity. The argument y is used for cases where we - * don’t want the leading 1 in the Taylor series. */ - real ret = real::R_1() - y, xn = x; - int i = 1; - - for (;;) - { - real newret = ret + xn; - if (newret == ret) - break; - ret = newret * ++i; - xn *= x; - } - - return ret / fast_fact(i); -} - -template<> real exp(real const &x) -{ - /* Strategy for exp(x): the Taylor series does not converge very fast - * with large positive or negative values. - * - * However, since the result is going to be in the form M*2^E, we first - * try to predict a value for E, which is approximately: - * E ≈ log2(exp(x)) = x / log(2) - * - * Let E0 be an integer close to x / log(2). We need to find a value x0 - * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - E0 log(2). - * - * Thus the final algorithm: - * int E0 = x / log(2) - * real x0 = x - E0 log(2) - * real x1 = exp(x0) - * return x1 * 2^E0 - */ - real::exponent_t e0 = x / real::R_LN2(); - real x0 = x - (real)e0 * real::R_LN2(); - real x1 = fast_exp_sub(x0, real::R_0()); - x1.m_exponent += e0; - return x1; -} - -template<> real exp2(real const &x) -{ - /* Strategy for exp2(x): see strategy in exp(). */ - real::exponent_t e0 = x; - real x0 = x - (real)e0; - real x1 = fast_exp_sub(x0 * real::R_LN2(), real::R_0()); - x1.m_exponent += e0; - return x1; -} - -template<> real erf(real const &x) -{ - /* Strategy for erf(x): - * - if x<0, erf(x) = -erf(-x) - * - if x<7, erf(x) = sum((-1)^n·x^(2n+1)/((2n+1)·n!))/sqrt(π/4) - * - if x≥7, erf(x) = 1+exp(-x²)/(x·sqrt(π))·sum((-1)^n·(2n-1)!!/(2x²)^n - * - * FIXME: do not compute factorials at each iteration, accumulate - * them instead (see fast_exp_sub). - * FIXME: For a potentially faster implementation, see “Expanding the - * Error Function erf(z)” at: - * http://www.mathematica-journal.com/2014/11/on-burmanns-theorem-and-its-application-to-problems-of-linear-and-nonlinear-heat-transfer-and-diffusion/#more-39602/ - */ - if (x.is_negative()) - return -erf(-x); - - real sum = real::R_0(); - real x2 = x * x; - - /* FIXME: this test is inefficient; the series converges slowly for x≥1 */ - if (x < real(7)) - { - real xn = x, xmul = x2; - for (int n = 0;; ++n, xn *= xmul) - { - real tmp = xn / (fast_fact(n) * (2 * n + 1)); - real newsum = (n & 1) ? sum - tmp : sum + tmp; - if (newsum == sum) - break; - sum = newsum; - } - return sum * real::R_2_SQRTPI(); - } - else - { - real xn = real::R_1(), xmul = inverse(x2 + x2); - /* FIXME: this does not converge well! We need to stop at 30 - * iterations and sacrifice some accuracy. */ - for (int n = 0; n < 30; ++n, xn *= xmul) - { - real tmp = xn * fast_fact(n * 2 - 1, 2); - real newsum = (n & 1) ? sum - tmp : sum + tmp; - if (newsum == sum) - break; - sum = newsum; - } - - return real::R_1() - exp(-x2) / (x * sqrt(real::R_PI())) * sum; - } -} - -template<> real sinh(real const &x) -{ - /* We cannot always use (exp(x)-exp(-x))/2 because we'll lose - * accuracy near zero. We only use this identity for |x|>0.5. If - * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */ - bool near_zero = (fabs(x) < real::R_1() / 2); - real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x); - real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x); - return (x1 - x2) / 2; -} - -template<> real tanh(real const &x) -{ - /* See sinh() for the strategy here */ - bool near_zero = (fabs(x) < real::R_1() / 2); - real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x); - real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x); - real x3 = near_zero ? x1 + x2 + real::R_2() : x1 + x2; - return (x1 - x2) / x3; -} - -template<> real cosh(real const &x) -{ - /* No need to worry about accuracy here; maybe the last bit is slightly - * off, but that's about it. */ - return (exp(x) + exp(-x)) / 2; -} - -template<> real frexp(real const &x, real::exponent_t *exp) -{ - if (!x) - { - *exp = 0; - return x; - } - - /* FIXME: check that this works */ - *exp = x.m_exponent; - - real ret = x; - ret.m_exponent = 0; - return ret; -} - -template<> real ldexp(real const &x, real::exponent_t exp) -{ - real ret = x; - if (ret) /* Only do something if non-zero */ - ret.m_exponent += exp; - return ret; -} - -template<> real modf(real const &x, real *iptr) -{ - real absx = fabs(x); - real tmp = floor(absx); - - *iptr = copysign(tmp, x); - return copysign(absx - tmp, x); -} - -template<> real nextafter(real const &x, real const &y) -{ - /* Linux manpage: “If x equals y, the functions return y.” */ - if (x == y) - return y; - - /* Ensure x is positive. */ - if (x.is_negative()) - return -nextafter(-x, -y); - - /* FIXME: broken for now */ - real ulp = ldexp(x, -x.total_bits()); - return x < y ? x + ulp : x - ulp; -} - -template<> real copysign(real const &x, real const &y) -{ - real ret = x; - ret.m_sign = y.m_sign; - return ret; -} - -template<> real floor(real const &x) -{ - /* Strategy for floor(x): - * - if negative, return -ceil(-x) - * - if zero or negative zero, return x - * - if less than one, return zero - * - otherwise, if e is the exponent, clear all bits except the - * first e. */ - if (x < -real::R_0()) - return -ceil(-x); - if (!x) - return x; - if (x < real::R_1()) - return real::R_0(); - - real ret = x; - real::exponent_t exponent = x.m_exponent; - - for (int i = 0; i < x.bigit_count(); ++i) - { - if (exponent <= 0) - ret.m_mantissa[i] = 0; - else if (exponent < real::bigit_bits()) - ret.m_mantissa[i] &= ~((1 << (real::bigit_bits() - exponent)) - 1); - - exponent -= real::bigit_bits(); - } - - return ret; -} - -template<> real ceil(real const &x) -{ - /* Strategy for ceil(x): - * - if negative, return -floor(-x) - * - if x == floor(x), return x - * - otherwise, return floor(x) + 1 */ - if (x < -real::R_0()) - return -floor(-x); - real ret = floor(x); - if (ret < x) - ret += real::R_1(); - return ret; -} - -template<> real round(real const &x) -{ - if (x < real::R_0()) - return -round(-x); - - return floor(x + (real::R_1() / 2)); -} - -template<> real fmod(real const &x, real const &y) -{ - if (!y) - return real::R_0(); /* FIXME: return NaN */ - - if (!x) - return x; - - real tmp = round(x / y); - return x - tmp * y; -} - -template<> real sin(real const &x) -{ - bool switch_sign = x.is_negative(); - - real absx = fmod(fabs(x), real::R_PI() * 2); - if (absx > real::R_PI()) - { - absx -= real::R_PI(); - switch_sign = !switch_sign; - } - - if (absx > real::R_PI_2()) - absx = real::R_PI() - absx; - - real ret = real::R_0(), fact = real::R_1(), xn = absx, mx2 = -absx * absx; - int i = 1; - for (;;) - { - real newret = ret + xn; - if (newret == ret) - break; - ret = newret * ((i + 1) * (i + 2)); - xn *= mx2; - i += 2; - } - ret /= fast_fact(i); - - /* Propagate sign */ - ret.m_sign ^= switch_sign; - return ret; -} - -template<> real cos(real const &x) -{ - return sin(real::R_PI_2() - x); -} - -template<> real tan(real const &x) -{ - /* Constrain input to [-π,π] */ - real y = fmod(x, real::R_PI()); - - /* Constrain input to [-π/2,π/2] */ - if (y < -real::R_PI_2()) - y += real::R_PI(); - else if (y > real::R_PI_2()) - y -= real::R_PI(); - - /* In [-π/4,π/4] return sin/cos */ - if (fabs(y) <= real::R_PI_4()) - return sin(y) / cos(y); - - /* Otherwise, return cos/sin */ - if (y > real::R_0()) - y = real::R_PI_2() - y; - else - y = -real::R_PI_2() - y; - - return cos(y) / sin(y); -} - -static inline real asinacos(real const &x, int is_asin) -{ - /* Strategy for asin(): in [-0.5..0.5], use a Taylor series around - * zero. In [0.5..1], use asin(x) = π/2 - 2*asin(sqrt((1-x)/2)), and - * in [-1..-0.5] just revert the sign. - * Strategy for acos(): use acos(x) = π/2 - asin(x) and try not to - * lose the precision around x=1. */ - real absx = fabs(x); - int around_zero = (absx < (real::R_1() / 2)); - - if (!around_zero) - absx = sqrt((real::R_1() - absx) / 2); - - real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1; - for (int i = 1; ; ++i) - { - xn *= x2; - real mul = (real)(2 * i + 1); - real newret = ret + ldexp(fact1 * xn / (mul * fact2), -2 * i); - if (newret == ret) - break; - ret = newret; - fact1 *= (real)((2 * i + 1) * (2 * i + 2)); - fact2 *= (real)((i + 1) * (i + 1)); - } - - if (x.is_negative()) - ret = -ret; - - if (around_zero) - ret = is_asin ? ret : real::R_PI_2() - ret; - else - { - real adjust = x.is_negative() ? real::R_PI() : real::R_0(); - if (is_asin) - ret = real::R_PI_2() - adjust - ret * 2; - else - ret = adjust + ret * 2; - } - - return ret; -} - -template<> real asin(real const &x) -{ - return asinacos(x, 1); -} - -template<> real acos(real const &x) -{ - return asinacos(x, 0); -} - -template<> real atan(real const &x) -{ - /* Computing atan(x): we choose a different Taylor series depending on - * the value of x to help with convergence. - * - * If |x| < 0.5 we evaluate atan(y) near 0: - * atan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 ... - * - * If 0.5 <= |x| < 1.5 we evaluate atan(1+y) near 0: - * atan(1+y) = π/4 + y/(1*2^1) - y^2/(2*2^1) + y^3/(3*2^2) - * - y^5/(5*2^3) + y^6/(6*2^3) - y^7/(7*2^4) - * + y^9/(9*2^5) - y^10/(10*2^5) + y^11/(11*2^6) ... - * - * If 1.5 <= |x| < 2 we evaluate atan(sqrt(3)+2y) near 0: - * atan(sqrt(3)+2y) = π/3 + 1/2 y - sqrt(3)/2 y^2/2 - * + y^3/3 - sqrt(3)/2 y^4/4 + 1/2 y^5/5 - * - 1/2 y^7/7 + sqrt(3)/2 y^8/8 - * - y^9/9 + sqrt(3)/2 y^10/10 - 1/2 y^11/11 - * + 1/2 y^13/13 - sqrt(3)/2 y^14/14 - * + y^15/15 - sqrt(3)/2 y^16/16 + 1/2 y^17/17 ... - * - * If |x| >= 2 we evaluate atan(y) near +∞: - * atan(y) = π/2 - y^-1 + y^-3/3 - y^-5/5 + y^-7/7 - y^-9/9 ... - */ - real absx = fabs(x); - - if (absx < (real::R_1() / 2)) - { - real ret = x, xn = x, mx2 = -x * x; - for (int i = 3; ; i += 2) - { - xn *= mx2; - real newret = ret + xn / (real)i; - if (newret == ret) - break; - ret = newret; - } - return ret; - } - - real ret = 0; - - if (absx < (real::R_3() / 2)) - { - real y = real::R_1() - absx; - real yn = y, my2 = -y * y; - for (int i = 0; ; i += 2) - { - real newret = ret + ldexp(yn / (real)(2 * i + 1), -i - 1); - yn *= y; - newret += ldexp(yn / (real)(2 * i + 2), -i - 1); - yn *= y; - newret += ldexp(yn / (real)(2 * i + 3), -i - 2); - if (newret == ret) - break; - ret = newret; - yn *= my2; - } - ret = real::R_PI_4() - ret; - } - else if (absx < real::R_2()) - { - real y = (absx - real::R_SQRT3()) / 2; - real yn = y, my2 = -y * y; - for (int i = 1; ; i += 6) - { - real newret = ret + ((yn / (real)i) / 2); - yn *= y; - newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 1); - yn *= y; - newret += yn / (real)(i + 2); - yn *= y; - newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 3); - yn *= y; - newret += (yn / (real)(i + 4)) / 2; - if (newret == ret) - break; - ret = newret; - yn *= my2; - } - ret = real::R_PI_3() + ret; - } - else - { - real y = inverse(absx); - real yn = y, my2 = -y * y; - ret = y; - for (int i = 3; ; i += 2) - { - yn *= my2; - real newret = ret + yn / (real)i; - if (newret == ret) - break; - ret = newret; - } - ret = real::R_PI_2() - ret; - } - - /* Propagate sign */ - ret.m_sign = x.m_sign; - return ret; -} - -template<> real atan2(real const &y, real const &x) -{ - if (!y) - { - if (!x.is_negative()) - return y; - return y.is_negative() ? -real::R_PI() : real::R_PI(); - } - - if (!x) - { - return y.is_negative() ? -real::R_PI() : real::R_PI(); - } - - /* FIXME: handle the Inf and NaN cases */ - real z = y / x; - real ret = atan(z); - if (x < real::R_0()) - ret += (y > real::R_0()) ? real::R_PI() : -real::R_PI(); - return ret; -} - -/* Franke’s function, used as a test for interpolation methods */ -template<> real franke(real const &x, real const &y) -{ - /* Compute 9x and 9y */ - real nx = x + x; nx += nx; nx += nx + x; - real ny = y + y; ny += ny; ny += ny + y; - - /* Temporary variables for the formula */ - real a = nx - real::R_2(); - real b = ny - real::R_2(); - real c = nx + real::R_1(); - real d = ny + real::R_1(); - real e = nx - real(7); - real f = ny - real::R_3(); - real g = nx - real(4); - real h = ny - real(7); - - return exp(-(a * a + b * b) * real(0.25)) * real(0.75) - + exp(-(c * c / real(49) + d * d / real::R_10())) * real(0.75) - + exp(-(e * e + f * f) * real(0.25)) * real(0.5) - - exp(-(g * g + h * h)) / real(5); -} - -/* The Peaks example function from Matlab */ -template<> real peaks(real const &x, real const &y) -{ - real x2 = x * x; - real y2 = y * y; - /* 3 * (1-x)^2 * exp(-x^2 - (y+1)^2) */ - real ret = real::R_3() - * (x2 - x - x + real::R_1()) - * exp(- x2 - y2 - y - y - real::R_1()); - /* -10 * (x/5 - x^3 - y^5) * exp(-x^2 - y^2) */ - ret -= (x + x - real::R_10() * (x2 * x + y2 * y2 * y)) * exp(-x2 - y2); - /* -1/3 * exp(-(x+1)^2 - y^2) */ - ret -= exp(-x2 - x - x - real::R_1() - y2) / real::R_3(); - return ret; -} - -template<> -std::ostream& operator <<(std::ostream &s, real const &x) -{ - bool hex = (s.flags() & std::ios_base::basefield) == std::ios_base::hex; - s << (hex ? x.xstr() : x.str((int)s.precision())); - return s; -} - -template<> std::string real::str(int ndigits) const -{ - std::stringstream ss; - real x = *this; - - if (x.is_negative()) - { - ss << '-'; - x = -x; - } - - if (!x) - { - ss << '0'; - return ss.str(); - } - - // Normalise x so that mantissa is in [1..9.999] - // FIXME: better use int64_t when the cast is implemented - // FIXME: does not work with R_MAX and probably R_MIN - int exponent = ceil(log10(x)); - x *= pow(R_10(), -(real)exponent); - - if (ndigits < 1) - ndigits = 1; - - // Add a bias to simulate some naive rounding - x += real(4.99f) * pow(R_10(), -(real)(ndigits + 1)); - - if (x < R_1()) - { - x *= R_10(); - exponent--; - } - - // Print digits - for (int i = 0; i < ndigits; ++i) - { - int digit = (int)floor(x); - ss << (char)('0' + digit); - if (i == 0) - ss << '.'; - x -= real(digit); - x *= R_10(); - } - - // Remove trailing zeroes - std::string ret = ss.str(); - ss.str(""); - size_t end = ret.find_last_not_of('0'); - if (end != std::string::npos) - ss << ret.substr(0, end + 1); - - // Print exponent information - if (exponent) - ss << 'e' << (exponent >= 0 ? '+' : '-') << lol::abs(exponent); - - return ss.str(); -} - -template<> std::string real::xstr() const -{ - std::stringstream ss; - if (is_negative()) - ss << '-'; - ss << "0x1." << std::hex << std::setfill('0'); - for (int i = 0; i < bigit_count(); ++i) - ss << std::setw(8) << m_mantissa[i]; - ss << std::dec; - - // Remove trailing zeroes - std::string ret = ss.str(); - ss.str(""); - size_t end = ret.find_last_not_of('0'); - if (end != std::string::npos) - ss << ret.substr(0, end + 1); - - ss << 'p' << m_exponent; - - return ss.str(); -} - -static real load_min() -{ - real ret = 1; - return ldexp(ret, std::numeric_limits::min()); -} - -static real load_max() -{ - /* FIXME: the last bits of the mantissa are not properly handled in this - * code! So we fallback to a slow but exact method. */ -#if 0 - real ret = 1; - ret = ldexp(ret, real::TOTAL_BITS - 1) - ret; - return ldexp(ret, real::EXPONENT_BIAS + 2 - real::TOTAL_BITS); -#endif - /* Generates 0x1.ffff..ffffp18446744073709551615 */ - char str[160]; - std::sprintf(str, "0x1.%llx%llx%llx%llx%llx%llx%llx%llxp%lld", - -1ll, -1ll, -1ll, -1ll, -1ll, -1ll, -1ll, -1ll, - (long long int)std::numeric_limits::max()); - return real(str); -} - -static real load_pi() -{ - /* Approximate π using Machin’s formula: 16*atan(1/5)-4*atan(1/239) */ - real ret = 0, x0 = 5, x1 = 239; - real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16, r4 = 4; - - for (int i = 1; ; i += 2) - { - real newret = ret + r16 / (x0 * (real)i) - r4 / (x1 * (real)i); - if (newret == ret) - break; - ret = newret; - x0 *= m0; - x1 *= m1; - } - - return ret; -} - -} /* namespace lol */ - diff --git a/src/scene.h b/src/scene.h index 46b9fd1f..42d2526c 100644 --- a/src/scene.h +++ b/src/scene.h @@ -25,7 +25,9 @@ #include "light.h" #include "camera.h" #include "mesh/mesh.h" + #include +#include #define LOL_MAX_LIGHT_COUNT 8 diff --git a/src/sys/getopt.cpp b/src/sys/getopt.cpp deleted file mode 100644 index 96221ea0..00000000 --- a/src/sys/getopt.cpp +++ /dev/null @@ -1,198 +0,0 @@ -// -// Lol Engine -// -// Copyright © 2002—2018 Sam Hocevar -// -// Lol Engine is free software. It comes without any warranty, to -// the extent permitted by applicable law. You can redistribute it -// and/or modify it under the terms of the Do What the Fuck You Want -// to Public License, Version 2, as published by the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#include - -#if HAVE_GETOPT_H -# include -#endif -#if HAVE_UNISTD_H -# include -#endif - -namespace lol -{ - -struct getopt_private -{ - getopt_private(int argc, char * const * _argv) - : m_argc(argc), - m_argv(_argv) - {} - -#if HAVE_GETOPT_LONG - typedef option optdesc; -#else - struct optdesc - { - char const *name; - int has_arg; - int *flag; - int val; - }; -#endif - - int m_argc; - char * const *m_argv; - - std::string m_optstring; - array m_opts; -}; - -getopt::getopt(int argc, char ** _argv) - : index(1), - arg(nullptr), - m_private(new getopt_private(argc, _argv)) -{ -} - -getopt::getopt(int argc, char * const * _argv) - : index(1), - arg(nullptr), - m_private(new getopt_private(argc, _argv)) -{ -} - -getopt::~getopt() -{ -} - -void getopt::add_opt(int short_opt, char const *long_opt, bool has_arg) -{ - getopt_private::optdesc o { long_opt, has_arg, nullptr, short_opt }; - m_private->m_opts.push(o); - - /* “The standards require that the argument [to isalnum()] is either - * EOF or a value that is representable in the type unsigned char.” */ - if ((int)(unsigned char)short_opt == short_opt && isalnum(short_opt)) - { - m_private->m_optstring += (char)short_opt; - if (has_arg) - m_private->m_optstring += ':'; - } -} - -int getopt::parse() -{ - int longindex = 0; // FIXME: what is this? - -#if HAVE_GETOPT_LONG - int ret; - optind = this->index; - optarg = this->arg; - m_private->m_opts.push(getopt_private::optdesc { nullptr, 0, nullptr, 0 }); - ret = getopt_long(m_private->m_argc, m_private->m_argv, m_private->m_optstring.c_str(), - (option const *)m_private->m_opts.data(), &longindex); - this->index = optind; - this->arg = optarg; - return ret; - -#else - /* XXX: this getopt_long implementation should not be trusted for other - * applications without any serious peer reviewing. It “just works” with - * zzuf and a few libcaca programs but may fail miserably in other - * programs. */ - char **argv = (char **)(uintptr_t)m_private->m_argv; - char *flag; - - if (this->index >= m_private->m_argc) - return -1; - - flag = argv[this->index]; - - if (flag[0] == '-' && flag[1] != '-') - { - char const *tmp; - int ret = flag[1]; - - if (ret == '\0') - return -1; - - tmp = strchr(m_private->m_optstring.c_str(), ret); - if (!tmp || ret == ':') - return '?'; - - ++this->index; - if (tmp[1] == ':') - { - if (flag[2] != '\0') - this->arg = flag + 2; - else - { - if (this->index >= m_private->m_argc) - goto too_few; - this->arg = argv[this->index++]; - } - return ret; - } - - if (flag[2] != '\0') - { - flag[1] = '-'; - --this->index; - ++argv[this->index]; - } - - return ret; - } - - if (flag[0] == '-' && flag[1] == '-') - { - if (flag[2] == '\0') - return -1; - - for (int i = 0; m_private->m_opts[i].name; ++i) - { - size_t l = strlen(m_private->m_opts[i].name); - - if (strncmp(flag + 2, m_private->m_opts[i].name, l)) - continue; - - switch (flag[2 + l]) - { - case '=': - if (!m_private->m_opts[i].has_arg) - goto bad_opt; - longindex = i; - ++this->index; - this->arg = flag + 2 + l + 1; - return m_private->m_opts[i].val; - case '\0': - longindex = i; - ++this->index; - if (m_private->m_opts[i].has_arg) - { - if (this->index >= m_private->m_argc) - goto too_few; - this->arg = argv[this->index++]; - } - return m_private->m_opts[i].val; - default: - break; - } - } - bad_opt: - fprintf(stderr, "%s: unrecognized option `%s'\n", argv[0], flag); - return '?'; - - too_few: - fprintf(stderr, "%s: option `%s' requires an argument\n", - argv[0], flag); - return '?'; - } - - return -1; -#endif -} - -} // namespace lol - diff --git a/src/t/base/string.cpp b/src/t/base/string.cpp index 6d8b279a..341b78d3 100644 --- a/src/t/base/string.cpp +++ b/src/t/base/string.cpp @@ -51,27 +51,27 @@ lolunit_declare_fixture(string_test) lolunit_declare_test(string_split) { auto l1 = split(std::string("abc")); - lolunit_assert(l1.count() == 1); + lolunit_assert(l1.size() == 1); lolunit_assert(l1[0] == "abc"); auto l2 = split(std::string("\nabc")); - lolunit_assert(l2.count() == 2); + lolunit_assert(l2.size() == 2); lolunit_assert(l2[0] == ""); lolunit_assert(l2[1] == "abc"); auto l3 = split(std::string("abc\n")); - lolunit_assert(l3.count() == 2); + lolunit_assert(l3.size() == 2); lolunit_assert(l3[0] == "abc"); lolunit_assert(l3[1] == ""); auto l4 = split(std::string("\n\n")); - lolunit_assert(l4.count() == 3); + lolunit_assert(l4.size() == 3); lolunit_assert(l4[0] == ""); lolunit_assert(l4[1] == ""); lolunit_assert(l4[2] == ""); auto l5 = split(std::string("abc\nde\n\nf\n")); - lolunit_assert(l5.count() == 5); + lolunit_assert(l5.size() == 5); lolunit_assert(l5[0] == "abc"); lolunit_assert(l5[1] == "de"); lolunit_assert(l5[2] == ""); diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp index 3ec3d661..0416263a 100644 --- a/src/t/math/polynomial.cpp +++ b/src/t/math/polynomial.cpp @@ -196,12 +196,12 @@ lolunit_declare_fixture(polynomial_test) * r(x) = 3 + x + x² * s(x) = 5 */ auto r = p / q; - lolunit_assert_equal(r.m1.degree(), 2); - lolunit_assert_doubles_equal(r.m1[0], 3.f, 1e-5f); - lolunit_assert_doubles_equal(r.m1[1], 1.f, 1e-5f); - lolunit_assert_doubles_equal(r.m1[2], 1.f, 1e-5f); - lolunit_assert_equal(r.m2.degree(), 0); - lolunit_assert_doubles_equal(r.m2[0], 5.f, 1e-5f); + lolunit_assert_equal(std::get<0>(r).degree(), 2); + lolunit_assert_doubles_equal(std::get<0>(r)[0], 3.f, 1e-5f); + lolunit_assert_doubles_equal(std::get<0>(r)[1], 1.f, 1e-5f); + lolunit_assert_doubles_equal(std::get<0>(r)[2], 1.f, 1e-5f); + lolunit_assert_equal(std::get<1>(r).degree(), 0); + lolunit_assert_doubles_equal(std::get<1>(r)[0], 5.f, 1e-5f); } lolunit_declare_test(composition_degree_2_2) @@ -241,7 +241,7 @@ lolunit_declare_fixture(polynomial_test) polynomial p { 42.f }; auto roots = p.roots(); - lolunit_assert_equal(roots.count(), 0); + lolunit_assert_equal(roots.size(), 0); } lolunit_declare_test(degree_1_root) @@ -250,7 +250,7 @@ lolunit_declare_fixture(polynomial_test) polynomial p { -6.f, 2.f }; auto roots = p.roots(); - lolunit_assert_equal(roots.count(), 1); + lolunit_assert_equal(roots.size(), 1); lolunit_assert_equal(roots[0], 3.f); } @@ -260,14 +260,14 @@ lolunit_declare_fixture(polynomial_test) polynomial p { 81.f, -18.f, 1.f }; auto roots1 = p.roots(); - lolunit_assert_equal(roots1.count(), 1); + lolunit_assert_equal(roots1.size(), 1); lolunit_assert_equal(roots1[0], 9.f); /* p(x) = 42 - 20x + 2x² */ polynomial q { 42.f, -20.f, 2.f }; auto roots2 = q.roots(); - lolunit_assert_equal(roots2.count(), 2); + lolunit_assert_equal(roots2.size(), 2); lolunit_assert_equal(roots2[0], 3.f); lolunit_assert_equal(roots2[1], 7.f); } @@ -277,7 +277,7 @@ lolunit_declare_fixture(polynomial_test) polynomial p { 1.f, 3.f, 3.f, 1.f }; auto roots1 = p.roots(); - lolunit_assert_equal(roots1.count(), 1); + lolunit_assert_equal(roots1.size(), 1); lolunit_assert_doubles_equal(roots1[0], -1, 0); } @@ -287,7 +287,7 @@ lolunit_declare_fixture(polynomial_test) auto roots1 = p.roots(); // Should have 2 solutions only, but precision leads to 3 solutions - lolunit_assert_equal(roots1.count(), 2); + lolunit_assert_equal(roots1.size(), 2); lolunit_assert_doubles_equal(roots1[0], -1, 1e-6); lolunit_assert_doubles_equal(roots1[1], -2, 1e-6); } @@ -297,7 +297,7 @@ lolunit_declare_fixture(polynomial_test) polynomial p { 6.f, 11.f, 6.f, 1.f }; auto roots1 = p.roots(); - lolunit_assert_equal(roots1.count(), 3); + lolunit_assert_equal(roots1.size(), 3); lolunit_assert_doubles_equal(roots1[0], -1, 1e-8); lolunit_assert_doubles_equal(roots1[1], -3, 1e-8); lolunit_assert_doubles_equal(roots1[2], -2, 1e-8); @@ -308,7 +308,7 @@ lolunit_declare_fixture(polynomial_test) polynomial p { -12000.f, 1000.f - 1200.f - 120.f, 100.f + 10.0f - 12.f, 1.f }; auto roots1 = p.roots(); - lolunit_assert_equal(roots1.count(), 3); + lolunit_assert_equal(roots1.size(), 3); lolunit_assert_doubles_equal(roots1[0], 12, 1e-5); lolunit_assert_doubles_equal(roots1[1], -100, 1e-5); diff --git a/src/t/math/quat.cpp b/src/t/math/quat.cpp index 0250a86f..acaaadd0 100644 --- a/src/t/math/quat.cpp +++ b/src/t/math/quat.cpp @@ -22,16 +22,16 @@ lolunit_declare_fixture(quaternion_test) void setup() { /* Generate identity quaternions */ - m_vectorpairs.push(vec3::axis_x, vec3::axis_x); - m_vectorpairs.push(2.f * vec3::axis_x, 3.f * vec3::axis_x); + push_vector_pair(vec3::axis_x, vec3::axis_x); + push_vector_pair(2.f * vec3::axis_x, 3.f * vec3::axis_x); /* Generate 90-degree rotations */ - m_vectorpairs.push(vec3::axis_x, vec3::axis_y); - m_vectorpairs.push(2.f * vec3::axis_x, 3.f * vec3::axis_y); + push_vector_pair(vec3::axis_x, vec3::axis_y); + push_vector_pair(2.f * vec3::axis_x, 3.f * vec3::axis_y); /* Generate 180-degree rotations */ - m_vectorpairs.push(vec3::axis_x, -vec3::axis_x); - m_vectorpairs.push(2.f * vec3::axis_x, -3.f * vec3::axis_x); + push_vector_pair(vec3::axis_x, -vec3::axis_x); + push_vector_pair(2.f * vec3::axis_x, -3.f * vec3::axis_x); /* Fill array with random test values */ for (int i = 0; i < 10000; ++i) @@ -40,7 +40,7 @@ lolunit_declare_fixture(quaternion_test) * vec3(rand(-1.f, 1.f), rand(-1.f, 1.f), rand(-1.f, 1.f)); vec3 v2 = lol::pow(10.f, rand(-5.f, 5.f)) * vec3(rand(-1.f, 1.f), rand(-1.f, 1.f), rand(-1.f, 1.f)); - m_vectorpairs.push(v1, v2); + push_vector_pair(v1, v2); } } @@ -203,41 +203,44 @@ lolunit_declare_fixture(quaternion_test) { for (auto pair : m_vectorpairs) { - vec3 a = pair.m1; - vec3 b = pair.m2; - vec3 da = normalize(a); - vec3 db = normalize(b); + vec3 a0 = std::get<0>(pair); + vec3 b0 = std::get<1>(pair); + vec3 a = normalize(a0); + vec3 b = normalize(b0); - quat q = quat::rotate(a, b); + quat q = quat::rotate(a0, b0); + + auto ctx = a0.tostring() + " " + b0.tostring(); + lolunit_set_context(ctx); /* Check that q is a unit quaternion */ lolunit_assert_doubles_equal(1.0, (double)norm(q), 1e-5); - /* Check that q transforms da into db */ - vec3 c = q.transform(da); + /* Check that q transforms a into b */ + vec3 c = q.transform(a); - lolunit_assert_doubles_equal(c.x, db.x, 1e-5); - lolunit_assert_doubles_equal(c.y, db.y, 1e-5); - lolunit_assert_doubles_equal(c.z, db.z, 1e-5); + lolunit_assert_doubles_equal(c.x, b.x, 2e-5); + lolunit_assert_doubles_equal(c.y, b.y, 2e-5); + lolunit_assert_doubles_equal(c.z, b.z, 2e-5); - /* Check that ~q transforms db into da */ - vec3 d = (~q).transform(db); + /* Check that ~q transforms b into a */ + vec3 d = (~q).transform(b); - lolunit_assert_doubles_equal(d.x, da.x, 1e-5); - lolunit_assert_doubles_equal(d.y, da.y, 1e-5); - lolunit_assert_doubles_equal(d.z, da.z, 1e-5); + lolunit_assert_doubles_equal(d.x, a.x, 2e-5); + lolunit_assert_doubles_equal(d.y, a.y, 2e-5); + lolunit_assert_doubles_equal(d.z, a.z, 2e-5); - if (distance(da, db) > 1e-6f) + if (distance(a, b) > 1e-6f) { - /* If da and db differ, check that the rotation axis is normal to both + /* If a and b differ, check that the rotation axis is normal to both * vectors, which is only true if the rotation uses the shortest path. */ vec3 axis = q.axis(); - lolunit_assert_doubles_equal(0.0, (double)dot(axis, da), 1e-5); - lolunit_assert_doubles_equal(0.0, (double)dot(axis, db), 1e-5); + lolunit_assert_doubles_equal(0.0, (double)dot(axis, a), 1e-5); + lolunit_assert_doubles_equal(0.0, (double)dot(axis, b), 1e-5); } else { - /* If da and db are roughly the same, check that the rotation angle + /* If a and b are roughly the same, check that the rotation angle * is zero. */ lolunit_assert_doubles_equal(0.0, (double)q.angle(), 1e-5); } @@ -533,7 +536,12 @@ lolunit_declare_fixture(quaternion_test) } private: - array m_vectorpairs; + void push_vector_pair(vec3 p, vec3 q) + { + m_vectorpairs.push_back(std::make_tuple(p, q)); + } + + std::vector> m_vectorpairs; }; } /* namespace lol */