Inclusion of MPFR for Turaev-Viro computation

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")

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)

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})

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

engine/maths/mfloat.h

@@ -0,0 +1,291 @@
+
+/**************************************************************************
+ *                                                                        *
+ *  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 <string.h>
+#include "regina-core.h"
+#include <mpfr.h>
+
+
+namespace regina {
+
+class MFloat {
+
+    private:
+
+        mpfr_t number_;
+        static unsigned long prec_;
+
+    public:
+
+        MFloat();
+
+        MFloat(const MFloat& other);
+
+        MFloat(MFloat&& other) noexcept;
+
+        MFloat& operator=(MFloat&& other);
+
+        explicit MFloat(double value);
+
+        explicit MFloat(unsigned long value);
+
+        static void setPrec(unsigned long prec) {
+            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() const;
+
+        void setPi();
+
+        ~MFloat();
+
+        double static extractDouble(MFloat&& rhs);
+
+        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();
+        }
+
+        friend std::ostream& operator<< ( std::ostream& out, const regina::MFloat& other);
+
+
+};
+
+unsigned long MFloat::prec_ = 32;
+
+inline MFloat::MFloat() {
+    mpfr_init2(number_,prec_);
+    mpfr_set_d(number_,0,MPFRRND);
+}
+
+inline MFloat::MFloat(const MFloat& other) {
+    mpfr_init2(number_,prec_);
+    mpfr_set(number_,other.number_,MPFRRND);
+}
+
+inline MFloat::MFloat(MFloat&& other) noexcept {
+    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_);
+    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) {
+    mpfr_set_d(number_,value,MPFRRND);
+}
+
+inline void MFloat::set(unsigned long value) {
+    mpfr_set_ui(number_,value,MPFRRND);
+}
+
+inline double MFloat::getDouble() const{
+    return mpfr_get_d(number_,MPFRRND);
+}
+
+inline MFloat::~MFloat() {
+    if(0 != number_->_mpfr_d)
+        mpfr_clear(number_);
+}
+
+inline double MFloat::extractDouble(MFloat&& rhs) {
+    return rhs.getDouble();
+}
+
+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);
+};
+
+std::ostream& operator << (std::ostream& out, const MFloat& value) {
+    out << std::to_string(value.getDouble());
+    return out;
+}
+
+} // namespace regina
+
+#endif

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.