diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt
index 61cc4be5c2fe3e428c540c32a1c30068e6586ec5..d5c833b5a22ad80a6d53d60cb49e85669bbd5465 100644
--- a/HySoP/CMakeLists.txt
+++ b/HySoP/CMakeLists.txt
@@ -104,7 +104,7 @@ set(PROJECT_LIBRARY_NAME ${PROJECT_NAME})
 #   - ${PROJECT_NAME}_SOURCE_DIR : where sources (.f and .h and this CMakeLists.txt) are located
 # Note that because of OutOfSourceBuild, binary_dir and source_dir must be different.
 
-set(LANGLIST)
+set(LANGLIST C)
 if(WITH_LIB_CXX)
   set(LANGLIST ${LANGLIST} CXX)
 endif()
@@ -326,8 +326,9 @@ if(WITH_LIB_CXX)
     endif()
   
     if(WITH_MAIN_CXX)
-        list(APPEND cxx_executable_sources "${CXX_MAIN_DIR}/poissonSolver.cpp")
+        list(APPEND cxx_executable_sources "${CXX_MAIN_DIR}/planner.cpp")
         list(APPEND cxx_executable_sources "${CXX_MAIN_DIR}/diffSolver.cpp")
+        list(APPEND cxx_executable_sources "${CXX_MAIN_DIR}/poissonSolver.cpp")
         foreach(cxx_main_source ${cxx_executable_sources})
             get_filename_component(cxx_exec_name "${cxx_main_source}" NAME_WE)
             add_executable(${cxx_exec_name} ${cxx_main_source})
diff --git a/HySoP/src/hysop++/main/diffSolver.cpp b/HySoP/src/hysop++/main/diffSolver.cpp
index 5cf98f3c9c4df4200ced9e09888325ac54f53b81..29548ea56da6d48c655a43ed58baeef64ff58146 100644
--- a/HySoP/src/hysop++/main/diffSolver.cpp
+++ b/HySoP/src/hysop++/main/diffSolver.cpp
@@ -24,18 +24,22 @@ static constexpr std::pair<fft::Extension,fft::Extension> pext[nExtensionsPair]
         std::make_pair(ext[1],ext[1]), //odd-odd
         std::make_pair(ext[0],ext[0]), //none-none
 };
-    
-static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+
+#ifdef HAS_QUADMATHS
+    static constexpr __float128  freqs[6] = { 1.0Q, 1.0Q, 0.75Q, 0.75Q, 0.50Q, 0.50Q };
+#else
+    static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#endif
 
 template <typename T>
 std::function<T(T)> func(std::size_t k) {
     switch(k) {
-        case 0: return [=](T x) {return cos(T(freqs[0])*x);};
-        case 1: return [=](T x) {return sin(T(freqs[1])*x);};
-        case 2: return [=](T x) {return cos(T(freqs[2])*x);};
-        case 3: return [=](T x) {return sin(T(freqs[3])*x);};
-        case 4: return [=](T x) {return cos(T(freqs[4])*x);};
-        case 5: return [=](T x) {return sin(T(freqs[5])*x);};
+        case 0: return [=](T x) {return std::cos(T(freqs[0])*x);};
+        case 1: return [=](T x) {return std::sin(T(freqs[1])*x);};
+        case 2: return [=](T x) {return std::cos(T(freqs[2])*x);};
+        case 3: return [=](T x) {return std::sin(T(freqs[3])*x);};
+        case 4: return [=](T x) {return std::cos(T(freqs[4])*x);};
+        case 5: return [=](T x) {return std::sin(T(freqs[5])*x);};
         default: return[=](T x) { return T(1); };
     }
 }
@@ -219,15 +223,15 @@ int main(int argc, const char *argv[]) {
     std::cout << std::endl;
 #endif
 
-//#ifdef FFTW_HAS_FFTW3Q
-    //std::cout << "== TEST 1D - __float128       ==" << std::endl;
-    //test<__float128,1,false>(5);
-    //std::cout << "== TEST 2D - __float128       ==" << std::endl;
-    //test<__float128,2,false>(3);
-    //std::cout << "== TEST 3D - __float128       ==" << std::endl;
-    //test<__float128,3,false>(1);
-    //std::cout << std::endl;
-//#endif
+#ifdef FFTW_HAS_FFTW3Q
+    std::cout << "== TEST 1D - __float128       ==" << std::endl;
+    test<__float128,1,false>(5);
+    std::cout << "== TEST 2D - __float128       ==" << std::endl;
+    test<__float128,2,false>(3);
+    std::cout << "== TEST 3D - __float128       ==" << std::endl;
+    test<__float128,3,false>(1);
+    std::cout << std::endl;
+#endif
 
     return EXIT_SUCCESS;
 }
