Some weird problem with quadrature related to types

Seems to be fixed now, not sure what the fix was

include/teqp/math/quadrature.hpp

@@ -12,92 +12,53 @@ namespace teqp{
  
  More coefficients here if needed: https://pomax.github.io/bezierinfo/legendre-gauss.html
 */
-template<int N, typename T, typename Function, typename Lim>
-inline auto quad(const Function& F, const Lim& a, const Lim& b){
-    std::common_type_t<Lim, T> summer = 0.0;
-    if constexpr(N==3){
-        static const std::valarray<double> x = {0, sqrt(3.0/5), -sqrt(3.0/5) };
-        static const std::valarray<double> w = {8/9, 5/9, 5/9};
-        for (auto i = 0; i < x.size(); ++i){
-            T arg = (b-a)/2.0*x[i] + (a+b)/2.0;
-            summer += F(arg);
+template<int N, typename T>
+inline auto quad(const std::function<T(double)>& F, const double& a, const double& b){
+    
+    // Locations x in [-1,1] where the function is to be evaluated, and the corresponding weight
+    static const std::map<int, std::tuple<std::vector<double>, std::vector<double>>> xw = {
+        {
+            3,
+            {{ 0, sqrt(3.0/5), -sqrt(3.0/5) },
+            {8.0/9.0, 5.0/9.0, 5.0/9.0}}
+        },
+        {
+            4,
+            {{sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)), -sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)), sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5)), -sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5))},
+            {(18.0+sqrt(30.0))/36.0, (18.0+sqrt(30.0))/36.0, (18.0-sqrt(30.0))/36.0, (18.0-sqrt(30.0))/36.0}}
+        },
+        {
+            5,
+            {{0, 1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), 1.0/3.0*sqrt(5+2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5+2*sqrt(10.0/7.0))},
+            {128/225.0, (322.0+13*sqrt(70))/900.0, (322.0+13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0}}
+        },
+        {
+            7,
+            {{0.0000000000000000,0.4058451513773972,-0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585},
+            { 0.4179591836734694,0.3818300505051189,0.3818300505051189,0.2797053914892766,0.2797053914892766,0.1294849661688697,0.1294849661688697}}
+        },
+        {
+            15,
+            {{0.0000000000000000,-0.2011940939974345,0.2011940939974345,-0.3941513470775634,0.3941513470775634,-0.5709721726085388,0.5709721726085388,-0.7244177313601701,0.7244177313601701,-0.8482065834104272,0.8482065834104272,-0.9372733924007060,0.9372733924007060,-0.9879925180204854,0.9879925180204854},
+            {0.2025782419255613,0.1984314853271116,0.1984314853271116,0.1861610000155622,0.1861610000155622,0.1662692058169939,0.1662692058169939,0.1395706779261543,0.1395706779261543,0.1071592204671719,0.1071592204671719,0.0703660474881081,0.0703660474881081,0.0307532419961173,0.0307532419961173}}
+        },
+        {
+            30,
+            {{-0.0514718425553177,0.0514718425553177,-0.1538699136085835,0.1538699136085835,-0.2546369261678899,0.2546369261678899,-0.3527047255308781,0.3527047255308781,-0.4470337695380892,0.4470337695380892,-0.5366241481420199,0.5366241481420199,-0.6205261829892429,0.6205261829892429,-0.6978504947933158,0.6978504947933158,-0.7677774321048262,0.7677774321048262,-0.8295657623827684,0.8295657623827684,-0.8825605357920527,0.8825605357920527,-0.9262000474292743,0.9262000474292743,-0.9600218649683075,0.9600218649683075,-0.9836681232797472,0.9836681232797472,-0.9968934840746495,0.9968934840746495
+            },
+            {0.1028526528935588,0.1028526528935588,0.1017623897484055,0.1017623897484055,0.0995934205867953,0.0995934205867953,0.0963687371746443,0.0963687371746443,0.0921225222377861,0.0921225222377861,0.0868997872010830,0.0868997872010830,0.0807558952294202,0.0807558952294202,0.0737559747377052,0.0737559747377052,0.0659742298821805,0.0659742298821805,0.0574931562176191,0.0574931562176191,0.0484026728305941,0.0484026728305941,0.0387991925696271,0.0387991925696271,0.0287847078833234,0.0287847078833234,0.0184664683110910,0.0184664683110910,0.0079681924961666,0.0079681924961666}}
         }
