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})
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
diff --git a/engine/maths/mfloat.h b/engine/maths/mfloat.h
new file mode 100644
index 0000000000..208290e05d
--- /dev/null
+++ b/engine/maths/mfloat.h
@@ -0,0 +1,327 @@
+
+/**************************************************************************
+ *                                                                        *
+ *  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_;
+        unsigned long prec_=defaultPrec_;
+        static unsigned long defaultPrec_;
+
+    public:
+
+        MFloat();
+
+        MFloat(const MFloat& other);
+
+        MFloat(MFloat&& other) noexcept;
+
+        MFloat& operator=(MFloat&& other);
+
+        explicit MFloat(double value);
+
+        explicit MFloat(double value, unsigned long prec);
+
+        explicit MFloat(unsigned long value);
+
+        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, unsigned long prec);
+
+        void set(unsigned long value, unsigned long prec);
+
+        double getDouble() const;
+
+        void setPi();
+        void setPi(unsigned long prec);
+
+        ~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::defaultPrec_ = sizeof(double);
+
+inline MFloat::MFloat() {
+    mpfr_init2(number_,prec_);
+    mpfr_set_d(number_,0,MPFRRND);
+}
+
+inline MFloat::MFloat(const MFloat& other) {
+    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;
+}
+
+inline MFloat::MFloat(double value) {
+    mpfr_init2(number_,prec_);
+    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 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, unsigned long prec) {
+    this->setPrec(prec);
+    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);
+}
+
+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;
+}
+
+inline bool MFloat::operator != (const MFloat& rhs) const {
+    return mpfr_equal_p(number_,rhs.number_)==0;
+}
+
+inline MFloat& MFloat::operator = (const MFloat& other) {
+    this->setPrec(other.prec_);
+    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
diff --git a/engine/triangulation/dim3/triangulation3.h b/engine/triangulation/dim3/triangulation3.h
index 3c74f7f6d5..56fcb439f8 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 = 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/engine/triangulation/dim3/turaevviro.cpp b/engine/triangulation/dim3/turaevviro.cpp
index 910cfe68fb..0586816fe9 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 <gmpxx.h>
 #include <map>
 
+
 // #define TV_BACKTRACK_DUMP_COLOURINGS
 // #define TV_IGNORE_CACHE
 
@@ -62,13 +64,15 @@
 namespace regina {
 
 namespace {
-    template <bool exact>
+	enum calcType {exact, floating, multiPrecision};
+    template <calcType method>
     struct TuraevViroDetails;
 
+
     template <>
-    struct TuraevViroDetails<true> {
-        using TVType = Cyclotomic;
-        using TVResult = Cyclotomic;
+    struct TuraevViroDetails<exact> {
+        typedef Cyclotomic TVType;
+        typedef Cyclotomic TVResult;
 
         static TVResult zero() {
             return Cyclotomic(1);
@@ -76,24 +80,37 @@ namespace {
     };
 
     template <>
-    struct TuraevViroDetails<false> {
-        using TVType = std::complex<double>;
-        using TVResult = double;
+    struct TuraevViroDetails<floating> {
+        typedef std::complex<double> TVType;
+        typedef double TVResult;
 
         static TVResult zero() {
             return 0;
         }
     };
 
+#ifdef __REGINA_MFLOAT
+    template <>
+    struct TuraevViroDetails<multiPrecision> {
+        typedef MFloat TVType;
+        typedef MFloat TVResult;
+
+        static TVResult zero() {
+            MFloat zero = MFloat(0.0);
+            return zero;
+        }
+    };
+#endif
+
     /**
      * Allows calculation of [n]! for arbitrary n.
      * Values are cached as they are calculated.
      */
