From dc60e18862b23282ec29974d4797e43e698fdadd Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Wed, 25 May 2022 15:39:58 +0200 Subject: [PATCH 01/13] base mfloat --- engine/maths/mfloat.h | 262 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 262 insertions(+) create mode 100644 engine/maths/mfloat.h diff --git a/engine/maths/mfloat.h b/engine/maths/mfloat.h new file mode 100644 index 0000000000..fca626a28e --- /dev/null +++ b/engine/maths/mfloat.h @@ -0,0 +1,262 @@ + +/************************************************************************** + * * + * Regina - A Normal Surface Theory Calculator * + * Computational Engine * + * * + * Copyright (c) 1999-2021, Ben Burton * + * For further details contact Ben Burton (bab@debian.org). * + * * + * This program is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License as * + * published by the Free Software Foundation; either version 2 of the * + * License, or (at your option) any later version. * + * * + * As an exception, when this program is distributed through (i) the * + * App Store by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or * + * (iii) Google Play by Google Inc., then that store may impose any * + * digital rights management, device limits and/or redistribution * + * restrictions that are required by its terms of service. * + * * + * This program is distributed in the hope that it will be useful, but * + * WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public * + * License along with this program; if not, write to the Free * + * Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, * + * MA 02110-1301, USA. * + * * + **************************************************************************/ + +#ifndef __REGINA_MFLOAT +#ifndef __DOXYGEN +#define __REGINA_MFLOAT +#endif + +#define MPFRRND MPFR_RNDN + +#include "regina-core.h" +#include + + +namespace regina { + +class MFloat { + + private: + + mpfr_t number_; + bool isInit_ = false; + + static unsigned long prec_; + + public: + + MFloat(); + + explicit MFloat(double value); + + explicit MFloat(unsigned long value); + + static void setPrec(unsigned long prec) { + prec_ = prec; + } + + void set(double value); + + void set(unsigned long value); + + double getDouble(); + + double getDoubleClear(); + + void setPi(); + + ~MFloat(); + + void clear(); + + bool operator== ( const MFloat& rhs ) const; + + bool operator != (const MFloat& rhs) const; + + MFloat& operator = (const MFloat& value); + + MFloat& operator = (const double value); + + void negate(); + + void invert(); + + MFloat inverse() const; + + MFloat& operator += (const double value); + + MFloat& operator -= (const double value); + + MFloat& operator *= (const double value); + + MFloat& operator /= (const double value); + + MFloat& operator += (const unsigned long value); + + MFloat& operator -= (const unsigned long value); + + MFloat& operator *= (const unsigned long value); + + MFloat& operator /= (const unsigned long value); + + MFloat& operator += (const MFloat& other); + + MFloat& operator -= (const MFloat& other); + + MFloat& operator *= (const MFloat& other); + + MFloat& operator /= (const MFloat& other); + + void sin (); + + void static freeCache() { + mpfr_free_cache(); + mpfr_mp_memory_cleanup(); + } + +}; + +unsigned long MFloat::prec_ = 32; + +inline MFloat::MFloat() { + mpfr_init2(number_,prec_); + isInit_ = true; +} + + +inline MFloat::MFloat(double value) { + mpfr_init2(number_,prec_); + isInit_ = true; + mpfr_set_d(number_,value,MPFRRND); +} + +inline MFloat::MFloat(unsigned long value) { + mpfr_init2(number_,prec_); + isInit_ = true; + mpfr_set_ui(number_,value,MPFRRND); +} + +inline void MFloat::set(double value) { + mpfr_set_d(number_,value,MPFRRND); +} + +inline void MFloat::set(unsigned long value) { + mpfr_set_ui(number_,value,MPFRRND); +} + +inline double MFloat::getDouble(){ + return mpfr_get_d(number_,MPFRRND); +} + +inline MFloat::~MFloat() { + if (!isInit_) { + mpfr_clear(number_); + } + isInit_ = false; +} + +inline void MFloat::setPi() { + mpfr_const_pi(number_,MPFRRND); +} + +bool MFloat::operator== ( const MFloat& rhs ) const { + return mpfr_equal_p(number_,rhs.number_)!=0; +} + +inline bool MFloat::operator != (const MFloat& rhs) const { + return mpfr_equal_p(number_,rhs.number_)==0; +} + +inline MFloat& MFloat::operator = (const MFloat& other) { + mpfr_set(number_,other.number_,MPFRRND); + return *this; +} + + +inline MFloat& MFloat::operator = (const double value) { + mpfr_set_d(number_,value,MPFRRND); + return *this; +} + +inline void MFloat::negate() { + mpfr_neg(number_,number_,MPFRRND); +} + +inline MFloat& MFloat::operator += (const double value) { + mpfr_add_d(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator -= (const double value) { + mpfr_sub_d(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator *= (const double value) { + mpfr_mul_d(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator /= (const double value) { + mpfr_div_d(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator += (const unsigned long value) { + mpfr_add_ui(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator -= (const unsigned long value) { + mpfr_sub_ui(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator *= (const unsigned long value) { + mpfr_mul_ui(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator /= (const unsigned long value) { + mpfr_div_ui(number_,number_,value,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator += (const MFloat& other) { + mpfr_add(number_,number_,other.number_,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator -= (const MFloat& other) { + mpfr_sub(number_,number_,other.number_,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator *= (const MFloat& other) { + mpfr_mul(number_,number_,other.number_,MPFRRND); + return *this; +} + +inline MFloat& MFloat::operator /= (const MFloat& other) { + mpfr_div(number_,number_,other.number_,MPFRRND); + return *this; +} + +inline void MFloat::sin () { + mpfr_sin(number_,number_,MPFRRND); +}; + + + +} // namespace regina + +#endif From 2adaef794a43d2c07a6e4a3522bcc690d34ec416 Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Sat, 28 May 2022 19:47:21 +0200 Subject: [PATCH 02/13] cmake update --- CMakeLists.txt | 9 +++++++ cmake/modules/FindMPFR.cmake | 47 ++++++++++++++++++++++++++++++++++++ engine/CMakeLists.txt | 2 +- 3 files changed, 57 insertions(+), 1 deletion(-) create mode 100644 cmake/modules/FindMPFR.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 5441fd335c..b3a75ea3fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -306,6 +306,15 @@ MACRO_LOG_FEATURE(JANSSON_FOUND "Essential: JSON support" "http://www.digip.org/jansson/" TRUE) + +FIND_PACKAGE(MPFR) +INCLUDE_DIRECTORIES(${MPFR_INCLUDE_DIR}) +MACRO_LOG_FEATURE(MPFR_FOUND + "MPFR" + "High precision computation" + "http://www.mpfr.org/" + TRUE) + IF (NOT DISABLE_PYTHON) MESSAGE (STATUS "Searching for a Python installation") diff --git a/cmake/modules/FindMPFR.cmake b/cmake/modules/FindMPFR.cmake new file mode 100644 index 0000000000..69e7468ac4 --- /dev/null +++ b/cmake/modules/FindMPFR.cmake @@ -0,0 +1,47 @@ +# Try to find the MPFR libraries +# MPFR_FOUND - system has MPFR lib +# MPFR_INCLUDE_DIR - the MPFR include directory +# MPFR_LIBRARIES_DIR - Directory where the MPFR libraries are located +# MPFR_LIBRARIES - the MPFR libraries + +include(FindPackageHandleStandardArgs) + +if(MPFR_INCLUDE_DIR) + set(MPFR_in_cache TRUE) +else() + set(MPFR_in_cache FALSE) +endif() +if(NOT MPFR_LIBRARIES) + set(MPFR_in_cache FALSE) +endif() + +# Is it already configured? +if (NOT MPFR_in_cache) + + find_path(MPFR_INCLUDE_DIR + NAMES mpfr.h + HINTS ENV MPFR_INC_DIR + ENV MPFR_DIR + PATH_SUFFIXES include + DOC "The directory containing the MPFR header files" + ) + + find_library(MPFR_LIBRARIES NAMES mpfr libmpfr-4 libmpfr-1 + HINTS ENV MPFR_LIB_DIR + ENV MPFR_DIR + PATH_SUFFIXES lib + DOC "Path to the MPFR library" + ) + + if ( MPFR_LIBRARIES ) + get_filename_component(MPFR_LIBRARIES_DIR ${MPFR_LIBRARIES} PATH CACHE ) + endif() + + # Attempt to load a user-defined configuration for MPFR if couldn't be found + if ( NOT MPFR_INCLUDE_DIR OR NOT MPFR_LIBRARIES_DIR ) + include( MPFRConfig OPTIONAL ) + endif() + +endif() + +find_package_handle_standard_args(MPFR "DEFAULT_MSG" MPFR_LIBRARIES MPFR_INCLUDE_DIR) diff --git a/engine/CMakeLists.txt b/engine/CMakeLists.txt index f5120ec05b..d6d99357d3 100644 --- a/engine/CMakeLists.txt +++ b/engine/CMakeLists.txt @@ -56,7 +56,7 @@ ADD_LIBRARY("regina-engine" SHARED ${SOURCES} ) # that build against libregina-engine. # Note: PUBLIC and PRIVATE only appeard in CMake 2.8.12. # For compatibility back to 2.8.7 we use LINK_PUBLIC and LINK_PRIVATE instead. -SET(ENGINE_LINK_LIBRARIES ${LIBXML2_LIBRARIES} ${GMP_LIBRARIES} ${GMPXX_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT}) +SET(ENGINE_LINK_LIBRARIES ${LIBXML2_LIBRARIES} ${GMP_LIBRARIES} ${GMPXX_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${MPFR_LIBRARIES}) TARGET_LINK_LIBRARIES(regina-engine LINK_PUBLIC ${ENGINE_LINK_LIBRARIES} LINK_PRIVATE ${STDFS_LIBRARY} "${ICONV_LIBRARY}" ${ZLIB_LIBRARIES} ${JANSSON_LIBRARIES} ${KVSTORE_LIBRARIES}) From b30c24a038e11d2a6418eb374cbe57a2403a487b Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Mon, 30 May 2022 04:48:07 +0200 Subject: [PATCH 03/13] More cmake --- engine/maths/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/engine/maths/CMakeLists.txt b/engine/maths/CMakeLists.txt index afcfaa0119..3a63555574 100644 --- a/engine/maths/CMakeLists.txt +++ b/engine/maths/CMakeLists.txt @@ -30,6 +30,7 @@ if (${REGINA_INSTALL_DEV}) matrix.h matrix2.h matrixops.h + mfloat.h numbertheory.h perm.h perm-impl.h From 1f6f3c19bc9204dc26012622ce2975ba4c149f58 Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Mon, 30 May 2022 04:48:40 +0200 Subject: [PATCH 04/13] basic tests --- testsuite/maths/CMakeLists.txt | 1 + testsuite/maths/mfloat.cpp | 197 +++++++++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 testsuite/maths/mfloat.cpp diff --git a/testsuite/maths/CMakeLists.txt b/testsuite/maths/CMakeLists.txt index d1b307183b..b58c662f4e 100644 --- a/testsuite/maths/CMakeLists.txt +++ b/testsuite/maths/CMakeLists.txt @@ -9,6 +9,7 @@ SET ( FILES laurent2.cpp matrix.cpp matrixops.cpp + mfloat.cpp numbertheory.cpp perm.cpp perm2.cpp diff --git a/testsuite/maths/mfloat.cpp b/testsuite/maths/mfloat.cpp new file mode 100644 index 0000000000..cbb311d9ab --- /dev/null +++ b/testsuite/maths/mfloat.cpp @@ -0,0 +1,197 @@ + +/************************************************************************** + * * + * Regina - A Normal Surface Theory Calculator * + * Test Suite * + * * + * Copyright (c) 1999-2021, Ben Burton * + * For further details contact Ben Burton (bab@debian.org). * + * * + * This program is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License as * + * published by the Free Software Foundation; either version 2 of the * + * License, or (at your option) any later version. * + * * + * As an exception, when this program is distributed through (i) the * + * App Store by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or * + * (iii) Google Play by Google Inc., then that store may impose any * + * digital rights management, device limits and/or redistribution * + * restrictions that are required by its terms of service. * + * * + * This program is distributed in the hope that it will be useful, but * + * WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public * + * License along with this program; if not, write to the Free * + * Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, * + * MA 02110-1301, USA. * + * * + **************************************************************************/ + +#include +#include +#include "maths/mfloat.h" +#include "utilities/exception.h" +#include "utilities/stringutils.h" +#include "testsuite/utilities/testutilities.h" + +using regina::MFloat; + +class MFloatTest : public CppUnit::TestFixture { + CPPUNIT_TEST_SUITE(MFloatTest); + + CPPUNIT_TEST(constructFromInteger); + CPPUNIT_TEST(constructFromDouble); + CPPUNIT_TEST(verifyAssignments); + CPPUNIT_TEST(basicArithmetic); + + CPPUNIT_TEST_SUITE_END(); + + public: + static constexpr double epsilon = 0.0000001; + /** Used for determining whether a number is "close enough" + to zero. This helps avoid the inaccuracies inherent in + using == with floating point numbers. */ + + public: + void setUp() override { + } + + void tearDown() override { + } + + template + std::string str(T x) { + std::ostringstream ans; + ans << x; + return ans.str(); + } + + bool inEpsilon(double d1,double d2) { + return d1-d2 Date: Mon, 30 May 2022 04:49:48 +0200 Subject: [PATCH 05/13] Memory leak-free backtrac --- engine/maths/mfloat.h | 68 +++- engine/triangulation/dim3/triangulation3.h | 4 +- engine/triangulation/dim3/turaevviro.cpp | 390 ++++++++++++++++----- python/dim3/triangulation3.cpp | 3 +- 4 files changed, 368 insertions(+), 97 deletions(-) diff --git a/engine/maths/mfloat.h b/engine/maths/mfloat.h index fca626a28e..c35d09fdef 100644 --- a/engine/maths/mfloat.h +++ b/engine/maths/mfloat.h @@ -37,6 +37,7 @@ #define MPFRRND MPFR_RNDN +#include #include "regina-core.h" #include @@ -56,6 +57,10 @@ class MFloat { MFloat(); + MFloat(const MFloat& other); + + MFloat(MFloat&& other) noexcept; + explicit MFloat(double value); explicit MFloat(unsigned long value); @@ -64,20 +69,26 @@ class MFloat { prec_ = prec; } + static std::string str(const MFloat& value){ + return std::to_string(value.getDouble()); + } + void set(double value); void set(unsigned long value); - double getDouble(); - - double getDoubleClear(); + double getDouble() const; void setPi(); ~MFloat(); void clear(); + + void extract(MFloat&& rhs); + double static extractDouble(MFloat&& rhs); + bool operator== ( const MFloat& rhs ) const; bool operator != (const MFloat& rhs) const; @@ -123,6 +134,10 @@ class MFloat { mpfr_mp_memory_cleanup(); } + friend std::ostream& operator<< ( std::ostream& out, const regina::MFloat& other); + + + }; unsigned long MFloat::prec_ = 32; @@ -132,6 +147,18 @@ inline MFloat::MFloat() { isInit_ = true; } +inline MFloat::MFloat(const MFloat& other) { + mpfr_init2(number_,prec_); + isInit_ = true; + mpfr_set(number_,other.number_,MPFRRND); +} + +inline MFloat::MFloat(MFloat&& other) noexcept { + if(other.isInit_) { + isInit_ = true; + std::swap (number_,other.number_); + } +} inline MFloat::MFloat(double value) { mpfr_init2(number_,prec_); @@ -146,24 +173,48 @@ inline MFloat::MFloat(unsigned long value) { } inline void MFloat::set(double value) { + if (!isInit_) { + mpfr_init2(number_,prec_); + } mpfr_set_d(number_,value,MPFRRND); } inline void MFloat::set(unsigned long value) { + if (!isInit_) { + mpfr_init2(number_,prec_); + } mpfr_set_ui(number_,value,MPFRRND); } -inline double MFloat::getDouble(){ +inline double MFloat::getDouble() const{ return mpfr_get_d(number_,MPFRRND); } -inline MFloat::~MFloat() { - if (!isInit_) { +inline MFloat::~MFloat() {} + +inline void MFloat::clear() { + if (isInit_) { mpfr_clear(number_); } isInit_ = false; } +inline void MFloat::extract(MFloat&& rhs) { + if (rhs.isInit_) { + mpfr_set(number_,rhs.number_,MPFRRND); + mpfr_clear(rhs.number_); + } +} + +inline double MFloat::extractDouble(MFloat&& rhs) { + double d = 0; + if (rhs.isInit_) { + d=rhs.getDouble(); + mpfr_clear(rhs.number_); + } + return d; +} + inline void MFloat::setPi() { mpfr_const_pi(number_,MPFRRND); } @@ -255,7 +306,10 @@ inline void MFloat::sin () { mpfr_sin(number_,number_,MPFRRND); }; - +std::ostream& operator << (std::ostream& out, const MFloat& value) { + out << std::to_string(value.getDouble()); + return out; +} } // namespace regina diff --git a/engine/triangulation/dim3/triangulation3.h b/engine/triangulation/dim3/triangulation3.h index 3c74f7f6d5..10d13b0144 100644 --- a/engine/triangulation/dim3/triangulation3.h +++ b/engine/triangulation/dim3/triangulation3.h @@ -847,12 +847,14 @@ class Triangulation<3> : public detail::TriangulationBase<3> { * as described above. * @param alg the algorithm with which to compute the invariant. If * you are not sure, the default value (ALG_DEFAULT) is a safe choice. + * @param acc is the number of bits to be used by the computation: + * under 64 included, doubles are used, otherwise mpfr is used. * @return the requested Turaev-Viro invariant. * * @see allCalculatedTuraevViro */ double turaevViroApprox(unsigned long r, unsigned long whichRoot = 1, - Algorithm alg = ALG_DEFAULT) const; + Algorithm alg = ALG_DEFAULT, unsigned long prec = 64) const; /** * Returns the cache of all Turaev-Viro state sum invariants that * have been calculated for this 3-manifold. diff --git a/engine/triangulation/dim3/turaevviro.cpp b/engine/triangulation/dim3/turaevviro.cpp index 910cfe68fb..75cf8a18ab 100644 --- a/engine/triangulation/dim3/turaevviro.cpp +++ b/engine/triangulation/dim3/turaevviro.cpp @@ -40,6 +40,7 @@ #include "libnormaliz/cone.h" #include "maths/cyclotomic.h" #include "maths/numbertheory.h" +#include "maths/mfloat.h" #include "progress/progresstracker.h" #include "treewidth/treedecomposition.h" #include "triangulation/dim3.h" @@ -47,6 +48,7 @@ #include #include + // #define TV_BACKTRACK_DUMP_COLOURINGS // #define TV_IGNORE_CACHE @@ -62,13 +64,14 @@ namespace regina { namespace { - template + enum calcType {exact, floating, multiPrecision}; + template struct TuraevViroDetails; template <> - struct TuraevViroDetails { - using TVType = Cyclotomic; - using TVResult = Cyclotomic; + struct TuraevViroDetails { + typedef Cyclotomic TVType; + typedef Cyclotomic TVResult; static TVResult zero() { return Cyclotomic(1); @@ -76,24 +79,35 @@ namespace { }; template <> - struct TuraevViroDetails { - using TVType = std::complex; - using TVResult = double; + struct TuraevViroDetails { + typedef std::complex TVType; + typedef double TVResult; static TVResult zero() { return 0; } }; + template <> + struct TuraevViroDetails { + typedef MFloat TVType; + typedef MFloat TVResult; + + static TVResult zero() { + MFloat zero = MFloat(0.0); + return zero; + } + }; + /** * Allows calculation of [n]! for arbitrary n. * Values are cached as they are calculated. */ - template + template class BracketFactorial { public: - using TVType = typename TuraevViroDetails::TVType; - using TVResult = typename TuraevViroDetails::TVResult; + using TVType = typename TuraevViroDetails::TVType; + using TVResult = typename TuraevViroDetails::TVResult; private: TVResult* bracket_; @@ -120,6 +134,11 @@ namespace { delete[] inv_; } + /** + * To clean up the content of the arrays of the destuctor. + */ + void clearMembers(unsigned long r); + /** * Returns the single value [index] (with no factorial symbol). * Requires index < r. @@ -146,7 +165,7 @@ namespace { }; template <> - BracketFactorial::BracketFactorial( + BracketFactorial::BracketFactorial( unsigned long r, unsigned long whichRoot) : bracket_(new TVResult[r]), fact_(new TVResult[r]), @@ -185,7 +204,7 @@ namespace { } template <> - BracketFactorial::BracketFactorial( + BracketFactorial::BracketFactorial( unsigned long r, unsigned long whichRoot) : bracket_(new TVResult[r]), fact_(new TVResult[r]), @@ -199,20 +218,70 @@ namespace { inv_[i] = inv_[i - 1] / bracket_[i]; } } + + template <> + BracketFactorial::BracketFactorial( + unsigned long r, unsigned long whichRoot) : + bracket_(new TVResult[r]), + fact_(new TVResult[r]), + inv_(new TVResult[r]) { + + TVResult angle,sinangle,siniangle,facti,invi; + + angle.setPi(); + angle *= whichRoot; + angle /= r; + + unsigned long one = 1; + + bracket_[0].set(one); + bracket_[1].set(one); + fact_[0].set(one); + fact_[1].set(one); + inv_[0].set(one); + inv_[1].set(one); + + sinangle = angle; + sinangle.sin(); + + siniangle = angle; + facti = fact_[1]; + invi = inv_[1]; + + for (unsigned long i = 2; i < r; i++) { + siniangle *= i; + siniangle.sin(); + bracket_[i] = siniangle; + bracket_[i] /= sinangle; + + facti *= bracket_[i]; + fact_[i] = facti; + invi /= bracket_[i]; + inv_[i] = invi; + + siniangle = angle; + } + + angle.clear(); + sinangle.clear(); + siniangle.clear(); + facti.clear(); + invi.clear(); + } /** * Represents the initial data as described in Section 7 of Turaev * and Viro's paper. */ - template + template struct InitialData { - using TVType = typename TuraevViroDetails::TVType; - using TVResult = typename TuraevViroDetails::TVResult; + using TVType = typename TuraevViroDetails::TVType; + using TVResult = typename TuraevViroDetails::TVResult; unsigned long r, whichRoot; /**< The Turaev-Viro parameters. */ bool halfField; - BracketFactorial fact; + BracketFactorial fact; /**< The cached values [n]!. */ TVType vertexContrib; /**< The vertex-based contribution to the Turaev-Viro invariant; @@ -220,6 +289,15 @@ namespace { InitialData(unsigned long newR, unsigned long newWhichRoot); + /** + * To clear BracketFactorial, to clear and delete an array of TVType of size n and to clear an element of type TVType. + */ + void clearData(); + static void clearVar(TVType& var); + static void clearArray(TVType* ar,unsigned long n); +// template + static void clearMap(std::map, TVType> *map); + static void negate(TVType& x); void initZero(TVType& x) const; @@ -310,6 +388,7 @@ namespace { ansToOverwrite -= term; } } + clearVar(term); } /** @@ -373,7 +452,7 @@ namespace { }; template <> - InitialData::InitialData( + InitialData::InitialData( unsigned long newR, unsigned long newWhichRoot) : r(newR), whichRoot(newWhichRoot), @@ -392,7 +471,7 @@ namespace { } template <> - InitialData::InitialData( + InitialData::InitialData( unsigned long newR, unsigned long newWhichRoot) : r(newR), whichRoot(newWhichRoot), @@ -401,48 +480,151 @@ namespace { double tmp = sin(M_PI * whichRoot / r); vertexContrib = 2.0 * tmp * tmp / r; } + + template <> + InitialData::InitialData( + unsigned long newR, unsigned long newWhichRoot) : + r(newR), + whichRoot(newWhichRoot), + halfField(r % 2 != 0 && whichRoot % 2 == 0), + fact(r, whichRoot) { + TVResult tmp; + tmp.setPi(); + tmp *= whichRoot; + tmp /= r; + tmp.sin(); + tmp *= tmp; + unsigned long two = 2; + tmp *= two; + tmp /= r; + vertexContrib = tmp; + tmp.clear(); + } + + template <> + inline void BracketFactorial::clearMembers(unsigned long r) { + for (unsigned long i = 0; i + inline void InitialData::clearData() { + fact.clearMembers(r); + vertexContrib.clear(); + } + + + + template <> + inline void InitialData::clearVar(TVType& var) {} + + template <> + inline void InitialData::clearVar(TVType& var) {} + + template <> + inline void InitialData::clearVar(TVType& var) { + var.clear(); + } + + template <> + inline void InitialData::clearArray(TVType* ar,unsigned long n) { + delete[] ar; + } + + template <> + inline void InitialData::clearArray(TVType* ar,unsigned long n) { + delete[] ar; + } + + template <> + inline void InitialData::clearArray(TVType* ar,unsigned long n) { + for (unsigned long i=0;i - inline void InitialData::negate(InitialData::TVType& x) { +// template + inline void InitialData::clearMap(std::map, TVType>* map) { + delete map; + } + + template <> +// template + inline void InitialData::clearMap(std::map, TVType>* map) { + delete[] map; + } + + template<> +// template + inline void InitialData::clearMap(std::map, TVType>* map) { + for(auto &x:*map) {x.second.clear();} + delete map; + } + + + template <> + inline void InitialData::negate(InitialData::TVType& x) { x.negate(); } template <> - inline void InitialData::negate(InitialData::TVType& x) { + inline void InitialData::negate(InitialData::TVType& x) { x = -x; } template <> - inline void InitialData::initZero(InitialData::TVType& x) + inline void InitialData::negate(InitialData::TVType& x) { + x.negate(); + } + + template <> + inline void InitialData::initZero(InitialData::TVType& x) const { x.init(halfField ? r : 2 * r); } template <> - inline void InitialData::initZero(InitialData::TVType& x) + inline void InitialData::initZero(InitialData::TVType& x) const { x = 0.0; } template <> - inline void InitialData::initOne(InitialData::TVType& x) + inline void InitialData::initZero(InitialData::TVType& x) + const { + x.set(0.0); + } + + template <> + inline void InitialData::initOne(InitialData::TVType& x) const { x.init(halfField ? r : 2 * r); x[0] = 1; } template <> - inline void InitialData::initOne(InitialData::TVType& x) + inline void InitialData::initOne(InitialData::TVType& x) const { x = 1.0; } - template - typename InitialData::TVType turaevViroBacktrack( + template <> + inline void InitialData::initOne(InitialData::TVType& x) + const { + x.set(1.0); + } + + template + typename InitialData::TVType turaevViroBacktrack( const Triangulation<3>& tri, - const InitialData& init, + const InitialData& init, ProgressTracker* tracker) { - using TVType = typename InitialData::TVType; + using TVType = typename InitialData::TVType; if (tracker) tracker->newStage("Enumerating colourings"); @@ -648,14 +830,17 @@ namespace { delete[] tetDone; delete[] tetDoneStart; - delete[] edgeCache; - delete[] triangleCache; - delete[] tetCache; + InitialData::clearArray(edgeCache,nEdges + 1); + InitialData::clearArray(triangleCache,nEdges + 1); + InitialData::clearArray(tetCache,nEdges + 1); + + InitialData::clearVar(valColour); + InitialData::clearVar(tmpTVType); if (tracker) { delete[] coeff; if (tracker->isCancelled()) - return TuraevViroDetails::zero(); + return TuraevViroDetails::zero(); } // Compute the vertex contributions separately, since these are @@ -666,12 +851,12 @@ namespace { return ans; } - template - typename InitialData::TVType turaevViroNaive( + template + typename InitialData::TVType turaevViroNaive( const Triangulation<3>& tri, - const InitialData& init, + const InitialData& init, ProgressTracker* tracker) { - using TVType = typename InitialData::TVType; + using TVType = typename InitialData::TVType; if (tracker) tracker->newStage("Enumerating colourings"); @@ -810,7 +995,7 @@ namespace { if (tracker) { delete[] coeff; if (tracker->isCancelled()) - return TuraevViroDetails::zero(); + return TuraevViroDetails::zero(); } // Compute the vertex contributions separately, since these are @@ -818,15 +1003,16 @@ namespace { for (i = 0; i < tri.countVertices(); i++) ans *= init.vertexContrib; + InitialData::clearVar(valColour); return ans; } - template - typename InitialData::TVType turaevViroTreewidth( + template + typename InitialData::TVType turaevViroTreewidth( const Triangulation<3>& tri, - InitialData& init, + InitialData& init, ProgressTracker* tracker) { - using TVType = typename InitialData::TVType; + using TVType = typename InitialData::TVType; // Progress: // - weight of forget/join bag processing is 0.9 @@ -850,7 +1036,7 @@ namespace { if (tracker) { if (tracker->isCancelled()) - return TuraevViroDetails::zero(); + return TuraevViroDetails::zero(); tracker->newStage("Analysing bags", 0.01); } @@ -1050,9 +1236,10 @@ namespace { // solutions if need be. existingSoln = partial[index]->try_emplace( std::move(seq), std::move(val)); - if (! existingSoln.second) + if (! existingSoln.second) { existingSoln.first->second += val; - + InitialData::clearVar(val); + } ++level; while (level < 6 && choiceType[level] != 0) ++level; @@ -1103,7 +1290,7 @@ namespace { } } - delete partial[child->index()]; + InitialData::clearMap(partial[child->index()]); partial[child->index()] = nullptr; } else { // Join bag. @@ -1195,7 +1382,8 @@ namespace { // We have two compatible solutions. // Combine them and store the corresponding // value, again aggregating if necessary. - TVType val = (*subit)->second * (*subit2)->second; + TVType val = (*subit)->second; + val *=(*subit2)->second; LightweightSequence seq(nEdges); for (i = 0; i < nEdges; ++i) @@ -1208,18 +1396,20 @@ namespace { existingSoln = partial[index]->try_emplace( std::move(seq), std::move(val)); - if (! existingSoln.second) + if (! existingSoln.second) { existingSoln.first->second += val; + InitialData::clearVar(val); + } } } delete[] leftIndexed; delete[] rightIndexed; - delete partial[child->index()]; + InitialData::clearMap(partial[child->index()]); partial[child->index()] = nullptr; - delete partial[sibling->index()]; + InitialData::clearMap(partial[sibling->index()]); partial[sibling->index()] = nullptr; } @@ -1245,7 +1435,7 @@ namespace { delete partial[i]; delete[] partial; - return TuraevViroDetails::zero(); + return TuraevViroDetails::zero(); } // We made it to the end. @@ -1255,7 +1445,7 @@ namespace { // only one colouring stored (in which all edge colours are aggregated). TVType ans = partial[nBags - 1]->begin()->second; - delete partial[nBags - 1]; + InitialData::clearMap(partial[nBags - 1]); delete[] partial; for (i = 0; i < tri.countVertices(); i++) @@ -1263,11 +1453,11 @@ namespace { return ans; } - template - typename InitialData::TVType turaevViroPolytope( + template + typename InitialData::TVType turaevViroPolytope( const Triangulation<3>& tri, - InitialData& init) { - using TVType = typename InitialData::TVType; + InitialData& init) { + using TVType = typename InitialData::TVType; std::vector > input; unsigned long nTri = tri.countTriangles(); @@ -1334,7 +1524,7 @@ namespace { } double Triangulation<3>::turaevViroApprox(unsigned long r, - unsigned long whichRoot, Algorithm alg) const { + unsigned long whichRoot, Algorithm alg, unsigned long prec) const { // Do some basic parameter checks. if (r < 3) return 0; @@ -1343,36 +1533,60 @@ double Triangulation<3>::turaevViroApprox(unsigned long r, if (gcd(r, whichRoot) > 1) return 0; - // Set up our initial data. - InitialData init(r, whichRoot); - - InitialData::TVType ans; - switch (alg) { - case ALG_BACKTRACK: - ans = turaevViroBacktrack(*this, init, nullptr); - break; - case ALG_NAIVE: - ans = turaevViroNaive(*this, init, nullptr); - break; - default: - ans = turaevViroTreewidth(*this, init, nullptr); - break; - } - /* - * Disable this check for now, since testing whether img(z) == 0 is - * error-prone due to floating-point approximation. - * - if (isNonZero(ans.imag())) { - // This should never happen, since the Turaev-Viro invariant is the - // square of the modulus of the Witten invariant for sl_2. - std::cerr << - "WARNING: The Turaev-Viro invariant has an imaginary component.\n" - " This should never happen.\n" - " Please report this (along with the 3-manifold that" - " was used) to Regina's authors." << std::endl; - } - */ - return ans.real(); + if (prec<65) { + // Set up our initial data. + InitialData init(r, whichRoot); + + InitialData::TVType ans; + switch (alg) { + case ALG_BACKTRACK: + ans = turaevViroBacktrack(*this, init, nullptr); + break; + case ALG_NAIVE: + ans = turaevViroNaive(*this, init, nullptr); + break; + default: + ans = turaevViroTreewidth(*this, init, nullptr); + break; + } + /* + * Disable this check for now, since testing whether img(z) == 0 is + * error-prone due to floating-point approximation. + * + if (isNonZero(ans.imag())) { + // This should never happen, since the Turaev-Viro invariant is the + // square of the modulus of the Witten invariant for sl_2. + std::cerr << + "WARNING: The Turaev-Viro invariant has an imaginary component.\n" + " This should never happen.\n" + " Please report this (along with the 3-manifold that" + " was used) to Regina's authors." << std::endl; + } + */ + return ans.real(); + } else { + // Set up our initial data. + MFloat::setPrec(prec); + InitialData init(r, whichRoot); + + double ans; + switch (alg) { + case ALG_BACKTRACK: + ans = MFloat::extractDouble(turaevViroBacktrack(*this, init, 0)); + break; + case ALG_NAIVE: + ans = MFloat::extractDouble(turaevViroNaive(*this, init, 0)); + break; + default: + ans = MFloat::extractDouble(turaevViroTreewidth(*this, init, 0)); + break; + } + init.clearData(); + MFloat::freeCache(); + std::cout.flush(); + return ans; + } + } Cyclotomic Triangulation<3>::turaevViro(unsigned long r, bool parity, @@ -1398,9 +1612,9 @@ Cyclotomic Triangulation<3>::turaevViro(unsigned long r, bool parity, #endif // Set up our initial data. - InitialData init(r, (parity ? 1 : 0)); + InitialData init(r, (parity ? 1 : 0)); - InitialData::TVType ans; + InitialData::TVType ans; switch (alg) { case ALG_BACKTRACK: ans = turaevViroBacktrack(*this, init, tracker); diff --git a/python/dim3/triangulation3.cpp b/python/dim3/triangulation3.cpp index d1c71a44a1..cdbca88f7a 100644 --- a/python/dim3/triangulation3.cpp +++ b/python/dim3/triangulation3.cpp @@ -273,7 +273,8 @@ void addTriangulation3(pybind11::module_& m) { pybind11::call_guard()) .def("turaevViroApprox", &Triangulation<3>::turaevViroApprox, pybind11::arg(), pybind11::arg("whichRoot") = 1, - pybind11::arg("alg") = regina::ALG_DEFAULT) + pybind11::arg("alg") = regina::ALG_DEFAULT, + pybind11::arg("prec") = 128) .def("allCalculatedTuraevViro", &Triangulation<3>::allCalculatedTuraevViro) .def("longitude", &Triangulation<3>::longitude, From b1afe93426ddad44d3e288e72e9ec656f2b2a63b Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Mon, 30 May 2022 18:36:14 +0200 Subject: [PATCH 06/13] all leaks fixed --- engine/triangulation/dim3/turaevviro.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/engine/triangulation/dim3/turaevviro.cpp b/engine/triangulation/dim3/turaevviro.cpp index 75cf8a18ab..c7680a24b3 100644 --- a/engine/triangulation/dim3/turaevviro.cpp +++ b/engine/triangulation/dim3/turaevviro.cpp @@ -295,7 +295,6 @@ namespace { void clearData(); static void clearVar(TVType& var); static void clearArray(TVType* ar,unsigned long n); -// template static void clearMap(std::map, TVType> *map); static void negate(TVType& x); @@ -410,6 +409,7 @@ namespace { tetContrib(colour0, colour1, colour3, colour5, colour4, colour2, tmp); ans *= tmp; + clearVar(tmp); int i; const Triangle<3>* triangle; @@ -548,19 +548,16 @@ namespace { } template <> -// template inline void InitialData::clearMap(std::map, TVType>* map) { delete map; } template <> -// template inline void InitialData::clearMap(std::map, TVType>* map) { - delete[] map; + delete map; } template<> -// template inline void InitialData::clearMap(std::map, TVType>* map) { for(auto &x:*map) {x.second.clear();} delete map; From ddc326f2ed701bd7d9750ffde41d41627a2acaca Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Mon, 30 May 2022 18:41:29 +0200 Subject: [PATCH 07/13] cancelled leaks fixed --- engine/triangulation/dim3/turaevviro.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/engine/triangulation/dim3/turaevviro.cpp b/engine/triangulation/dim3/turaevviro.cpp index c7680a24b3..b231718eba 100644 --- a/engine/triangulation/dim3/turaevviro.cpp +++ b/engine/triangulation/dim3/turaevviro.cpp @@ -836,8 +836,10 @@ namespace { if (tracker) { delete[] coeff; - if (tracker->isCancelled()) + if (tracker->isCancelled()) { + InitialData::clearVar(ans); return TuraevViroDetails::zero(); + } } // Compute the vertex contributions separately, since these are @@ -988,11 +990,14 @@ namespace { delete[] colour; delete[] sortedEdges; delete[] edgePos; + InitialData::clearVar(valColour); if (tracker) { delete[] coeff; - if (tracker->isCancelled()) + if (tracker->isCancelled()) { + InitialData::clearVar(ans); return TuraevViroDetails::zero(); + } } // Compute the vertex contributions separately, since these are @@ -1000,7 +1005,7 @@ namespace { for (i = 0; i < tri.countVertices(); i++) ans *= init.vertexContrib; - InitialData::clearVar(valColour); + return ans; } @@ -1429,7 +1434,9 @@ namespace { // We don't know which elements of partial[] have been // deallocated, so check them all. for (i = 0; i < nBags; ++i) - delete partial[i]; + if(partial[i]!=nullptr) { + InitialData::clearMap(partial[i]); + } delete[] partial; return TuraevViroDetails::zero(); From 322755f4e15d90a43865f6236d83fe741f1a7626 Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Mon, 30 May 2022 19:15:50 +0200 Subject: [PATCH 08/13] fixed default precision for python --- python/dim3/triangulation3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dim3/triangulation3.cpp b/python/dim3/triangulation3.cpp index cdbca88f7a..f5ec1552c3 100644 --- a/python/dim3/triangulation3.cpp +++ b/python/dim3/triangulation3.cpp @@ -274,7 +274,7 @@ void addTriangulation3(pybind11::module_& m) { .def("turaevViroApprox", &Triangulation<3>::turaevViroApprox, pybind11::arg(), pybind11::arg("whichRoot") = 1, pybind11::arg("alg") = regina::ALG_DEFAULT, - pybind11::arg("prec") = 128) + pybind11::arg("prec") = 64) .def("allCalculatedTuraevViro", &Triangulation<3>::allCalculatedTuraevViro) .def("longitude", &Triangulation<3>::longitude, From 661f7b017d4e18b43de97408ff8ec3954fbc810c Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Tue, 7 Jun 2022 17:34:25 +0200 Subject: [PATCH 09/13] removed clear mess --- engine/maths/mfloat.h | 54 +++------- engine/triangulation/dim3/turaevviro.cpp | 124 +++-------------------- 2 files changed, 30 insertions(+), 148 deletions(-) diff --git a/engine/maths/mfloat.h b/engine/maths/mfloat.h index c35d09fdef..013c8a7a8d 100644 --- a/engine/maths/mfloat.h +++ b/engine/maths/mfloat.h @@ -49,8 +49,7 @@ class MFloat { private: mpfr_t number_; - bool isInit_ = false; - + int flag = 0; static unsigned long prec_; public: @@ -61,6 +60,8 @@ class MFloat { MFloat(MFloat&& other) noexcept; + MFloat& operator=(MFloat&& other); + explicit MFloat(double value); explicit MFloat(unsigned long value); @@ -83,10 +84,6 @@ class MFloat { ~MFloat(); - void clear(); - - void extract(MFloat&& rhs); - double static extractDouble(MFloat&& rhs); bool operator== ( const MFloat& rhs ) const; @@ -137,52 +134,47 @@ class MFloat { friend std::ostream& operator<< ( std::ostream& out, const regina::MFloat& other); - }; unsigned long MFloat::prec_ = 32; inline MFloat::MFloat() { mpfr_init2(number_,prec_); - isInit_ = true; + mpfr_set_d(number_,0,MPFRRND); } inline MFloat::MFloat(const MFloat& other) { mpfr_init2(number_,prec_); - isInit_ = true; mpfr_set(number_,other.number_,MPFRRND); } inline MFloat::MFloat(MFloat&& other) noexcept { - if(other.isInit_) { - isInit_ = true; - std::swap (number_,other.number_); + number_->_mpfr_d = 0; + mpfr_swap (number_,other.number_); +} + +inline MFloat& MFloat::operator=(MFloat&& other) { + if (this != &other) { + mpfr_swap(number_,other.number_); } + return *this; } inline MFloat::MFloat(double value) { mpfr_init2(number_,prec_); - isInit_ = true; mpfr_set_d(number_,value,MPFRRND); } inline MFloat::MFloat(unsigned long value) { mpfr_init2(number_,prec_); - isInit_ = true; mpfr_set_ui(number_,value,MPFRRND); } inline void MFloat::set(double value) { - if (!isInit_) { - mpfr_init2(number_,prec_); - } mpfr_set_d(number_,value,MPFRRND); } inline void MFloat::set(unsigned long value) { - if (!isInit_) { - mpfr_init2(number_,prec_); - } mpfr_set_ui(number_,value,MPFRRND); } @@ -190,29 +182,13 @@ inline double MFloat::getDouble() const{ return mpfr_get_d(number_,MPFRRND); } -inline MFloat::~MFloat() {} - -inline void MFloat::clear() { - if (isInit_) { +inline MFloat::~MFloat() { + if(0 != number_->_mpfr_d) mpfr_clear(number_); - } - isInit_ = false; -} - -inline void MFloat::extract(MFloat&& rhs) { - if (rhs.isInit_) { - mpfr_set(number_,rhs.number_,MPFRRND); - mpfr_clear(rhs.number_); - } } inline double MFloat::extractDouble(MFloat&& rhs) { - double d = 0; - if (rhs.isInit_) { - d=rhs.getDouble(); - mpfr_clear(rhs.number_); - } - return d; + return rhs.getDouble(); } inline void MFloat::setPi() { diff --git a/engine/triangulation/dim3/turaevviro.cpp b/engine/triangulation/dim3/turaevviro.cpp index b231718eba..9fb95063bb 100644 --- a/engine/triangulation/dim3/turaevviro.cpp +++ b/engine/triangulation/dim3/turaevviro.cpp @@ -134,11 +134,6 @@ namespace { delete[] inv_; } - /** - * To clean up the content of the arrays of the destuctor. - */ - void clearMembers(unsigned long r); - /** * Returns the single value [index] (with no factorial symbol). * Requires index < r. @@ -261,12 +256,6 @@ namespace { siniangle = angle; } - - angle.clear(); - sinangle.clear(); - siniangle.clear(); - facti.clear(); - invi.clear(); } /** @@ -289,14 +278,6 @@ namespace { InitialData(unsigned long newR, unsigned long newWhichRoot); - /** - * To clear BracketFactorial, to clear and delete an array of TVType of size n and to clear an element of type TVType. - */ - void clearData(); - static void clearVar(TVType& var); - static void clearArray(TVType* ar,unsigned long n); - static void clearMap(std::map, TVType> *map); - static void negate(TVType& x); void initZero(TVType& x) const; @@ -387,7 +368,6 @@ namespace { ansToOverwrite -= term; } } - clearVar(term); } /** @@ -409,7 +389,6 @@ namespace { tetContrib(colour0, colour1, colour3, colour5, colour4, colour2, tmp); ans *= tmp; - clearVar(tmp); int i; const Triangle<3>* triangle; @@ -498,69 +477,6 @@ namespace { tmp *= two; tmp /= r; vertexContrib = tmp; - tmp.clear(); - } - - template <> - inline void BracketFactorial::clearMembers(unsigned long r) { - for (unsigned long i = 0; i - inline void InitialData::clearData() { - fact.clearMembers(r); - vertexContrib.clear(); - } - - - - template <> - inline void InitialData::clearVar(TVType& var) {} - - template <> - inline void InitialData::clearVar(TVType& var) {} - - template <> - inline void InitialData::clearVar(TVType& var) { - var.clear(); - } - - template <> - inline void InitialData::clearArray(TVType* ar,unsigned long n) { - delete[] ar; - } - - template <> - inline void InitialData::clearArray(TVType* ar,unsigned long n) { - delete[] ar; - } - - template <> - inline void InitialData::clearArray(TVType* ar,unsigned long n) { - for (unsigned long i=0;i - inline void InitialData::clearMap(std::map, TVType>* map) { - delete map; - } - - template <> - inline void InitialData::clearMap(std::map, TVType>* map) { - delete map; - } - - template<> - inline void InitialData::clearMap(std::map, TVType>* map) { - for(auto &x:*map) {x.second.clear();} - delete map; } @@ -827,17 +743,14 @@ namespace { delete[] tetDone; delete[] tetDoneStart; - InitialData::clearArray(edgeCache,nEdges + 1); - InitialData::clearArray(triangleCache,nEdges + 1); - InitialData::clearArray(tetCache,nEdges + 1); - InitialData::clearVar(valColour); - InitialData::clearVar(tmpTVType); + delete[] edgeCache; + delete[] triangleCache; + delete[] tetCache;; if (tracker) { delete[] coeff; if (tracker->isCancelled()) { - InitialData::clearVar(ans); return TuraevViroDetails::zero(); } } @@ -990,12 +903,10 @@ namespace { delete[] colour; delete[] sortedEdges; delete[] edgePos; - InitialData::clearVar(valColour); if (tracker) { delete[] coeff; if (tracker->isCancelled()) { - InitialData::clearVar(ans); return TuraevViroDetails::zero(); } } @@ -1240,7 +1151,6 @@ namespace { std::move(seq), std::move(val)); if (! existingSoln.second) { existingSoln.first->second += val; - InitialData::clearVar(val); } ++level; while (level < 6 && choiceType[level] != 0) @@ -1292,7 +1202,7 @@ namespace { } } - InitialData::clearMap(partial[child->index()]); + delete partial[child->index()]; partial[child->index()] = nullptr; } else { // Join bag. @@ -1400,7 +1310,6 @@ namespace { std::move(seq), std::move(val)); if (! existingSoln.second) { existingSoln.first->second += val; - InitialData::clearVar(val); } } } @@ -1408,10 +1317,10 @@ namespace { delete[] leftIndexed; delete[] rightIndexed; - InitialData::clearMap(partial[child->index()]); + delete partial[child->index()]; partial[child->index()] = nullptr; - InitialData::clearMap(partial[sibling->index()]); + delete partial[child->index()]; partial[sibling->index()] = nullptr; } @@ -1434,9 +1343,8 @@ namespace { // We don't know which elements of partial[] have been // deallocated, so check them all. for (i = 0; i < nBags; ++i) - if(partial[i]!=nullptr) { - InitialData::clearMap(partial[i]); - } + if(partial[i]!=nullptr) + delete partial[i]; delete[] partial; return TuraevViroDetails::zero(); @@ -1449,7 +1357,7 @@ namespace { // only one colouring stored (in which all edge colours are aggregated). TVType ans = partial[nBags - 1]->begin()->second; - InitialData::clearMap(partial[nBags - 1]); + delete partial[nBags - 1]; delete[] partial; for (i = 0; i < tri.countVertices(); i++) @@ -1573,22 +1481,20 @@ double Triangulation<3>::turaevViroApprox(unsigned long r, MFloat::setPrec(prec); InitialData init(r, whichRoot); - double ans; + InitialData::TVType ans; switch (alg) { case ALG_BACKTRACK: - ans = MFloat::extractDouble(turaevViroBacktrack(*this, init, 0)); + ans = turaevViroBacktrack(*this, init, 0); break; case ALG_NAIVE: - ans = MFloat::extractDouble(turaevViroNaive(*this, init, 0)); + ans = turaevViroNaive(*this, init, 0); break; default: - ans = MFloat::extractDouble(turaevViroTreewidth(*this, init, 0)); + ans = turaevViroTreewidth(*this, init, 0); break; } - init.clearData(); - MFloat::freeCache(); - std::cout.flush(); - return ans; + + return ans.getDouble(); } } From 4c71d127517740151958f59c9b8c16ac39064098 Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Tue, 7 Jun 2022 17:35:17 +0200 Subject: [PATCH 10/13] removed flag too --- engine/maths/mfloat.h | 1 - 1 file changed, 1 deletion(-) diff --git a/engine/maths/mfloat.h b/engine/maths/mfloat.h index 013c8a7a8d..a36f18495c 100644 --- a/engine/maths/mfloat.h +++ b/engine/maths/mfloat.h @@ -49,7 +49,6 @@ class MFloat { private: mpfr_t number_; - int flag = 0; static unsigned long prec_; public: From 3955c2ed69d7cf8840a94c968a35ca95c262082b Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Tue, 7 Jun 2022 17:40:03 +0200 Subject: [PATCH 11/13] fixed test suite --- testsuite/maths/mfloat.cpp | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/testsuite/maths/mfloat.cpp b/testsuite/maths/mfloat.cpp index cbb311d9ab..29b59195fb 100644 --- a/testsuite/maths/mfloat.cpp +++ b/testsuite/maths/mfloat.cpp @@ -81,7 +81,6 @@ class MFloatTest : public CppUnit::TestFixture { " incorrectly assigned from ul."; CPPUNIT_FAIL(msg.str()); } - num.clear(); } void verifyConstructd(double val) { @@ -92,7 +91,6 @@ class MFloatTest : public CppUnit::TestFixture { " incorrectly assigned from ul."; CPPUNIT_FAIL(msg.str()); } - num.clear(); } void constructFromInteger() { @@ -126,9 +124,6 @@ class MFloatTest : public CppUnit::TestFixture { msg << "Incorrect MFloat assignment."; CPPUNIT_FAIL(msg.str()); } - a.clear(); - b.clear(); - c.clear(); } void verifyArithmetic(double v1, double v2) { @@ -176,14 +171,6 @@ class MFloatTest : public CppUnit::TestFixture { "."; CPPUNIT_FAIL(msg.str()); } - //MFloat::toStr(num1); - num1.clear(); - num2.clear(); - tmp.clear(); - sum.clear(); - diff.clear(); - prod.clear(); - diff.clear(); } void basicArithmetic() { From 745ca1208bd92b1c71a6521e71a35d21e045ad6f Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Wed, 8 Jun 2022 18:41:32 +0200 Subject: [PATCH 12/13] default precision set to sizeof(double), ifdef around mpfr in TV, precision not set as global variable --- engine/maths/mfloat.h | 56 +++++++++-- engine/triangulation/dim3/triangulation3.h | 2 +- engine/triangulation/dim3/turaevviro.cpp | 102 +++++++++++++++------ python/dim3/triangulation3.cpp | 2 +- testsuite/maths/mfloat.cpp | 3 +- 5 files changed, 124 insertions(+), 41 deletions(-) diff --git a/engine/maths/mfloat.h b/engine/maths/mfloat.h index a36f18495c..208290e05d 100644 --- a/engine/maths/mfloat.h +++ b/engine/maths/mfloat.h @@ -49,7 +49,8 @@ class MFloat { private: mpfr_t number_; - static unsigned long prec_; + unsigned long prec_=defaultPrec_; + static unsigned long defaultPrec_; public: @@ -63,23 +64,34 @@ class MFloat { explicit MFloat(double value); + explicit MFloat(double value, unsigned long prec); + explicit MFloat(unsigned long value); - static void setPrec(unsigned long prec) { - prec_ = prec; + explicit MFloat(unsigned long value, unsigned long prec); + + void setPrec(unsigned long prec) { + if (prec_ != prec) + prec_ = prec; + if (mpfr_get_prec(number_) != prec) + mpfr_set_prec(number_,prec); + } + static void setDefaultPrec(unsigned long prec) { + defaultPrec_ = prec; } static std::string str(const MFloat& value){ return std::to_string(value.getDouble()); } - void set(double value); + void set(double value, unsigned long prec); - void set(unsigned long value); + void set(unsigned long value, unsigned long prec); double getDouble() const; void setPi(); + void setPi(unsigned long prec); ~MFloat(); @@ -135,7 +147,7 @@ class MFloat { }; -unsigned long MFloat::prec_ = 32; +unsigned long MFloat::defaultPrec_ = sizeof(double); inline MFloat::MFloat() { mpfr_init2(number_,prec_); @@ -143,17 +155,20 @@ inline MFloat::MFloat() { } inline MFloat::MFloat(const MFloat& other) { - mpfr_init2(number_,prec_); + mpfr_init2(number_,other.prec_); + prec_ = other.prec_; mpfr_set(number_,other.number_,MPFRRND); } inline MFloat::MFloat(MFloat&& other) noexcept { + prec_ = other.prec_; number_->_mpfr_d = 0; mpfr_swap (number_,other.number_); } inline MFloat& MFloat::operator=(MFloat&& other) { if (this != &other) { + prec_ = other.prec_; mpfr_swap(number_,other.number_); } return *this; @@ -164,16 +179,30 @@ inline MFloat::MFloat(double value) { mpfr_set_d(number_,value,MPFRRND); } +inline MFloat::MFloat(double value, unsigned long prec) { + mpfr_init2(number_,prec); + this->setPrec(prec); + mpfr_set_d(number_,value,MPFRRND); +} + inline MFloat::MFloat(unsigned long value) { mpfr_init2(number_,prec_); mpfr_set_ui(number_,value,MPFRRND); } -inline void MFloat::set(double value) { +inline MFloat::MFloat(unsigned long value, unsigned long prec) { + mpfr_init2(number_,prec); + this->setPrec(prec); + mpfr_set_ui(number_,value,MPFRRND); +} + +inline void MFloat::set(double value, unsigned long prec) { + this->setPrec(prec); mpfr_set_d(number_,value,MPFRRND); } -inline void MFloat::set(unsigned long value) { +inline void MFloat::set(unsigned long value, unsigned long prec) { + this->setPrec(prec); mpfr_set_ui(number_,value,MPFRRND); } @@ -182,8 +211,9 @@ inline double MFloat::getDouble() const{ } inline MFloat::~MFloat() { - if(0 != number_->_mpfr_d) + if(0 != number_->_mpfr_d) { mpfr_clear(number_); + } } inline double MFloat::extractDouble(MFloat&& rhs) { @@ -194,6 +224,11 @@ inline void MFloat::setPi() { mpfr_const_pi(number_,MPFRRND); } +inline void MFloat::setPi(unsigned long prec) { + this->setPrec(prec); + mpfr_const_pi(number_,MPFRRND); +} + bool MFloat::operator== ( const MFloat& rhs ) const { return mpfr_equal_p(number_,rhs.number_)!=0; } @@ -203,6 +238,7 @@ inline bool MFloat::operator != (const MFloat& rhs) const { } inline MFloat& MFloat::operator = (const MFloat& other) { + this->setPrec(other.prec_); mpfr_set(number_,other.number_,MPFRRND); return *this; } diff --git a/engine/triangulation/dim3/triangulation3.h b/engine/triangulation/dim3/triangulation3.h index 10d13b0144..1de4aed50a 100644 --- a/engine/triangulation/dim3/triangulation3.h +++ b/engine/triangulation/dim3/triangulation3.h @@ -854,7 +854,7 @@ class Triangulation<3> : public detail::TriangulationBase<3> { * @see allCalculatedTuraevViro */ double turaevViroApprox(unsigned long r, unsigned long whichRoot = 1, - Algorithm alg = ALG_DEFAULT, unsigned long prec = 64) const; + Algorithm alg = ALG_DEFAULT, unsigned long prec = sizeof(double)) const; /** * Returns the cache of all Turaev-Viro state sum invariants that * have been calculated for this 3-manifold. diff --git a/engine/triangulation/dim3/turaevviro.cpp b/engine/triangulation/dim3/turaevviro.cpp index 9fb95063bb..0586816fe9 100644 --- a/engine/triangulation/dim3/turaevviro.cpp +++ b/engine/triangulation/dim3/turaevviro.cpp @@ -68,6 +68,7 @@ namespace { template struct TuraevViroDetails; + template <> struct TuraevViroDetails { typedef Cyclotomic TVType; @@ -88,6 +89,7 @@ namespace { } }; +#ifdef __REGINA_MFLOAT template <> struct TuraevViroDetails { typedef MFloat TVType; @@ -98,6 +100,7 @@ namespace { return zero; } }; +#endif /** * Allows calculation of [n]! for arbitrary n. @@ -124,6 +127,7 @@ namespace { * Requires r >= 3. */ BracketFactorial(unsigned long r, unsigned long whichRoot); + BracketFactorial(unsigned long r, unsigned long whichRoot, unsigned long prec); /** * Clean up memory. @@ -213,28 +217,29 @@ namespace { inv_[i] = inv_[i - 1] / bracket_[i]; } } - + +#ifdef __REGINA_MFLOAT template <> BracketFactorial::BracketFactorial( - unsigned long r, unsigned long whichRoot) : + unsigned long r, unsigned long whichRoot, unsigned long prec) : bracket_(new TVResult[r]), fact_(new TVResult[r]), inv_(new TVResult[r]) { TVResult angle,sinangle,siniangle,facti,invi; - - angle.setPi(); + + angle.setPi(prec); angle *= whichRoot; angle /= r; unsigned long one = 1; - bracket_[0].set(one); - bracket_[1].set(one); - fact_[0].set(one); - fact_[1].set(one); - inv_[0].set(one); - inv_[1].set(one); + bracket_[0].set(one,prec); + bracket_[1].set(one,prec); + fact_[0].set(one,prec); + fact_[1].set(one,prec); + inv_[0].set(one,prec); + inv_[1].set(one,prec); sinangle = angle; sinangle.sin(); @@ -257,6 +262,7 @@ namespace { siniangle = angle; } } +#endif /** * Represents the initial data as described in Section 7 of Turaev @@ -267,8 +273,8 @@ namespace { using TVType = typename TuraevViroDetails::TVType; using TVResult = typename TuraevViroDetails::TVResult; - unsigned long r, whichRoot; - /**< The Turaev-Viro parameters. */ + unsigned long r, whichRoot, prec; + /**< The Turaev-Viro parameters. Prec is used only for multiPrecision. */ bool halfField; BracketFactorial fact; /**< The cached values [n]!. */ @@ -276,12 +282,17 @@ namespace { /**< The vertex-based contribution to the Turaev-Viro invariant; this is the inverse square of the distinguished value w. */ + /** + * Constructor, prec is for multiPrecision. + */ InitialData(unsigned long newR, unsigned long newWhichRoot); + InitialData(unsigned long newR, unsigned long newWhichRoot, unsigned long prec); static void negate(TVType& x); void initZero(TVType& x) const; void initOne(TVType& x) const; + void initLongInt(TVType& x, unsigned long i) const; /** * Determines whether (i/2, j/2, k/2) is an admissible triple. @@ -385,7 +396,8 @@ namespace { unsigned long colour2, unsigned long colour3, unsigned long colour4, unsigned long colour5, TVType& ans) const { - TVType tmp(halfField ? r : 2 * r); + TVType tmp; + initLongInt(tmp,halfField ? r : 2 * r); tetContrib(colour0, colour1, colour3, colour5, colour4, colour2, tmp); ans *= tmp; @@ -459,16 +471,18 @@ namespace { double tmp = sin(M_PI * whichRoot / r); vertexContrib = 2.0 * tmp * tmp / r; } - + +#ifdef __REGINA_MFLOAT template <> InitialData::InitialData( - unsigned long newR, unsigned long newWhichRoot) : + unsigned long newR, unsigned long newWhichRoot, unsigned long prec) : r(newR), whichRoot(newWhichRoot), + prec(prec), halfField(r % 2 != 0 && whichRoot % 2 == 0), - fact(r, whichRoot) { + fact(r, whichRoot,prec) { TVResult tmp; - tmp.setPi(); + tmp.setPi(prec); tmp *= whichRoot; tmp /= r; tmp.sin(); @@ -478,7 +492,7 @@ namespace { tmp /= r; vertexContrib = tmp; } - +#endif template <> inline void InitialData::negate(InitialData::TVType& x) { @@ -490,10 +504,12 @@ namespace { x = -x; } +#ifdef __REGINA_MFLOAT template <> inline void InitialData::negate(InitialData::TVType& x) { x.negate(); } +#endif template <> inline void InitialData::initZero(InitialData::TVType& x) @@ -507,11 +523,13 @@ namespace { x = 0.0; } +#ifdef __REGINA_MFLOAT template <> inline void InitialData::initZero(InitialData::TVType& x) const { - x.set(0.0); + x.set(0.0,prec); } +#endif template <> inline void InitialData::initOne(InitialData::TVType& x) @@ -526,12 +544,35 @@ namespace { x = 1.0; } +#ifdef __REGINA_MFLOAT template <> inline void InitialData::initOne(InitialData::TVType& x) const { - x.set(1.0); + x.set(1.0,prec); + } +#endif + + template <> + inline void InitialData::initLongInt(InitialData::TVType& x, unsigned long i) + const { + x.init(halfField ? r : 2 * r); + x[0] = i; } + template <> + inline void InitialData::initLongInt(InitialData::TVType& x, unsigned long i) + const { + x = (double) i; + } + +#ifdef __REGINA_MFLOAT + template <> + inline void InitialData::initLongInt(InitialData::TVType& x, unsigned long i) + const { + x.set(i,prec); + } +#endif + template typename InitialData::TVType turaevViroBacktrack( const Triangulation<3>& tri, @@ -621,8 +662,10 @@ namespace { std::fill(colour, colour + nEdges, 0); long curr = 0; - TVType valColour(init.halfField ? init.r : 2 * init.r); - TVType tmpTVType(init.halfField ? init.r : 2 * init.r); + TVType valColour; + init.initLongInt(valColour,init.halfField ? init.r : 2 * init.r); + TVType tmpTVType; + init.initLongInt(tmpTVType,init.halfField ? init.r : 2 * init.r); bool admissible; const Tetrahedron<3>* tet; const Triangle<3>* triangle; @@ -800,7 +843,8 @@ namespace { std::fill(colour, colour + nEdges, 0); long curr = 0; - TVType valColour(init.halfField ? init.r : 2 * init.r); + TVType valColour; + init.initLongInt(valColour,init.halfField ? init.r : 2 * init.r); bool admissible; long index1, index2; const Tetrahedron<3>* tet; @@ -1444,8 +1488,9 @@ double Triangulation<3>::turaevViroApprox(unsigned long r, return 0; if (gcd(r, whichRoot) > 1) return 0; - - if (prec<65) { +#ifdef __REGINA_MFLOAT + if (prec<=sizeof(double)) { +#endif // Set up our initial data. InitialData init(r, whichRoot); @@ -1476,10 +1521,11 @@ double Triangulation<3>::turaevViroApprox(unsigned long r, } */ return ans.real(); +#ifdef __REGINA_MFLOAT } else { // Set up our initial data. - MFloat::setPrec(prec); - InitialData init(r, whichRoot); + InitialData init(r, whichRoot,prec); + InitialData::TVType ans; switch (alg) { @@ -1496,7 +1542,7 @@ double Triangulation<3>::turaevViroApprox(unsigned long r, return ans.getDouble(); } - +#endif } Cyclotomic Triangulation<3>::turaevViro(unsigned long r, bool parity, diff --git a/python/dim3/triangulation3.cpp b/python/dim3/triangulation3.cpp index f5ec1552c3..697db25c35 100644 --- a/python/dim3/triangulation3.cpp +++ b/python/dim3/triangulation3.cpp @@ -274,7 +274,7 @@ void addTriangulation3(pybind11::module_& m) { .def("turaevViroApprox", &Triangulation<3>::turaevViroApprox, pybind11::arg(), pybind11::arg("whichRoot") = 1, pybind11::arg("alg") = regina::ALG_DEFAULT, - pybind11::arg("prec") = 64) + pybind11::arg("prec") = sizeof(double)) .def("allCalculatedTuraevViro", &Triangulation<3>::allCalculatedTuraevViro) .def("longitude", &Triangulation<3>::longitude, diff --git a/testsuite/maths/mfloat.cpp b/testsuite/maths/mfloat.cpp index 29b59195fb..bbf1bf5ff5 100644 --- a/testsuite/maths/mfloat.cpp +++ b/testsuite/maths/mfloat.cpp @@ -29,6 +29,7 @@ * MA 02110-1301, USA. * * * **************************************************************************/ +#ifdef __REGINA_MFLOAT #include #include @@ -181,4 +182,4 @@ class MFloatTest : public CppUnit::TestFixture { } }; - +#endif From f729baaedb600554eb9312c94f480a2fb90ade80 Mon Sep 17 00:00:00 2001 From: Owen Rouille Date: Wed, 8 Jun 2022 19:12:50 +0200 Subject: [PATCH 13/13] sizeof * 8 --- engine/triangulation/dim3/triangulation3.h | 2 +- python/dim3/triangulation3.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/engine/triangulation/dim3/triangulation3.h b/engine/triangulation/dim3/triangulation3.h index 1de4aed50a..56fcb439f8 100644 --- a/engine/triangulation/dim3/triangulation3.h +++ b/engine/triangulation/dim3/triangulation3.h @@ -854,7 +854,7 @@ class Triangulation<3> : public detail::TriangulationBase<3> { * @see allCalculatedTuraevViro */ double turaevViroApprox(unsigned long r, unsigned long whichRoot = 1, - Algorithm alg = ALG_DEFAULT, unsigned long prec = sizeof(double)) const; + Algorithm alg = ALG_DEFAULT, unsigned long prec = sizeof(double)*8) const; /** * Returns the cache of all Turaev-Viro state sum invariants that * have been calculated for this 3-manifold. diff --git a/python/dim3/triangulation3.cpp b/python/dim3/triangulation3.cpp index 697db25c35..57106d9991 100644 --- a/python/dim3/triangulation3.cpp +++ b/python/dim3/triangulation3.cpp @@ -274,7 +274,7 @@ void addTriangulation3(pybind11::module_& m) { .def("turaevViroApprox", &Triangulation<3>::turaevViroApprox, pybind11::arg(), pybind11::arg("whichRoot") = 1, pybind11::arg("alg") = regina::ALG_DEFAULT, - pybind11::arg("prec") = sizeof(double)) + pybind11::arg("prec") = sizeof(double)*8) .def("allCalculatedTuraevViro", &Triangulation<3>::allCalculatedTuraevViro) .def("longitude", &Triangulation<3>::longitude,