diff --git a/DiFfRG/Mathematica/DiFfRG/CodeTools/MakeKernel.m b/DiFfRG/Mathematica/DiFfRG/CodeTools/MakeKernel.m
index 40f6e76a9..690a54dd9 100644
--- a/DiFfRG/Mathematica/DiFfRG/CodeTools/MakeKernel.m
+++ b/DiFfRG/Mathematica/DiFfRG/CodeTools/MakeKernel.m
@@ -12,6 +12,7 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in
 These are prepended to the respective methods of the integration kernel, allowing one to e.g. define specific angles one needs for the flow code.
 
 The options \"KernelReturnTransform\" and \"ConstantReturnTransform\" (default Identity) accept a Mathematica function applied to the optimized expression before code generation, which lets you wrap the return value, e.g. \"KernelReturnTransform\" -> Re renders the kernel return as real(...).
+The option \"ComputeType\" (default \"double\") is the value type the kernel is evaluated and integrated in: \"double\", \"float\", \"DiFfRG::complex<double>\" or \"DiFfRG::complex<float>\". With a float type the kernel, its literals and the integrator run in single precision, while map()/get() still write double results.
 The option \"KernelTraits\" declares integrator traits on the emitted kernel class, as a list of names or name -> Boolean rules, e.g. \"KernelTraits\" -> {\"matsubara_finite_extent\"} or {\"matsubara_split\" -> True}. Each becomes a `static constexpr bool <name> = <value>;` member. \"MatsubaraEven\" stays a separate option because it carries a symbolic evenness check that a generic mechanism cannot.";
 
 MakeKernel::Invalid = "The given arguments are invalid. See MakeKernel::usage";
@@ -30,6 +31,8 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in
 
 MakeKernel::notEven = "MatsubaraEven requested for kernel \"`1`\" but it is not even in \"`2`\"; emitting the standard kernel (the integrator keeps the explicit kernel(+f0)+kernel(-f0) form). This is expected when the loop contains fermionic dressings evaluated at f0-shifted arguments.";
 
+MakeKernel::ctypedeprecated = "The option \"ctype\" is deprecated; use \"ComputeType\" -> `1` instead.";
+
 MakeKernel::InvalidTrait = "KernelTraits entry `1` is not a C++ identifier optionally followed by -> True|False.";
 
 Begin["`Private`"]