-    template <bool exact>
+    template <calcType method>
     class BracketFactorial {
         public:
-            using TVType = typename TuraevViroDetails<exact>::TVType;
-            using TVResult = typename TuraevViroDetails<exact>::TVResult;
+            using TVType = typename TuraevViroDetails<method>::TVType;
+            using TVResult = typename TuraevViroDetails<method>::TVResult;
 
         private:
             TVResult* bracket_;
@@ -110,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.
@@ -146,7 +164,7 @@ namespace {
     };
 
     template <>
-    BracketFactorial<true>::BracketFactorial(
+    BracketFactorial<exact>::BracketFactorial(
             unsigned long r, unsigned long whichRoot) :
             bracket_(new TVResult[r]),
             fact_(new TVResult[r]),
@@ -185,7 +203,7 @@ namespace {
     }
 
     template <>
-    BracketFactorial<false>::BracketFactorial(
+    BracketFactorial<floating>::BracketFactorial(
             unsigned long r, unsigned long whichRoot) :
             bracket_(new TVResult[r]),
             fact_(new TVResult[r]),
@@ -200,30 +218,81 @@ namespace {
         }
     }
 
+#ifdef __REGINA_MFLOAT
+    template <>
+    BracketFactorial<multiPrecision>::BracketFactorial(
+            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(prec);
+        angle *= whichRoot;
+        angle /= r;
+
+        unsigned long one = 1;
+
+        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();
+
+        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;
+        }
+    }
+#endif
+
     /**
      * Represents the initial data as described in Section 7 of Turaev
      * and Viro's paper.
      */
-    template <bool exact>
+    template <calcType method>
     struct InitialData {
-        using TVType = typename TuraevViroDetails<exact>::TVType;
-        using TVResult = typename TuraevViroDetails<exact>::TVResult;
+        using TVType = typename TuraevViroDetails<method>::TVType;
+        using TVResult = typename TuraevViroDetails<method>::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<exact> fact;
+        BracketFactorial<method> fact;
             /**< The cached values [n]!. */
         TVType vertexContrib;
             /**< 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.
@@ -327,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;
@@ -373,7 +443,7 @@ namespace {
     };
 
     template <>
-    InitialData<true>::InitialData(
+    InitialData<exact>::InitialData(
             unsigned long newR, unsigned long newWhichRoot) :
             r(newR),
             whichRoot(newWhichRoot),
@@ -392,7 +462,7 @@ namespace {
     }
 
     template <>
-    InitialData<false>::InitialData(
+    InitialData<floating>::InitialData(
             unsigned long newR, unsigned long newWhichRoot) :
             r(newR),
             whichRoot(newWhichRoot),
@@ -402,47 +472,113 @@ namespace {
         vertexContrib = 2.0 * tmp * tmp / r;
     }
 
+#ifdef __REGINA_MFLOAT
     template <>
-    inline void InitialData<true>::negate(InitialData<true>::TVType& x) {
+    InitialData<multiPrecision>::InitialData(
+            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,prec) {
+        TVResult tmp;
+        tmp.setPi(prec);
+        tmp *= whichRoot;
+        tmp /= r;
+        tmp.sin();
+        tmp *= tmp;
+        unsigned long two = 2;
+        tmp *= two;
+        tmp /= r;
+        vertexContrib = tmp;
+    }
+#endif
+
+    template <>
+    inline void InitialData<exact>::negate(InitialData<exact>::TVType& x) {
         x.negate();
     }
 
     template <>
-    inline void InitialData<false>::negate(InitialData<false>::TVType& x) {
+    inline void InitialData<floating>::negate(InitialData<floating>::TVType& x) {
         x = -x;
     }
 
+#ifdef __REGINA_MFLOAT
     template <>
-    inline void InitialData<true>::initZero(InitialData<true>::TVType& x)
+    inline void InitialData<multiPrecision>::negate(InitialData<multiPrecision>::TVType& x) {
+        x.negate();
+    }
+#endif
+
+    template <>
+    inline void InitialData<exact>::initZero(InitialData<exact>::TVType& x)
             const {
         x.init(halfField ? r : 2 * r);
     }
 
     template <>
-    inline void InitialData<false>::initZero(InitialData<false>::TVType& x)
+    inline void InitialData<floating>::initZero(InitialData<floating>::TVType& x)
             const {
         x = 0.0;
     }
 
+#ifdef __REGINA_MFLOAT
     template <>
-    inline void InitialData<true>::initOne(InitialData<true>::TVType& x)
+    inline void InitialData<multiPrecision>::initZero(InitialData<multiPrecision>::TVType& x)
+            const {
+        x.set(0.0,prec);
+    }
+#endif
+
+    template <>
+    inline void InitialData<exact>::initOne(InitialData<exact>::TVType& x)
             const {
         x.init(halfField ? r : 2 * r);
         x[0] = 1;
     }
 
     template <>
-    inline void InitialData<false>::initOne(InitialData<false>::TVType& x)
+    inline void InitialData<floating>::initOne(InitialData<floating>::TVType& x)
             const {
         x = 1.0;
     }
 
-    template <bool exact>
-    typename InitialData<exact>::TVType turaevViroBacktrack(
+#ifdef __REGINA_MFLOAT
+    template <>
+    inline void InitialData<multiPrecision>::initOne(InitialData<multiPrecision>::TVType& x)
+            const {
+        x.set(1.0,prec);
+    }
+#endif
+
+    template <>
+    inline void InitialData<exact>::initLongInt(InitialData<exact>::TVType& x, unsigned long i)
+            const {
+        x.init(halfField ? r : 2 * r);
+        x[0] = i;
+    }
+
+    template <>
+    inline void InitialData<floating>::initLongInt(InitialData<floating>::TVType& x, unsigned long i)
+            const {
+        x = (double) i;
+    }
+
+#ifdef __REGINA_MFLOAT
+    template <>
+    inline void InitialData<multiPrecision>::initLongInt(InitialData<multiPrecision>::TVType& x, unsigned long i)
+            const {
+        x.set(i,prec);
+    }
+#endif
+
+    template <calcType method>
+    typename InitialData<method>::TVType turaevViroBacktrack(
             const Triangulation<3>& tri,
-            const InitialData<exact>& init,
+            const InitialData<method>& init,
             ProgressTracker* tracker) {
-        using TVType = typename InitialData<exact>::TVType;
+        using TVType = typename InitialData<method>::TVType;
 
         if (tracker)
             tracker->newStage("Enumerating colourings");
@@ -526,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;
@@ -648,14 +786,16 @@ namespace {
         delete[] tetDone;
         delete[] tetDoneStart;
 
+
         delete[] edgeCache;
         delete[] triangleCache;
-        delete[] tetCache;
+        delete[] tetCache;;
 
         if (tracker) {
             delete[] coeff;
-            if (tracker->isCancelled())
-                return TuraevViroDetails<exact>::zero();
+            if (tracker->isCancelled()) {
+                return TuraevViroDetails<method>::zero();
+            }
         }
 
         // Compute the vertex contributions separately, since these are
@@ -666,12 +806,12 @@ namespace {
         return ans;
     }
 
-    template <bool exact>
-    typename InitialData<exact>::TVType turaevViroNaive(
+    template <calcType method>
+    typename InitialData<method>::TVType turaevViroNaive(
             const Triangulation<3>& tri,
-            const InitialData<exact>& init,
+            const InitialData<method>& init,
             ProgressTracker* tracker) {
-        using TVType = typename InitialData<exact>::TVType;
+        using TVType = typename InitialData<method>::TVType;
 
         if (tracker)
             tracker->newStage("Enumerating colourings");
@@ -703,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;
@@ -809,8 +950,9 @@ namespace {
 
         if (tracker) {
             delete[] coeff;
-            if (tracker->isCancelled())
-                return TuraevViroDetails<exact>::zero();
+            if (tracker->isCancelled()) {
+                return TuraevViroDetails<method>::zero();
+            }
         }
 
         // Compute the vertex contributions separately, since these are
@@ -818,15 +960,16 @@ namespace {
         for (i = 0; i < tri.countVertices(); i++)
             ans *= init.vertexContrib;
 
+
         return ans;
     }
 
-    template <bool exact>
-    typename InitialData<exact>::TVType turaevViroTreewidth(
+    template <calcType method>
+    typename InitialData<method>::TVType turaevViroTreewidth(
             const Triangulation<3>& tri,
-            InitialData<exact>& init,
+            InitialData<method>& init,
             ProgressTracker* tracker) {
-        using TVType = typename InitialData<exact>::TVType;
+        using TVType = typename InitialData<method>::TVType;
 
         // Progress:
         // - weight of forget/join bag processing is 0.9
@@ -850,7 +993,7 @@ namespace {
 
         if (tracker) {
             if (tracker->isCancelled())
-                return TuraevViroDetails<exact>::zero();
+                return TuraevViroDetails<method>::zero();
             tracker->newStage("Analysing bags", 0.01);
         }
 
@@ -1050,9 +1193,9 @@ 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;
-
+                            }
                             ++level;
                             while (level < 6 && choiceType[level] != 0)
                                 ++level;
@@ -1195,7 +1338,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<int> seq(nEdges);
                             for (i = 0; i < nEdges; ++i)
@@ -1208,8 +1352,9 @@ namespace {
 
                             existingSoln = partial[index]->try_emplace(
                                 std::move(seq), std::move(val));
-                            if (! existingSoln.second)
+                            if (! existingSoln.second) {
                                 existingSoln.first->second += val;
+                            }
                         }
                 }
 
@@ -1219,7 +1364,7 @@ namespace {
                 delete partial[child->index()];
                 partial[child->index()] = nullptr;
 
-                delete partial[sibling->index()];
+                delete partial[child->index()];
                 partial[sibling->index()] = nullptr;
             }
 
@@ -1242,10 +1387,11 @@ 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)
+                    delete partial[i];
             delete[] partial;
 
-            return TuraevViroDetails<exact>::zero();
+            return TuraevViroDetails<method>::zero();
         }
 
         // We made it to the end.
@@ -1263,11 +1409,11 @@ namespace {
         return ans;
     }
 
-    template <bool exact>
-    typename InitialData<exact>::TVType turaevViroPolytope(
+    template <calcType method>
+    typename InitialData<method>::TVType turaevViroPolytope(
             const Triangulation<3>& tri,
-            InitialData<exact>& init) {
-        using TVType = typename InitialData<exact>::TVType;
+            InitialData<method>& init) {
+        using TVType = typename InitialData<method>::TVType;
 
         std::vector<std::vector<mpz_class> > input;
         unsigned long nTri = tri.countTriangles();
@@ -1334,7 +1480,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;
@@ -1342,37 +1488,61 @@ double Triangulation<3>::turaevViroApprox(unsigned long r,
         return 0;
     if (gcd(r, whichRoot) > 1)
         return 0;
-
-    // Set up our initial data.
-    InitialData<false> init(r, whichRoot);
-
-    InitialData<false>::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();
+#ifdef __REGINA_MFLOAT
+    if (prec<=sizeof(double)) {
+#endif
+    	// Set up our initial data.
+		InitialData<floating> init(r, whichRoot);
+
+		InitialData<floating>::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();
+#ifdef __REGINA_MFLOAT
+	} else {
+		// Set up our initial data.
+		InitialData<multiPrecision> init(r, whichRoot,prec);
+
+
+		InitialData<multiPrecision>::TVType ans;
+		switch (alg) {
+		    case ALG_BACKTRACK:
+		        ans = turaevViroBacktrack(*this, init, 0);
+		        break;
+		    case ALG_NAIVE:
+		        ans = turaevViroNaive(*this, init, 0);
+		        break;
+		    default:
+		        ans = turaevViroTreewidth(*this, init, 0);
+		        break;
+		}
+
+		return ans.getDouble();
+	}
+#endif
 }
 
 Cyclotomic Triangulation<3>::turaevViro(unsigned long r, bool parity,
@@ -1398,9 +1568,9 @@ Cyclotomic Triangulation<3>::turaevViro(unsigned long r, bool parity,
 #endif
 
     // Set up our initial data.
-    InitialData<true> init(r, (parity ? 1 : 0));
+    InitialData<exact> init(r, (parity ? 1 : 0));
 
-    InitialData<true>::TVType ans;
+    InitialData<exact>::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..57106d9991 100644
--- a/python/dim3/triangulation3.cpp
+++ b/python/dim3/triangulation3.cpp
@@ -273,7 +273,8 @@ void addTriangulation3(pybind11::module_& m) {
             pybind11::call_guard<pybind11::gil_scoped_release>())
         .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") = sizeof(double)*8)
         .def("allCalculatedTuraevViro",
             &Triangulation<3>::allCalculatedTuraevViro)
         .def("longitude", &Triangulation<3>::longitude,
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..bbf1bf5ff5
--- /dev/null
+++ b/testsuite/maths/mfloat.cpp
@@ -0,0 +1,185 @@
+
+/**************************************************************************
+ *                                                                        *
+ *  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.                                                   *
+ *                                                                        *
+ **************************************************************************/
+#ifdef __REGINA_MFLOAT
+
+#include <sstream>
+#include <cppunit/extensions/HelperMacros.h>
+#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 <typename T>
+        std::string str(T x) {
+            std::ostringstream ans;
+            ans << x;
+            return ans.str();
+        }
+
+        bool inEpsilon(double d1,double d2) {
+			return d1-d2<epsilon && d2-d1 <epsilon;
+		}
+
+        void verifyConstructi(unsigned long val) {
+            MFloat num(val);
+            if (!inEpsilon(num.getDouble(),val)) {
+                std::ostringstream msg;
+                msg << "MFloat " << val <<
+                    " incorrectly assigned from ul.";
+                CPPUNIT_FAIL(msg.str());
+            }
+        }
+
+        void verifyConstructd(double val) {
+            MFloat num(val);
+            if (!inEpsilon(num.getDouble(),val)) {
+                std::ostringstream msg;
+                msg << "MFloat " << val <<
+                    " incorrectly assigned from ul.";
+                CPPUNIT_FAIL(msg.str());
+            }
+        }
+
+        void constructFromInteger() {
+            verifyConstructi(-1);
+            verifyConstructi(0);
+            verifyConstructi(1);
+        }
+
+        void constructFromDouble() {
+            verifyConstructd(-1);
+            verifyConstructd(0);
+            verifyConstructd(1);
+        }
+
+		void verifyAssignments() {
+			MFloat a;
+			a = (double) 1;
+			MFloat b = a;
+			b = (double) 2;
+			MFloat c(double(3));
+			MFloat d(c);
+			d = (double) 4;
+
+
+			if (!inEpsilon(a.getDouble(),1)
+				&& !inEpsilon(b.getDouble(),2)
+				&& !inEpsilon(c.getDouble(),4)
+				&& !inEpsilon(d.getDouble(),4)
+			) {
+                std::ostringstream msg;
+                msg << "Incorrect MFloat assignment.";
+                CPPUNIT_FAIL(msg.str());
+            }
+		}
+
+        void verifyArithmetic(double v1, double v2) {
+            MFloat num1(v1);
+			MFloat num2(v2);
+
+			MFloat sum(v1+v2);
+			MFloat diff(v1-v2);
+			MFloat prod(v1*v2);
+			MFloat div(v1/v2);
+
+			MFloat tmp((double)0);
+			tmp = num1;
+			tmp += num2;
+            if (!inEpsilon(sum.getDouble(),tmp.getDouble())) {
+                std::ostringstream msg;
+                msg << "Incorrect MFloat sum " << sum.getDouble() <<
+                    ".";
+                CPPUNIT_FAIL(msg.str());
+            }
+
+            tmp = num1;
+			tmp -= num2;
+            if (!inEpsilon(diff.getDouble(),tmp.getDouble())) {
+                std::ostringstream msg;
+                msg << "Incorrect MFloat diff " << diff.getDouble() <<
+                    ".";
+                CPPUNIT_FAIL(msg.str());
+            }
+
+            tmp = num1;
+			tmp *= num2;
+            if (!inEpsilon(prod.getDouble(),tmp.getDouble())) {
+                std::ostringstream msg;
+                msg << "Incorrect MFloat prod " << prod.getDouble() <<
+                    ".";
+                CPPUNIT_FAIL(msg.str());
+            }
+
+            tmp = num1;
+			tmp /= num2;
+            if (!inEpsilon(div.getDouble(),tmp.getDouble())) {
+                std::ostringstream msg;
+                msg << "Incorrect MFloat div " << div.getDouble() <<
+                    ".";
+                CPPUNIT_FAIL(msg.str());
+            }
+        }
+
+        void basicArithmetic() {
+			verifyArithmetic(1,1);
+			verifyArithmetic(-1,1);
+			verifyArithmetic(1,-1);
+			verifyArithmetic(0,1);
+		}
+
+};
+#endif