diff --git a/HySoP/src/hysop++/main/poissonSolver.cpp b/HySoP/src/hysop++/main/poissonSolver.cpp
index 6806499bd90fa7b4d06684f9d3ee6f867513a5d1..dfb829a49d914ac168ecf7a068c9170e29bf3394 100644
--- a/HySoP/src/hysop++/main/poissonSolver.cpp
+++ b/HySoP/src/hysop++/main/poissonSolver.cpp
@@ -24,17 +24,21 @@ static constexpr std::pair<domain::Boundary,domain::Boundary> pbds[nBoundaryPair
         std::make_pair(bds[0],bds[0]), //none-none
 };
     
-static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#ifdef HAS_QUADMATHS
+    static constexpr __float128  freqs[6] = { 1.0Q, 1.0Q, 0.75Q, 0.75Q, 0.50Q, 0.50Q };
+#else
+    static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#endif
 
 template <typename T>
 std::function<T(T)> func(std::size_t k) {
     switch(k) {
-        case 0: return [=](T x) {return cos(T(freqs[0])*x);};
-        case 1: return [=](T x) {return sin(T(freqs[1])*x);};
-        case 2: return [=](T x) {return cos(T(freqs[2])*x);};
-        case 3: return [=](T x) {return sin(T(freqs[3])*x);};
-        case 4: return [=](T x) {return cos(T(freqs[4])*x);};
-        case 5: return [=](T x) {return sin(T(freqs[5])*x);};
+        case 0: return [=](T x) {return std::cos(T(freqs[0])*x);};
+        case 1: return [=](T x) {return std::sin(T(freqs[1])*x);};
+        case 2: return [=](T x) {return std::cos(T(freqs[2])*x);};
+        case 3: return [=](T x) {return std::sin(T(freqs[3])*x);};
+        case 4: return [=](T x) {return std::cos(T(freqs[4])*x);};
+        case 5: return [=](T x) {return std::sin(T(freqs[5])*x);};
         default: return[=](T x) { return T(1); };
     }
 }