@@ -69,13 +72,44 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in
 
 KernelSpecQ[spec_Association] :=
     Module[{validKeys, validKeyTypes},
-        validKeys = CheckKey[spec, "Name", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"] && CheckKey[spec, "Integrator", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"] && CheckKey[spec, "d", IntegerQ[#] && # >= 0&, "Must be an Integer >= 0"] && CheckKey[spec, "AD", BooleanQ, "Must be a Boolean"] && CheckKey[spec, "Device", MemberQ[{"Threads", "TBB", "GPU"}, #]&, "Must be Threads, TBB or GPU."] && CheckKey[spec, "Type", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"];
+        validKeys = CheckKey[spec, "Name", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"] && CheckKey[spec, "Integrator", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"] && CheckKey[spec, "d", IntegerQ[#] && # >= 0&, "Must be an Integer >= 0"] && CheckKey[spec, "AD", BooleanQ, "Must be a Boolean"] && CheckKey[spec, "Device", MemberQ[{"Threads", "TBB", "GPU"}, #]&, "Must be Threads, TBB or GPU."] && CheckKey[spec, "Type", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"] && CheckKey[spec, "ComputeType", StringQ[#] && StringLength[#] > 0&, "Cannot be empty"];
         Return[validKeys];
     ];
 
 GetStandardKernelDefinitions[] :=
     $StandardKernelDefinitions
 
+(* "ctype" is the old name of "ComputeType"; an explicit ctype still wins so old scripts keep working. *)
+resolveComputeType[spec_Association] :=
+    Module[{s = spec},
+        If[s["ctype"] =!= Automatic,
+            Message[MakeKernel::ctypedeprecated, ToString[s["ctype"], InputForm]];
+            s["ComputeType"] = s["ctype"]
+        ];
+        KeyDrop[s, "ctype"]
+    ];
+
+singlePrecisionQ[spec_Association] :=
+    StringContainsQ[spec["ComputeType"], "float"];
+
+(* Real type of the kernel's arguments. *)
+realType[spec_Association] :=
+    If[singlePrecisionQ[spec], "float", "double"];
+
+(* map()/get() hand results back in double even when the integrator runs in float, so the
+   caller's double state vectors need no staging. *)
+resultType[spec_Association] :=
+    StringReplace[spec["ComputeType"], "float" -> "double"];
+
+(* Print the kernel's literals as float (0.5f, complex<float>) in single precision; a double
+   literal would silently promote the whole expression back to double. *)
+SetAttributes[withComputePrecision, HoldRest];
+
+withComputePrecision[spec_Association, body_] :=
+    Block[{FunKit`Private`$codePrecision = If[singlePrecisionQ[spec], "single", FunKit`Private`$codePrecision]},
+        body
+    ];
+
 (* Emit arbitrary integrator traits as class members. The traits are what the finite-T
    integrator branches on (matsubara_even, matsubara_finite_extent, matsubara_split), and a
    `requires { requires K::trait; }` detector reads FALSE for anything it cannot see -- a
@@ -114,7 +148,7 @@ integrator branches on (matsubara_even, matsubara_finite_extent, matsubara_split
 
 (* Internal functions added here with Internal`*::usage *)
 
-Options[MakeKernel] = {"Coordinates" -> {}, "CoordinateArguments" -> {}, "IntegrationVariables" -> {}, "KernelDefinitions" -> $StandardKernelDefinitions, "Regulator" -> "DiFfRG::PolynomialExpRegulator", "RegulatorOpts" -> {"", ""}, "KernelBody" -> "", "KernelReturnType" -> "auto", "KernelReturnTransform" -> Identity, "ConstantBody" -> "", "ConstantReturnType" -> "auto", "ConstantReturnTransform" -> Identity, "Parameters" -> {}, "Name" -> "", "d" -> -1, "Integrator" -> "", "AD" -> False, "ctype" -> "double", "Device" -> "TBB", "Type" -> "double", "SplitKernel" -> False, "SeparateLookups" -> False, "Decorator" -> "static KOKKOS_FUNCTION", "MatsubaraEven" -> False, "KernelTraits" -> {}};
+Options[MakeKernel] = {"Coordinates" -> {}, "CoordinateArguments" -> {}, "IntegrationVariables" -> {}, "KernelDefinitions" -> $StandardKernelDefinitions, "Regulator" -> "DiFfRG::PolynomialExpRegulator", "RegulatorOpts" -> {"", ""}, "KernelBody" -> "", "KernelReturnType" -> "auto", "KernelReturnTransform" -> Identity, "ConstantBody" -> "", "ConstantReturnType" -> "auto", "ConstantReturnTransform" -> Identity, "Parameters" -> {}, "Name" -> "", "d" -> -1, "Integrator" -> "", "AD" -> False, "ComputeType" -> "double", "ctype" -> Automatic, "Device" -> "TBB", "Type" -> "double", "SplitKernel" -> False, "SeparateLookups" -> False, "Decorator" -> "static KOKKOS_FUNCTION", "MatsubaraEven" -> False, "KernelTraits" -> {}};
 
 MakeKernel[__] :=
     (
@@ -128,6 +162,7 @@ integrator branches on (matsubara_even, matsubara_finite_extent, matsubara_split
 MakeKernel[kernelExpr_, constExpr_, OptionsPattern[]] :=
     Module[{expr, const, exec, kernel, constant, kernelClass, kernelHeader, integratorHeader, integratorCpp, integratorTemplateParams, tparams = <|"Name" -> "...t", "Type" -> "auto&&", "Reference" -> False, "Const" -> False|>, kernelDefs = OptionValue["KernelDefinitions"], coordinates = OptionValue["Coordinates"], getArgs = OptionValue["CoordinateArguments"], intVariables = OptionValue["IntegrationVariables"], preArguments, regulator, params, adSpecs, explParamAD, arguments, outputPath, sources, returnType, returnTypePointer, spec, parameters, parametersKernel, matsubaraEvenTrait, kernelTraits},
         spec = Association @@ Thread[Rule @@ {#, OptionValue[MakeKernel, #]}]& @ Keys[Options[MakeKernel]];
+        spec = resolveComputeType[spec];
         If[Not @ KernelSpecQ[spec],
             Message[MakeKernel::InvalidSpec];
             Abort[]
@@ -166,9 +201,9 @@ GENUINELY even in that variable (kernel(+f0) == kernel(-f0)), the integrator
         const = constExpr;
         While[ListQ[const], const = Plus @@ const];
         intVariables = FunKit`Private`prepParam /@ intVariables;
-        intVariables = Map[Append[#, "Type" -> "double"]&, intVariables];
+        intVariables = Map[Append[#, "Type" -> realType[spec]]&, intVariables];
         getArgs = FunKit`Private`prepParam /@ getArgs;
-        getArgs = Map[Append[#, "Type" -> "double"]&, getArgs];
+        getArgs = Map[Append[#, "Type" -> realType[spec]]&, getArgs];
         (********************************************************************)
         (* First, the kernel itself *)
         (********************************************************************)
@@ -188,13 +223,13 @@ GENUINELY even in that variable (kernel(+f0) == kernel(-f0)), the integrator
                 ,
                 spec["Parameters"]
             ];
-        kernel =
+        kernel = withComputePrecision[spec,
             If[TrueQ[OptionValue["SplitKernel"]] || TrueQ[OptionValue["SeparateLookups"]],
                 FunKit`MakeCppFunctionSplit[expr, "Name" -> "kernel", "Return" -> OptionValue["KernelReturnType"], "Suffix" -> "", "Prefix" -> "static KOKKOS_INLINE_FUNCTION", "Decorator" -> OptionValue["Decorator"], "SeparateLookups" -> OptionValue["SeparateLookups"], "Parameters" -> Join[intVariables, getArgs, parametersKernel], "Body" -> StringTemplate["using namespace DiFfRG;using namespace DiFfRG::compute;\n`1`"][OptionValue["KernelBody"]], "ReturnTransform" -> OptionValue["KernelReturnTransform"]]
                 ,
                 FunKit`MakeCppFunction[expr, "Name" -> "kernel", "Return" -> OptionValue["KernelReturnType"], "Suffix" -> "", "Prefix" -> "static KOKKOS_INLINE_FUNCTION", "Parameters" -> Join[intVariables, getArgs, parametersKernel], "Body" -> StringTemplate["using namespace DiFfRG;using namespace DiFfRG::compute;\n`1`"][OptionValue["KernelBody"]], "ReturnTransform" -> OptionValue["KernelReturnTransform"]]
-            ];
-        constant = FunKit`MakeCppFunction[constExpr, "Name" -> "constant", "Return" -> OptionValue["ConstantReturnType"], "Suffix" -> "", "Prefix" -> "static KOKKOS_INLINE_FUNCTION", "Parameters" -> Join[getArgs, parametersKernel], "Body" -> StringTemplate["using namespace DiFfRG;using namespace DiFfRG::compute;\n`1`"][OptionValue["ConstantBody"]], "ReturnTransform" -> OptionValue["ConstantReturnTransform"]];
+            ]];
+        constant = withComputePrecision[spec, FunKit`MakeCppFunction[constExpr, "Name" -> "constant", "Return" -> OptionValue["ConstantReturnType"], "Suffix" -> "", "Prefix" -> "static KOKKOS_INLINE_FUNCTION", "Parameters" -> Join[getArgs, parametersKernel], "Body" -> StringTemplate["using namespace DiFfRG;using namespace DiFfRG::compute;\n`1`"][OptionValue["ConstantBody"]], "ReturnTransform" -> OptionValue["ConstantReturnTransform"]]];
         kernelTraits = kernelTraitMembers[OptionValue["KernelTraits"]];
         kernelClass = FunKit`MakeCppClass["TemplateTypes" -> {"_Regulator"}, "Name" -> OptionValue["Name"] <> "_kernel", "MembersPublic" -> Join[{"using Regulator = _Regulator;"}, matsubaraEvenTrait, kernelTraits, {kernel, constant}], "MembersPrivate" -> kernelDefs];
         kernelHeader = FunKit`MakeCppHeader["Includes" -> {"DiFfRG/physics/interpolation.hh", "DiFfRG/physics/physics.hh"}, "Body" -> {"namespace DiFfRG {", kernelClass, StringTemplate["} using DiFfRG::`1`_kernel;"][spec["Name"]]}];
@@ -222,11 +257,11 @@ GENUINELY even in that variable (kernel(+f0) == kernel(-f0)), the integrator
             ];
         integratorTemplateParams = TemplateParameterGeneration[spec];
         integratorTemplateParams = StringRiffle[integratorTemplateParams, ", "];
-        returnType = spec["ctype"];
+        returnType = resultType[spec];
         returnTypePointer = StringTemplate["`1`*"][returnType];
         adSpecs =
             If[spec["AD"],
-                Map[Merge[{#, <|"IntegratorTemplateParams" -> StringRiffle[TemplateParameterGeneration[spec, #["Replacements"]], ", "], "ReturnType" -> spec["ctype"] /. #["Replacements"], "Params" -> Last @ processParameters[params, #["Replacements"]]|>}, Last]&, $ADSpecializations]
+                Map[Merge[{#, <|"IntegratorTemplateParams" -> StringRiffle[TemplateParameterGeneration[spec, #["Replacements"]], ", "], "ReturnType" -> resultType[spec] /. #["Replacements"], "Params" -> Last @ processParameters[params, #["Replacements"]]|>}, Last]&, $ADSpecializations]
                 ,
                 {}
             ];
diff --git a/DiFfRG/Mathematica/DiFfRG/CodeTools/TemplateParameterGeneration.m b/DiFfRG/Mathematica/DiFfRG/CodeTools/TemplateParameterGeneration.m
index b2c41688f..dfe849624 100644
--- a/DiFfRG/Mathematica/DiFfRG/CodeTools/TemplateParameterGeneration.m
+++ b/DiFfRG/Mathematica/DiFfRG/CodeTools/TemplateParameterGeneration.m
@@ -40,16 +40,16 @@
 
 appendKeyType[templateParameter_List, params_Association, {}] :=
     Module[{},
-        If[KeyExistsQ[params, "ctype"],
-            Append[templateParameter, ToString[params["ctype"]]],
+        If[KeyExistsQ[params, "ComputeType"],
+            Append[templateParameter, ToString[params["ComputeType"]]],
             Append[templateParameter, "double"]
         ]
     ]
 
 appendKeyType[templateParameter_List, params_Association, ADReplacements_] :=
     Module[{},
-        If[KeyExistsQ[params, "ctype"],
-            Append[templateParameter, ToString[params["ctype"]] /. ADReplacements],
+        If[KeyExistsQ[params, "ComputeType"],
+            Append[templateParameter, ToString[params["ComputeType"]] /. ADReplacements],
             Append[templateParameter, "autodiff::real"]
         ]
     ]
diff --git a/DiFfRG/Mathematica/DiFfRG/Tests/MakeKernelSecondOrderADTests.m b/DiFfRG/Mathematica/DiFfRG/Tests/MakeKernelSecondOrderADTests.m
index 9172027b8..5abb7ed7a 100644
--- a/DiFfRG/Mathematica/DiFfRG/Tests/MakeKernelSecondOrderADTests.m
+++ b/DiFfRG/Mathematica/DiFfRG/Tests/MakeKernelSecondOrderADTests.m
@@ -33,7 +33,7 @@
             "Name" -> "pion",
             "Integrator" -> "Integrator_p2",
             "d" -> 3,
-            "ctype" -> "DiFfRG::complex<double>",
+            "ComputeType" -> "DiFfRG::complex<double>",
             "Parameters" -> {
                 <|"Name" -> "k", "Type" -> "double", "Const" -> True, "AD" -> False|>,
                 <|"Name" -> "T", "Type" -> "double", "Const" -> True, "AD" -> False|>,
@@ -80,3 +80,53 @@
     Block[{Print}, Get["FunKit`"]];
     AUMPCHECK[Quiet[generatesSecondOrderComplexAD[], OptionValue::nodef]];
 ];
+
+(* A float ComputeType must leave no double in the kernel (it would promote the arithmetic back),
+   while map()/get() keep double destinations for the caller's state vectors. *)
+generatesFloatKernel[typeOption_] :=
+    Module[{tmp, header, kernelHeader},
+        tmp = FileNameJoin[{AUMPTestTempDirectory[], "generated"}];
+        CreateDirectory[tmp];
+        SetFlowDirectory[tmp <> "/"];
+        CreateDirectory[FileNameJoin[{tmp, "flows", "fl", "src"}], CreateIntermediateDirectories -> True];
+        MakeKernel[
+            0.5 l1 + 2 k + I l1,
+            "IntegrationVariables" -> {"l1"},
+            "Coordinates" -> {"LinearCoordinates1D<float>"},
+            "Name" -> "fl",
+            "Integrator" -> "Integrator_p2",
+            "d" -> 3,
+            typeOption -> "DiFfRG::complex<float>",
+            "Parameters" -> {<|"Name" -> "k", "Type" -> "float", "Const" -> True|>}
+        ];
+        header = Import[FileNameJoin[{tmp, "flows", "fl", "fl.hh"}], "Text"];
+        kernelHeader = Import[FileNameJoin[{tmp, "flows", "fl", "kernel.hh"}], "Text"];
+        containsAll[
+            header,
+            {
+                "Integrator_p2<3, DiFfRG::complex<float>, fl_kernel<Regulator>, DiFfRG::TBB_exec> integrator;",
+                "map(DiFfRG::complex<double>* dest, const LinearCoordinates1D<float>& coordinates",
+                "void get(DiFfRG::complex<double>& dest"
+            }
+        ] &&
+            containsAll[kernelHeader, {"const float& l1", "0.5f", "complex<float>("}] &&
+            StringFreeQ[kernelHeader, "double"]
+    ];
+
+AUMPTestCase["MakeKernel with a float ComputeType emits a pure-float kernel", {"make-kernel", "precision", "funkit", "form"},
+    AUMPAssume[Length @ PacletFind["FunKit"] > 0, "FunKit is not installed"];
+    AUMPAssume[formAvailableQ[], "FORM is not installed"];
+    Needs["FormTracer`"];
+    FormTracer`DefineFormExecutable[formExecutablePath[]];
+    Block[{Print}, Get["FunKit`"]];
+    AUMPCHECK[Quiet[generatesFloatKernel["ComputeType"], OptionValue::nodef]];
+];
+
+AUMPTestCase["MakeKernel still accepts the deprecated ctype option", {"make-kernel", "precision", "funkit", "form"},
+    AUMPAssume[Length @ PacletFind["FunKit"] > 0, "FunKit is not installed"];
+    AUMPAssume[formAvailableQ[], "FORM is not installed"];
+    Needs["FormTracer`"];
+    FormTracer`DefineFormExecutable[formExecutablePath[]];
+    Block[{Print}, Get["FunKit`"]];
+    AUMPCHECK[Quiet[generatesFloatKernel["ctype"], {OptionValue::nodef, MakeKernel::ctypedeprecated}]];
+];
diff --git a/DiFfRG/Mathematica/DiFfRG/Tests/TemplateParameterGenerationTests.m b/DiFfRG/Mathematica/DiFfRG/Tests/TemplateParameterGenerationTests.m
index b79ba5a12..4ff7cc544 100644
--- a/DiFfRG/Mathematica/DiFfRG/Tests/TemplateParameterGenerationTests.m
+++ b/DiFfRG/Mathematica/DiFfRG/Tests/TemplateParameterGenerationTests.m
@@ -3,14 +3,14 @@
 
 AUMPTestCase["TemplateParameterGeneration supports GPU float kernels", {"template-parameters"},
     AUMPCHECKEqual[
-        TemplateParameterGeneration[<|"d" -> 2, "Name" -> "MyKernel", "ctype" -> "float", "Device" -> "GPU"|>],
+        TemplateParameterGeneration[<|"d" -> 2, "Name" -> "MyKernel", "ComputeType" -> "float", "Device" -> "GPU"|>],
         {"2", "float", "MyKernel_kernel<Regulator>", "DiFfRG::GPU_exec"}
     ];
 ];
 
 AUMPTestCase["TemplateParameterGeneration supports TBB double kernels", {"template-parameters"},
     AUMPCHECKEqual[
-        TemplateParameterGeneration[<|"d" -> 3, "Name" -> "Integrator", "ctype" -> "double", "Device" -> "TBB"|>],
+        TemplateParameterGeneration[<|"d" -> 3, "Name" -> "Integrator", "ComputeType" -> "double", "Device" -> "TBB"|>],
         {"3", "double", "Integrator_kernel<Regulator>", "DiFfRG::TBB_exec"}
     ];
 ];
@@ -18,7 +18,7 @@
 AUMPTestCase["TemplateParameterGeneration applies real AD replacements", {"template-parameters", "ad"},
     AUMPCHECKEqual[
         TemplateParameterGeneration[
-            <|"d" -> 3, "Name" -> "Integrator", "ctype" -> "double", "Device" -> "TBB"|>,
+            <|"d" -> 3, "Name" -> "Integrator", "ComputeType" -> "double", "Device" -> "TBB"|>,
             {"double" -> "autodiff::real"}
         ],
         {"3", "autodiff::real", "Integrator_kernel<Regulator>", "DiFfRG::TBB_exec"}
@@ -28,7 +28,7 @@
 AUMPTestCase["TemplateParameterGeneration applies second-order complex AD replacements", {"template-parameters", "ad"},
     AUMPCHECKEqual[
         TemplateParameterGeneration[
-            <|"d" -> 3, "Name" -> "pion", "ctype" -> "DiFfRG::complex<double>", "Device" -> "TBB"|>,
+            <|"d" -> 3, "Name" -> "pion", "ComputeType" -> "DiFfRG::complex<double>", "Device" -> "TBB"|>,
             {"DiFfRG::complex<double>" -> "cxReal<2, double>"}
         ],
         {"3", "cxReal<2, double>", "pion_kernel<Regulator>", "DiFfRG::TBB_exec"}
@@ -37,12 +37,12 @@
 
 AUMPTestCase["TemplateParameterGeneration supports Threads execution", {"template-parameters"},
     AUMPCHECKEqual[
-        TemplateParameterGeneration[<|"d" -> 1, "Name" -> "Test", "ctype" -> "float", "Device" -> "Threads"|>],
+        TemplateParameterGeneration[<|"d" -> 1, "Name" -> "Test", "ComputeType" -> "float", "Device" -> "Threads"|>],
         {"1", "float", "Test_kernel<Regulator>", "DiFfRG::KokkosHost_exec"}
     ];
 ];
 
-AUMPTestCase["TemplateParameterGeneration defaults ctype to double", {"template-parameters", "defaults"},
+AUMPTestCase["TemplateParameterGeneration defaults ComputeType to double", {"template-parameters", "defaults"},
     AUMPCHECKEqual[
         TemplateParameterGeneration[<|"d" -> 4, "Name" -> "DefaultType", "Device" -> "TBB"|>],
         {"4", "double", "DefaultType_kernel<Regulator>", "DiFfRG::TBB_exec"}
@@ -51,7 +51,7 @@
 
 AUMPTestCase["TemplateParameterGeneration defaults Device to TBB", {"template-parameters", "defaults"},
     AUMPCHECKEqual[
-        TemplateParameterGeneration[<|"d" -> 2, "Name" -> "DefaultDevice", "ctype" -> "double"|>],
+        TemplateParameterGeneration[<|"d" -> 2, "Name" -> "DefaultDevice", "ComputeType" -> "double"|>],
         {"2", "double", "DefaultDevice_kernel<Regulator>", "DiFfRG::TBB_exec"}
     ];
 ];
diff --git a/DiFfRG/include/DiFfRG/common/types.hh b/DiFfRG/include/DiFfRG/common/types.hh
index 36b06978c..666e55404 100644
--- a/DiFfRG/include/DiFfRG/common/types.hh
+++ b/DiFfRG/include/DiFfRG/common/types.hh
@@ -75,6 +75,25 @@ namespace DiFfRG
 
     template <typename CT> using ctype = typename internal::_ctype<CT>::value;
 
+    namespace internal
+    {
+      template <typename T> struct _double_precision {
+        using value = T;
+      };
+
+      template <> struct _double_precision<float> {
+        using value = double;
+      };
+
+      template <> struct _double_precision<complex<float>> {
+        using value = complex<double>;
+      };
+    } // namespace internal
+
+    /// The double-precision counterpart of a single-precision type (float -> double,
+    /// complex<float> -> complex<double>); every other type maps to itself.
+    template <typename T> using double_precision = typename internal::_double_precision<T>::value;
+
     template <typename T> inline constexpr bool is_autodiff = internal::_is_autodiff<T>;
   } // namespace get_type
 } // namespace DiFfRG
diff --git a/DiFfRG/include/DiFfRG/physics/integration/finiteT/quadrature_integrator_fT.hh b/DiFfRG/include/DiFfRG/physics/integration/finiteT/quadrature_integrator_fT.hh
index fbfa7c21e..8a092c270 100644
--- a/DiFfRG/include/DiFfRG/physics/integration/finiteT/quadrature_integrator_fT.hh
+++ b/DiFfRG/include/DiFfRG/physics/integration/finiteT/quadrature_integrator_fT.hh
@@ -284,15 +284,25 @@ namespace DiFfRG
       dest = m_result_host();
     }
 
+    /// Single-precision integration handing back a double result.
     template <typename OT, typename... T>
-      requires(!std::is_same_v<OT, NT>)
+      requires internal::is_widened_result<OT, NT>
+    void get(OT &dest, const T &...t) const
+    {
+      NT result;
+      get(result, t...);
+      dest = OT(result);
+    }
+
+    template <typename OT, typename... T>
+      requires(!std::is_same_v<OT, NT> && !internal::is_widened_result<OT, NT>)
     void get(OT &dest, const T &...t) const
     {
       get(space, dest, t...);
     }
 
     template <typename OT, typename... Args>
-      requires(!std::is_same_v<OT, NT>)
+      requires(!std::is_same_v<OT, NT> && !internal::is_widened_result<OT, NT>)
     void get(ExecutionSpace &space, OT &dest, const Args &...t) const
     {
       const auto args = device::make_tuple(t...);
@@ -535,8 +545,9 @@ namespace DiFfRG
       return volume;
     }
 
-    template <typename Coordinates, typename... Args>
-    auto map(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+      requires internal::is_map_result<OT, NT>
+    auto map(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       auto &scheduler = MapScheduler::instance();
 
@@ -544,7 +555,7 @@ namespace DiFfRG
       if (scheduler.active() && scheduler.plan_contains(integrator_id())) MapCompletion::flush();
 
       const MapSlice slice =
-          scheduler.schedule(integrator_id(), dest, sizeof(NT), coordinates.size(), quadrature_volume(),
+          scheduler.schedule(integrator_id(), dest, sizeof(OT), coordinates.size(), quadrature_volume(),
                              /* splittable */ true, map_target<ExecutionSpace>());
 
       if (slice.count == 0) {
@@ -556,10 +567,11 @@ namespace DiFfRG
       return map_dist(dest + slice.offset, SubCoordinates(coordinates, slice.offset, slice.count), args...);
     }
 
-    template <typename Coordinates, typename... Args>
-    auto map_dist(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+    auto map_dist(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       const size_t n = coordinates.size();
+      auto &stage = m_staging.template get<OT>();
 
       if constexpr (std::is_same_v<typename ExecutionSpace::memory_space, CPU_memory>) {
         // Host backend: "device" memory is host memory, so there is nothing to stage. The work is
@@ -573,16 +585,10 @@ namespace DiFfRG
         // One staging buffer per integrator, so a second map() before a flush would clobber the
         // first result. Land the outstanding one first; in the normal call pattern (each flow
         // mapped once per flush interval) this never triggers.
-        if (m_dest_pinned_size > 0 && MapCompletion::has_pending(m_dest_pinned.data())) MapCompletion::flush();
-
-        auto dest_device_view = device_scratch(n);
+        if (stage.pinned_size > 0 && MapCompletion::has_pending(stage.pinned.data())) MapCompletion::flush();
 
-        if (m_dest_pinned_size < n) {
-          m_dest_pinned = Kokkos::View<NT *, PinnedHost_memory>(
-              Kokkos::view_alloc(Kokkos::WithoutInitializing, "MapIntegrators_fT_pinned_view"), n);
-          m_dest_pinned_size = n;
-        }
-        auto pinned_view = Kokkos::View<NT *, PinnedHost_memory>(m_dest_pinned, Kokkos::make_pair(size_t(0), n));
+        auto dest_device_view = stage.device_view(space, n);
+        auto pinned_view = stage.pinned_view(n);
 
         map(space, dest_device_view, coordinates, args...);
 
@@ -594,7 +600,7 @@ namespace DiFfRG
 
         // Caller has promised not to read `dest` until its DeferredMaps scope closes, so leave the
         // result in staging and keep the host running ahead of the device.
-        MapCompletion::record(dest, m_dest_pinned.data(), n * sizeof(NT));
+        MapCompletion::record(dest, stage.pinned.data(), n * sizeof(OT));
         // Original contract outside such a scope: `dest` is valid on return. flush() fences, lands
         // the staged copy and -- under MPI -- exchanges this batch's slices, all of which must
         // happen before the caller looks at `dest`.
@@ -745,20 +751,10 @@ namespace DiFfRG
       grid_size[dim - 1] = matsubara_nodes.size();
     }
 
-    /// Grow-only scratch in the integrator's own execution space, reused across calls.
-    Kokkos::View<NT *, ExecutionSpace> device_scratch(const size_t n)
-    {
-      if (m_dest_device_size < n) {
-        m_dest_device = Kokkos::View<NT *, ExecutionSpace>(Kokkos::view_alloc(space, "MapIntegrators_device_view"), n);
-        m_dest_device_size = n;
-      }
-      return Kokkos::View<NT *, ExecutionSpace>(m_dest_device, Kokkos::make_pair(size_t(0), n));
-    }
-
     /// Run a host-backend map now, or queue it for flush time if a deferral scope is open. See
     /// MapCompletion::record_work. Compiled out in a CUDA-less build.
-    template <typename Coordinates, typename... Args>
-    void run_or_queue_host(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+    void run_or_queue_host(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       if constexpr (internal::has_device_backend) {
         if (MapCompletion::deferral_enabled()) {
@@ -771,14 +767,14 @@ namespace DiFfRG
     }
 
     /// The host-backend body of map_dist(). Queued jobs run one after another, so the shared
-    /// `m_dest_device` scratch is written and drained before the next job touches it.
-    template <typename Coordinates, typename... Args>
-    void run_host(NT *dest, const Coordinates &coordinates, const Args &...args)
+    /// staging device scratch is written and drained before the next job touches it.
+    template <typename OT, typename Coordinates, typename... Args>
+    void run_host(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       const size_t n = coordinates.size();
-      auto dest_device_view = device_scratch(n);
+      auto dest_device_view = m_staging.template get<OT>().device_view(space, n);
       // create unmanaged host view for dest
-      auto dest_view = Kokkos::View<NT *, CPU_memory, Kokkos::MemoryUnmanaged>(dest, n);
+      auto dest_view = Kokkos::View<OT *, CPU_memory, Kokkos::MemoryUnmanaged>(dest, n);
 
       // run the map function
       map(space, dest_device_view, coordinates, args...);
@@ -827,11 +823,7 @@ namespace DiFfRG
     // Keyed on the coordinates' to_string() identity; flat layout [grid_point * cdim + d].
     mutable Kokkos::View<ctype *, typename ExecutionSpace::memory_space> m_positions;
     mutable std::string m_positions_key;
-    mutable Kokkos::View<NT *, ExecutionSpace> m_dest_device;
-    mutable size_t m_dest_device_size = 0;
-    /// Page-locked staging for the device path, so the result copy is genuinely asynchronous.
-    mutable Kokkos::View<NT *, PinnedHost_memory> m_dest_pinned;
-    mutable size_t m_dest_pinned_size = 0;
+    mutable internal::MapStagingSet<NT, ExecutionSpace> m_staging;
     mutable Kokkos::View<NT, typename ExecutionSpace::memory_space> m_result_view;
     mutable typename Kokkos::View<NT, typename ExecutionSpace::memory_space>::host_mirror_type m_result_host;
     mutable bool m_result_views_initialized = false;
@@ -892,8 +884,18 @@ namespace DiFfRG
       dest = KERNEL::constant(t...) + TBBReduction<dim, NT, decltype(functor)>(grid_size, functor);
     }
 
-    template <typename Coordinates, typename... Args>
-    void map(execution_space &, NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename... Args>
+      requires(internal::is_widened_result<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, Args...>)
+    void get(OT &dest, const Args &...t) const
+    {
+      NT result;
+      get(result, t...);
+      dest = OT(result);
+    }
+
+    template <typename OT, typename Coordinates, typename... Args>
+      requires internal::is_map_result<OT, NT>
+    void map(execution_space &, OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       const auto m_args = device::tie(args...);
 
@@ -908,8 +910,9 @@ namespace DiFfRG
       });
     }
 
-    template <typename Coordinates, typename... Args>
-    auto map(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+      requires internal::is_map_result<OT, NT>
+    auto map(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       auto space = execution_space();
       auto &scheduler = MapScheduler::instance();
@@ -917,7 +920,7 @@ namespace DiFfRG
       if (scheduler.active() && scheduler.plan_contains(this->integrator_id())) MapCompletion::flush();
 
       const MapSlice slice =
-          scheduler.schedule(this->integrator_id(), dest, sizeof(NT), coordinates.size(), Base::quadrature_volume(),
+          scheduler.schedule(this->integrator_id(), dest, sizeof(OT), coordinates.size(), Base::quadrature_volume(),
                              /* splittable */ true, map_target<execution_space>());
 
       if (slice.count == 0) {
@@ -936,8 +939,8 @@ namespace DiFfRG
   private:
     /// Run now, or queue for flush time inside a deferral scope -- see MapCompletion::record_work.
     /// Compiled out in a CUDA-less build.
-    template <typename Coordinates, typename... Args>
-    void run_or_queue(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+    void run_or_queue(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       if constexpr (internal::has_device_backend) {
         if (MapCompletion::deferral_enabled()) {
diff --git a/DiFfRG/include/DiFfRG/physics/integration/quadrature_integrator.hh b/DiFfRG/include/DiFfRG/physics/integration/quadrature_integrator.hh
index 6620bbb4f..e22b3aff6 100644
--- a/DiFfRG/include/DiFfRG/physics/integration/quadrature_integrator.hh
+++ b/DiFfRG/include/DiFfRG/physics/integration/quadrature_integrator.hh
@@ -30,6 +30,62 @@ namespace DiFfRG
     c.forward(c.from_linear_index(size_t(0)));
   };
 
+  namespace internal
+  {
+    /// A single-precision integrator may hand its results back in double precision.
+    template <typename OT, typename NT>
+    inline constexpr bool is_widened_result =
+        !std::is_same_v<OT, NT> && std::is_same_v<OT, get_type::double_precision<NT>>;
+
+    /// Element types an integrator with value type NT can map() into.
+    template <typename OT, typename NT>
+    inline constexpr bool is_map_result = std::is_same_v<OT, NT> || is_widened_result<OT, NT>;
+
+    /// Grow-only result buffers of map(): device scratch, plus page-locked staging so the copy back
+    /// is genuinely asynchronous (see MapCompletion).
+    template <typename T, typename ExecutionSpace> struct MapStaging {
+      Kokkos::View<T *, ExecutionSpace> device;
+      size_t device_size = 0;
+      Kokkos::View<T *, PinnedHost_memory> pinned;
+      size_t pinned_size = 0;
+
+      auto device_view(const ExecutionSpace &space, const size_t n)
+      {
+        if (device_size < n) {
+          device = Kokkos::View<T *, ExecutionSpace>(Kokkos::view_alloc(space, "MapIntegrators_device_view"), n);
+          device_size = n;
+        }
+        return Kokkos::View<T *, ExecutionSpace>(device, Kokkos::make_pair(size_t(0), n));
+      }
+
+      auto pinned_view(const size_t n)
+      {
+        if (pinned_size < n) {
+          pinned = Kokkos::View<T *, PinnedHost_memory>(
+              Kokkos::view_alloc(Kokkos::WithoutInitializing, "MapIntegrators_pinned_view"), n);
+          pinned_size = n;
+        }
+        return Kokkos::View<T *, PinnedHost_memory>(pinned, Kokkos::make_pair(size_t(0), n));
+      }
+    };
+
+    /// map() result buffers for value type NT and, when that is single precision, for its double
+    /// counterpart. The kernel writes each result straight into the destination-typed buffer, so a
+    /// float integrator mapping into double converts on the device.
+    template <typename NT, typename ExecutionSpace> struct MapStagingSet {
+      MapStaging<NT, ExecutionSpace> native;
+      MapStaging<get_type::double_precision<NT>, ExecutionSpace> widened;
+
+      template <typename OT> auto &get()
+      {
+        if constexpr (std::is_same_v<OT, NT>)
+          return native;
+        else
+          return widened;
+      }
+    };
+  } // namespace internal
+
   /**
    * @brief This class performs numerical integration over a d-dimensional hypercube using quadrature rules.
    *
@@ -98,8 +154,18 @@ namespace DiFfRG
       dest = m_result_host();
     }
 
+    /// Single-precision integration handing back a double result.
+    template <typename OT, typename... T>
+      requires(internal::is_widened_result<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, T...>)
+    void get(OT &dest, const T &...t) const
+    {
+      NT result;
+      get(result, t...);
+      dest = OT(result);
+    }
+
     template <typename OT, typename... T>
-      requires(!std::is_same_v<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, T...>)
+      requires(!std::is_same_v<OT, NT> && !internal::is_widened_result<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, T...>)
     void get(OT &dest, const T &...t) const
     {
       ExecutionSpace space;
@@ -107,7 +173,7 @@ namespace DiFfRG
     }
 
     template <typename OT, typename... T>
-      requires(!std::is_same_v<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, T...>)
+      requires(!std::is_same_v<OT, NT> && !internal::is_widened_result<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, T...>)
     void get(ExecutionSpace &space, OT &dest, const T &...t) const
     {
       const auto args = device::make_tuple(t...);
@@ -384,8 +450,9 @@ namespace DiFfRG
       return volume;
     }
 
-    template <typename Coordinates, typename... Args>
-    auto map(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+      requires internal::is_map_result<OT, NT>
+    auto map(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       auto &scheduler = MapScheduler::instance();
 
@@ -395,7 +462,7 @@ namespace DiFfRG
       // other ranks perform, leaving the collectives mismatched.
       if (scheduler.active() && scheduler.plan_contains(integrator_id())) MapCompletion::flush();
 
-      const MapSlice slice = scheduler.schedule(integrator_id(), dest, sizeof(NT), coordinates.size(),
+      const MapSlice slice = scheduler.schedule(integrator_id(), dest, sizeof(OT), coordinates.size(),
                                                 quadrature_volume(), /* splittable */ true,
                                                 map_target<ExecutionSpace>());
 
@@ -410,18 +477,12 @@ namespace DiFfRG
       return map_dist(dest + slice.offset, SubCoordinates(coordinates, slice.offset, slice.count), args...);
     }
 
-    template <typename Coordinates, typename... Args>
-    auto map_dist(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+    auto map_dist(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       const size_t n = coordinates.size();
-
-      // Reuse cached device view if large enough, otherwise reallocate (grow-only)
-      if (m_dest_device_size < n) {
-        m_dest_device =
-            Kokkos::View<NT *, ExecutionSpace>(Kokkos::view_alloc(space, "MapIntegrators_device_view"), n);
-        m_dest_device_size = n;
-      }
-      auto dest_device_view = Kokkos::View<NT *, ExecutionSpace>(m_dest_device, Kokkos::make_pair(size_t(0), n));
+      auto &stage = m_staging.template get<OT>();
+      auto dest_device_view = stage.device_view(space, n);
 
       if constexpr (std::is_same_v<typename ExecutionSpace::memory_space, CPU_memory>) {
         // Host backend: "device" memory is host memory, so there is nothing to stage. The work is
@@ -435,15 +496,8 @@ namespace DiFfRG
         // One staging buffer per integrator, so a second map() before a flush would clobber the
         // first result. Land the outstanding one first; in the normal call pattern (each flow
         // mapped once per flush interval) this never triggers.
-        if (m_dest_pinned_size > 0 && MapCompletion::has_pending(m_dest_pinned.data())) MapCompletion::flush();
-
-        if (m_dest_pinned_size < n) {
-          m_dest_pinned = Kokkos::View<NT *, PinnedHost_memory>(
-              Kokkos::view_alloc(Kokkos::WithoutInitializing, "MapIntegrators_pinned_view"), n);
-          m_dest_pinned_size = n;
-        }
-        auto pinned_view =
-            Kokkos::View<NT *, PinnedHost_memory>(m_dest_pinned, Kokkos::make_pair(size_t(0), n));
+        if (stage.pinned_size > 0 && MapCompletion::has_pending(stage.pinned.data())) MapCompletion::flush();
+        auto pinned_view = stage.pinned_view(n);
 
         map(space, dest_device_view, coordinates, args...);
 
@@ -452,7 +506,7 @@ namespace DiFfRG
 
         // Caller has promised not to read `dest` until its DeferredMaps scope closes, so leave the
         // result in staging and keep the host running ahead of the device.
-        MapCompletion::record(dest, m_dest_pinned.data(), n * sizeof(NT));
+        MapCompletion::record(dest, stage.pinned.data(), n * sizeof(OT));
         // Original contract outside such a scope: `dest` is valid on return. flush() fences, lands
         // the staged copy and -- under MPI -- exchanges this batch's slices, all of which must
         // happen before the caller looks at `dest`.
@@ -473,8 +527,8 @@ namespace DiFfRG
      *
      * Compiled out entirely in a CUDA-less build, where there is no device to overlap with.
      */
-    template <typename Coordinates, typename... Args>
-    void run_or_queue_host(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+    void run_or_queue_host(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       if constexpr (internal::has_device_backend) {
         if (MapCompletion::deferral_enabled()) {
@@ -489,18 +543,13 @@ namespace DiFfRG
     }
 
     /// The host-backend body of map_dist(). Queued jobs run one after another, so the shared
-    /// `m_dest_device` scratch is written and drained before the next job touches it.
-    template <typename Coordinates, typename... Args>
-    void run_host(NT *dest, const Coordinates &coordinates, const Args &...args)
+    /// staging device scratch is written and drained before the next job touches it.
+    template <typename OT, typename Coordinates, typename... Args>
+    void run_host(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       const size_t n = coordinates.size();
-      if (m_dest_device_size < n) {
-        m_dest_device =
-            Kokkos::View<NT *, ExecutionSpace>(Kokkos::view_alloc(space, "MapIntegrators_device_view"), n);
-        m_dest_device_size = n;
-      }
-      auto dest_device_view = Kokkos::View<NT *, ExecutionSpace>(m_dest_device, Kokkos::make_pair(size_t(0), n));
-      auto dest_view = Kokkos::View<NT *, CPU_memory, Kokkos::MemoryUnmanaged>(dest, n);
+      auto dest_device_view = m_staging.template get<OT>().device_view(space, n);
+      auto dest_view = Kokkos::View<OT *, CPU_memory, Kokkos::MemoryUnmanaged>(dest, n);
       map(space, dest_device_view, coordinates, args...);
       Kokkos::deep_copy(space, dest_view, dest_device_view);
     }
@@ -524,11 +573,7 @@ namespace DiFfRG
     // Keyed on the coordinates' to_string() identity; flat layout [grid_point * cdim + d].
     mutable Kokkos::View<ctype *, typename ExecutionSpace::memory_space> m_positions;
     mutable std::string m_positions_key;
-    mutable Kokkos::View<NT *, ExecutionSpace> m_dest_device;
-    mutable size_t m_dest_device_size = 0;
-    // Page-locked staging for the result copy; see map_dist() and MapCompletion.
-    mutable Kokkos::View<NT *, PinnedHost_memory> m_dest_pinned;
-    mutable size_t m_dest_pinned_size = 0;
+    mutable internal::MapStagingSet<NT, ExecutionSpace> m_staging;
     mutable Kokkos::View<NT, typename ExecutionSpace::memory_space> m_result_view;
     mutable typename Kokkos::View<NT, typename ExecutionSpace::memory_space>::host_mirror_type m_result_host;
     mutable bool m_result_views_initialized = false;
@@ -581,8 +626,18 @@ namespace DiFfRG
       dest = KERNEL::constant(t...) + TBBReduction<dim, NT, decltype(functor)>(grid_size, functor);
     }
 
-    template <typename Coordinates, typename... Args>
-    void map(execution_space &, NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename... T>
+      requires(internal::is_widened_result<OT, NT> && is_valid_kernel<NT, KERNEL, ctype, dim, T...>)
+    void get(OT &dest, const T &...t) const
+    {
+      NT result;
+      get(result, t...);
+      dest = OT(result);
+    }
+
+    template <typename OT, typename Coordinates, typename... Args>
+      requires internal::is_map_result<OT, NT>
+    void map(execution_space &, OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       const auto m_args = device::tie(args...);
 
@@ -597,8 +652,9 @@ namespace DiFfRG
       });
     }
 
-    template <typename Coordinates, typename... Args>
-    auto map(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+      requires internal::is_map_result<OT, NT>
+    auto map(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       auto space = execution_space();
       auto &scheduler = MapScheduler::instance();
@@ -608,7 +664,7 @@ namespace DiFfRG
       // overload for why that distinction is what keeps the collectives matched.
       if (scheduler.active() && scheduler.plan_contains(this->integrator_id())) MapCompletion::flush();
 
-      const MapSlice slice = scheduler.schedule(this->integrator_id(), dest, sizeof(NT), coordinates.size(),
+      const MapSlice slice = scheduler.schedule(this->integrator_id(), dest, sizeof(OT), coordinates.size(),
                                                 Base::quadrature_volume(), /* splittable */ true,
                                                 map_target<execution_space>());
 
@@ -632,8 +688,8 @@ namespace DiFfRG
     /// Run now, or queue for flush time inside a deferral scope. `tbb::parallel_for` is
     /// synchronous, so running it inline would stall the launch of any device map issued after it;
     /// see MapCompletion::record_work. Compiled out in a CUDA-less build.
-    template <typename Coordinates, typename... Args>
-    void run_or_queue(NT *dest, const Coordinates &coordinates, const Args &...args)
+    template <typename OT, typename Coordinates, typename... Args>
+    void run_or_queue(OT *dest, const Coordinates &coordinates, const Args &...args)
     {
       if constexpr (internal::has_device_backend) {
         if (MapCompletion::deferral_enabled()) {
diff --git a/DiFfRG/include/DiFfRG/physics/regulators.hh b/DiFfRG/include/DiFfRG/physics/regulators.hh
index abe061825..ffa81660a 100644
--- a/DiFfRG/include/DiFfRG/physics/regulators.hh
+++ b/DiFfRG/include/DiFfRG/physics/regulators.hh
@@ -2,6 +2,7 @@
 
 // standard library
 #include <cmath>
+#include <type_traits>
 
 // DiFfRG
 #include <DiFfRG/common/kokkos.hh>
@@ -29,7 +30,8 @@ namespace DiFfRG
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RBdot(const NT1 k2, const NT2 q2)
     {
-      return 2. * k2 * (k2 > q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return T(2) * k2 * (k2 > q2);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RF(const NT1 k2, const NT2 q2)
@@ -41,7 +43,8 @@ namespace DiFfRG
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RFdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
-      return 0.5 * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return T(0.5) * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
     }
 
     // Give an explicit compile error if a call to dq2RB is made
@@ -49,7 +52,7 @@ namespace DiFfRG
   };
 
   struct BosonicRegulatorOpts {
-    static constexpr double b = 2.;
+    static constexpr int b = 2;
   };
   /**
    * @brief Implements one of the standard exponential regulators, i.e.
@@ -68,7 +71,8 @@ namespace DiFfRG
    * @tparam b The exponent in the regulator.
    */
   template <class OPTS = BosonicRegulatorOpts> struct BosonicRegulator {
-    static constexpr double b = OPTS::b;
+    // INTEGER, as in ExponentialRegulator: b is the exponent of powr<b>, a non-type int parameter.
+    static constexpr int b = OPTS::b;
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RB(const NT1 k2, const NT2 q2)
     {
@@ -83,7 +87,8 @@ namespace DiFfRG
       const auto xb = powr<b>(q2 / k2);
       const auto mexp = exp(xb);
       const auto mexpm1 = expm1(xb);
-      return b * mexp * powr<b - 1>(q2 / k2) * (xb - mexpm1) / powr<2>(mexpm1);
+      using T = std::decay_t<decltype(xb)>;
+      return T(b) * mexp * powr<b - 1>(q2 / k2) * (xb - mexpm1) / powr<2>(mexpm1);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RBdot(const NT1 k2, const NT2 q2)
@@ -93,7 +98,8 @@ namespace DiFfRG
       const auto xb = powr<b>(q2 / k2);
       const auto mexp = exp(xb);
       const auto mexpm1 = expm1(xb);
-      return 2. * k2 * xb * (mexpm1 * (1. - b) + b * mexp * xb) / powr<2>(mexpm1);
+      using T = std::decay_t<decltype(xb)>;
+      return T(2) * k2 * xb * (mexpm1 * T(1 - b) + T(b) * mexp * xb) / powr<2>(mexpm1);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RF(const NT1 k2, const NT2 q2)
@@ -106,13 +112,15 @@ namespace DiFfRG
     {
       using Kokkos::sqrt;
       const auto q = sqrt(q2);
-      return (-1. + q * (1. + dq2RB(k2, q2)) / (q + RF(k2, q2))) / (2. * q);
+      using T = std::decay_t<decltype(q)>;
+      return (T(-1) + q * (T(1) + dq2RB(k2, q2)) / (q + RF(k2, q2))) / (T(2) * q);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RFdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
-      return 0.5 * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return T(0.5) * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
     }
   };
 
@@ -153,14 +161,16 @@ namespace DiFfRG
       using Kokkos::exp;
       const auto xbm1 = powr<b - 1>(q2 / k2);
       const auto xb = xbm1 * (q2 / k2);
-      return -b * xbm1 * exp(-xb);
+      using T = std::decay_t<decltype(xb)>;
+      return T(-b) * xbm1 * exp(-xb);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RBdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::exp;
       const auto xb = powr<b>(q2 / k2);
-      return 2. * exp(-xb) * k2 * (1. + b * xb);
+      using T = std::decay_t<decltype(xb)>;
+      return T(2) * exp(-xb) * k2 * (T(1) + T(b) * xb);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RF(const NT1 k2, const NT2 q2)
@@ -173,13 +183,15 @@ namespace DiFfRG
     {
       using Kokkos::sqrt;
       const auto q = sqrt(q2);
-      return (-1. + q * (1. + dq2RB(k2, q2)) / (q + RF(k2, q2))) / (2. * q);
+      using T = std::decay_t<decltype(q)>;
+      return (T(-1) + q * (T(1) + dq2RB(k2, q2)) / (q + RF(k2, q2))) / (T(2) * q);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RFdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
-      return 0.5 * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return T(0.5) * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
     }
   };
 
@@ -208,15 +220,22 @@ namespace DiFfRG
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RB(const NT1 k2, const NT2 q2)
     {
       using Kokkos::exp;
-      return (k2 - q2) / (1. + exp((q2 / k2 - 1.) / alpha));
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return (k2 - q2) / (T(1) + exp((q2 / k2 - T(1)) / T(alpha)));
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RBdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::exp;
       const auto x = q2 / k2;
-      const auto mexp = exp((x - 1.) / alpha);
-      return mexp > 1e50 ? 0. : 2. * k2 * (1. + mexp * (x - powr<2>(x) + alpha) / alpha) / powr<2>(1. + mexp);
+      using T = std::decay_t<decltype(x)>;
+      const auto mexp = exp((x - T(1)) / T(alpha));
+      // Past the cutoff the result underflows; return 0 before mexp^2 overflows to inf/inf. 1e50 is
+      // beyond float range, so single precision cuts off at 1e30 instead.
+      constexpr double cutoff = std::is_same_v<T, float> ? 1e30 : 1e50;
+      return mexp > T(cutoff)
+                 ? T(0)
+                 : T(2) * k2 * (T(1) + mexp * (x - powr<2>(x) + T(alpha)) / T(alpha)) / powr<2>(T(1) + mexp);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RF(const NT1 k2, const NT2 q2)
@@ -228,7 +247,8 @@ namespace DiFfRG
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RFdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
-      return 0.5 * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return T(0.5) * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
     }
   };
 
@@ -260,289 +280,292 @@ namespace DiFfRG
     constexpr static double c = OPTS::c;
     constexpr static double b0 = OPTS::b0;
 
+    // Every literal and coefficient is converted to x's type, so a float x stays float throughout.
     static KOKKOS_INLINE_FUNCTION auto get_f(const auto x)
     {
-      constexpr double b0 = OPTS::b0;
+      using T = std::decay_t<decltype(x)>;
+      const T b0 = T(OPTS::b0), c = T(OPTS::c);
       if constexpr (order == 1) return x;
       if constexpr (order == 2)
-        return x * (-2. + c * (2. + x + 2. * b0 * x)) * powr<-1>(-2. + x + 2. * c * (1. + b0 * x));
+        return x * (-T(2.) + c * (T(2.) + x + T(2.) * b0 * x)) * powr<-1>(-T(2.) + x + T(2.) * c * (T(1.) + b0 * x));
       if constexpr (order == 3)
-        return x * (-3. * (2. + x + 2. * b0 * x) + c * (6. + (3. + 6. * b0) * x + (2. + 3. * b0) * powr<2>(x))) *
-               powr<-1>(3. * b0 * (-2. + x) * x + 6. * c * (1. + b0 * x) + 2. * (-3. + powr<2>(x)));
+        return x * (-T(3.) * (T(2.) + x + T(2.) * b0 * x) + c * (T(6.) + (T(3.) + T(6.) * b0) * x + (T(2.) + T(3.) * b0) * powr<2>(x))) *
+               powr<-1>(T(3.) * b0 * (-T(2.) + x) * x + T(6.) * c * (T(1.) + b0 * x) + T(2.) * (-T(3.) + powr<2>(x)));
       if constexpr (order == 4)
-        return 0.08333333333333333 * x *
-               (12. + 6. * x + 4. * (1. + 3. * b0) * powr<2>(x) +
-                3. * (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<3>(x)) *
-               powr<-1>(1. + b0 * powr<2>(x) + 0.25 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<3>(x));
+        return T(0.08333333333333333) * x *
+               (T(12.) + T(6.) * x + T(4.) * (T(1.) + T(3.) * b0) * powr<2>(x) +
+                T(3.) * (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<3>(x)) *
+               powr<-1>(T(1.) + b0 * powr<2>(x) + T(0.25) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<3>(x));
       if constexpr (order == 5)
-        return 0.016666666666666666 * x *
-               (60. + 30. * x + 20. * (1. + 3. * b0) * powr<2>(x) + 15. * (1. + 2. * b0) * powr<3>(x) +
-                4. * (3. + 5. * b0) * c * powr<-1>(-1. + c) * powr<4>(x)) *
-               powr<-1>(1. + b0 * powr<2>(x) + 0.06666666666666667 * (3. + 5. * b0) * powr<-1>(-1. + c) * powr<4>(x));
+        return T(0.016666666666666666) * x *
+               (T(60.) + T(30.) * x + T(20.) * (T(1.) + T(3.) * b0) * powr<2>(x) + T(15.) * (T(1.) + T(2.) * b0) * powr<3>(x) +
+                T(4.) * (T(3.) + T(5.) * b0) * c * powr<-1>(-T(1.) + c) * powr<4>(x)) *
+               powr<-1>(T(1.) + b0 * powr<2>(x) + T(0.06666666666666667) * (T(3.) + T(5.) * b0) * powr<-1>(-T(1.) + c) * powr<4>(x));
       if constexpr (order == 6)
-        return 0.016666666666666666 * x *
-               (60. + 30. * x + 20. * powr<2>(x) + 15. * (1. + 4. * b0) * powr<3>(x) +
-                6. * (2. + 5. * b0) * powr<4>(x) + 10. * (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<5>(x)) *
-               powr<-1>(1. + b0 * powr<3>(x) + 0.16666666666666666 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<5>(x));
+        return T(0.016666666666666666) * x *
+               (T(60.) + T(30.) * x + T(20.) * powr<2>(x) + T(15.) * (T(1.) + T(4.) * b0) * powr<3>(x) +
+                T(6.) * (T(2.) + T(5.) * b0) * powr<4>(x) + T(10.) * (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<5>(x)) *
+               powr<-1>(T(1.) + b0 * powr<3>(x) + T(0.16666666666666666) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<5>(x));
       if constexpr (order == 7)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + (0.25 + b0) * powr<4>(x) +
-                0.1 * (2. + 5. * b0) * powr<5>(x) + 0.16666666666666666 * (1. + 2. * b0) * powr<6>(x) +
-                0.03571428571428571 * (4. + 7. * b0) * c * powr<-1>(-1. + c) * powr<7>(x)) *
-               powr<-1>(1. + b0 * powr<3>(x) + 0.03571428571428571 * (4. + 7. * b0) * powr<-1>(-1. + c) * powr<6>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + (T(0.25) + b0) * powr<4>(x) +
+                T(0.1) * (T(2.) + T(5.) * b0) * powr<5>(x) + T(0.16666666666666666) * (T(1.) + T(2.) * b0) * powr<6>(x) +
+                T(0.03571428571428571) * (T(4.) + T(7.) * b0) * c * powr<-1>(-T(1.) + c) * powr<7>(x)) *
+               powr<-1>(T(1.) + b0 * powr<3>(x) + T(0.03571428571428571) * (T(4.) + T(7.) * b0) * powr<-1>(-T(1.) + c) * powr<6>(x));
       if constexpr (order == 8)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + (0.2 + b0) * powr<5>(x) +
-                0.16666666666666666 * (1. + 3. * b0) * powr<6>(x) + 0.047619047619047616 * (3. + 7. * b0) * powr<7>(x) +
-                0.125 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<8>(x)) *
-               powr<-1>(1. + b0 * powr<4>(x) + 0.125 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<7>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + (T(0.2) + b0) * powr<5>(x) +
+                T(0.16666666666666666) * (T(1.) + T(3.) * b0) * powr<6>(x) + T(0.047619047619047616) * (T(3.) + T(7.) * b0) * powr<7>(x) +
+                T(0.125) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<8>(x)) *
+               powr<-1>(T(1.) + b0 * powr<4>(x) + T(0.125) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<7>(x));
       if constexpr (order == 9)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + (0.2 + b0) * powr<5>(x) +
-                0.16666666666666666 * (1. + 3. * b0) * powr<6>(x) + 0.047619047619047616 * (3. + 7. * b0) * powr<7>(x) +
-                0.125 * (1. + 2. * b0) * powr<8>(x) +
-                0.022222222222222223 * (5. + 9. * b0) * c * powr<-1>(-1. + c) * powr<9>(x)) *
-               powr<-1>(1. + b0 * powr<4>(x) + 0.022222222222222223 * (5. + 9. * b0) * powr<-1>(-1. + c) * powr<8>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + (T(0.2) + b0) * powr<5>(x) +
+                T(0.16666666666666666) * (T(1.) + T(3.) * b0) * powr<6>(x) + T(0.047619047619047616) * (T(3.) + T(7.) * b0) * powr<7>(x) +
+                T(0.125) * (T(1.) + T(2.) * b0) * powr<8>(x) +
+                T(0.022222222222222223) * (T(5.) + T(9.) * b0) * c * powr<-1>(-T(1.) + c) * powr<9>(x)) *
+               powr<-1>(T(1.) + b0 * powr<4>(x) + T(0.022222222222222223) * (T(5.) + T(9.) * b0) * powr<-1>(-T(1.) + c) * powr<8>(x));
       if constexpr (order == 10)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                (0.16666666666666666 + b0) * powr<6>(x) + 0.07142857142857142 * (2. + 7. * b0) * powr<7>(x) +
-                0.041666666666666664 * (3. + 8. * b0) * powr<8>(x) +
-                0.027777777777777776 * (4. + 9. * b0) * powr<9>(x) +
-                0.1 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<10>(x)) *
-               powr<-1>(1. + b0 * powr<5>(x) + 0.1 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<9>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                (T(0.16666666666666666) + b0) * powr<6>(x) + T(0.07142857142857142) * (T(2.) + T(7.) * b0) * powr<7>(x) +
+                T(0.041666666666666664) * (T(3.) + T(8.) * b0) * powr<8>(x) +
+                T(0.027777777777777776) * (T(4.) + T(9.) * b0) * powr<9>(x) +
+                T(0.1) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<10>(x)) *
+               powr<-1>(T(1.) + b0 * powr<5>(x) + T(0.1) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<9>(x));
       if constexpr (order == 11)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                (0.16666666666666666 + b0) * powr<6>(x) + 0.07142857142857142 * (2. + 7. * b0) * powr<7>(x) +
-                0.041666666666666664 * (3. + 8. * b0) * powr<8>(x) +
-                0.027777777777777776 * (4. + 9. * b0) * powr<9>(x) + 0.1 * (1. + 2. * b0) * powr<10>(x) +
-                0.015151515151515152 * (6. + 11. * b0) * c * powr<-1>(-1. + c) * powr<11>(x)) *
-               powr<-1>(1. + b0 * powr<5>(x) +
-                        0.015151515151515152 * (6. + 11. * b0) * powr<-1>(-1. + c) * powr<10>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                (T(0.16666666666666666) + b0) * powr<6>(x) + T(0.07142857142857142) * (T(2.) + T(7.) * b0) * powr<7>(x) +
+                T(0.041666666666666664) * (T(3.) + T(8.) * b0) * powr<8>(x) +
+                T(0.027777777777777776) * (T(4.) + T(9.) * b0) * powr<9>(x) + T(0.1) * (T(1.) + T(2.) * b0) * powr<10>(x) +
+                T(0.015151515151515152) * (T(6.) + T(11.) * b0) * c * powr<-1>(-T(1.) + c) * powr<11>(x)) *
+               powr<-1>(T(1.) + b0 * powr<5>(x) +
+                        T(0.015151515151515152) * (T(6.) + T(11.) * b0) * powr<-1>(-T(1.) + c) * powr<10>(x));
       if constexpr (order == 12)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                0.16666666666666666 * powr<6>(x) + (0.14285714285714285 + b0) * powr<7>(x) +
-                0.125 * (1. + 4. * b0) * powr<8>(x) + 0.1111111111111111 * (1. + 3. * b0) * powr<9>(x) +
-                0.05 * (2. + 5. * b0) * powr<10>(x) + 0.01818181818181818 * (5. + 11. * b0) * powr<11>(x) +
-                0.08333333333333333 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<12>(x)) *
-               powr<-1>(1. + b0 * powr<6>(x) + 0.08333333333333333 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<11>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                T(0.16666666666666666) * powr<6>(x) + (T(0.14285714285714285) + b0) * powr<7>(x) +
+                T(0.125) * (T(1.) + T(4.) * b0) * powr<8>(x) + T(0.1111111111111111) * (T(1.) + T(3.) * b0) * powr<9>(x) +
+                T(0.05) * (T(2.) + T(5.) * b0) * powr<10>(x) + T(0.01818181818181818) * (T(5.) + T(11.) * b0) * powr<11>(x) +
+                T(0.08333333333333333) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<12>(x)) *
+               powr<-1>(T(1.) + b0 * powr<6>(x) + T(0.08333333333333333) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<11>(x));
       if constexpr (order == 13)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                0.16666666666666666 * powr<6>(x) + (0.14285714285714285 + b0) * powr<7>(x) +
-                0.125 * (1. + 4. * b0) * powr<8>(x) + 0.1111111111111111 * (1. + 3. * b0) * powr<9>(x) +
-                0.05 * (2. + 5. * b0) * powr<10>(x) + 0.01818181818181818 * (5. + 11. * b0) * powr<11>(x) +
-                0.08333333333333333 * (1. + 2. * b0) * powr<12>(x) +
-                0.01098901098901099 * (7. + 13. * b0) * c * powr<-1>(-1. + c) * powr<13>(x)) *
-               powr<-1>(1. + b0 * powr<6>(x) + 0.01098901098901099 * (7. + 13. * b0) * powr<-1>(-1. + c) * powr<12>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                T(0.16666666666666666) * powr<6>(x) + (T(0.14285714285714285) + b0) * powr<7>(x) +
+                T(0.125) * (T(1.) + T(4.) * b0) * powr<8>(x) + T(0.1111111111111111) * (T(1.) + T(3.) * b0) * powr<9>(x) +
+                T(0.05) * (T(2.) + T(5.) * b0) * powr<10>(x) + T(0.01818181818181818) * (T(5.) + T(11.) * b0) * powr<11>(x) +
+                T(0.08333333333333333) * (T(1.) + T(2.) * b0) * powr<12>(x) +
+                T(0.01098901098901099) * (T(7.) + T(13.) * b0) * c * powr<-1>(-T(1.) + c) * powr<13>(x)) *
+               powr<-1>(T(1.) + b0 * powr<6>(x) + T(0.01098901098901099) * (T(7.) + T(13.) * b0) * powr<-1>(-T(1.) + c) * powr<12>(x));
       if constexpr (order == 14)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                0.16666666666666666 * powr<6>(x) + 0.14285714285714285 * powr<7>(x) + (0.125 + b0) * powr<8>(x) +
-                0.05555555555555555 * (2. + 9. * b0) * powr<9>(x) +
-                0.03333333333333333 * (3. + 10. * b0) * powr<10>(x) +
-                0.022727272727272728 * (4. + 11. * b0) * powr<11>(x) +
-                0.016666666666666666 * (5. + 12. * b0) * powr<12>(x) +
-                0.01282051282051282 * (6. + 13. * b0) * powr<13>(x) +
-                0.07142857142857142 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<14>(x)) *
-               powr<-1>(1. + b0 * powr<7>(x) + 0.07142857142857142 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<13>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                T(0.16666666666666666) * powr<6>(x) + T(0.14285714285714285) * powr<7>(x) + (T(0.125) + b0) * powr<8>(x) +
+                T(0.05555555555555555) * (T(2.) + T(9.) * b0) * powr<9>(x) +
+                T(0.03333333333333333) * (T(3.) + T(10.) * b0) * powr<10>(x) +
+                T(0.022727272727272728) * (T(4.) + T(11.) * b0) * powr<11>(x) +
+                T(0.016666666666666666) * (T(5.) + T(12.) * b0) * powr<12>(x) +
+                T(0.01282051282051282) * (T(6.) + T(13.) * b0) * powr<13>(x) +
+                T(0.07142857142857142) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<14>(x)) *
+               powr<-1>(T(1.) + b0 * powr<7>(x) + T(0.07142857142857142) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<13>(x));
       if constexpr (order == 15)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                0.16666666666666666 * powr<6>(x) + 0.14285714285714285 * powr<7>(x) + (0.125 + b0) * powr<8>(x) +
-                0.05555555555555555 * (2. + 9. * b0) * powr<9>(x) +
-                0.03333333333333333 * (3. + 10. * b0) * powr<10>(x) +
-                0.022727272727272728 * (4. + 11. * b0) * powr<11>(x) +
-                0.016666666666666666 * (5. + 12. * b0) * powr<12>(x) +
-                0.01282051282051282 * (6. + 13. * b0) * powr<13>(x) +
-                0.07142857142857142 * (1. + 2. * b0) * powr<14>(x) +
-                0.008333333333333333 * (8. + 15. * b0) * c * powr<-1>(-1. + c) * powr<15>(x)) *
-               powr<-1>(1. + b0 * powr<7>(x) +
-                        0.008333333333333333 * (8. + 15. * b0) * powr<-1>(-1. + c) * powr<14>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                T(0.16666666666666666) * powr<6>(x) + T(0.14285714285714285) * powr<7>(x) + (T(0.125) + b0) * powr<8>(x) +
+                T(0.05555555555555555) * (T(2.) + T(9.) * b0) * powr<9>(x) +
+                T(0.03333333333333333) * (T(3.) + T(10.) * b0) * powr<10>(x) +
+                T(0.022727272727272728) * (T(4.) + T(11.) * b0) * powr<11>(x) +
+                T(0.016666666666666666) * (T(5.) + T(12.) * b0) * powr<12>(x) +
+                T(0.01282051282051282) * (T(6.) + T(13.) * b0) * powr<13>(x) +
+                T(0.07142857142857142) * (T(1.) + T(2.) * b0) * powr<14>(x) +
+                T(0.008333333333333333) * (T(8.) + T(15.) * b0) * c * powr<-1>(-T(1.) + c) * powr<15>(x)) *
+               powr<-1>(T(1.) + b0 * powr<7>(x) +
+                        T(0.008333333333333333) * (T(8.) + T(15.) * b0) * powr<-1>(-T(1.) + c) * powr<14>(x));
       if constexpr (order == 16)
-        return (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                0.16666666666666666 * powr<6>(x) + 0.14285714285714285 * powr<7>(x) + 0.125 * powr<8>(x) +
-                (0.1111111111111111 + b0) * powr<9>(x) + 0.1 * (1. + 5. * b0) * powr<10>(x) +
-                0.030303030303030304 * (3. + 11. * b0) * powr<11>(x) +
-                0.08333333333333333 * (1. + 3. * b0) * powr<12>(x) +
-                0.015384615384615385 * (5. + 13. * b0) * powr<13>(x) +
-                0.023809523809523808 * (3. + 7. * b0) * powr<14>(x) +
-                (0.06666666666666667 + 0.14285714285714285 * b0) * powr<15>(x) +
-                0.0625 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<16>(x)) *
-               powr<-1>(1. + b0 * powr<8>(x) + 0.0625 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<15>(x));
+        return (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                T(0.16666666666666666) * powr<6>(x) + T(0.14285714285714285) * powr<7>(x) + T(0.125) * powr<8>(x) +
+                (T(0.1111111111111111) + b0) * powr<9>(x) + T(0.1) * (T(1.) + T(5.) * b0) * powr<10>(x) +
+                T(0.030303030303030304) * (T(3.) + T(11.) * b0) * powr<11>(x) +
+                T(0.08333333333333333) * (T(1.) + T(3.) * b0) * powr<12>(x) +
+                T(0.015384615384615385) * (T(5.) + T(13.) * b0) * powr<13>(x) +
+                T(0.023809523809523808) * (T(3.) + T(7.) * b0) * powr<14>(x) +
+                (T(0.06666666666666667) + T(0.14285714285714285) * b0) * powr<15>(x) +
+                T(0.0625) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<16>(x)) *
+               powr<-1>(T(1.) + b0 * powr<8>(x) + T(0.0625) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<15>(x));
       static_assert(order <= 16, "Rational regulator of order > 16 not implemented");
     }
     static KOKKOS_INLINE_FUNCTION auto get_df(const auto x)
     {
-      constexpr double b0 = OPTS::b0;
-      if constexpr (order == 1) return 1.;
+      using T = std::decay_t<decltype(x)>;
+      const T b0 = T(OPTS::b0), c = T(OPTS::c);
+      if constexpr (order == 1) return T(1.);
       if constexpr (order == 2)
-        return (4. + c * (-8. - 4. * (1. + 2. * b0) * x + (1. + 2. * b0) * powr<2>(x)) +
-                2. * powr<2>(c) * (2. + (2. + 4. * b0) * x + b0 * (1. + 2. * b0) * powr<2>(x))) *
-               powr<-2>(-2. + x + 2. * c * (1. + b0 * x));
+        return (T(4.) + c * (-T(8.) - T(4.) * (T(1.) + T(2.) * b0) * x + (T(1.) + T(2.) * b0) * powr<2>(x)) +
+                T(2.) * powr<2>(c) * (T(2.) + (T(2.) + T(4.) * b0) * x + b0 * (T(1.) + T(2.) * b0) * powr<2>(x))) *
+               powr<-2>(-T(2.) + x + T(2.) * c * (T(1.) + b0 * x));
       if constexpr (order == 3)
-        return (1. + x + 2. * b0 * x +
-                0.3333333333333333 * (-1. + 3. * c + b0 * (-3. + 6. * c) + 3. * (-1. + c) * powr<2>(b0)) *
-                    powr<-1>(-1. + c) * powr<2>(x) +
-                0.3333333333333333 * b0 * (2. + 3. * b0) * c * powr<-1>(-1. + c) * powr<3>(x) +
-                0.027777777777777776 * c * powr<2>(2. + 3. * b0) * powr<-2>(-1. + c) * powr<4>(x)) *
-               powr<-2>(1. + b0 * x + 0.16666666666666666 * (2. + 3. * b0) * powr<-1>(-1. + c) * powr<2>(x));
+        return (T(1.) + x + T(2.) * b0 * x +
+                T(0.3333333333333333) * (-T(1.) + T(3.) * c + b0 * (-T(3.) + T(6.) * c) + T(3.) * (-T(1.) + c) * powr<2>(b0)) *
+                    powr<-1>(-T(1.) + c) * powr<2>(x) +
+                T(0.3333333333333333) * b0 * (T(2.) + T(3.) * b0) * c * powr<-1>(-T(1.) + c) * powr<3>(x) +
+                T(0.027777777777777776) * c * powr<2>(T(2.) + T(3.) * b0) * powr<-2>(-T(1.) + c) * powr<4>(x)) *
+               powr<-2>(T(1.) + b0 * x + T(0.16666666666666666) * (T(2.) + T(3.) * b0) * powr<-1>(-T(1.) + c) * powr<2>(x));
       if constexpr (order == 4)
-        return (1. + x + (1. + 2. * b0) * powr<2>(x) +
-                0.5 * (1. + 2. * b0) * (-1. + 2. * c) * powr<-1>(-1. + c) * powr<3>(x) +
-                0.041666666666666664 * (-3. + 2. * b0 * (-7. + 4. * c) + 24. * (-1. + c) * powr<2>(b0)) *
-                    powr<-1>(-1. + c) * powr<4>(x) +
-                0.5 * b0 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<5>(x) +
-                0.0625 * c * powr<2>(1. + 2. * b0) * powr<-2>(-1. + c) * powr<6>(x)) *
-               powr<-2>(1. + b0 * powr<2>(x) + 0.25 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<3>(x));
+        return (T(1.) + x + (T(1.) + T(2.) * b0) * powr<2>(x) +
+                T(0.5) * (T(1.) + T(2.) * b0) * (-T(1.) + T(2.) * c) * powr<-1>(-T(1.) + c) * powr<3>(x) +
+                T(0.041666666666666664) * (-T(3.) + T(2.) * b0 * (-T(7.) + T(4.) * c) + T(24.) * (-T(1.) + c) * powr<2>(b0)) *
+                    powr<-1>(-T(1.) + c) * powr<4>(x) +
+                T(0.5) * b0 * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<5>(x) +
+                T(0.0625) * c * powr<2>(T(1.) + T(2.) * b0) * powr<-2>(-T(1.) + c) * powr<6>(x)) *
+               powr<-2>(T(1.) + b0 * powr<2>(x) + T(0.25) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<3>(x));
       if constexpr (order == 5)
-        return ((1. + b0 * powr<2>(x) + 0.06666666666666667 * (3. + 5. * b0) * powr<-1>(-1. + c) * powr<4>(x)) *
-                    (1. + x + (1. + 3. * b0) * powr<2>(x) + (1. + 2. * b0) * powr<3>(x) +
-                     0.3333333333333333 * (3. + 5. * b0) * c * powr<-1>(-1. + c) * powr<4>(x)) -
-                0.016666666666666666 * x *
-                    (2. * b0 * x + 0.26666666666666666 * (3. + 5. * b0) * powr<-1>(-1. + c) * powr<3>(x)) *
-                    (60. + 30. * x + 20. * (1. + 3. * b0) * powr<2>(x) + 15. * (1. + 2. * b0) * powr<3>(x) +
-                     4. * (3. + 5. * b0) * c * powr<-1>(-1. + c) * powr<4>(x))) *
-               powr<-2>(1. + b0 * powr<2>(x) + 0.06666666666666667 * (3. + 5. * b0) * powr<-1>(-1. + c) * powr<4>(x));
+        return ((T(1.) + b0 * powr<2>(x) + T(0.06666666666666667) * (T(3.) + T(5.) * b0) * powr<-1>(-T(1.) + c) * powr<4>(x)) *
+                    (T(1.) + x + (T(1.) + T(3.) * b0) * powr<2>(x) + (T(1.) + T(2.) * b0) * powr<3>(x) +
+                     T(0.3333333333333333) * (T(3.) + T(5.) * b0) * c * powr<-1>(-T(1.) + c) * powr<4>(x)) -
+                T(0.016666666666666666) * x *
+                    (T(2.) * b0 * x + T(0.26666666666666666) * (T(3.) + T(5.) * b0) * powr<-1>(-T(1.) + c) * powr<3>(x)) *
+                    (T(60.) + T(30.) * x + T(20.) * (T(1.) + T(3.) * b0) * powr<2>(x) + T(15.) * (T(1.) + T(2.) * b0) * powr<3>(x) +
+                     T(4.) * (T(3.) + T(5.) * b0) * c * powr<-1>(-T(1.) + c) * powr<4>(x))) *
+               powr<-2>(T(1.) + b0 * powr<2>(x) + T(0.06666666666666667) * (T(3.) + T(5.) * b0) * powr<-1>(-T(1.) + c) * powr<4>(x));
       if constexpr (order == 6)
-        return ((1. + b0 * powr<3>(x) + 0.16666666666666666 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<5>(x)) *
-                    (1. + x + powr<2>(x) + (1. + 4. * b0) * powr<3>(x) + (1. + 2.5 * b0) * powr<4>(x) +
-                     (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<5>(x)) -
-                0.016666666666666666 * x *
-                    (3. * b0 * powr<2>(x) + 0.8333333333333334 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<4>(x)) *
-                    (60. + 30. * x + 20. * powr<2>(x) + 15. * (1. + 4. * b0) * powr<3>(x) +
-                     6. * (2. + 5. * b0) * powr<4>(x) + 10. * (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<5>(x))) *
-               powr<-2>(1. + b0 * powr<3>(x) + 0.16666666666666666 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<5>(x));
+        return ((T(1.) + b0 * powr<3>(x) + T(0.16666666666666666) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<5>(x)) *
+                    (T(1.) + x + powr<2>(x) + (T(1.) + T(4.) * b0) * powr<3>(x) + (T(1.) + T(2.5) * b0) * powr<4>(x) +
+                     (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<5>(x)) -
+                T(0.016666666666666666) * x *
+                    (T(3.) * b0 * powr<2>(x) + T(0.8333333333333334) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<4>(x)) *
+                    (T(60.) + T(30.) * x + T(20.) * powr<2>(x) + T(15.) * (T(1.) + T(4.) * b0) * powr<3>(x) +
+                     T(6.) * (T(2.) + T(5.) * b0) * powr<4>(x) + T(10.) * (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<5>(x))) *
+               powr<-2>(T(1.) + b0 * powr<3>(x) + T(0.16666666666666666) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<5>(x));
       if constexpr (order == 7)
-        return ((1. + b0 * powr<3>(x) + 0.03571428571428571 * (4. + 7. * b0) * powr<-1>(-1. + c) * powr<6>(x)) *
-                    (1. + x + powr<2>(x) + (1. + 4. * b0) * powr<3>(x) + (1. + 2.5 * b0) * powr<4>(x) +
-                     (1. + 2. * b0) * powr<5>(x) + 0.25 * (4. + 7. * b0) * c * powr<-1>(-1. + c) * powr<6>(x)) -
-                1. * (3. * b0 * powr<2>(x) + 0.21428571428571427 * (4. + 7. * b0) * powr<-1>(-1. + c) * powr<5>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + (0.25 + b0) * powr<4>(x) +
-                     0.1 * (2. + 5. * b0) * powr<5>(x) + 0.16666666666666666 * (1. + 2. * b0) * powr<6>(x) +
-                     0.03571428571428571 * (4. + 7. * b0) * c * powr<-1>(-1. + c) * powr<7>(x))) *
-               powr<-2>(1. + b0 * powr<3>(x) + 0.03571428571428571 * (4. + 7. * b0) * powr<-1>(-1. + c) * powr<6>(x));
+        return ((T(1.) + b0 * powr<3>(x) + T(0.03571428571428571) * (T(4.) + T(7.) * b0) * powr<-1>(-T(1.) + c) * powr<6>(x)) *
+                    (T(1.) + x + powr<2>(x) + (T(1.) + T(4.) * b0) * powr<3>(x) + (T(1.) + T(2.5) * b0) * powr<4>(x) +
+                     (T(1.) + T(2.) * b0) * powr<5>(x) + T(0.25) * (T(4.) + T(7.) * b0) * c * powr<-1>(-T(1.) + c) * powr<6>(x)) -
+                T(1.) * (T(3.) * b0 * powr<2>(x) + T(0.21428571428571427) * (T(4.) + T(7.) * b0) * powr<-1>(-T(1.) + c) * powr<5>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + (T(0.25) + b0) * powr<4>(x) +
+                     T(0.1) * (T(2.) + T(5.) * b0) * powr<5>(x) + T(0.16666666666666666) * (T(1.) + T(2.) * b0) * powr<6>(x) +
+                     T(0.03571428571428571) * (T(4.) + T(7.) * b0) * c * powr<-1>(-T(1.) + c) * powr<7>(x))) *
+               powr<-2>(T(1.) + b0 * powr<3>(x) + T(0.03571428571428571) * (T(4.) + T(7.) * b0) * powr<-1>(-T(1.) + c) * powr<6>(x));
       if constexpr (order == 8)
-        return ((1. + b0 * powr<4>(x) + 0.125 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<7>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + (1. + 5. * b0) * powr<4>(x) + (1. + 3. * b0) * powr<5>(x) +
-                     (1. + 2.3333333333333335 * b0) * powr<6>(x) + (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<7>(x)) -
-                1. * (4. * b0 * powr<3>(x) + 0.875 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<6>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) +
-                     (0.2 + b0) * powr<5>(x) + 0.16666666666666666 * (1. + 3. * b0) * powr<6>(x) +
-                     0.047619047619047616 * (3. + 7. * b0) * powr<7>(x) +
-                     0.125 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<8>(x))) *
-               powr<-2>(1. + b0 * powr<4>(x) + 0.125 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<7>(x));
+        return ((T(1.) + b0 * powr<4>(x) + T(0.125) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<7>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + (T(1.) + T(5.) * b0) * powr<4>(x) + (T(1.) + T(3.) * b0) * powr<5>(x) +
+                     (T(1.) + T(2.3333333333333335) * b0) * powr<6>(x) + (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<7>(x)) -
+                T(1.) * (T(4.) * b0 * powr<3>(x) + T(0.875) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<6>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) +
+                     (T(0.2) + b0) * powr<5>(x) + T(0.16666666666666666) * (T(1.) + T(3.) * b0) * powr<6>(x) +
+                     T(0.047619047619047616) * (T(3.) + T(7.) * b0) * powr<7>(x) +
+                     T(0.125) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<8>(x))) *
+               powr<-2>(T(1.) + b0 * powr<4>(x) + T(0.125) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<7>(x));
       if constexpr (order == 9)
-        return ((1. + b0 * powr<4>(x) + 0.022222222222222223 * (5. + 9. * b0) * powr<-1>(-1. + c) * powr<8>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + (1. + 5. * b0) * powr<4>(x) + (1. + 3. * b0) * powr<5>(x) +
-                     (1. + 2.3333333333333335 * b0) * powr<6>(x) + (1. + 2. * b0) * powr<7>(x) +
-                     0.2 * (5. + 9. * b0) * c * powr<-1>(-1. + c) * powr<8>(x)) -
-                1. * (4. * b0 * powr<3>(x) + 0.17777777777777778 * (5. + 9. * b0) * powr<-1>(-1. + c) * powr<7>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) +
-                     (0.2 + b0) * powr<5>(x) + 0.16666666666666666 * (1. + 3. * b0) * powr<6>(x) +
-                     0.047619047619047616 * (3. + 7. * b0) * powr<7>(x) + 0.125 * (1. + 2. * b0) * powr<8>(x) +
-                     0.022222222222222223 * (5. + 9. * b0) * c * powr<-1>(-1. + c) * powr<9>(x))) *
-               powr<-2>(1. + b0 * powr<4>(x) + 0.022222222222222223 * (5. + 9. * b0) * powr<-1>(-1. + c) * powr<8>(x));
+        return ((T(1.) + b0 * powr<4>(x) + T(0.022222222222222223) * (T(5.) + T(9.) * b0) * powr<-1>(-T(1.) + c) * powr<8>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + (T(1.) + T(5.) * b0) * powr<4>(x) + (T(1.) + T(3.) * b0) * powr<5>(x) +
+                     (T(1.) + T(2.3333333333333335) * b0) * powr<6>(x) + (T(1.) + T(2.) * b0) * powr<7>(x) +
+                     T(0.2) * (T(5.) + T(9.) * b0) * c * powr<-1>(-T(1.) + c) * powr<8>(x)) -
+                T(1.) * (T(4.) * b0 * powr<3>(x) + T(0.17777777777777778) * (T(5.) + T(9.) * b0) * powr<-1>(-T(1.) + c) * powr<7>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) +
+                     (T(0.2) + b0) * powr<5>(x) + T(0.16666666666666666) * (T(1.) + T(3.) * b0) * powr<6>(x) +
+                     T(0.047619047619047616) * (T(3.) + T(7.) * b0) * powr<7>(x) + T(0.125) * (T(1.) + T(2.) * b0) * powr<8>(x) +
+                     T(0.022222222222222223) * (T(5.) + T(9.) * b0) * c * powr<-1>(-T(1.) + c) * powr<9>(x))) *
+               powr<-2>(T(1.) + b0 * powr<4>(x) + T(0.022222222222222223) * (T(5.) + T(9.) * b0) * powr<-1>(-T(1.) + c) * powr<8>(x));
       if constexpr (order == 10)
-        return ((1. + b0 * powr<5>(x) + 0.1 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<9>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + (1. + 6. * b0) * powr<5>(x) +
-                     (1. + 3.5 * b0) * powr<6>(x) + (1. + 2.6666666666666665 * b0) * powr<7>(x) +
-                     (1. + 2.25 * b0) * powr<8>(x) + (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<9>(x)) -
-                1. * (5. * b0 * powr<4>(x) + 0.9 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<8>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     (0.16666666666666666 + b0) * powr<6>(x) + 0.07142857142857142 * (2. + 7. * b0) * powr<7>(x) +
-                     0.041666666666666664 * (3. + 8. * b0) * powr<8>(x) +
-                     0.027777777777777776 * (4. + 9. * b0) * powr<9>(x) +
-                     0.1 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<10>(x))) *
-               powr<-2>(1. + b0 * powr<5>(x) + 0.1 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<9>(x));
+        return ((T(1.) + b0 * powr<5>(x) + T(0.1) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<9>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + (T(1.) + T(6.) * b0) * powr<5>(x) +
+                     (T(1.) + T(3.5) * b0) * powr<6>(x) + (T(1.) + T(2.6666666666666665) * b0) * powr<7>(x) +
+                     (T(1.) + T(2.25) * b0) * powr<8>(x) + (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<9>(x)) -
+                T(1.) * (T(5.) * b0 * powr<4>(x) + T(0.9) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<8>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     (T(0.16666666666666666) + b0) * powr<6>(x) + T(0.07142857142857142) * (T(2.) + T(7.) * b0) * powr<7>(x) +
+                     T(0.041666666666666664) * (T(3.) + T(8.) * b0) * powr<8>(x) +
+                     T(0.027777777777777776) * (T(4.) + T(9.) * b0) * powr<9>(x) +
+                     T(0.1) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<10>(x))) *
+               powr<-2>(T(1.) + b0 * powr<5>(x) + T(0.1) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<9>(x));
       if constexpr (order == 11)
-        return ((1. + b0 * powr<5>(x) + 0.015151515151515152 * (6. + 11. * b0) * powr<-1>(-1. + c) * powr<10>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + (1. + 6. * b0) * powr<5>(x) +
-                     (1. + 3.5 * b0) * powr<6>(x) + (1. + 2.6666666666666665 * b0) * powr<7>(x) +
-                     (1. + 2.25 * b0) * powr<8>(x) + (1. + 2. * b0) * powr<9>(x) +
-                     0.16666666666666666 * (6. + 11. * b0) * c * powr<-1>(-1. + c) * powr<10>(x)) -
-                1. * (5. * b0 * powr<4>(x) + 0.15151515151515152 * (6. + 11. * b0) * powr<-1>(-1. + c) * powr<9>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     (0.16666666666666666 + b0) * powr<6>(x) + 0.07142857142857142 * (2. + 7. * b0) * powr<7>(x) +
-                     0.041666666666666664 * (3. + 8. * b0) * powr<8>(x) +
-                     0.027777777777777776 * (4. + 9. * b0) * powr<9>(x) + 0.1 * (1. + 2. * b0) * powr<10>(x) +
-                     0.015151515151515152 * (6. + 11. * b0) * c * powr<-1>(-1. + c) * powr<11>(x))) *
-               powr<-2>(1. + b0 * powr<5>(x) +
-                        0.015151515151515152 * (6. + 11. * b0) * powr<-1>(-1. + c) * powr<10>(x));
+        return ((T(1.) + b0 * powr<5>(x) + T(0.015151515151515152) * (T(6.) + T(11.) * b0) * powr<-1>(-T(1.) + c) * powr<10>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + (T(1.) + T(6.) * b0) * powr<5>(x) +
+                     (T(1.) + T(3.5) * b0) * powr<6>(x) + (T(1.) + T(2.6666666666666665) * b0) * powr<7>(x) +
+                     (T(1.) + T(2.25) * b0) * powr<8>(x) + (T(1.) + T(2.) * b0) * powr<9>(x) +
+                     T(0.16666666666666666) * (T(6.) + T(11.) * b0) * c * powr<-1>(-T(1.) + c) * powr<10>(x)) -
+                T(1.) * (T(5.) * b0 * powr<4>(x) + T(0.15151515151515152) * (T(6.) + T(11.) * b0) * powr<-1>(-T(1.) + c) * powr<9>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     (T(0.16666666666666666) + b0) * powr<6>(x) + T(0.07142857142857142) * (T(2.) + T(7.) * b0) * powr<7>(x) +
+                     T(0.041666666666666664) * (T(3.) + T(8.) * b0) * powr<8>(x) +
+                     T(0.027777777777777776) * (T(4.) + T(9.) * b0) * powr<9>(x) + T(0.1) * (T(1.) + T(2.) * b0) * powr<10>(x) +
+                     T(0.015151515151515152) * (T(6.) + T(11.) * b0) * c * powr<-1>(-T(1.) + c) * powr<11>(x))) *
+               powr<-2>(T(1.) + b0 * powr<5>(x) +
+                        T(0.015151515151515152) * (T(6.) + T(11.) * b0) * powr<-1>(-T(1.) + c) * powr<10>(x));
       if constexpr (order == 12)
-        return ((1. + b0 * powr<6>(x) + 0.08333333333333333 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<11>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + (1. + 7. * b0) * powr<6>(x) +
-                     (1. + 4. * b0) * powr<7>(x) + (1. + 3. * b0) * powr<8>(x) + (1. + 2.5 * b0) * powr<9>(x) +
-                     (1. + 2.2 * b0) * powr<10>(x) + (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<11>(x)) -
-                1. * (6. * b0 * powr<5>(x) + 0.9166666666666666 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<10>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     0.16666666666666666 * powr<6>(x) + (0.14285714285714285 + b0) * powr<7>(x) +
-                     0.125 * (1. + 4. * b0) * powr<8>(x) + 0.1111111111111111 * (1. + 3. * b0) * powr<9>(x) +
-                     0.05 * (2. + 5. * b0) * powr<10>(x) + 0.01818181818181818 * (5. + 11. * b0) * powr<11>(x) +
-                     0.08333333333333333 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<12>(x))) *
-               powr<-2>(1. + b0 * powr<6>(x) + 0.08333333333333333 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<11>(x));
+        return ((T(1.) + b0 * powr<6>(x) + T(0.08333333333333333) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<11>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + (T(1.) + T(7.) * b0) * powr<6>(x) +
+                     (T(1.) + T(4.) * b0) * powr<7>(x) + (T(1.) + T(3.) * b0) * powr<8>(x) + (T(1.) + T(2.5) * b0) * powr<9>(x) +
+                     (T(1.) + T(2.2) * b0) * powr<10>(x) + (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<11>(x)) -
+                T(1.) * (T(6.) * b0 * powr<5>(x) + T(0.9166666666666666) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<10>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     T(0.16666666666666666) * powr<6>(x) + (T(0.14285714285714285) + b0) * powr<7>(x) +
+                     T(0.125) * (T(1.) + T(4.) * b0) * powr<8>(x) + T(0.1111111111111111) * (T(1.) + T(3.) * b0) * powr<9>(x) +
+                     T(0.05) * (T(2.) + T(5.) * b0) * powr<10>(x) + T(0.01818181818181818) * (T(5.) + T(11.) * b0) * powr<11>(x) +
+                     T(0.08333333333333333) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<12>(x))) *
+               powr<-2>(T(1.) + b0 * powr<6>(x) + T(0.08333333333333333) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<11>(x));
       if constexpr (order == 13)
-        return ((1. + b0 * powr<6>(x) + 0.01098901098901099 * (7. + 13. * b0) * powr<-1>(-1. + c) * powr<12>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + (1. + 7. * b0) * powr<6>(x) +
-                     (1. + 4. * b0) * powr<7>(x) + (1. + 3. * b0) * powr<8>(x) + (1. + 2.5 * b0) * powr<9>(x) +
-                     (1. + 2.2 * b0) * powr<10>(x) + (1. + 2. * b0) * powr<11>(x) +
-                     0.14285714285714285 * (7. + 13. * b0) * c * powr<-1>(-1. + c) * powr<12>(x)) -
-                1. * (6. * b0 * powr<5>(x) + 0.13186813186813187 * (7. + 13. * b0) * powr<-1>(-1. + c) * powr<11>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     0.16666666666666666 * powr<6>(x) + (0.14285714285714285 + b0) * powr<7>(x) +
-                     0.125 * (1. + 4. * b0) * powr<8>(x) + 0.1111111111111111 * (1. + 3. * b0) * powr<9>(x) +
-                     0.05 * (2. + 5. * b0) * powr<10>(x) + 0.01818181818181818 * (5. + 11. * b0) * powr<11>(x) +
-                     0.08333333333333333 * (1. + 2. * b0) * powr<12>(x) +
-                     0.01098901098901099 * (7. + 13. * b0) * c * powr<-1>(-1. + c) * powr<13>(x))) *
-               powr<-2>(1. + b0 * powr<6>(x) + 0.01098901098901099 * (7. + 13. * b0) * powr<-1>(-1. + c) * powr<12>(x));
+        return ((T(1.) + b0 * powr<6>(x) + T(0.01098901098901099) * (T(7.) + T(13.) * b0) * powr<-1>(-T(1.) + c) * powr<12>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + (T(1.) + T(7.) * b0) * powr<6>(x) +
+                     (T(1.) + T(4.) * b0) * powr<7>(x) + (T(1.) + T(3.) * b0) * powr<8>(x) + (T(1.) + T(2.5) * b0) * powr<9>(x) +
+                     (T(1.) + T(2.2) * b0) * powr<10>(x) + (T(1.) + T(2.) * b0) * powr<11>(x) +
+                     T(0.14285714285714285) * (T(7.) + T(13.) * b0) * c * powr<-1>(-T(1.) + c) * powr<12>(x)) -
+                T(1.) * (T(6.) * b0 * powr<5>(x) + T(0.13186813186813187) * (T(7.) + T(13.) * b0) * powr<-1>(-T(1.) + c) * powr<11>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     T(0.16666666666666666) * powr<6>(x) + (T(0.14285714285714285) + b0) * powr<7>(x) +
+                     T(0.125) * (T(1.) + T(4.) * b0) * powr<8>(x) + T(0.1111111111111111) * (T(1.) + T(3.) * b0) * powr<9>(x) +
+                     T(0.05) * (T(2.) + T(5.) * b0) * powr<10>(x) + T(0.01818181818181818) * (T(5.) + T(11.) * b0) * powr<11>(x) +
+                     T(0.08333333333333333) * (T(1.) + T(2.) * b0) * powr<12>(x) +
+                     T(0.01098901098901099) * (T(7.) + T(13.) * b0) * c * powr<-1>(-T(1.) + c) * powr<13>(x))) *
+               powr<-2>(T(1.) + b0 * powr<6>(x) + T(0.01098901098901099) * (T(7.) + T(13.) * b0) * powr<-1>(-T(1.) + c) * powr<12>(x));
       if constexpr (order == 14)
-        return ((1. + b0 * powr<7>(x) + 0.07142857142857142 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<13>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + powr<6>(x) +
-                     (1. + 8. * b0) * powr<7>(x) + (1. + 4.5 * b0) * powr<8>(x) +
-                     (1. + 3.3333333333333335 * b0) * powr<9>(x) + (1. + 2.75 * b0) * powr<10>(x) +
-                     (1. + 2.4 * b0) * powr<11>(x) + (1. + 2.1666666666666665 * b0) * powr<12>(x) +
-                     (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<13>(x)) -
-                1. * (7. * b0 * powr<6>(x) + 0.9285714285714286 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<12>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     0.16666666666666666 * powr<6>(x) + 0.14285714285714285 * powr<7>(x) + (0.125 + b0) * powr<8>(x) +
-                     0.05555555555555555 * (2. + 9. * b0) * powr<9>(x) +
-                     0.03333333333333333 * (3. + 10. * b0) * powr<10>(x) +
-                     0.022727272727272728 * (4. + 11. * b0) * powr<11>(x) +
-                     0.016666666666666666 * (5. + 12. * b0) * powr<12>(x) +
-                     0.01282051282051282 * (6. + 13. * b0) * powr<13>(x) +
-                     0.07142857142857142 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<14>(x))) *
-               powr<-2>(1. + b0 * powr<7>(x) + 0.07142857142857142 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<13>(x));
+        return ((T(1.) + b0 * powr<7>(x) + T(0.07142857142857142) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<13>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + powr<6>(x) +
+                     (T(1.) + T(8.) * b0) * powr<7>(x) + (T(1.) + T(4.5) * b0) * powr<8>(x) +
+                     (T(1.) + T(3.3333333333333335) * b0) * powr<9>(x) + (T(1.) + T(2.75) * b0) * powr<10>(x) +
+                     (T(1.) + T(2.4) * b0) * powr<11>(x) + (T(1.) + T(2.1666666666666665) * b0) * powr<12>(x) +
+                     (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<13>(x)) -
+                T(1.) * (T(7.) * b0 * powr<6>(x) + T(0.9285714285714286) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<12>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     T(0.16666666666666666) * powr<6>(x) + T(0.14285714285714285) * powr<7>(x) + (T(0.125) + b0) * powr<8>(x) +
+                     T(0.05555555555555555) * (T(2.) + T(9.) * b0) * powr<9>(x) +
+                     T(0.03333333333333333) * (T(3.) + T(10.) * b0) * powr<10>(x) +
+                     T(0.022727272727272728) * (T(4.) + T(11.) * b0) * powr<11>(x) +
+                     T(0.016666666666666666) * (T(5.) + T(12.) * b0) * powr<12>(x) +
+                     T(0.01282051282051282) * (T(6.) + T(13.) * b0) * powr<13>(x) +
+                     T(0.07142857142857142) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<14>(x))) *
+               powr<-2>(T(1.) + b0 * powr<7>(x) + T(0.07142857142857142) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<13>(x));
       if constexpr (order == 15)
-        return ((1. + b0 * powr<7>(x) + 0.008333333333333333 * (8. + 15. * b0) * powr<-1>(-1. + c) * powr<14>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + powr<6>(x) +
-                     (1. + 8. * b0) * powr<7>(x) + (1. + 4.5 * b0) * powr<8>(x) +
-                     (1. + 3.3333333333333335 * b0) * powr<9>(x) + (1. + 2.75 * b0) * powr<10>(x) +
-                     (1. + 2.4 * b0) * powr<11>(x) + (1. + 2.1666666666666665 * b0) * powr<12>(x) +
-                     (1. + 2. * b0) * powr<13>(x) + 0.125 * (8. + 15. * b0) * c * powr<-1>(-1. + c) * powr<14>(x)) -
-                1. * (7. * b0 * powr<6>(x) + 0.11666666666666667 * (8. + 15. * b0) * powr<-1>(-1. + c) * powr<13>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     0.16666666666666666 * powr<6>(x) + 0.14285714285714285 * powr<7>(x) + (0.125 + b0) * powr<8>(x) +
-                     0.05555555555555555 * (2. + 9. * b0) * powr<9>(x) +
-                     0.03333333333333333 * (3. + 10. * b0) * powr<10>(x) +
-                     0.022727272727272728 * (4. + 11. * b0) * powr<11>(x) +
-                     0.016666666666666666 * (5. + 12. * b0) * powr<12>(x) +
-                     0.01282051282051282 * (6. + 13. * b0) * powr<13>(x) +
-                     0.07142857142857142 * (1. + 2. * b0) * powr<14>(x) +
-                     0.008333333333333333 * (8. + 15. * b0) * c * powr<-1>(-1. + c) * powr<15>(x))) *
-               powr<-2>(1. + b0 * powr<7>(x) +
-                        0.008333333333333333 * (8. + 15. * b0) * powr<-1>(-1. + c) * powr<14>(x));
+        return ((T(1.) + b0 * powr<7>(x) + T(0.008333333333333333) * (T(8.) + T(15.) * b0) * powr<-1>(-T(1.) + c) * powr<14>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + powr<6>(x) +
+                     (T(1.) + T(8.) * b0) * powr<7>(x) + (T(1.) + T(4.5) * b0) * powr<8>(x) +
+                     (T(1.) + T(3.3333333333333335) * b0) * powr<9>(x) + (T(1.) + T(2.75) * b0) * powr<10>(x) +
+                     (T(1.) + T(2.4) * b0) * powr<11>(x) + (T(1.) + T(2.1666666666666665) * b0) * powr<12>(x) +
+                     (T(1.) + T(2.) * b0) * powr<13>(x) + T(0.125) * (T(8.) + T(15.) * b0) * c * powr<-1>(-T(1.) + c) * powr<14>(x)) -
+                T(1.) * (T(7.) * b0 * powr<6>(x) + T(0.11666666666666667) * (T(8.) + T(15.) * b0) * powr<-1>(-T(1.) + c) * powr<13>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     T(0.16666666666666666) * powr<6>(x) + T(0.14285714285714285) * powr<7>(x) + (T(0.125) + b0) * powr<8>(x) +
+                     T(0.05555555555555555) * (T(2.) + T(9.) * b0) * powr<9>(x) +
+                     T(0.03333333333333333) * (T(3.) + T(10.) * b0) * powr<10>(x) +
+                     T(0.022727272727272728) * (T(4.) + T(11.) * b0) * powr<11>(x) +
+                     T(0.016666666666666666) * (T(5.) + T(12.) * b0) * powr<12>(x) +
+                     T(0.01282051282051282) * (T(6.) + T(13.) * b0) * powr<13>(x) +
+                     T(0.07142857142857142) * (T(1.) + T(2.) * b0) * powr<14>(x) +
+                     T(0.008333333333333333) * (T(8.) + T(15.) * b0) * c * powr<-1>(-T(1.) + c) * powr<15>(x))) *
+               powr<-2>(T(1.) + b0 * powr<7>(x) +
+                        T(0.008333333333333333) * (T(8.) + T(15.) * b0) * powr<-1>(-T(1.) + c) * powr<14>(x));
       if constexpr (order == 16)
-        return ((1. + b0 * powr<8>(x) + 0.0625 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<15>(x)) *
-                    (1. + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + powr<6>(x) + powr<7>(x) +
-                     (1. + 9. * b0) * powr<8>(x) + (1. + 5. * b0) * powr<9>(x) +
-                     (1. + 3.6666666666666665 * b0) * powr<10>(x) + (1. + 3. * b0) * powr<11>(x) +
-                     (1. + 2.6 * b0) * powr<12>(x) + (1. + 2.3333333333333335 * b0) * powr<13>(x) +
-                     (1. + 2.142857142857143 * b0) * powr<14>(x) +
-                     (c + 2. * b0 * c) * powr<-1>(-1. + c) * powr<15>(x)) -
-                1. * (8. * b0 * powr<7>(x) + 0.9375 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<14>(x)) *
-                    (x + 0.5 * powr<2>(x) + 0.3333333333333333 * powr<3>(x) + 0.25 * powr<4>(x) + 0.2 * powr<5>(x) +
-                     0.16666666666666666 * powr<6>(x) + 0.14285714285714285 * powr<7>(x) + 0.125 * powr<8>(x) +
-                     (0.1111111111111111 + b0) * powr<9>(x) + 0.1 * (1. + 5. * b0) * powr<10>(x) +
-                     0.030303030303030304 * (3. + 11. * b0) * powr<11>(x) +
-                     0.08333333333333333 * (1. + 3. * b0) * powr<12>(x) +
-                     0.015384615384615385 * (5. + 13. * b0) * powr<13>(x) +
-                     0.023809523809523808 * (3. + 7. * b0) * powr<14>(x) +
-                     (0.06666666666666667 + 0.14285714285714285 * b0) * powr<15>(x) +
-                     0.0625 * (1. + 2. * b0) * c * powr<-1>(-1. + c) * powr<16>(x))) *
-               powr<-2>(1. + b0 * powr<8>(x) + 0.0625 * (1. + 2. * b0) * powr<-1>(-1. + c) * powr<15>(x));
+        return ((T(1.) + b0 * powr<8>(x) + T(0.0625) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<15>(x)) *
+                    (T(1.) + x + powr<2>(x) + powr<3>(x) + powr<4>(x) + powr<5>(x) + powr<6>(x) + powr<7>(x) +
+                     (T(1.) + T(9.) * b0) * powr<8>(x) + (T(1.) + T(5.) * b0) * powr<9>(x) +
+                     (T(1.) + T(3.6666666666666665) * b0) * powr<10>(x) + (T(1.) + T(3.) * b0) * powr<11>(x) +
+                     (T(1.) + T(2.6) * b0) * powr<12>(x) + (T(1.) + T(2.3333333333333335) * b0) * powr<13>(x) +
+                     (T(1.) + T(2.142857142857143) * b0) * powr<14>(x) +
+                     (c + T(2.) * b0 * c) * powr<-1>(-T(1.) + c) * powr<15>(x)) -
+                T(1.) * (T(8.) * b0 * powr<7>(x) + T(0.9375) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<14>(x)) *
+                    (x + T(0.5) * powr<2>(x) + T(0.3333333333333333) * powr<3>(x) + T(0.25) * powr<4>(x) + T(0.2) * powr<5>(x) +
+                     T(0.16666666666666666) * powr<6>(x) + T(0.14285714285714285) * powr<7>(x) + T(0.125) * powr<8>(x) +
+                     (T(0.1111111111111111) + b0) * powr<9>(x) + T(0.1) * (T(1.) + T(5.) * b0) * powr<10>(x) +
+                     T(0.030303030303030304) * (T(3.) + T(11.) * b0) * powr<11>(x) +
+                     T(0.08333333333333333) * (T(1.) + T(3.) * b0) * powr<12>(x) +
+                     T(0.015384615384615385) * (T(5.) + T(13.) * b0) * powr<13>(x) +
+                     T(0.023809523809523808) * (T(3.) + T(7.) * b0) * powr<14>(x) +
+                     (T(0.06666666666666667) + T(0.14285714285714285) * b0) * powr<15>(x) +
+                     T(0.0625) * (T(1.) + T(2.) * b0) * c * powr<-1>(-T(1.) + c) * powr<16>(x))) *
+               powr<-2>(T(1.) + b0 * powr<8>(x) + T(0.0625) * (T(1.) + T(2.) * b0) * powr<-1>(-T(1.) + c) * powr<15>(x));
       static_assert(order <= 16, "Rational regulator of order > 16 not implemented");
     }
 
@@ -560,7 +583,8 @@ namespace DiFfRG
       const auto x = q2 / k2;
       const auto f = get_f(x);
       const auto df = get_df(x);
-      return 2. * k2 * exp(-f) * (1. + df * x);
+      using T = std::decay_t<decltype(x)>;
+      return T(2) * k2 * exp(-f) * (T(1) + df * x);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto dq2RB(const NT1 k2, const NT2 q2)
@@ -581,14 +605,16 @@ namespace DiFfRG
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RFdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
-      return 0.5 * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
+      return T(0.5) * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
     }
 
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto dq2RF(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
       const auto q = sqrt(q2);
-      return (-1. + q * (1. + dq2RB(k2, q2)) / (q + RF(k2, q2))) / (2. * q);
+      using T = std::decay_t<decltype(q)>;
+      return (T(-1) + q * (T(1) + dq2RB(k2, q2)) / (q + RF(k2, q2))) / (T(2) * q);
     }
   };
 
@@ -665,7 +691,7 @@ namespace DiFfRG
     template <typename NT1, typename NT2> static KOKKOS_INLINE_FUNCTION auto RFdot(const NT1 k2, const NT2 q2)
     {
       using Kokkos::sqrt;
-      using T = decltype(k2 * q2);
+      using T = std::decay_t<decltype(k2 * q2)>;
       return T(0.5) * RBdot(k2, q2) / sqrt(RB(k2, q2) + q2);
     }
 
diff --git a/DiFfRG/tests/physics/CMakeLists.txt b/DiFfRG/tests/physics/CMakeLists.txt
index 50600c286..8906c884a 100644
--- a/DiFfRG/tests/physics/CMakeLists.txt
+++ b/DiFfRG/tests/physics/CMakeLists.txt
@@ -9,3 +9,9 @@ add_subdirectory(interpolation)
 # Test loop integrals
 # ##############################################################################
 
+
+# ##############################################################################
+# Regulators
+# ##############################################################################
+
+setup_test(regulators.cc)
diff --git a/DiFfRG/tests/physics/integration/finiteT/integrator_fT_p2.cc b/DiFfRG/tests/physics/integration/finiteT/integrator_fT_p2.cc
index 9579108c6..23dcf22a1 100644
--- a/DiFfRG/tests/physics/integration/finiteT/integrator_fT_p2.cc
+++ b/DiFfRG/tests/physics/integration/finiteT/integrator_fT_p2.cc
@@ -132,6 +132,49 @@ TEMPLATE_TEST_CASE_SIG("Test finite T momentum integrals", "[integration][quadra
   SECTION("GPU") { check(GPU_exec(), (double)0); }
 }
 
+TEST_CASE("Integrator_fT_p2 in single precision writes double results", "[integration][quadrature][float]")
+{
+  DiFfRG::Init();
+
+  // The float integrator must agree with the double one to float accuracy, whether it hands its
+  // results back through get() or map() into double destinations.
+  auto check = [](auto execution_space) {
+    using ExecutionSpace = std::decay_t<decltype(execution_space)>;
+    constexpr int dim = 4;
+    const float m2 = 0.5f;
+
+    QuadratureProvider quadrature_provider;
+    Integrator_fT_p2<dim, float, PolyIntegrand<2, float, -1>, ExecutionSpace> integrator_f(quadrature_provider,
+                                                                                           {32}, 2.f);
+    Integrator_fT_p2<dim, double, PolyIntegrand<2, double, -1>, ExecutionSpace> integrator_d(quadrature_provider,
+                                                                                             {32}, 2.);
+    integrator_f.set_T(0.1f);
+    integrator_d.set_T(0.1);
+    integrator_f.set_k(0.8f);
+    integrator_d.set_k(0.8);
+    integrator_f.set_typical_E(1.f);
+    integrator_d.set_typical_E(1.);
+
+    double integral_f = 0., integral_d = 0.;
+    integrator_f.get(integral_f, 0.25f, 1.f, 0.f, 0.f, 0.f, m2, 0.f, 1.f, 0.f);
+    integrator_d.get(integral_d, 0.25, 1., 0., 0., 0., m2, 0., 1., 0.);
+    CHECK(integral_f == Catch::Approx(integral_d).epsilon(1e-5));
+
+    // map() passes the grid position as the kernel's constant term.
+    const LinearCoordinates1D<float> coordinates_f(16, 0.f, 1.f);
+    const LinearCoordinates1D<double> coordinates_d(16, 0., 1.);
+    std::vector<double> mapped_f(coordinates_f.size()), mapped_d(coordinates_d.size());
+    integrator_f.map(mapped_f.data(), coordinates_f, 1.f, 0.f, 0.f, 0.f, m2, 0.f, 1.f, 0.f);
+    integrator_d.map(mapped_d.data(), coordinates_d, 1., 0., 0., 0., m2, 0., 1., 0.);
+    for (size_t i = 0; i < mapped_f.size(); ++i)
+      CHECK(mapped_f[i] == Catch::Approx(mapped_d[i]).epsilon(1e-5));
+  };
+
+  SECTION("TBB") { check(TBB_exec()); }
+  SECTION("Threads") { check(KokkosHost_exec()); }
+  SECTION("GPU") { check(GPU_exec()); }
+}
+
 TEST_CASE("Test integrator_fT_p2 bug", "[integration][quadrature]")
 {
   using NT = double;
diff --git a/DiFfRG/tests/physics/integration/vacuum/integrator_p2.cc b/DiFfRG/tests/physics/integration/vacuum/integrator_p2.cc
index 8c8f6d14d..4d9de929b 100644
--- a/DiFfRG/tests/physics/integration/vacuum/integrator_p2.cc
+++ b/DiFfRG/tests/physics/integration/vacuum/integrator_p2.cc
@@ -114,6 +114,43 @@ TEMPLATE_TEST_CASE_SIG("Test 1D momentum integrals", "[integration][quadrature]"
   SECTION("GPU") { check(GPU_exec(), (double)0); }
 }
 
+TEST_CASE("Integrator_p2 in single precision writes double results", "[integration][quadrature][float]")
+{
+  DiFfRG::Init();
+
+  // The float integrator must agree with the double one to float accuracy, whether it hands its
+  // results back through get() or map() into double destinations.
+  auto check = [](auto execution_space) {
+    using ExecutionSpace = std::decay_t<decltype(execution_space)>;
+    constexpr int dim = 4;
+    const std::array<float, 4> x{0.3f, -0.7f, 0.2f, 0.9f};
+
+    QuadratureProvider quadrature_provider;
+    Integrator_p2<dim, float, PolyIntegrand<1, float>, ExecutionSpace> integrator_f(quadrature_provider, {32}, 2.f);
+    Integrator_p2<dim, double, PolyIntegrand<1, double>, ExecutionSpace> integrator_d(quadrature_provider, {32}, 2.);
+    integrator_f.set_k(0.8f);
+    integrator_d.set_k(0.8);
+
+    double integral_f = 0., integral_d = 0.;
+    integrator_f.get(integral_f, 0.25f, x[0], x[1], x[2], x[3]);
+    integrator_d.get(integral_d, 0.25, x[0], x[1], x[2], x[3]);
+    CHECK(integral_f == Catch::Approx(integral_d).epsilon(1e-5));
+
+    // map() passes the grid position as the kernel's constant term.
+    const LinearCoordinates1D<float> coordinates_f(16, 0.f, 1.f);
+    const LinearCoordinates1D<double> coordinates_d(16, 0., 1.);
+    std::vector<double> mapped_f(coordinates_f.size()), mapped_d(coordinates_d.size());
+    integrator_f.map(mapped_f.data(), coordinates_f, x[0], x[1], x[2], x[3]);
+    integrator_d.map(mapped_d.data(), coordinates_d, x[0], x[1], x[2], x[3]);
+    for (size_t i = 0; i < mapped_f.size(); ++i)
+      CHECK(mapped_f[i] == Catch::Approx(mapped_d[i]).epsilon(1e-5));
+  };
+
+  SECTION("TBB") { check(TBB_exec()); }
+  SECTION("Threads") { check(KokkosHost_exec()); }
+  SECTION("GPU") { check(GPU_exec()); }
+}
+
 TEST_CASE("Integrator_p2 propagates second-order complex autodiff through TBB", "[integration][quadrature][autodiff]")
 {
   DiFfRG::Init();
diff --git a/DiFfRG/tests/physics/regulators.cc b/DiFfRG/tests/physics/regulators.cc
new file mode 100644
index 000000000..d2fd1d47a
--- /dev/null
+++ b/DiFfRG/tests/physics/regulators.cc
@@ -0,0 +1,62 @@
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch_all.hpp>
+
+#include <DiFfRG/physics/regulators.hh>
+
+#include <algorithm>
+#include <cmath>
+#include <type_traits>
+#include <vector>
+
+using namespace DiFfRG;
+
+namespace
+{
+  template <int n> struct RationalExpOrder : RationalExpRegulatorOpts {
+    static constexpr int order = n;
+  };
+  template <int n> struct PolynomialExpOrder : PolynomialExpRegulatorOpts {
+    static constexpr int order = n;
+  };
+
+  // Every regulator function evaluated in float must stay float -- a double constant anywhere in
+  // its body would promote the result -- and agree with the double evaluation to float accuracy.
+  // The error is taken against the function's scale over the sample, not pointwise: several of
+  // these functions pass through zero.
+  template <typename Reg, typename F> void check_function(const char *name, F &&f)
+  {
+    static_assert(std::is_same_v<decltype(f(1.f, 1.f)), float>, "regulator function promotes float");
+    std::vector<double> err, ref;
+    const double k2 = 1.3;
+    for (int i = 0; i <= 400; ++i) {
+      const double q2 = k2 * (1e-3 + 4. * i / 400.);
+      const double d = f(k2, q2);
+      const float s = f(float(k2), float(q2));
+      if (!std::isfinite(d)) continue;
+      ref.push_back(std::abs(d));
+      err.push_back(std::abs(double(s) - d));
+    }
+    REQUIRE(!ref.empty());
+    const double scale = *std::max_element(ref.begin(), ref.end());
+    const double max_err = *std::max_element(err.begin(), err.end());
+    INFO(name << ": max|float - double| = " << max_err << ", scale = " << scale);
+    CHECK(max_err <= 1e-4 * scale);
+  }
+} // namespace
+
+TEMPLATE_TEST_CASE("Regulators evaluate in single precision", "[physics][regulators][float]", LitimRegulator<>,
+                   BosonicRegulator<>, ExponentialRegulator<>, SmoothedLitimRegulator<>, RationalExpRegulator<>,
+                   RationalExpRegulator<RationalExpOrder<2>>, RationalExpRegulator<RationalExpOrder<5>>,
+                   RationalExpRegulator<RationalExpOrder<12>>, RationalExpRegulator<RationalExpOrder<16>>,
+                   PolynomialExpRegulator<>, PolynomialExpRegulator<PolynomialExpOrder<3>>)
+{
+  using R = TestType;
+  check_function<R>("RB", [](auto k2, auto q2) { return R::RB(k2, q2); });
+  check_function<R>("RBdot", [](auto k2, auto q2) { return R::RBdot(k2, q2); });
+  check_function<R>("RF", [](auto k2, auto q2) { return R::RF(k2, q2); });
+  check_function<R>("RFdot", [](auto k2, auto q2) { return R::RFdot(k2, q2); });
+  if constexpr (requires { R::dq2RB(1.f, 1.f); })
+    check_function<R>("dq2RB", [](auto k2, auto q2) { return R::dq2RB(k2, q2); });
+  if constexpr (requires { R::dq2RF(1.f, 1.f); })
+    check_function<R>("dq2RF", [](auto k2, auto q2) { return R::dq2RF(k2, q2); });
+}
diff --git a/Examples/FourFermi/Four-Fermion.m b/Examples/FourFermi/Four-Fermion.m
index de3edf4ac..a622ae2a5 100644
--- a/Examples/FourFermi/Four-Fermion.m
+++ b/Examples/FourFermi/Four-Fermion.m
@@ -183,16 +183,16 @@
 Finalize[expr_]:=SimplifyAllMomenta[q,ExpandScalarProductsFiniteT[expr]]
 
 
-kernelLambda1=<|"Path"->"qbqqbq","Name"->"lambda1","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda2=<|"Path"->"qbqqbq","Name"->"lambda2","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda3=<|"Path"->"qbqqbq","Name"->"lambda3","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda4=<|"Path"->"qbqqbq","Name"->"lambda4","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda5=<|"Path"->"qbqqbq","Name"->"lambda5","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda6=<|"Path"->"qbqqbq","Name"->"lambda6","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda7=<|"Path"->"qbqqbq","Name"->"lambda7","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda8=<|"Path"->"qbqqbq","Name"->"lambda8","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda9=<|"Path"->"qbqqbq","Name"->"lambda9","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
-kernelLambda10=<|"Path"->"qbqqbq","Name"->"lambda10","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
+kernelLambda1=<|"Path"->"qbqqbq","Name"->"lambda1","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda2=<|"Path"->"qbqqbq","Name"->"lambda2","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda3=<|"Path"->"qbqqbq","Name"->"lambda3","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda4=<|"Path"->"qbqqbq","Name"->"lambda4","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda5=<|"Path"->"qbqqbq","Name"->"lambda5","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda6=<|"Path"->"qbqqbq","Name"->"lambda6","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda7=<|"Path"->"qbqqbq","Name"->"lambda7","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda8=<|"Path"->"qbqqbq","Name"->"lambda8","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda9=<|"Path"->"qbqqbq","Name"->"lambda9","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
+kernelLambda10=<|"Path"->"qbqqbq","Name"->"lambda10","Type"->"Quadratureq0","Angles"->0,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
 
 kernels={kernelLambda1,kernelLambda2,kernelLambda3,kernelLambda4,kernelLambda5,kernelLambda6,kernelLambda7,kernelLambda8,kernelLambda9,kernelLambda10};
 
diff --git a/Examples/QuarkMesonLPAprime/QuarkMesonLPAprime.m b/Examples/QuarkMesonLPAprime/QuarkMesonLPAprime.m
index 266b38ff8..a4fda2819 100644
--- a/Examples/QuarkMesonLPAprime/QuarkMesonLPAprime.m
+++ b/Examples/QuarkMesonLPAprime/QuarkMesonLPAprime.m
@@ -243,12 +243,12 @@
 
 (* ::Input::Initialization:: *)
 (*Anomalous dimensions*)
-kerneletaQ=<|"Path"->"qbq","Name"->"etaQ","Type"->"Quadratureq0","Angles"->1,"d"->4,"AD"->True,"ctype"->"double","Device"->"CPU"|>;
-kerneletaPhi=<|"Path"->"PhiPhi","Name"->"etaPhi","Type"->"Quadratureq0","Angles"->1,"d"->4,"AD"->True,"ctype"->"double","Device"->"CPU"|>;
+kerneletaQ=<|"Path"->"qbq","Name"->"etaQ","Type"->"Quadratureq0","Angles"->1,"d"->4,"AD"->True,"ComputeType"->"double","Device"->"CPU"|>;
+kerneletaPhi=<|"Path"->"PhiPhi","Name"->"etaPhi","Type"->"Quadratureq0","Angles"->1,"d"->4,"AD"->True,"ComputeType"->"double","Device"->"CPU"|>;
 
 (*Mesons*)
-kernelV=<|"Path"->"V","Name"->"V","Type"->"Quadrature","Angles"->0,"d"->3,"AD"->True,"ctype"->"double","Device"->"CPU"|>;
-kernelhPhi=<|"Path"->"Phiqbq","Name"->"hPhi0","Type"->"Quadratureq0","Angles"->1,"d"->4,"AD"->False,"ctype"->"double","Device"->"CPU"|>;
+kernelV=<|"Path"->"V","Name"->"V","Type"->"Quadrature","Angles"->0,"d"->3,"AD"->True,"ComputeType"->"double","Device"->"CPU"|>;
+kernelhPhi=<|"Path"->"Phiqbq","Name"->"hPhi0","Type"->"Quadratureq0","Angles"->1,"d"->4,"AD"->False,"ComputeType"->"double","Device"->"CPU"|>;
 
 kernels={
 kerneletaQ,kerneletaPhi,
diff --git a/Examples/YangMills/Full/Yang-Mills.wl b/Examples/YangMills/Full/Yang-Mills.wl
index 79ae3d68c..3bbb41bb5 100644
--- a/Examples/YangMills/Full/Yang-Mills.wl
+++ b/Examples/YangMills/Full/Yang-Mills.wl
@@ -219,7 +219,7 @@ kernelParameterList = {
 commonKernelOptions = {
   "d" -> 4,
   "AD" -> False,
-  "ctype" -> "double",
+  "ComputeType" -> "double",
   "Device" -> "GPU",
   "Type" -> "double",
   "Parameters" -> kernelParameterList
diff --git a/Examples/YangMills/SP/Yang-Mills.wl b/Examples/YangMills/SP/Yang-Mills.wl
index 58965692b..480052073 100644
--- a/Examples/YangMills/SP/Yang-Mills.wl
+++ b/Examples/YangMills/SP/Yang-Mills.wl
@@ -171,7 +171,7 @@ kernelParameterList = Join[
 commonKernelOptions = {
   "d" -> 4,
   "AD" -> False,
-  "ctype" -> "double",
+  "ComputeType" -> "double",
   "Device" -> "GPU",
   "Type" -> "double",
   "Parameters" -> kernelParameterList,
diff --git a/Examples/YangMills/SP/parameter.toml b/Examples/YangMills/SP/parameter.toml
index afe57e791..bc584956c 100644
--- a/Examples/YangMills/SP/parameter.toml
+++ b/Examples/YangMills/SP/parameter.toml
@@ -17,7 +17,7 @@
   u_upper = -1e-3
 
 [tuning.STI]
-  tune      = false
+  tune      = true
   tol       = 0.0001
   step_size = 0.01
 
@@ -33,7 +33,7 @@
   coarse_factor     = 10
 
 [physical]
-  Lambda       = 100
+  Lambda       = 200
   Nc           = 3
   m2A          = -71.98873562456
   alphaA4      = 0.1