-        return (b-a)/2.0*summer;
-    }
-    else if constexpr(N==4){
-        static const std::valarray<double> x = {
-             sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)),
-            -sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)),
-             sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5)),
-            -sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5))
-        };
-        static const std::valarray<double> w = {
-            (18.0+sqrt(30.0))/36.0,
-            (18.0+sqrt(30.0))/36.0,
-            (18.0-sqrt(30.0))/36.0,
-            (18.0-sqrt(30.0))/36.0
-        };
-        for (auto i = 0; i < x.size(); ++i){
-            T arg = (b-a)/2.0*x[i] + (a+b)/2.0;
-            summer += w[i]*F(arg);
-        }
-        return (b-a)/2.0*summer;
-    }
-    else if constexpr(N==5){
-        static const std::vector<double> x = {0,
-             1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)),
-            -1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)),
-             1.0/3.0*sqrt(5+2*sqrt(10.0/7.0)),
-            -1.0/3.0*sqrt(5+2*sqrt(10.0/7.0))
-        };
-        static const std::vector<double> w = {
-            128/225.0,
-            (322.0+13*sqrt(70))/900.0,
-            (322.0+13*sqrt(70))/900.0,
-            (322.0-13*sqrt(70))/900.0,
-            (322.0-13*sqrt(70))/900.0
-        };
-        for (auto i = 0; i < x.size(); ++i){
-            T arg = (b-a)/2.0*x[i] + (a+b)/2.0;
-            summer += w[i]*F(arg);
-        }
-        return (b-a)/2.0*summer;
-    }
-    else if constexpr(N==7){
-        static const std::vector<double> x = {0.0000000000000000,0.4058451513773972,-0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585
-        };
-        static const std::vector<double> w = { 0.4179591836734694,0.3818300505051189,0.3818300505051189,0.2797053914892766,0.2797053914892766,0.1294849661688697,0.1294849661688697
-        };
-        for (auto i = 0; i < x.size(); ++i){
-            T arg = (b-a)/2.0*x[i] + (a+b)/2.0;
-            summer += w[i]*F(arg);
-        }
-        return (b-a)/2.0*summer;
-    }
-    else if constexpr(N==15){
-        static const std::vector<double> x = {
-            0.0000000000000000,-0.2011940939974345,0.2011940939974345,-0.3941513470775634,0.3941513470775634,-0.5709721726085388,0.5709721726085388,-0.7244177313601701,0.7244177313601701,-0.8482065834104272,0.8482065834104272,-0.9372733924007060,0.9372733924007060,-0.9879925180204854,0.9879925180204854
-        };
-        static const std::vector<double> w = { 0.2025782419255613,0.1984314853271116,0.1984314853271116,0.1861610000155622,0.1861610000155622,0.1662692058169939,0.1662692058169939,0.1395706779261543,0.1395706779261543,0.1071592204671719,0.1071592204671719,0.0703660474881081,0.0703660474881081,0.0307532419961173,0.0307532419961173
-        };
-        for (auto i = 0; i < x.size(); ++i){
-            T arg = (b-a)/2.0*x[i] + (a+b)/2.0;
-            summer += w[i]*F(arg);
-        }
-        return (b-a)/2.0*summer;
-    }
-    else if constexpr(N==30){
-        static const std::vector<double> x = {
-            -0.0514718425553177,0.0514718425553177,-0.1538699136085835,0.1538699136085835,-0.2546369261678899,0.2546369261678899,-0.3527047255308781,0.3527047255308781,-0.4470337695380892,0.4470337695380892,-0.5366241481420199,0.5366241481420199,-0.6205261829892429,0.6205261829892429,-0.6978504947933158,0.6978504947933158,-0.7677774321048262,0.7677774321048262,-0.8295657623827684,0.8295657623827684,-0.8825605357920527,0.8825605357920527,-0.9262000474292743,0.9262000474292743,-0.9600218649683075,0.9600218649683075,-0.9836681232797472,0.9836681232797472,-0.9968934840746495,0.9968934840746495
-        };
-        static const std::vector<double> w = { 0.1028526528935588,0.1028526528935588,0.1017623897484055,0.1017623897484055,0.0995934205867953,0.0995934205867953,0.0963687371746443,0.0963687371746443,0.0921225222377861,0.0921225222377861,0.0868997872010830,0.0868997872010830,0.0807558952294202,0.0807558952294202,0.0737559747377052,0.0737559747377052,0.0659742298821805,0.0659742298821805,0.0574931562176191,0.0574931562176191,0.0484026728305941,0.0484026728305941,0.0387991925696271,0.0387991925696271,0.0287847078833234,0.0287847078833234,0.0184664683110910,0.0184664683110910,0.0079681924961666,0.0079681924961666
-        };
-        for (auto i = 0; i < x.size(); ++i){
-            T arg = (b-a)/2.0*x[i] + (a+b)/2.0;
-            summer += w[i]*F(arg);
-        }
-        return (b-a)/2.0*summer;
+    };
+    
+    T summer = 0.0;
+    // Lookup the coefficients without making a copy or re-allocating
+    const std::tuple<std::vector<double>, std::vector<double>> &pair = xw.at(N);
+    const std::vector<double> &x = std::get<0>(pair);
+    const std::vector<double> &w = std::get<1>(pair);
+    for (auto i = 0; i < N; ++i){
+        double arg = (b-a)/2.0*x[i] + (a+b)/2.0;
+        summer += w[i]*F(arg);
     }
+    return (b-a)/2.0*summer;
 }
-
 }