@@ -166,14 +170,17 @@ void test(bool includePeriodicBds=false) {
     
 int main(int argc, const char *argv[]) {
 
+#ifdef FFTW_HAS_FFTW3F
     std::cout << "== TEST 1D - float       ==";
     test<float,1,true>();
     std::cout << "== TEST 2D - float       ==";
     test<float,2,true>();
     std::cout << "== TEST 3D - float       ==";
-    test<float,3,true>();
+    test<float,3,false>();
     std::cout << std::endl;
+#endif
     
+#ifdef FFTW_HAS_FFTW3D
     std::cout << "== TEST 1D - double      ==";
     test<double,1,false>();
     std::cout << "== TEST 2D - double      ==";
@@ -181,7 +188,9 @@ int main(int argc, const char *argv[]) {
     std::cout << "== TEST 3D - double      ==";
     test<double,3,false>();
     std::cout << std::endl;
+#endif
     
+#ifdef FFTW_HAS_FFTW3L
     std::cout << "== TEST 1D - long double ==";
     test<long double,1,false>();
     std::cout << "== TEST 2D - long double ==";
@@ -189,6 +198,16 @@ int main(int argc, const char *argv[]) {
     std::cout << "== TEST 3D - long double ==";
     test<long double,3,false>();
     std::cout << std::endl;
-
+#endif
+
+#ifdef FFTW_HAS_FFTW3Q
+    std::cout << "== TEST 1D - __float128 ==";
+    test<__float128,1,false>();
+    std::cout << "== TEST 2D - __float128 ==";
+    test<__float128,2,false>();
+    std::cout << "== TEST 3D - __float128 ==";
+    test<__float128,3,false>();
+    std::cout << std::endl;
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/HySoP/src/hysop++/src/domain/domain.h b/HySoP/src/hysop++/src/domain/domain.h
index ff228c3d1b738fb183fbebd6bbf47ca9e71a65f9..9a9b0d27c189e77bf36545b3c0e1492bd662cf5d 100644
--- a/HySoP/src/hysop++/src/domain/domain.h
+++ b/HySoP/src/hysop++/src/domain/domain.h
@@ -26,12 +26,12 @@ namespace hysop {
                 public:
                     Domain() :
                         m_shape{0}, m_dataShape{0}, m_leftDataOffset{0},
-                        m_domainConfig(), m_domainSize{0.0L}, m_spaceStep{0.0L}, m_data() {}
+                        m_domainConfig(), m_domainSize{0}, m_spaceStep{0}, m_data() {}
 
                     Domain(const Shape<Dim>& p_shape, const DomainConfiguration<Dim>& p_domainConfig, const DomainSize& p_domainSize) :
                         m_shape{0}, m_dataShape{0}, m_leftDataOffset{0},
                         m_domainConfig(p_domainConfig), m_domainSize(p_domainSize), 
-                        m_spaceStep{0.0L}, m_data() {
+                        m_spaceStep{0}, m_data() {
                             this->reshape(p_shape);
                         }
     
@@ -124,7 +124,7 @@ namespace hysop {
         template <typename Functor, typename... Args>
             Domain<T,Dim>& Domain<T,Dim>::apply(const Functor& f,  Args&&... fargs) {
                 hysop::Index<Dim> idx(m_dataShape);
-                std::array<T,Dim> X{0.0L};
+                std::array<T,Dim> X{0};
                 T* data = m_data.origin();
                 for (std::size_t k = 0; k < m_data.num_elements(); k++) {
                     data[k] = static_cast<T>(f(X, std::forward<Args>(fargs)...));
diff --git a/HySoP/src/hysop++/src/fft/transform.h b/HySoP/src/hysop++/src/fft/transform.h
index 06cc243339bc5dce32f79754a28d5d1a9caf3740..c300483a247472025b79a2509265090107908c8e 100644
--- a/HySoP/src/hysop++/src/fft/transform.h
+++ b/HySoP/src/hysop++/src/fft/transform.h
@@ -56,7 +56,7 @@ namespace hysop {
                 }
 
                 template <typename T>
-                std::complex<T> omega(std::size_t k, std::size_t N, T L = T(1.0L), bool lastDim=false) const {
+                std::complex<T> omega(std::size_t k, std::size_t N, T L = T(1), bool lastDim=false) const {
                     using namespace hysop::constants;
                     switch(kind) {
                         case(FFTW_FORWARD):
diff --git a/HySoP/src/hysop++/src/maths/quad_maths.h b/HySoP/src/hysop++/src/maths/quad_maths.h
index 9991df6a91e929831453170e8bde2cf487fdc107..a411d804162045cb911d18c49faf7a5ac2dd18ed 100644
--- a/HySoP/src/hysop++/src/maths/quad_maths.h
+++ b/HySoP/src/hysop++/src/maths/quad_maths.h
@@ -5,12 +5,13 @@
 #define HYSOP_QUAD_MATHS_H
 
 #include <cfloat>
+#include <cmath>
 #include <limits>
+#include <quadmath.h>
+
 #include <iostream>
 #include <iomanip>
 
-#include <quadmath.h>
-
 /* missing gcc defines */
 #define FLT128_RADIX FLT_RADIX
 #define FLT128_HAS_DENORM      true
@@ -63,6 +64,7 @@ namespace std {
         static constexpr float_round_style round_style = round_to_nearest;
     };
 
+    inline int fpclassify(__float128 arg) { return std::fpclassify(static_cast<long double>(arg)); }
 
     inline __float128 abs(__float128 x) { return cabsq(__complex128{x,0.0Q}); }
     inline __float128 acos(__float128 x) { return acosq(x); }
diff --git a/HySoP/src/hysop++/src/solver/fftPoissonSolver.h b/HySoP/src/hysop++/src/solver/fftPoissonSolver.h
index 4662efb836c38244d262334766cb4543c39174fa..5e9ac41c019be13e10763f40a40c74912001b6d4 100644
--- a/HySoP/src/hysop++/src/solver/fftPoissonSolver.h
+++ b/HySoP/src/hysop++/src/solver/fftPoissonSolver.h
@@ -4,6 +4,7 @@
 
 #include <cmath>
 
+#include "maths/quad_maths.h"
 #include "solver/poissonSolver.h"
 #include "fft/planner.h"
 #include "data/accumulatorIndex.h"
diff --git a/HySoP/src/hysop++/src/utils/constants.h b/HySoP/src/hysop++/src/utils/constants.h
index 344a143e3a4eed241e69f0c08b7b2a1476531167..07ef13458cd91f6df1e9b7ef42f2b9ba2bf44a46 100644
--- a/HySoP/src/hysop++/src/utils/constants.h
+++ b/HySoP/src/hysop++/src/utils/constants.h
@@ -3,13 +3,19 @@
 #define CONSTANTS_H
 
 #include <cmath>
+#include "maths/quad_maths.h"
 #include "utils/types.h"
 
 namespace hysop {
     namespace constants {
-        static constexpr hysop::types::complex I = hysop::types::complex(0.0L,1.0L);
-        static constexpr hysop::types::complex Z = hysop::types::complex(0.0L,0.0L);
+        static constexpr hysop::types::complex I = hysop::types::complex(0,1);
+        static constexpr hysop::types::complex Z = hysop::types::complex(0,0);
+
+#ifdef HAS_QUADMATHS
+        static const __float128 pi = acosq(-1.0Q);
+#else
         static constexpr long double pi = acosl(-1.0L);
+#endif
     }
 }
 
diff --git a/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp
index 7e66b594cc9ddd717bade09d8357624ef5b2a5d4..9144b1fdcdf32470922e522d1f180a52b327c08b 100644
--- a/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp
+++ b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp
@@ -26,17 +26,21 @@ static constexpr std::pair<fft::Extension,fft::Extension> pext[nExtensionsPair]
         std::make_pair(ext[0],ext[0]), //none-none
 };
 
-static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#ifdef HAS_QUADMATHS
+    static constexpr __float128  freqs[6] = { 1.0Q, 1.0Q, 0.75Q, 0.75Q, 0.50Q, 0.50Q };
+#else
+    static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#endif
 
 template <typename T>
 std::function<T(T)> func(std::size_t k) {
     switch(k) {
-        case 0: return [=](T x) {return cos(T(freqs[0])*x);};
-        case 1: return [=](T x) {return sin(T(freqs[1])*x);};
-        case 2: return [=](T x) {return cos(T(freqs[2])*x);};
-        case 3: return [=](T x) {return sin(T(freqs[3])*x);};
-        case 4: return [=](T x) {return cos(T(freqs[4])*x);};
-        case 5: return [=](T x) {return sin(T(freqs[5])*x);};
+        case 0: return [=](T x) {return std::cos(T(freqs[0])*x);};
+        case 1: return [=](T x) {return std::sin(T(freqs[1])*x);};
+        case 2: return [=](T x) {return std::cos(T(freqs[2])*x);};
+        case 3: return [=](T x) {return std::sin(T(freqs[3])*x);};
+        case 4: return [=](T x) {return std::cos(T(freqs[4])*x);};
+        case 5: return [=](T x) {return std::sin(T(freqs[5])*x);};
         default: return[=](T x) { return T(1); };
     }
 }
@@ -189,6 +193,7 @@ void test(std::size_t p_maxOrder, bool includePeriodicBds=false) {
     }
 }
     
+#ifdef FFTW_HAS_FFTW3F
 TEST_F(DiffSolverTest, FloatDerivatives) {
     std::cout << std::endl;
     std::cout << "== TEST 1D - float       ==" << std::endl;
@@ -198,7 +203,9 @@ TEST_F(DiffSolverTest, FloatDerivatives) {
     std::cout << "== TEST 3D - float       ==" << std::endl;
     test<float,3,false>(1);
 }
+#endif
 
+#ifdef FFTW_HAS_FFTW3D
 TEST_F(DiffSolverTest, DoubleDerivatives) {
     std::cout << std::endl;
     std::cout << "== TEST 1D - double       ==" << std::endl;
@@ -208,7 +215,9 @@ TEST_F(DiffSolverTest, DoubleDerivatives) {
     std::cout << "== TEST 3D - double       ==" << std::endl;
     test<double,3,false>(1);
 }
+#endif
     
+#ifdef FFTW_HAS_FFTW3L
 TEST_F(DiffSolverTest, LongDoubleDerivatives) {
     std::cout << std::endl;
     std::cout << "== TEST 1D - long double       ==" << std::endl;
@@ -218,4 +227,17 @@ TEST_F(DiffSolverTest, LongDoubleDerivatives) {
     std::cout << "== TEST 3D - long double       ==" << std::endl;
     test<long double,3,false>(1);
 }
+#endif
+
+#ifdef FFTW_HAS_FFTW3Q
+TEST_F(DiffSolverTest, QuadFloatDerivatives) {
+    std::cout << std::endl;
+    std::cout << "== TEST 1D - __float128       ==" << std::endl;
+    test<__float128,1,false>(5);
+    std::cout << "== TEST 2D - __float128       ==" << std::endl;
+    test<__float128,2,false>(3);
+    std::cout << "== TEST 3D - __float128       ==" << std::endl;
+    test<__float128,3,false>(1);
+}
+#endif
 
diff --git a/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp
index e144e167fb729e80d1b2f6b8ca96cec2712a575f..9711a4cbac573c3a286dbadcb85bc3d4b7304c90 100644
--- a/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp
+++ b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp
@@ -26,17 +26,21 @@ static constexpr std::pair<domain::Boundary,domain::Boundary> pbds[nBoundaryPair
         std::make_pair(bds[0],bds[0]), //none-none
 };
     
-static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#ifdef HAS_QUADMATHS
+    static constexpr __float128  freqs[6] = { 1.0Q, 1.0Q, 0.75Q, 0.75Q, 0.50Q, 0.50Q };
+#else
+    static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L };
+#endif
 
 template <typename T>
 std::function<T(T)> func(std::size_t k) {
     switch(k) {
-        case 0: return [=](T x) {return cos(T(freqs[0])*x);};
-        case 1: return [=](T x) {return sin(T(freqs[1])*x);};
-        case 2: return [=](T x) {return cos(T(freqs[2])*x);};
-        case 3: return [=](T x) {return sin(T(freqs[3])*x);};
-        case 4: return [=](T x) {return cos(T(freqs[4])*x);};
-        case 5: return [=](T x) {return sin(T(freqs[5])*x);};
+        case 0: return [=](T x) {return std::cos(T(freqs[0])*x);};
+        case 1: return [=](T x) {return std::sin(T(freqs[1])*x);};
+        case 2: return [=](T x) {return std::cos(T(freqs[2])*x);};
+        case 3: return [=](T x) {return std::sin(T(freqs[3])*x);};
+        case 4: return [=](T x) {return std::cos(T(freqs[4])*x);};
+        case 5: return [=](T x) {return std::sin(T(freqs[5])*x);};
         default: return[=](T x) { return T(1); };
     }
 }
@@ -154,7 +158,8 @@ void test(bool includePeriodicBds=false) {
         std::cout << " ~= " <<  std::fixed << std::setprecision(0) << meanDists << " eps" << std::endl; 
 }
 
-    
+   
+#ifdef FFTW_HAS_FFTW3F
 TEST_F(PoissonSolverTest, FloatPoissonSolver) {
     std::cout << std::endl;
     std::cout << "== TEST 1D - float       ==";
@@ -164,7 +169,9 @@ TEST_F(PoissonSolverTest, FloatPoissonSolver) {
     std::cout << "== TEST 3D - float       ==";
     test<float,3,false>();
 }
+#endif
 
+#ifdef FFTW_HAS_FFTW3D
 TEST_F(PoissonSolverTest, DoublePoissonSolver) {
     std::cout << std::endl;
     std::cout << "== TEST 1D - double      ==";
@@ -174,7 +181,9 @@ TEST_F(PoissonSolverTest, DoublePoissonSolver) {
     std::cout << "== TEST 3D - double      ==";
     test<double,3,false>();
 }
+#endif
     
+#ifdef FFTW_HAS_FFTW3L
 TEST_F(PoissonSolverTest, LongDoublePoissonSolver) {
     std::cout << std::endl;
     std::cout << "== TEST 1D - long double ==";
@@ -184,3 +193,16 @@ TEST_F(PoissonSolverTest, LongDoublePoissonSolver) {
     std::cout << "== TEST 3D - long double ==";
     test<long double,3,false>();
 }
+#endif
+
+#ifdef FFTW_HAS_FFTW3Q
+TEST_F(PoissonSolverTest, QuadFloatPoissonSolver) {
+    std::cout << std::endl;
+    std::cout << "== TEST 1D - __float128 ==";
+    test<__float128,1,false>();
+    std::cout << "== TEST 2D - __float128 ==";
+    test<__float128,2,false>();
+    std::cout << "== TEST 3D - __float128 ==";
+    test<__float128,3,false>();
+}
+#endif
diff --git a/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt b/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt
index 19cbe1e7480bc0b20e7cf62c96f01cffbd4d98dc..ae4b35dd2c7c115b8e73a278ec049d0bd4336748 100644
--- a/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt
+++ b/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt
@@ -6,8 +6,9 @@ get_filename_component(test_name ${CMAKE_CURRENT_SOURCE_DIR} NAME)
 add_definitions(${CXX_EXTRA_DEFINES})
 add_executable(${test_name} ${SRCS})
 add_dependencies(${test_name} ${HYSOP_CXX_LIBRARY_DEP})
+
 target_link_libraries(${test_name} ${HYSOP_LIBRARY})
-target_link_libraries(${test_name} ${GTEST_LIBRARIES} ${EXT_LIBRARIES})
+target_link_libraries(${test_name} ${GTEST_LIBRARIES} ${CXX_EXT_LIBS})
 
 add_test("${test_name}" "${test_name}")
 
diff --git a/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp
index 4dd35271a16453de9905eb9e78abb32e8001d411..e4001fa854c5f182df4004faa498fee1c408cb39 100644
--- a/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp
+++ b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp
@@ -8,7 +8,7 @@ void evalTest(std::size_t p_size, std::size_t p_samples);
     
 TEST_F(PolynomialTest, EvalTest1D) {
     const std::size_t order=10;
-    const std::size_t samples=1024;
+    const std::size_t samples=4096;
 
     evalTest<bool,1>(order,samples); 
 
@@ -27,6 +27,9 @@ TEST_F(PolynomialTest, EvalTest1D) {
     evalTest<float,1>(order,samples); 
     evalTest<double,1>(order,samples); 
     evalTest<long double,1>(order,samples); 
+#ifdef HAS_QUADMATHS
+    evalTest<__float128,1>(order,samples); 
+#endif
 
     evalTest<std::size_t,1>(order,samples);
     evalTest<std::ptrdiff_t,1>(order,samples);
@@ -34,29 +37,38 @@ TEST_F(PolynomialTest, EvalTest1D) {
 
 TEST_F(PolynomialTest, EvalTest2D) {
     const std::size_t order=10;
-    const std::size_t samples=128;
+    const std::size_t samples=32;
 
     evalTest<float,2>(order,samples); 
     evalTest<double,2>(order,samples); 
     evalTest<long double,2>(order,samples); 
+#ifdef HAS_QUADMATHS
+    evalTest<__float128,2>(order,samples); 
+#endif
 }
 
 TEST_F(PolynomialTest, EvalTest3D) {
     const std::size_t order=10;
-    const std::size_t samples=16;
+    const std::size_t samples=4;
 
     evalTest<float,3>(order,samples); 
     evalTest<double,3>(order,samples); 
     evalTest<long double,3>(order,samples); 
+#ifdef HAS_QUADMATHS
+    evalTest<__float128,3>(order,samples); 
+#endif
 }
 
 TEST_F(PolynomialTest, EvalTest4D) {
     const std::size_t order=10;
-    const std::size_t samples=4;
+    const std::size_t samples=2;
 
     evalTest<float,4>(order,samples); 
     evalTest<double,4>(order,samples); 
     evalTest<long double,4>(order,samples); 
+#ifdef HAS_QUADMATHS
+    evalTest<__float128,4>(order,samples); 
+#endif
 }
 
 template <typename T, std::size_t Dim>