Symbolic SPH kernels & c++ tests

from __future__ import division

import matplotlib.pyplot as plt
import mpmath
import numpy as np
from sympy import *

import shamrock

# Set precision to 32 digits
mpmath.mp.dps = 32

q = symbols("q")

Utilities

def getnorm(w, R):
    c1D = sympify(1) / (2 * integrate(w, (q, 0, R)))
    c2D = sympify(1) / (integrate(2 * pi * q * w, (q, 0, R)))
    c3D = sympify(1) / (integrate(4 * pi * q * q * w, (q, 0, R)))
    return (c1D, c2D, c3D)


def replace_pi_constants(cpp_code):
    """Replace M_PI and M_1_PI with shambase::constants::pi<Tscal>"""
    cpp_code = cpp_code.replace("M_1_PI", "1/shambase::constants::pi<Tscal>")
    cpp_code = cpp_code.replace("M_PI", "shambase::constants::pi<Tscal>")
    return cpp_code


def sympy_to_cpp(expr):
    """Convert a sympy expression to C++ code with proper formatting"""

    import re

    # Expand the expression first
    # expr = expand(expr)
    cpp_code = ccode(expr)

    # Replace pow(q, n) with q*q*q... for powers up to 11
    # for power in range(11, 1, -1):  # From 11 down to 2
    #     q_mult = '*'.join(['q'] * power)
    #     cpp_code = cpp_code.replace(f'pow(q, {power})', q_mult)

    # Replace pow(base, n) with sham::pow_constexpr<n>(base)
    # Use regex to match pow(anything, number) pattern
    def replace_pow_with_constexpr(match):
        base = match.group(1)
        exponent = match.group(2)
        return f"sham::pow_constexpr<{exponent}>({base})"

    cpp_code = re.sub(r"pow\(([^,]+),\s*(\d+)\)", replace_pow_with_constexpr, cpp_code)

    # Wrap divisions in parentheses for clarity
    # Pattern: number/number followed by * (but not already in parentheses)
    # Match something like: 1155.0/512.0*q  -> (1155.0/512.0)*q
    cpp_code = re.sub(r"(\d+\.?\d*)/(\d+\.?\d*)(\*)", r"(\1/\2)\3", cpp_code)

    # Replace pi constants
    cpp_code = replace_pi_constants(cpp_code)

    return cpp_code


def process_power_variables(cpp_exprs, indent=8):
    """
    Self-contained power variable processing:
    1. Extracts sham::pow_constexpr calls from expressions
    2. Replaces them with temporary variable names (t1, t2, ...)
    3. Generates header declarations

    Returns: (modified_exprs, power_header_string)
    """
    spaces = " " * indent

    # Find all sham::pow_constexpr<n>(...) calls
    # Need to handle nested parentheses properly
    def find_pow_constexpr_calls(text):
        """Find all sham::pow_constexpr calls with proper parenthesis matching"""
        calls = []
        pattern = "sham::pow_constexpr<"
        pos = 0

        while True:
            pos = text.find(pattern, pos)
            if pos == -1:
                break

            # Find the power number
            power_start = pos + len(pattern)
            power_end = text.find(">", power_start)
            power = text[power_start:power_end]

            # Find the matching closing parenthesis
            paren_start = text.find("(", power_end)
            if paren_start == -1:
                break

            # Count parentheses to find the matching closing one
            paren_count = 1
            i = paren_start + 1
            while i < len(text) and paren_count > 0:
                if text[i] == "(":
                    paren_count += 1
                elif text[i] == ")":
                    paren_count -= 1
                i += 1

            if paren_count == 0:
                # Found matching closing parenthesis
                full_match = text[pos:i]
                base = text[paren_start + 1 : i - 1]
                calls.append((full_match, power, base))
                pos = i
            else:
                pos += 1

        return calls

    # Dictionary to store unique pow_constexpr calls
    # Key: full match string, Value: (var_name, power, base)
    pow_calls = {}
    var_counter = 1

    vars = {}

    for expr in cpp_exprs:
        for full_match, power, base in find_pow_constexpr_calls(expr):
            vars[full_match] = (power, base)

    # print(vars)

    # sort them lexicographically by match (from the end of the key)
    sorted_vars = sorted(vars.items(), key=lambda x: x[0])
    for full_match, (power, base) in sorted_vars:
        if full_match not in pow_calls:
            var_name = f"t{var_counter}"
            pow_calls[full_match] = (var_name, power, base)
            var_counter += 1

    # Replace pow_constexpr calls with variable names in expressions
    # Sort by length (longest first) to avoid partial replacements
    sorted_pow_calls = sorted(pow_calls.items(), key=lambda x: len(x[0]), reverse=True)
    for full_match, (var_name, power, base) in sorted_pow_calls:
        cpp_exprs = [expr.replace(full_match, var_name) for expr in cpp_exprs]

    # Generate header declarations
    power_header = ""
    for full_match, (var_name, power, base) in pow_calls.items():
        power_header += f"{spaces}    Tscal {var_name} = sham::pow_constexpr<{power}>({base});\n"

    return cpp_exprs, power_header


def generate_function_body(pieces, cpp_exprs, indent=8):
    """Generate the body section with if-else conditions and return statements"""
    spaces = " " * indent
    body = ""

    for i, ((expr, cond), cpp_expr) in enumerate(zip(pieces, cpp_exprs)):
        if cond == True:
            body += f"{spaces}    else\n"
            body += f"{spaces}        return {cpp_expr};\n"
        else:
            if_keyword = "if" if i == 0 else "else if"
            # Convert condition to string
            cond_str = ccode(cond)

            body += f"{spaces}    {if_keyword} ({cond_str}) {{\n"
            body += f"{spaces}        return {cpp_expr};\n"
            body += f"{spaces}    }}"
            if i < len(pieces) - 1 and pieces[i + 1][1] != True:
                body += " "
            else:
                body += "\n"

    return body


def generate_piecewise_function(func_expr, func_name, indent=8):
    """
    Generate C++ code for a piecewise function.
    Works with separate header (declarations) and body (logic) sections.
    """
    spaces = " " * indent

    # Extract pieces from Piecewise object
    pieces = []
    for arg in func_expr.args:
        # Each arg is an ExprCondPair (expr, cond) tuple
        if hasattr(arg, "expr") and hasattr(arg, "cond"):
            pieces.append((arg.expr, arg.cond))
        elif isinstance(arg, tuple) and len(arg) == 2:
            pieces.append(arg)

    # Phase 1: Generate raw C++ expressions for the body
    cpp_exprs = []
    for expr, _ in pieces:
        cpp_exprs.append(sympy_to_cpp(expr))

    # Phase 2: Process power variables (self-contained: detect, replace, generate header)
    cpp_exprs, power_header = process_power_variables(cpp_exprs, indent)

    # Phase 3: Extract constants and replace in body
    # constants = extract_constants(cpp_exprs)
    # cpp_exprs = replace_constants_in_body(cpp_exprs, constants)

    # Phase 4: Generate constant declarations header
    # const_header = generate_constant_header(constants, indent)

    # Phase 5: Generate body (if-else structure)
    body = generate_function_body(pieces, cpp_exprs, indent)

    # Phase 6: Combine everything (constants first, then powers, then body)
    code = f"{spaces}inline static Tscal {func_name}(Tscal q) {{\n"
    # if const_header:
    #    code += const_header
    if power_header:
        code += power_header
    code += body
    code += f"{spaces}}}\n"

    return code


def kernel_to_shamrock(kernel_gen):
    """
    Generate complete C++ kernel definition from SymPy expression

    Args:
        f_func: Function that returns (R, f_expr, name) tuple
        generate_phi: Whether to generate phi_tilde_3d (requires manual derivation)
    """
    R, f_expr, name = kernel_gen()

    # Compute normalization constants
    c_norm_1d, c_norm_2d, c_norm_3d = getnorm(f_expr, R)

    class text_body:
        def __init__(self):
            self.text = ""

        def __call__(self, text=""):
            self.text += text + "\n"

    text = text_body()

    # Generate class header
    text("template<class Tscal>")
    text(f"class KernelDef{name} {{")
    text("    public:")
    text(
        f"    inline static constexpr Tscal Rkern = {ccode(R)}; ///< Compact support radius of the kernel"
    )
    text("    /// default hfact to be used for this kernel")
    text("    inline static constexpr Tscal hfactd = 1.2;")
    text()

    # Normalize constants with pi handling
    norm_1d_str = replace_pi_constants(ccode(c_norm_1d))
    norm_2d_str = replace_pi_constants(ccode(c_norm_2d))
    norm_3d_str = replace_pi_constants(ccode(c_norm_3d))

    text("    /// 1D norm of the kernel")
    text(f"    inline static constexpr Tscal norm_1d = {norm_1d_str};")
    text("    /// 2D norm of the kernel")
    text(f"    inline static constexpr Tscal norm_2d = {norm_2d_str};")
    text("    /// 3D norm of the kernel")
    text(f"    inline static constexpr Tscal norm_3d = {norm_3d_str};")
    text()

    from sympy.polys.polyfuncs import horner

    # Generate f function
    text(generate_piecewise_function(f_expr, "f", indent=4))
    text()

    # Generate df function
    df_expr = diff(f_expr, q)
    text(generate_piecewise_function(df_expr, "df", indent=4))

    # Generate ddf function
    ddf_expr = diff(df_expr, q)
    text(generate_piecewise_function(ddf_expr, "ddf", indent=4))

    # phi_expr
    phi_tilde_3d_expr = integrate(integrate(f_expr * 4 * pi * q * q, q) * (1 / q**2), q).simplify()

    # correct phi to have 0 at infinity
    filter_cond = next(expr for expr, cond in phi_tilde_3d_expr.args if cond == True)
    lim_phi = limit(filter_cond, q, oo)
    phi_tilde_3d_expr -= lim_phi
    phi_tilde_3d_expr = phi_tilde_3d_expr.simplify()

    text(generate_piecewise_function(phi_tilde_3d_expr, "phi_tilde_3d", indent=4))

    # phi_tilde_3d_prime
    phi_tilde_3d_prime_expr = diff(phi_tilde_3d_expr, q)

    text(generate_piecewise_function(phi_tilde_3d_prime_expr, "phi_tilde_3d_prime", indent=4))

    text("};")
    text()

    return {
        "R": R,
        "name": name,
        "norm_1d": c_norm_1d.evalf(),
        "norm_2d": c_norm_2d.evalf(),
        "norm_3d": c_norm_3d.evalf(),
        "Rkern": R,
        "f": f_expr,
        "df": df_expr,
        "ddf": ddf_expr,
        "phi_tilde_3d": phi_tilde_3d_expr,
        "phi_tilde_3d_prime": phi_tilde_3d_prime_expr,
        "text": text.text,
    }

Testing the kernels

def test_kernel(ret, tolerance=1e-12):

    R = ret["R"]
    name = ret["name"]
    norm_1d = ret["norm_1d"]
    norm_2d = ret["norm_2d"]
    norm_3d = ret["norm_3d"]
    Rkern = ret["Rkern"]
    f_expr = ret["f"]
    df_expr = ret["df"]
    ddf_expr = ret["ddf"]
    phi_tilde_3d_expr = ret["phi_tilde_3d"]
    phi_tilde_3d_prime_expr = ret["phi_tilde_3d_prime"]
    text = ret["text"]

    print("------------------------------------------")
    print(f"Testing kernel {name} matching with Shamrock code")
    print("------------------------------------------")

    f = lambdify((q), f_expr, modules=["mpmath"])
    df = lambdify((q), df_expr, modules=["mpmath"])
    ddf = lambdify((q), ddf_expr, modules=["mpmath"])
    phi_tilde_3d = lambdify((q), phi_tilde_3d_expr, modules=["mpmath"])
    phi_tilde_3d_prime = lambdify((q), phi_tilde_3d_prime_expr, modules=["mpmath"])

    print("Testing norms:")
    shamrock_norm_1d = getattr(shamrock.math.sphkernel, f"{name}_norm_1d")()
    shamrock_norm_2d = getattr(shamrock.math.sphkernel, f"{name}_norm_2d")()
    shamrock_norm_3d = getattr(shamrock.math.sphkernel, f"{name}_norm_3d")()

    print(
        f"Shamrock norm_1d = {shamrock_norm_1d}, expected={norm_1d}, delta={abs(shamrock_norm_1d - norm_1d)}"
    )
    print(
        f"Shamrock norm_2d = {shamrock_norm_2d}, expected={norm_2d}, delta={abs(shamrock_norm_2d - norm_2d)}"
    )
    print(
        f"Shamrock norm_3d = {shamrock_norm_3d}, expected={norm_3d}, delta={abs(shamrock_norm_3d - norm_3d)}"
    )

    # test the norm constants down to 1e-12
    assert abs(shamrock_norm_1d - norm_1d) < tolerance
    assert abs(shamrock_norm_2d - norm_2d) < tolerance
    assert abs(shamrock_norm_3d - norm_3d) < tolerance

    print()
    print("Testing kernel radius:")
    shamrock_Rkern = getattr(shamrock.math.sphkernel, f"{name}_Rkern")()

    print(f"Shamrock Rkern = {shamrock_Rkern}, delta={abs(shamrock_Rkern - Rkern)}")
    assert abs(shamrock_Rkern - Rkern) < tolerance

    print()
    print("Testing kernel functions:")
    shamrock_f = getattr(shamrock.math.sphkernel, f"{name}_f")
    shamrock_df = getattr(shamrock.math.sphkernel, f"{name}_df")
    shamrock_ddf = getattr(shamrock.math.sphkernel, f"{name}_ddf")
    shamrock_phi_tilde_3d = getattr(shamrock.math.sphkernel, f"{name}_phi_tilde_3d")
    shamrock_phi_tilde_3d_prime = getattr(shamrock.math.sphkernel, f"{name}_phi_tilde_3d_prime")

    print(f"Shamrock f(q) = {shamrock_f}")
    print(f"Shamrock df(q) = {shamrock_df}")
    print(f"Shamrock ddf(q) = {shamrock_ddf}")
    print(f"Shamrock phi_tilde_3d(q) = {shamrock_phi_tilde_3d}")
    print(f"Shamrock phi_tilde_3d_prime(q) = {shamrock_phi_tilde_3d_prime}")

    q_arr = np.linspace(0, max(1.1 * float(Rkern), 5.0), 1000)
    shamrock_f = [shamrock_f(x) for x in q_arr]
    shamrock_df = [shamrock_df(x) for x in q_arr]
    shamrock_ddf = [shamrock_ddf(x) for x in q_arr]
    shamrock_phi_tilde_3d = [shamrock_phi_tilde_3d(x) for x in q_arr]
    shamrock_phi_tilde_3d_prime = [shamrock_phi_tilde_3d_prime(x) for x in q_arr]

    sympy_f = [f(x) for x in q_arr]
    sympy_df = [df(x) for x in q_arr]
    sympy_ddf = [ddf(x) for x in q_arr]
    sympy_phi_tilde_3d = [phi_tilde_3d(x) for x in q_arr]
    sympy_phi_tilde_3d_prime = [phi_tilde_3d_prime(x) for x in q_arr]

    # compute the absolute error
    abs_err_f = np.max(np.abs(np.array(shamrock_f) - np.array(sympy_f)))
    abs_err_df = np.max(np.abs(np.array(shamrock_df) - np.array(sympy_df)))
    abs_err_ddf = np.max(np.abs(np.array(shamrock_ddf) - np.array(sympy_ddf)))
    abs_err_phi_tilde_3d = np.max(
        np.abs(np.array(shamrock_phi_tilde_3d) - np.array(sympy_phi_tilde_3d))
    )
    abs_err_phi_tilde_3d_prime = np.max(
        np.abs(np.array(shamrock_phi_tilde_3d_prime) - np.array(sympy_phi_tilde_3d_prime))
    )

    print(f"Absolute error f(q) = {abs_err_f}")
    print(f"Absolute error df(q) = {abs_err_df}")
    print(f"Absolute error ddf(q) = {abs_err_ddf}")
    print(f"Absolute error phi_tilde_3d(q) = {abs_err_phi_tilde_3d}")
    print(f"Absolute error phi_tilde_3d_prime(q) = {abs_err_phi_tilde_3d_prime}")

    assert abs_err_f < tolerance
    assert abs_err_df < tolerance * 10
    assert abs_err_ddf < tolerance * 100
    assert abs_err_phi_tilde_3d < tolerance * 100
    assert abs_err_phi_tilde_3d_prime < tolerance * 1000

    print("------------------------------------------")
    print("")

    return {
        "q_arr": q_arr,
        "shamrock_Cf": np.array(shamrock_f) * norm_3d,
        "shamrock_Cdf": np.array(shamrock_df) * norm_3d,
        "shamrock_Cddf": np.array(shamrock_ddf) * norm_3d,
        "shamrock_Cphi_tilde_3d": np.array(shamrock_phi_tilde_3d) * norm_3d,
        "shamrock_Cphi_tilde_3d_prime": np.array(shamrock_phi_tilde_3d_prime) * norm_3d,
    }


def print_kernel_info(ret):
    print("------------------------------------------")
    print(f"Sympy expression for kernel {ret['name']}")
    print("------------------------------------------")
    print(f"f(q)   = {ret['f']}")
    print(f"df(q)  = {ret['df']}")
    print(f"ddf(q) = {ret['ddf']}")
    print(f"phi_tilde_3d(q) = {ret['phi_tilde_3d']}")
    print(f"phi_tilde_3d_prime(q) = {ret['phi_tilde_3d_prime']}")
    print("------------------------------------------")
    print("")


def print_kernel_cpp_code(ret):
    print("------------------------------------------")
    print(f"C++ generated code for kernel {ret['name']}")
    print("------------------------------------------")
    print(ret["text"])
    print("------------------------------------------")
    print("")

All the SPH kernels

Cubic splines kernels (M-series)

def m4():
    Rkern = 2
    R = sympify(Rkern)
    f = Piecewise(
        (sympify(1) / 4 * (R - q) ** 3 - (R / 2 - q) ** 3, q < R / 2),
        (sympify(1) / 4 * (R - q) ** 3, q < R),
        (0, True),
    )
    return (R, f, "M4")


def m5():
    Rkern = sympify(5) / 2
    R = Rkern
    term1 = (R - q) ** 4
    term2 = -5 * (sympify(3) / 5 * R - q) ** 4
    term3 = 10 * (sympify(1) / 5 * R - q) ** 4
    f = Piecewise(
        (term1 + term2 + term3, q < sympify(1) / 5 * R),
        (term1 + term2, q < sympify(3) / 5 * R),
        (term1, q < R),
        (0, True),
    )
    return (R, f, "M5")


def m6():
    Rkern = 3
    R = sympify(Rkern)
    term1 = (R - q) ** 5
    term2 = -6 * (sympify(2) / 3 * R - q) ** 5
    term3 = 15 * (sympify(1) / 3 * R - q) ** 5
    f = Piecewise(
        (term1 + term2 + term3, q < sympify(1) / 3 * R),
        (term1 + term2, q < sympify(2) / 3 * R),
        (term1, q < R),
        (0, True),
    )
    return (R, f, "M6")


def m7():
    Rkern = sympify(7) / 2
    R = Rkern
    term1 = (R - q) ** 6
    term2 = -7 * (sympify(5) / 7 * R - q) ** 6
    term3 = 21 * (sympify(3) / 7 * R - q) ** 6
    term4 = -35 * (sympify(1) / 7 * R - q) ** 6
    f = Piecewise(
        (term1 + term2 + term3 + term4, q < sympify(1) / 7 * R),
        (term1 + term2 + term3, q < sympify(3) / 7 * R),
        (term1 + term2, q < sympify(5) / 7 * R),
        (term1, q < R),
        (0, True),
    )
    return (R, f, "M7")


def m8():
    Rkern = 4
    R = sympify(Rkern)
    term1 = (4 - q) ** 7
    term2 = -8 * (3 - q) ** 7
    term3 = 28 * (2 - q) ** 7
    term4 = -56 * (1 - q) ** 7
    f = Piecewise(
        (term1 + term2 + term3 + term4, q < 1),
        (term1 + term2 + term3, q < 2),
        (term1 + term2, q < 3),
        (term1, q < R),
        (0, True),
    )
    return (R, f, "M8")


def m9():
    Rkern = sympify(9) / 2
    R = Rkern
    term1 = (R - q) ** 8
    term2 = -9 * (sympify(7) / 9 * R - q) ** 8
    term3 = 36 * (sympify(5) / 9 * R - q) ** 8
    term4 = -84 * (sympify(3) / 9 * R - q) ** 8
    term5 = 126 * (sympify(1) / 9 * R - q) ** 8
    f = Piecewise(
        (term1 + term2 + term3 + term4 + term5, q < sympify(1) / 9 * R),
        (term1 + term2 + term3 + term4, q < sympify(3) / 9 * R),
        (term1 + term2 + term3, q < sympify(5) / 9 * R),
        (term1 + term2, q < sympify(7) / 9 * R),
        (term1, q < R),
        (0, True),
    )
    return (R, f, "M9")


def m10():
    Rkern = 5
    R = sympify(Rkern)
    term1 = (R - q) ** 9
    term2 = -10 * (sympify(4) / 5 * R - q) ** 9
    term3 = 45 * (sympify(3) / 5 * R - q) ** 9
    term4 = -120 * (sympify(2) / 5 * R - q) ** 9
    term5 = 210 * (sympify(1) / 5 * R - q) ** 9
    f = Piecewise(
        (term1 + term2 + term3 + term4 + term5, q < sympify(1) / 5 * R),
        (term1 + term2 + term3 + term4, q < sympify(2) / 5 * R),
        (term1 + term2 + term3, q < sympify(3) / 5 * R),
        (term1 + term2, q < sympify(4) / 5 * R),
        (term1, q < R),
        (0, True),
    )
    return (R, f, "M10")

Wendland kernels (C-series)

def c2():
    Rkern = 2
    R = sympify(Rkern)
    f = Piecewise(((1 - q / 2) ** 4 * (1 + 2 * q), q < R), (0, True))
    return (R, f, "C2")


def c4():
    Rkern = 2
    R = sympify(Rkern)
    f = Piecewise(((1 - q / 2) ** 6 * (1 + 3 * q + sympify(35) / 12 * q**2), q < R), (0, True))
    return (R, f, "C4")


def c6():
    Rkern = 2
    R = sympify(Rkern)
    f = Piecewise(
        ((1 - q / 2) ** 8 * (1 + 4 * q + sympify(25) / 4 * q**2 + 4 * q**3), q < R), (0, True)
    )
    return (R, f, "C6")

Helper function to multiply into piecewise without expanding

def multiply_piecewise(piecewise_expr, factor):
    """Multiply a factor into each piece of a Piecewise expression without expanding"""
    new_args = []
    for arg in piecewise_expr.args:
        if hasattr(arg, "expr") and hasattr(arg, "cond"):
            # Multiply the expression part, keep condition unchanged
            new_expr = factor * arg.expr
            new_args.append((new_expr, arg.cond))
        elif isinstance(arg, tuple) and len(arg) == 2:
            new_expr = factor * arg[0]
            new_args.append((new_expr, arg[1]))

    return Piecewise(*new_args)

Helper function to shift and scale a kernel

def shift_scale_kernel(piecewise_expr, shift_val, scale_val):
    return piecewise_fold(
        Piecewise(
            (1, q < shift_val),
            (piecewise_expr.subs({q: (q - shift_val) * scale_val}), q >= shift_val),
        )
    )

Double hump

def m4dh():
    R, f, _ = m4()
    f = multiply_piecewise(f, q * q)
    return (R, f, "M4DH")


def m4dh3():
    R, f, _ = m4()
    f = multiply_piecewise(f, q * q * q)
    return (R, f, "M4DH3")


def m4dh5():
    R, f, _ = m4()
    f = multiply_piecewise(f, q * q * q * q * q)
    return (R, f, "M4DH5")


def m4dh7():
    R, f, _ = m4()
    f = multiply_piecewise(f, q * q * q * q * q * q * q)
    return (R, f, "M4DH7")

M4Shift kernels

def m4shift2():
    R, f, _ = m4()
    # For q < 1: return 1
    # For q >= 1: return M4((q - 1) * 2)
    f_shifted = shift_scale_kernel(f, shift_val=1, scale_val=2)
    return (R, f_shifted, "M4Shift2")


def m4shift4():
    R, f, _ = m4()
    # For q < 1.5: return 1
    # For q >= 1.5: return M4((q - 1.5) * 4)
    f_shifted = shift_scale_kernel(f, shift_val=sympify(3) / 2, scale_val=4)
    return (R, f_shifted, "M4Shift4")


def m4shift8():
    R, f, _ = m4()
    # For q < 1.75: return 1
    # For q >= 1.75: return M4((q - 1.75) * 8)
    f_shifted = shift_scale_kernel(f, shift_val=sympify(7) / 4, scale_val=8)
    return (R, f_shifted, "M4Shift8")


def m4shift16():
    R, f, _ = m4()
    # For q < 1.875: return 1
    # For q >= 1.875: return M4((q - 1.875) * 16)
    f_shifted = shift_scale_kernel(f, shift_val=sympify(15) / 8, scale_val=16)
    return (R, f_shifted, "M4Shift16")

M4 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4
------------------------------------------
f(q)   = Piecewise((-(1 - q)**3 + (2 - q)**3/4, q < 1), ((2 - q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((3*(1 - q)**2 - 3*(2 - q)**2/4, q < 1), (-3*(2 - q)**2/4, q < 2), (0, True))
ddf(q) = Piecewise((9*q/2 - 3, q < 1), (3 - 3*q/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(3*q**3 - 9*q**2 + 20) - 42)/30, q < 1), (pi*(-q**6 + 9*q**5 - 30*q**4 + 40*q**3 - 48*q + 2)/(30*q), q < 2), (-pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(9*q**2 - 18*q) + 2*q*(3*q**3 - 9*q**2 + 20))/30, q < 1), (pi*(-6*q**5 + 45*q**4 - 120*q**3 + 120*q**2 - 48)/(30*q) - pi*(-q**6 + 9*q**5 - 30*q**4 + 40*q**3 - 48*q + 2)/(30*q**2), q < 2), (pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4
------------------------------------------
template<class Tscal>
class KernelDefM4 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 2.0/3.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (10.0/7.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = 1/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1 - q);
        Tscal t2 = sham::pow_constexpr<3>(2 - q);
        if (q < 1) {
            return -t1 + (1.0/4.0)*t2;
        }         else if (q < 2) {
            return (1.0/4.0)*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        if (q < 1) {
            return 3*t1 - (3.0/4.0)*t2;
        }         else if (q < 2) {
            return -(3.0/4.0)*t2;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        if (q < 1) {
            return (9.0/2.0)*q - 3;
        }         else if (q < 2) {
            return 3 - (3.0/2.0)*q;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 1) {
            return (1.0/30.0)*shambase::constants::pi<Tscal>*(t1*(3*t2 - 9*t1 + 20) - 42);
        }         else if (q < 2) {
            return (1.0/30.0)*shambase::constants::pi<Tscal>*(-t5 + 9*t4 - 30*t3 + 40*t2 - 48*q + 2)/q;
        }
        else
            return -shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 1) {
            return (1.0/30.0)*shambase::constants::pi<Tscal>*(t1*(9*t1 - 18*q) + 2*q*(3*t2 - 9*t1 + 20));
        }         else if (q < 2) {
            return (1.0/30.0)*shambase::constants::pi<Tscal>*(-6*t4 + 45*t3 - 120*t2 + 120*t1 - 48)/q - (1.0/30.0)*shambase::constants::pi<Tscal>*(-t5 + 9*t4 - 30*t3 + 40*t2 - 48*q + 2)/t1;
        }
        else
            return shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4 = test_kernel(ret)

------------------------------------------
Testing kernel M4 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.6666666666666666, expected=0.666666666666667, delta=0
Shamrock norm_2d = 0.4547284088339867, expected=0.454728408833987, delta=0
Shamrock norm_3d = 0.3183098861837907, expected=0.318309886183791, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4_f of PyCapsule object at 0x7f5b0176fc30>
Shamrock df(q) = <built-in method M4_df of PyCapsule object at 0x7f5b0176fcc0>
Shamrock ddf(q) = <built-in method M4_ddf of PyCapsule object at 0x7f5b0176fd50>
Shamrock phi_tilde_3d(q) = <built-in method M4_phi_tilde_3d of PyCapsule object at 0x7f5b0176fdb0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4_phi_tilde_3d_prime of PyCapsule object at 0x7f5b0176fde0>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 6.6613381477509392425417900085449e-16
Absolute error ddf(q) = 1.9428902930940239457413554191589e-16
Absolute error phi_tilde_3d(q) = 5.4873597759449205760603571505981e-15
Absolute error phi_tilde_3d_prime(q) = 9.4052927422763077193482903413079e-15
------------------------------------------

M5 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m5)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M5
------------------------------------------
f(q)   = Piecewise((10*(1/2 - q)**4 - 5*(3/2 - q)**4 + (5/2 - q)**4, q < 1/2), (-5*(3/2 - q)**4 + (5/2 - q)**4, q < 3/2), ((5/2 - q)**4, q < 5/2), (0, True))
df(q)  = Piecewise((-40*(1/2 - q)**3 + 20*(3/2 - q)**3 - 4*(5/2 - q)**3, q < 1/2), (20*(3/2 - q)**3 - 4*(5/2 - q)**3, q < 3/2), (-4*(5/2 - q)**3, q < 5/2), (0, True))
ddf(q) = Piecewise((120*(1/2 - q)**2 - 60*(3/2 - q)**2 + 12*(5/2 - q)**2, q < 1/2), (-60*(3/2 - q)**2 + 12*(5/2 - q)**2, q < 3/2), (12*(5/2 - q)**2, q < 5/2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(192*q**6 - 1008*q**4 + 3220*q**2 - 8393)/336, q < 1/2), (pi*(-128*q**7 + 896*q**6 - 2016*q**5 + 560*q**4 + 3080*q**3 - 8386*q - 1)/(336*q), q < 3/2), (pi*(64*q**7 - 896*q**6 + 5040*q**5 - 14000*q**4 + 17500*q**3 - 21875*q + 2185)/(672*q), q < 5/2), (-20*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(1152*q**5 - 4032*q**3 + 6440*q)/336, q < 1/2), (pi*(-896*q**6 + 5376*q**5 - 10080*q**4 + 2240*q**3 + 9240*q**2 - 8386)/(336*q) - pi*(-128*q**7 + 896*q**6 - 2016*q**5 + 560*q**4 + 3080*q**3 - 8386*q - 1)/(336*q**2), q < 3/2), (pi*(448*q**6 - 5376*q**5 + 25200*q**4 - 56000*q**3 + 52500*q**2 - 21875)/(672*q) - pi*(64*q**7 - 896*q**6 + 5040*q**5 - 14000*q**4 + 17500*q**3 - 21875*q + 2185)/(672*q**2), q < 5/2), (20*pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M5
------------------------------------------
template<class Tscal>
class KernelDefM5 {
    public:
    inline static constexpr Tscal Rkern = 5.0/2.0; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 1.0/24.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (96.0/1199.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1.0/20.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<4>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<4>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<4>(5.0/2.0 - q);
        if (q < 1.0/2.0) {
            return 10*t1 - 5*t2 + t3;
        }         else if (q < 3.0/2.0) {
            return -5*t2 + t3;
        }         else if (q < 5.0/2.0) {
            return t3;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<3>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<3>(5.0/2.0 - q);
        if (q < 1.0/2.0) {
            return -40*t1 + 20*t2 - 4*t3;
        }         else if (q < 3.0/2.0) {
            return 20*t2 - 4*t3;
        }         else if (q < 5.0/2.0) {
            return -4*t3;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<2>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<2>(5.0/2.0 - q);
        if (q < 1.0/2.0) {
            return 120*t1 - 60*t2 + 12*t3;
        }         else if (q < 3.0/2.0) {
            return -60*t2 + 12*t3;
        }         else if (q < 5.0/2.0) {
            return 12*t3;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        if (q < 1.0/2.0) {
            return (1.0/336.0)*shambase::constants::pi<Tscal>*(192*t5 - 1008*t3 + 3220*t1 - 8393);
        }         else if (q < 3.0/2.0) {
            return (1.0/336.0)*shambase::constants::pi<Tscal>*(-128*t6 + 896*t5 - 2016*t4 + 560*t3 + 3080*t2 - 8386*q - 1)/q;
        }         else if (q < 5.0/2.0) {
            return (1.0/672.0)*shambase::constants::pi<Tscal>*(64*t6 - 896*t5 + 5040*t4 - 14000*t3 + 17500*t2 - 21875*q + 2185)/q;
        }
        else
            return -20*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        if (q < 1.0/2.0) {
            return (1.0/336.0)*shambase::constants::pi<Tscal>*(1152*t4 - 4032*t2 + 6440*q);
        }         else if (q < 3.0/2.0) {
            return (1.0/336.0)*shambase::constants::pi<Tscal>*(-896*t5 + 5376*t4 - 10080*t3 + 2240*t2 + 9240*t1 - 8386)/q - (1.0/336.0)*shambase::constants::pi<Tscal>*(-128*t6 + 896*t5 - 2016*t4 + 560*t3 + 3080*t2 - 8386*q - 1)/t1;
        }         else if (q < 5.0/2.0) {
            return (1.0/672.0)*shambase::constants::pi<Tscal>*(448*t5 - 5376*t4 + 25200*t3 - 56000*t2 + 52500*t1 - 21875)/q - (1.0/672.0)*shambase::constants::pi<Tscal>*(64*t6 - 896*t5 + 5040*t4 - 14000*t3 + 17500*t2 - 21875*q + 2185)/t1;
        }
        else
            return 20*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m5 = test_kernel(ret)

------------------------------------------
Testing kernel M5 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.041666666666666664, expected=0.0416666666666667, delta=0
Shamrock norm_2d = 0.025486029252413597, expected=0.0254860292524136, delta=0
Shamrock norm_3d = 0.015915494309189534, expected=0.0159154943091895, delta=0

Testing kernel radius:
Shamrock Rkern = 2.5, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M5_f of PyCapsule object at 0x7f5b01780210>
Shamrock df(q) = <built-in method M5_df of PyCapsule object at 0x7f5b017802a0>
Shamrock ddf(q) = <built-in method M5_ddf of PyCapsule object at 0x7f5b01780330>
Shamrock phi_tilde_3d(q) = <built-in method M5_phi_tilde_3d of PyCapsule object at 0x7f5b01780390>
Shamrock phi_tilde_3d_prime(q) = <built-in method M5_phi_tilde_3d_prime of PyCapsule object at 0x7f5b017803c0>
Absolute error f(q) = 1.1020960353947708389243550387522e-14
Absolute error df(q) = 1.6987646535410197641121131341036e-14
Absolute error ddf(q) = 2.7343555911051482541505584418919e-14
Absolute error phi_tilde_3d(q) = 2.1792635850890083577813759428242e-13
Absolute error phi_tilde_3d_prime(q) = 5.1144045415903564046637237745082e-13
------------------------------------------

M6 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m6)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M6
------------------------------------------
f(q)   = Piecewise((15*(1 - q)**5 - 6*(2 - q)**5 + (3 - q)**5, q < 1), (-6*(2 - q)**5 + (3 - q)**5, q < 2), ((3 - q)**5, q < 3), (0, True))
df(q)  = Piecewise((-75*(1 - q)**4 + 30*(2 - q)**4 - 5*(3 - q)**4, q < 1), (30*(2 - q)**4 - 5*(3 - q)**4, q < 2), (-5*(3 - q)**4, q < 3), (0, True))
ddf(q) = Piecewise((300*(1 - q)**3 - 120*(2 - q)**3 + 20*(3 - q)**3, q < 1), (-120*(2 - q)**3 + 20*(3 - q)**3, q < 2), (20*(3 - q)**3, q < 3), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(-5*q**7 + 20*q**6 - 84*q**4 + 308*q**2 - 956)/7, q < 1), (pi*(5*q**8 - 60*q**7 + 280*q**6 - 588*q**5 + 350*q**4 + 476*q**3 - 1892*q - 5)/(14*q), q < 2), (pi*(-q**8 + 20*q**7 - 168*q**6 + 756*q**5 - 1890*q**4 + 2268*q**3 - 2916*q + 507)/(14*q), q < 3), (-120*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(-35*q**6 + 120*q**5 - 336*q**3 + 616*q)/7, q < 1), (pi*(40*q**7 - 420*q**6 + 1680*q**5 - 2940*q**4 + 1400*q**3 + 1428*q**2 - 1892)/(14*q) - pi*(5*q**8 - 60*q**7 + 280*q**6 - 588*q**5 + 350*q**4 + 476*q**3 - 1892*q - 5)/(14*q**2), q < 2), (pi*(-8*q**7 + 140*q**6 - 1008*q**5 + 3780*q**4 - 7560*q**3 + 6804*q**2 - 2916)/(14*q) - pi*(-q**8 + 20*q**7 - 168*q**6 + 756*q**5 - 1890*q**4 + 2268*q**3 - 2916*q + 507)/(14*q**2), q < 3), (120*pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M6
------------------------------------------
template<class Tscal>
class KernelDefM6 {
    public:
    inline static constexpr Tscal Rkern = 3; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 1.0/120.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (7.0/478.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1.0/120.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<5>(1 - q);
        Tscal t2 = sham::pow_constexpr<5>(2 - q);
        Tscal t3 = sham::pow_constexpr<5>(3 - q);
        if (q < 1) {
            return 15*t1 - 6*t2 + t3;
        }         else if (q < 2) {
            return -6*t2 + t3;
        }         else if (q < 3) {
            return t3;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<4>(1 - q);
        Tscal t2 = sham::pow_constexpr<4>(2 - q);
        Tscal t3 = sham::pow_constexpr<4>(3 - q);
        if (q < 1) {
            return -75*t1 + 30*t2 - 5*t3;
        }         else if (q < 2) {
            return 30*t2 - 5*t3;
        }         else if (q < 3) {
            return -5*t3;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1 - q);
        Tscal t2 = sham::pow_constexpr<3>(2 - q);
        Tscal t3 = sham::pow_constexpr<3>(3 - q);
        if (q < 1) {
            return 300*t1 - 120*t2 + 20*t3;
        }         else if (q < 2) {
            return -120*t2 + 20*t3;
        }         else if (q < 3) {
            return 20*t3;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        if (q < 1) {
            return (1.0/7.0)*shambase::constants::pi<Tscal>*(-5*t6 + 20*t5 - 84*t3 + 308*t1 - 956);
        }         else if (q < 2) {
            return (1.0/14.0)*shambase::constants::pi<Tscal>*(5*t7 - 60*t6 + 280*t5 - 588*t4 + 350*t3 + 476*t2 - 1892*q - 5)/q;
        }         else if (q < 3) {
            return (1.0/14.0)*shambase::constants::pi<Tscal>*(-t7 + 20*t6 - 168*t5 + 756*t4 - 1890*t3 + 2268*t2 - 2916*q + 507)/q;
        }
        else
            return -120*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        if (q < 1) {
            return (1.0/7.0)*shambase::constants::pi<Tscal>*(-35*t5 + 120*t4 - 336*t2 + 616*q);
        }         else if (q < 2) {
            return (1.0/14.0)*shambase::constants::pi<Tscal>*(40*t6 - 420*t5 + 1680*t4 - 2940*t3 + 1400*t2 + 1428*t1 - 1892)/q - (1.0/14.0)*shambase::constants::pi<Tscal>*(5*t7 - 60*t6 + 280*t5 - 588*t4 + 350*t3 + 476*t2 - 1892*q - 5)/t1;
        }         else if (q < 3) {
            return (1.0/14.0)*shambase::constants::pi<Tscal>*(-8*t6 + 140*t5 - 1008*t4 + 3780*t3 - 7560*t2 + 6804*t1 - 2916)/q - (1.0/14.0)*shambase::constants::pi<Tscal>*(-t7 + 20*t6 - 168*t5 + 756*t4 - 1890*t3 + 2268*t2 - 2916*q + 507)/t1;
        }
        else
            return 120*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m6 = test_kernel(ret)

------------------------------------------
Testing kernel M6 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.008333333333333333, expected=0.00833333333333333, delta=0
Shamrock norm_2d = 0.00466144184787978, expected=0.00466144184787978, delta=0
Shamrock norm_3d = 0.0026525823848649226, expected=0.00265258238486492, delta=4.33680868994202E-19

Testing kernel radius:
Shamrock Rkern = 3.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M6_f of PyCapsule object at 0x7f5b017807b0>
Shamrock df(q) = <built-in method M6_df of PyCapsule object at 0x7f5b01780840>
Shamrock ddf(q) = <built-in method M6_ddf of PyCapsule object at 0x7f5b017808d0>
Shamrock phi_tilde_3d(q) = <built-in method M6_phi_tilde_3d of PyCapsule object at 0x7f5b01780930>
Shamrock phi_tilde_3d_prime(q) = <built-in method M6_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01780960>
Absolute error f(q) = 7.105427357601002e-14
Absolute error df(q) = 1.4921397450962104e-13
Absolute error ddf(q) = 2.2737367544323206e-13
Absolute error phi_tilde_3d(q) = 3.7781705394188353281542744566612e-12
Absolute error phi_tilde_3d_prime(q) = 4.3695323848306376217663142727157e-12
------------------------------------------

M7 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m7)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M7
------------------------------------------
f(q)   = Piecewise((-35*(1/2 - q)**6 + 21*(3/2 - q)**6 - 7*(5/2 - q)**6 + (7/2 - q)**6, q < 1/2), (21*(3/2 - q)**6 - 7*(5/2 - q)**6 + (7/2 - q)**6, q < 3/2), (-7*(5/2 - q)**6 + (7/2 - q)**6, q < 5/2), ((7/2 - q)**6, q < 7/2), (0, True))
df(q)  = Piecewise((210*(1/2 - q)**5 - 126*(3/2 - q)**5 + 42*(5/2 - q)**5 - 6*(7/2 - q)**5, q < 1/2), (-126*(3/2 - q)**5 + 42*(5/2 - q)**5 - 6*(7/2 - q)**5, q < 3/2), (42*(5/2 - q)**5 - 6*(7/2 - q)**5, q < 5/2), (-6*(7/2 - q)**5, q < 7/2), (0, True))
ddf(q) = Piecewise((-1050*(1/2 - q)**4 + 630*(3/2 - q)**4 - 210*(5/2 - q)**4 + 30*(7/2 - q)**4, q < 1/2), (630*(3/2 - q)**4 - 210*(5/2 - q)**4 + 30*(7/2 - q)**4, q < 3/2), (-210*(5/2 - q)**4 + 30*(7/2 - q)**4, q < 5/2), (30*(7/2 - q)**4, q < 7/2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(-1280*q**8 + 11520*q**6 - 66528*q**4 + 282576*q**2 - 1018341)/1152, q < 1/2), (pi*(3840*q**9 - 34560*q**8 + 103680*q**7 - 53760*q**6 - 235872*q**5 - 10080*q**4 + 1131984*q**3 - 4073409*q + 5)/(4608*q), q < 3/2), (pi*(-768*q**9 + 13824*q**8 - 103680*q**7 + 408576*q**6 - 852768*q**5 + 729792*q**4 + 198576*q**3 - 1948131*q - 29522)/(2304*q), q < 5/2), (pi*(256*q**9 - 6912*q**8 + 80640*q**7 - 526848*q**6 + 2074464*q**5 - 4840416*q**4 + 5647152*q**3 - 7411887*q + 1894081)/(4608*q), q < 7/2), (-840*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(-10240*q**7 + 69120*q**5 - 266112*q**3 + 565152*q)/1152, q < 1/2), (pi*(34560*q**8 - 276480*q**7 + 725760*q**6 - 322560*q**5 - 1179360*q**4 - 40320*q**3 + 3395952*q**2 - 4073409)/(4608*q) - pi*(3840*q**9 - 34560*q**8 + 103680*q**7 - 53760*q**6 - 235872*q**5 - 10080*q**4 + 1131984*q**3 - 4073409*q + 5)/(4608*q**2), q < 3/2), (pi*(-6912*q**8 + 110592*q**7 - 725760*q**6 + 2451456*q**5 - 4263840*q**4 + 2919168*q**3 + 595728*q**2 - 1948131)/(2304*q) - pi*(-768*q**9 + 13824*q**8 - 103680*q**7 + 408576*q**6 - 852768*q**5 + 729792*q**4 + 198576*q**3 - 1948131*q - 29522)/(2304*q**2), q < 5/2), (pi*(2304*q**8 - 55296*q**7 + 564480*q**6 - 3161088*q**5 + 10372320*q**4 - 19361664*q**3 + 16941456*q**2 - 7411887)/(4608*q) - pi*(256*q**9 - 6912*q**8 + 80640*q**7 - 526848*q**6 + 2074464*q**5 - 4840416*q**4 + 5647152*q**3 - 7411887*q + 1894081)/(4608*q**2), q < 7/2), (840*pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M7
------------------------------------------
template<class Tscal>
class KernelDefM7 {
    public:
    inline static constexpr Tscal Rkern = 7.0/2.0; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 1.0/720.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (256.0/113149.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1.0/840.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<6>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<6>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<6>(5.0/2.0 - q);
        Tscal t4 = sham::pow_constexpr<6>(7.0/2.0 - q);
        if (q < 1.0/2.0) {
            return -35*t1 + 21*t2 - 7*t3 + t4;
        }         else if (q < 3.0/2.0) {
            return 21*t2 - 7*t3 + t4;
        }         else if (q < 5.0/2.0) {
            return -7*t3 + t4;
        }         else if (q < 7.0/2.0) {
            return t4;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<5>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<5>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<5>(5.0/2.0 - q);
        Tscal t4 = sham::pow_constexpr<5>(7.0/2.0 - q);
        if (q < 1.0/2.0) {
            return 210*t1 - 126*t2 + 42*t3 - 6*t4;
        }         else if (q < 3.0/2.0) {
            return -126*t2 + 42*t3 - 6*t4;
        }         else if (q < 5.0/2.0) {
            return 42*t3 - 6*t4;
        }         else if (q < 7.0/2.0) {
            return -6*t4;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<4>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<4>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<4>(5.0/2.0 - q);
        Tscal t4 = sham::pow_constexpr<4>(7.0/2.0 - q);
        if (q < 1.0/2.0) {
            return -1050*t1 + 630*t2 - 210*t3 + 30*t4;
        }         else if (q < 3.0/2.0) {
            return 630*t2 - 210*t3 + 30*t4;
        }         else if (q < 5.0/2.0) {
            return -210*t3 + 30*t4;
        }         else if (q < 7.0/2.0) {
            return 30*t4;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        Tscal t8 = sham::pow_constexpr<9>(q);
        if (q < 1.0/2.0) {
            return (1.0/1152.0)*shambase::constants::pi<Tscal>*(-1280*t7 + 11520*t5 - 66528*t3 + 282576*t1 - 1018341);
        }         else if (q < 3.0/2.0) {
            return (1.0/4608.0)*shambase::constants::pi<Tscal>*(3840*t8 - 34560*t7 + 103680*t6 - 53760*t5 - 235872*t4 - 10080*t3 + 1131984*t2 - 4073409*q + 5)/q;
        }         else if (q < 5.0/2.0) {
            return (1.0/2304.0)*shambase::constants::pi<Tscal>*(-768*t8 + 13824*t7 - 103680*t6 + 408576*t5 - 852768*t4 + 729792*t3 + 198576*t2 - 1948131*q - 29522)/q;
        }         else if (q < 7.0/2.0) {
            return (1.0/4608.0)*shambase::constants::pi<Tscal>*(256*t8 - 6912*t7 + 80640*t6 - 526848*t5 + 2074464*t4 - 4840416*t3 + 5647152*t2 - 7411887*q + 1894081)/q;
        }
        else
            return -840*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        Tscal t8 = sham::pow_constexpr<9>(q);
        if (q < 1.0/2.0) {
            return (1.0/1152.0)*shambase::constants::pi<Tscal>*(-10240*t6 + 69120*t4 - 266112*t2 + 565152*q);
        }         else if (q < 3.0/2.0) {
            return (1.0/4608.0)*shambase::constants::pi<Tscal>*(34560*t7 - 276480*t6 + 725760*t5 - 322560*t4 - 1179360*t3 - 40320*t2 + 3395952*t1 - 4073409)/q - (1.0/4608.0)*shambase::constants::pi<Tscal>*(3840*t8 - 34560*t7 + 103680*t6 - 53760*t5 - 235872*t4 - 10080*t3 + 1131984*t2 - 4073409*q + 5)/t1;
        }         else if (q < 5.0/2.0) {
            return (1.0/2304.0)*shambase::constants::pi<Tscal>*(-6912*t7 + 110592*t6 - 725760*t5 + 2451456*t4 - 4263840*t3 + 2919168*t2 + 595728*t1 - 1948131)/q - (1.0/2304.0)*shambase::constants::pi<Tscal>*(-768*t8 + 13824*t7 - 103680*t6 + 408576*t5 - 852768*t4 + 729792*t3 + 198576*t2 - 1948131*q - 29522)/t1;
        }         else if (q < 7.0/2.0) {
            return (1.0/4608.0)*shambase::constants::pi<Tscal>*(2304*t7 - 55296*t6 + 564480*t5 - 3161088*t4 + 10372320*t3 - 19361664*t2 + 16941456*t1 - 7411887)/q - (1.0/4608.0)*shambase::constants::pi<Tscal>*(256*t8 - 6912*t7 + 80640*t6 - 526848*t5 + 2074464*t4 - 4840416*t3 + 5647152*t2 - 7411887*q + 1894081)/t1;
        }
        else
            return 840*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m7 = test_kernel(ret, tolerance=1e-11)

------------------------------------------
Testing kernel M7 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.001388888888888889, expected=0.00138888888888889, delta=0
Shamrock norm_2d = 0.0007201772076028106, expected=0.000720177207602811, delta=0
Shamrock norm_3d = 0.0003789403406949889, expected=0.000378940340694989, delta=0

Testing kernel radius:
Shamrock Rkern = 3.5, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M7_f of PyCapsule object at 0x7f5b01780d20>
Shamrock df(q) = <built-in method M7_df of PyCapsule object at 0x7f5b01780db0>
Shamrock ddf(q) = <built-in method M7_ddf of PyCapsule object at 0x7f5b01780e40>
Shamrock phi_tilde_3d(q) = <built-in method M7_phi_tilde_3d of PyCapsule object at 0x7f5b01780ea0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M7_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01780ed0>
Absolute error f(q) = 7.1169918592659551989348027572597e-13
Absolute error df(q) = 1.5356443651156504336746027078175e-12
Absolute error ddf(q) = 2.5422845301193941132174144420942e-12
Absolute error phi_tilde_3d(q) = 5.0618446198491799591642322417426e-11
Absolute error phi_tilde_3d_prime(q) = 9.1132909954461041619386612291864e-11
------------------------------------------

M8 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m8)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M8
------------------------------------------
f(q)   = Piecewise((-56*(1 - q)**7 + 28*(2 - q)**7 - 8*(3 - q)**7 + (4 - q)**7, q < 1), (28*(2 - q)**7 - 8*(3 - q)**7 + (4 - q)**7, q < 2), (-8*(3 - q)**7 + (4 - q)**7, q < 3), ((4 - q)**7, q < 4), (0, True))
df(q)  = Piecewise((392*(1 - q)**6 - 196*(2 - q)**6 + 56*(3 - q)**6 - 7*(4 - q)**6, q < 1), (-196*(2 - q)**6 + 56*(3 - q)**6 - 7*(4 - q)**6, q < 2), (56*(3 - q)**6 - 7*(4 - q)**6, q < 3), (-7*(4 - q)**6, q < 4), (0, True))
ddf(q) = Piecewise((-2352*(1 - q)**5 + 1176*(2 - q)**5 - 336*(3 - q)**5 + 42*(4 - q)**5, q < 1), (1176*(2 - q)**5 - 336*(3 - q)**5 + 42*(4 - q)**5, q < 2), (-336*(3 - q)**5 + 42*(4 - q)**5, q < 3), (42*(4 - q)**5, q < 4), (0, True))
phi_tilde_3d(q) = Piecewise((2*pi*(q**2*(7*q**7 - 35*q**6 + 240*q**4 - 1512*q**2 + 7248) - 29740)/9, q < 1), (2*pi*(-21*q**10 + 315*q**9 - 1890*q**8 + 5400*q**7 - 5880*q**6 - 2268*q**5 - 2940*q**4 + 37080*q**3 - 148770*q + 14)/(45*q), q < 2), (2*pi*(7*q**10 - 175*q**9 + 1890*q**8 - 11400*q**7 + 41160*q**6 - 86940*q**5 + 91140*q**4 - 16680*q**3 - 130850*q - 7154)/(45*q), q < 3), (2*pi*(-q**10 + 35*q**9 - 540*q**8 + 4800*q**7 - 26880*q**6 + 96768*q**5 - 215040*q**4 + 245760*q**3 - 327680*q + 110944)/(45*q), q < 4), (-6720*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((2*pi*(q**2*(49*q**6 - 210*q**5 + 960*q**3 - 3024*q) + 2*q*(7*q**7 - 35*q**6 + 240*q**4 - 1512*q**2 + 7248))/9, q < 1), (2*pi*(-210*q**9 + 2835*q**8 - 15120*q**7 + 37800*q**6 - 35280*q**5 - 11340*q**4 - 11760*q**3 + 111240*q**2 - 148770)/(45*q) - 2*pi*(-21*q**10 + 315*q**9 - 1890*q**8 + 5400*q**7 - 5880*q**6 - 2268*q**5 - 2940*q**4 + 37080*q**3 - 148770*q + 14)/(45*q**2), q < 2), (2*pi*(70*q**9 - 1575*q**8 + 15120*q**7 - 79800*q**6 + 246960*q**5 - 434700*q**4 + 364560*q**3 - 50040*q**2 - 130850)/(45*q) - 2*pi*(7*q**10 - 175*q**9 + 1890*q**8 - 11400*q**7 + 41160*q**6 - 86940*q**5 + 91140*q**4 - 16680*q**3 - 130850*q - 7154)/(45*q**2), q < 3), (2*pi*(-10*q**9 + 315*q**8 - 4320*q**7 + 33600*q**6 - 161280*q**5 + 483840*q**4 - 860160*q**3 + 737280*q**2 - 327680)/(45*q) - 2*pi*(-q**10 + 35*q**9 - 540*q**8 + 4800*q**7 - 26880*q**6 + 96768*q**5 - 215040*q**4 + 245760*q**3 - 327680*q + 110944)/(45*q**2), q < 4), (6720*pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M8
------------------------------------------
template<class Tscal>
class KernelDefM8 {
    public:
    inline static constexpr Tscal Rkern = 4; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 1.0/5040.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (9.0/29740.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1.0/6720.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<7>(1 - q);
        Tscal t2 = sham::pow_constexpr<7>(2 - q);
        Tscal t3 = sham::pow_constexpr<7>(3 - q);
        Tscal t4 = sham::pow_constexpr<7>(4 - q);
        if (q < 1) {
            return -56*t1 + 28*t2 - 8*t3 + t4;
        }         else if (q < 2) {
            return 28*t2 - 8*t3 + t4;
        }         else if (q < 3) {
            return -8*t3 + t4;
        }         else if (q < 4) {
            return t4;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<6>(1 - q);
        Tscal t2 = sham::pow_constexpr<6>(2 - q);
        Tscal t3 = sham::pow_constexpr<6>(3 - q);
        Tscal t4 = sham::pow_constexpr<6>(4 - q);
        if (q < 1) {
            return 392*t1 - 196*t2 + 56*t3 - 7*t4;
        }         else if (q < 2) {
            return -196*t2 + 56*t3 - 7*t4;
        }         else if (q < 3) {
            return 56*t3 - 7*t4;
        }         else if (q < 4) {
            return -7*t4;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<5>(1 - q);
        Tscal t2 = sham::pow_constexpr<5>(2 - q);
        Tscal t3 = sham::pow_constexpr<5>(3 - q);
        Tscal t4 = sham::pow_constexpr<5>(4 - q);
        if (q < 1) {
            return -2352*t1 + 1176*t2 - 336*t3 + 42*t4;
        }         else if (q < 2) {
            return 1176*t2 - 336*t3 + 42*t4;
        }         else if (q < 3) {
            return -336*t3 + 42*t4;
        }         else if (q < 4) {
            return 42*t4;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<2>(q);
        Tscal t3 = sham::pow_constexpr<3>(q);
        Tscal t4 = sham::pow_constexpr<4>(q);
        Tscal t5 = sham::pow_constexpr<5>(q);
        Tscal t6 = sham::pow_constexpr<6>(q);
        Tscal t7 = sham::pow_constexpr<7>(q);
        Tscal t8 = sham::pow_constexpr<8>(q);
        Tscal t9 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (2.0/9.0)*shambase::constants::pi<Tscal>*(t2*(7*t7 - 35*t6 + 240*t4 - 1512*t2 + 7248) - 29740);
        }         else if (q < 2) {
            return (2.0/45.0)*shambase::constants::pi<Tscal>*(-21*t1 + 315*t9 - 1890*t8 + 5400*t7 - 5880*t6 - 2268*t5 - 2940*t4 + 37080*t3 - 148770*q + 14)/q;
        }         else if (q < 3) {
            return (2.0/45.0)*shambase::constants::pi<Tscal>*(7*t1 - 175*t9 + 1890*t8 - 11400*t7 + 41160*t6 - 86940*t5 + 91140*t4 - 16680*t3 - 130850*q - 7154)/q;
        }         else if (q < 4) {
            return (2.0/45.0)*shambase::constants::pi<Tscal>*(-t1 + 35*t9 - 540*t8 + 4800*t7 - 26880*t6 + 96768*t5 - 215040*t4 + 245760*t3 - 327680*q + 110944)/q;
        }
        else
            return -6720*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<2>(q);
        Tscal t3 = sham::pow_constexpr<3>(q);
        Tscal t4 = sham::pow_constexpr<4>(q);
        Tscal t5 = sham::pow_constexpr<5>(q);
        Tscal t6 = sham::pow_constexpr<6>(q);
        Tscal t7 = sham::pow_constexpr<7>(q);
        Tscal t8 = sham::pow_constexpr<8>(q);
        Tscal t9 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (2.0/9.0)*shambase::constants::pi<Tscal>*(t2*(49*t6 - 210*t5 + 960*t3 - 3024*q) + 2*q*(7*t7 - 35*t6 + 240*t4 - 1512*t2 + 7248));
        }         else if (q < 2) {
            return (2.0/45.0)*shambase::constants::pi<Tscal>*(-210*t9 + 2835*t8 - 15120*t7 + 37800*t6 - 35280*t5 - 11340*t4 - 11760*t3 + 111240*t2 - 148770)/q - (2.0/45.0)*shambase::constants::pi<Tscal>*(-21*t1 + 315*t9 - 1890*t8 + 5400*t7 - 5880*t6 - 2268*t5 - 2940*t4 + 37080*t3 - 148770*q + 14)/t2;
        }         else if (q < 3) {
            return (2.0/45.0)*shambase::constants::pi<Tscal>*(70*t9 - 1575*t8 + 15120*t7 - 79800*t6 + 246960*t5 - 434700*t4 + 364560*t3 - 50040*t2 - 130850)/q - (2.0/45.0)*shambase::constants::pi<Tscal>*(7*t1 - 175*t9 + 1890*t8 - 11400*t7 + 41160*t6 - 86940*t5 + 91140*t4 - 16680*t3 - 130850*q - 7154)/t2;
        }         else if (q < 4) {
            return (2.0/45.0)*shambase::constants::pi<Tscal>*(-10*t9 + 315*t8 - 4320*t7 + 33600*t6 - 161280*t5 + 483840*t4 - 860160*t3 + 737280*t2 - 327680)/q - (2.0/45.0)*shambase::constants::pi<Tscal>*(-t1 + 35*t9 - 540*t8 + 4800*t7 - 26880*t6 + 96768*t5 - 215040*t4 + 245760*t3 - 327680*q + 110944)/t2;
        }
        else
            return 6720*shambase::constants::pi<Tscal>/t2;
    }

};


------------------------------------------

Test the kernel

test_result_m8 = test_kernel(ret, tolerance=1e-10)

------------------------------------------
Testing kernel M8 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.0001984126984126984, expected=0.000198412698412698, delta=0
Shamrock norm_2d = 9.63278068478183e-05, expected=0.0000963278068478183, delta=0
Shamrock norm_3d = 4.736754258687362e-05, expected=0.0000473675425868736, delta=6.77626357803440E-21

Testing kernel radius:
Shamrock Rkern = 4.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M8_f of PyCapsule object at 0x7f5b017812c0>
Shamrock df(q) = <built-in method M8_df of PyCapsule object at 0x7f5b01781350>
Shamrock ddf(q) = <built-in method M8_ddf of PyCapsule object at 0x7f5b017813e0>
Shamrock phi_tilde_3d(q) = <built-in method M8_phi_tilde_3d of PyCapsule object at 0x7f5b01781440>
Shamrock phi_tilde_3d_prime(q) = <built-in method M8_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01781470>
Absolute error f(q) = 5.002220859751105e-12
Absolute error df(q) = 1.1482370609883219e-11
Absolute error ddf(q) = 2.9103830456733704e-11
Absolute error phi_tilde_3d(q) = 0.000000001022101002292260246261593670214
Absolute error phi_tilde_3d_prime(q) = 0.0000000018193685467461683891335635391754
------------------------------------------

M9 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m9)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M9
------------------------------------------
f(q)   = Piecewise((126*(1/2 - q)**8 - 84*(3/2 - q)**8 + 36*(5/2 - q)**8 - 9*(7/2 - q)**8 + (9/2 - q)**8, q < 1/2), (-84*(3/2 - q)**8 + 36*(5/2 - q)**8 - 9*(7/2 - q)**8 + (9/2 - q)**8, q < 3/2), (36*(5/2 - q)**8 - 9*(7/2 - q)**8 + (9/2 - q)**8, q < 5/2), (-9*(7/2 - q)**8 + (9/2 - q)**8, q < 7/2), ((9/2 - q)**8, q < 9/2), (0, True))
df(q)  = Piecewise((-1008*(1/2 - q)**7 + 672*(3/2 - q)**7 - 288*(5/2 - q)**7 + 72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 1/2), (672*(3/2 - q)**7 - 288*(5/2 - q)**7 + 72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 3/2), (-288*(5/2 - q)**7 + 72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 5/2), (72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 7/2), (-8*(9/2 - q)**7, q < 9/2), (0, True))
ddf(q) = Piecewise((7056*(1/2 - q)**6 - 4704*(3/2 - q)**6 + 2016*(5/2 - q)**6 - 504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 1/2), (-4704*(3/2 - q)**6 + 2016*(5/2 - q)**6 - 504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 3/2), (2016*(5/2 - q)**6 - 504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 5/2), (-504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 7/2), (56*(9/2 - q)**6, q < 9/2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(7168*q**10 - 98560*q**8 + 908160*q**6 - 6408864*q**4 + 34283436*q**2 - 157802293)/2816, q < 1/2), (pi*(-28672*q**11 + 315392*q**10 - 1182720*q**9 + 887040*q**8 + 3801600*q**7 + 413952*q**6 - 32199552*q**5 + 36960*q**4 + 171412560*q**3 - 789011388*q - 7)/(14080*q), q < 3/2), (pi*(14336*q**11 - 315392*q**10 + 2956800*q**9 - 15079680*q**8 + 43718400*q**7 - 66646272*q**6 + 43243200*q**5 - 53850720*q**4 + 191620440*q**3 - 792042570*q + 826679)/(14080*q), q < 5/2), (pi*(-4096*q**11 + 135168*q**10 - 1971200*q**9 + 16600320*q**8 - 88281600*q**7 + 302953728*q**6 - 649756800*q**5 + 771149280*q**4 - 324004560*q**3 - 577198820*q - 96829571)/(14080*q), q < 7/2), (pi*(1024*q**11 - 45056*q**10 + 887040*q**9 - 10264320*q**8 + 76982400*q**7 - 387991296*q**6 + 1309470624*q**5 - 2806008480*q**4 + 3156759540*q**3 - 4261625379*q + 1783667601)/(28160*q), q < 9/2), (-60480*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(71680*q**9 - 788480*q**7 + 5448960*q**5 - 25635456*q**3 + 68566872*q)/2816, q < 1/2), (pi*(-315392*q**10 + 3153920*q**9 - 10644480*q**8 + 7096320*q**7 + 26611200*q**6 + 2483712*q**5 - 160997760*q**4 + 147840*q**3 + 514237680*q**2 - 789011388)/(14080*q) - pi*(-28672*q**11 + 315392*q**10 - 1182720*q**9 + 887040*q**8 + 3801600*q**7 + 413952*q**6 - 32199552*q**5 + 36960*q**4 + 171412560*q**3 - 789011388*q - 7)/(14080*q**2), q < 3/2), (pi*(157696*q**10 - 3153920*q**9 + 26611200*q**8 - 120637440*q**7 + 306028800*q**6 - 399877632*q**5 + 216216000*q**4 - 215402880*q**3 + 574861320*q**2 - 792042570)/(14080*q) - pi*(14336*q**11 - 315392*q**10 + 2956800*q**9 - 15079680*q**8 + 43718400*q**7 - 66646272*q**6 + 43243200*q**5 - 53850720*q**4 + 191620440*q**3 - 792042570*q + 826679)/(14080*q**2), q < 5/2), (pi*(-45056*q**10 + 1351680*q**9 - 17740800*q**8 + 132802560*q**7 - 617971200*q**6 + 1817722368*q**5 - 3248784000*q**4 + 3084597120*q**3 - 972013680*q**2 - 577198820)/(14080*q) - pi*(-4096*q**11 + 135168*q**10 - 1971200*q**9 + 16600320*q**8 - 88281600*q**7 + 302953728*q**6 - 649756800*q**5 + 771149280*q**4 - 324004560*q**3 - 577198820*q - 96829571)/(14080*q**2), q < 7/2), (pi*(11264*q**10 - 450560*q**9 + 7983360*q**8 - 82114560*q**7 + 538876800*q**6 - 2327947776*q**5 + 6547353120*q**4 - 11224033920*q**3 + 9470278620*q**2 - 4261625379)/(28160*q) - pi*(1024*q**11 - 45056*q**10 + 887040*q**9 - 10264320*q**8 + 76982400*q**7 - 387991296*q**6 + 1309470624*q**5 - 2806008480*q**4 + 3156759540*q**3 - 4261625379*q + 1783667601)/(28160*q**2), q < 9/2), (60480*pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M9
------------------------------------------
template<class Tscal>
class KernelDefM9 {
    public:
    inline static constexpr Tscal Rkern = 9.0/2.0; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 1.0/40320.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (512.0/14345663.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1.0/60480.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<8>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<8>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<8>(5.0/2.0 - q);
        Tscal t4 = sham::pow_constexpr<8>(7.0/2.0 - q);
        Tscal t5 = sham::pow_constexpr<8>(9.0/2.0 - q);
        if (q < 1.0/2.0) {
            return 126*t1 - 84*t2 + 36*t3 - 9*t4 + t5;
        }         else if (q < 3.0/2.0) {
            return -84*t2 + 36*t3 - 9*t4 + t5;
        }         else if (q < 5.0/2.0) {
            return 36*t3 - 9*t4 + t5;
        }         else if (q < 7.0/2.0) {
            return -9*t4 + t5;
        }         else if (q < 9.0/2.0) {
            return t5;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<7>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<7>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<7>(5.0/2.0 - q);
        Tscal t4 = sham::pow_constexpr<7>(7.0/2.0 - q);
        Tscal t5 = sham::pow_constexpr<7>(9.0/2.0 - q);
        if (q < 1.0/2.0) {
            return -1008*t1 + 672*t2 - 288*t3 + 72*t4 - 8*t5;
        }         else if (q < 3.0/2.0) {
            return 672*t2 - 288*t3 + 72*t4 - 8*t5;
        }         else if (q < 5.0/2.0) {
            return -288*t3 + 72*t4 - 8*t5;
        }         else if (q < 7.0/2.0) {
            return 72*t4 - 8*t5;
        }         else if (q < 9.0/2.0) {
            return -8*t5;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<6>(1.0/2.0 - q);
        Tscal t2 = sham::pow_constexpr<6>(3.0/2.0 - q);
        Tscal t3 = sham::pow_constexpr<6>(5.0/2.0 - q);
        Tscal t4 = sham::pow_constexpr<6>(7.0/2.0 - q);
        Tscal t5 = sham::pow_constexpr<6>(9.0/2.0 - q);
        if (q < 1.0/2.0) {
            return 7056*t1 - 4704*t2 + 2016*t3 - 504*t4 + 56*t5;
        }         else if (q < 3.0/2.0) {
            return -4704*t2 + 2016*t3 - 504*t4 + 56*t5;
        }         else if (q < 5.0/2.0) {
            return 2016*t3 - 504*t4 + 56*t5;
        }         else if (q < 7.0/2.0) {
            return -504*t4 + 56*t5;
        }         else if (q < 9.0/2.0) {
            return 56*t5;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(q);
        Tscal t5 = sham::pow_constexpr<4>(q);
        Tscal t6 = sham::pow_constexpr<5>(q);
        Tscal t7 = sham::pow_constexpr<6>(q);
        Tscal t8 = sham::pow_constexpr<7>(q);
        Tscal t9 = sham::pow_constexpr<8>(q);
        Tscal t10 = sham::pow_constexpr<9>(q);
        if (q < 1.0/2.0) {
            return (1.0/2816.0)*shambase::constants::pi<Tscal>*(7168*t1 - 98560*t9 + 908160*t7 - 6408864*t5 + 34283436*t3 - 157802293);
        }         else if (q < 3.0/2.0) {
            return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-28672*t2 + 315392*t1 - 1182720*t10 + 887040*t9 + 3801600*t8 + 413952*t7 - 32199552*t6 + 36960*t5 + 171412560*t4 - 789011388*q - 7)/q;
        }         else if (q < 5.0/2.0) {
            return (1.0/14080.0)*shambase::constants::pi<Tscal>*(14336*t2 - 315392*t1 + 2956800*t10 - 15079680*t9 + 43718400*t8 - 66646272*t7 + 43243200*t6 - 53850720*t5 + 191620440*t4 - 792042570*q + 826679)/q;
        }         else if (q < 7.0/2.0) {
            return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-4096*t2 + 135168*t1 - 1971200*t10 + 16600320*t9 - 88281600*t8 + 302953728*t7 - 649756800*t6 + 771149280*t5 - 324004560*t4 - 577198820*q - 96829571)/q;
        }         else if (q < 9.0/2.0) {
            return (1.0/28160.0)*shambase::constants::pi<Tscal>*(1024*t2 - 45056*t1 + 887040*t10 - 10264320*t9 + 76982400*t8 - 387991296*t7 + 1309470624*t6 - 2806008480*t5 + 3156759540*t4 - 4261625379*q + 1783667601)/q;
        }
        else
            return -60480*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(q);
        Tscal t5 = sham::pow_constexpr<4>(q);
        Tscal t6 = sham::pow_constexpr<5>(q);
        Tscal t7 = sham::pow_constexpr<6>(q);
        Tscal t8 = sham::pow_constexpr<7>(q);
        Tscal t9 = sham::pow_constexpr<8>(q);
        Tscal t10 = sham::pow_constexpr<9>(q);
        if (q < 1.0/2.0) {
            return (1.0/2816.0)*shambase::constants::pi<Tscal>*(71680*t10 - 788480*t8 + 5448960*t6 - 25635456*t4 + 68566872*q);
        }         else if (q < 3.0/2.0) {
            return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-315392*t1 + 3153920*t10 - 10644480*t9 + 7096320*t8 + 26611200*t7 + 2483712*t6 - 160997760*t5 + 147840*t4 + 514237680*t3 - 789011388)/q - (1.0/14080.0)*shambase::constants::pi<Tscal>*(-28672*t2 + 315392*t1 - 1182720*t10 + 887040*t9 + 3801600*t8 + 413952*t7 - 32199552*t6 + 36960*t5 + 171412560*t4 - 789011388*q - 7)/t3;
        }         else if (q < 5.0/2.0) {
            return (1.0/14080.0)*shambase::constants::pi<Tscal>*(157696*t1 - 3153920*t10 + 26611200*t9 - 120637440*t8 + 306028800*t7 - 399877632*t6 + 216216000*t5 - 215402880*t4 + 574861320*t3 - 792042570)/q - (1.0/14080.0)*shambase::constants::pi<Tscal>*(14336*t2 - 315392*t1 + 2956800*t10 - 15079680*t9 + 43718400*t8 - 66646272*t7 + 43243200*t6 - 53850720*t5 + 191620440*t4 - 792042570*q + 826679)/t3;
        }         else if (q < 7.0/2.0) {
            return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-45056*t1 + 1351680*t10 - 17740800*t9 + 132802560*t8 - 617971200*t7 + 1817722368*t6 - 3248784000*t5 + 3084597120*t4 - 972013680*t3 - 577198820)/q - (1.0/14080.0)*shambase::constants::pi<Tscal>*(-4096*t2 + 135168*t1 - 1971200*t10 + 16600320*t9 - 88281600*t8 + 302953728*t7 - 649756800*t6 + 771149280*t5 - 324004560*t4 - 577198820*q - 96829571)/t3;
        }         else if (q < 9.0/2.0) {
            return (1.0/28160.0)*shambase::constants::pi<Tscal>*(11264*t1 - 450560*t10 + 7983360*t9 - 82114560*t8 + 538876800*t7 - 2327947776*t6 + 6547353120*t5 - 11224033920*t4 + 9470278620*t3 - 4261625379)/q - (1.0/28160.0)*shambase::constants::pi<Tscal>*(1024*t2 - 45056*t1 + 887040*t10 - 10264320*t9 + 76982400*t8 - 387991296*t7 + 1309470624*t6 - 2806008480*t5 + 3156759540*t4 - 4261625379*q + 1783667601)/t3;
        }
        else
            return 60480*shambase::constants::pi<Tscal>/t3;
    }

};


------------------------------------------

Test the kernel

test_result_m9 = test_kernel(ret, tolerance=1e-9)

------------------------------------------
Testing kernel M9 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.48015873015873e-05, expected=0.0000248015873015873, delta=0
Shamrock norm_2d = 1.1360552783520763e-05, expected=0.0000113605527835208, delta=0
Shamrock norm_3d = 5.263060287430401e-06, expected=0.00000526306028743040, delta=0

Testing kernel radius:
Shamrock Rkern = 4.5, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M9_f of PyCapsule object at 0x7f5b01781800>
Shamrock df(q) = <built-in method M9_df of PyCapsule object at 0x7f5b01781890>
Shamrock ddf(q) = <built-in method M9_ddf of PyCapsule object at 0x7f5b01781920>
Shamrock phi_tilde_3d(q) = <built-in method M9_phi_tilde_3d of PyCapsule object at 0x7f5b01781980>
Shamrock phi_tilde_3d_prime(q) = <built-in method M9_phi_tilde_3d_prime of PyCapsule object at 0x7f5b017819b0>
Absolute error f(q) = 1.3116584471113417810615033274653e-10
Absolute error df(q) = 1.6072615378338547344603540902465e-10
Absolute error ddf(q) = 2.873446049366578385242747112273e-10
Absolute error phi_tilde_3d(q) = 0.000000022175857131947026943452188448116
Absolute error phi_tilde_3d_prime(q) = 0.000000039433226222822340380178224869235
------------------------------------------

M10 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m10)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M10
------------------------------------------
f(q)   = Piecewise((210*(1 - q)**9 - 120*(2 - q)**9 + 45*(3 - q)**9 - 10*(4 - q)**9 + (5 - q)**9, q < 1), (-120*(2 - q)**9 + 45*(3 - q)**9 - 10*(4 - q)**9 + (5 - q)**9, q < 2), (45*(3 - q)**9 - 10*(4 - q)**9 + (5 - q)**9, q < 3), (-10*(4 - q)**9 + (5 - q)**9, q < 4), ((5 - q)**9, q < 5), (0, True))
df(q)  = Piecewise((-1890*(1 - q)**8 + 1080*(2 - q)**8 - 405*(3 - q)**8 + 90*(4 - q)**8 - 9*(5 - q)**8, q < 1), (1080*(2 - q)**8 - 405*(3 - q)**8 + 90*(4 - q)**8 - 9*(5 - q)**8, q < 2), (-405*(3 - q)**8 + 90*(4 - q)**8 - 9*(5 - q)**8, q < 3), (90*(4 - q)**8 - 9*(5 - q)**8, q < 4), (-9*(5 - q)**8, q < 5), (0, True))
ddf(q) = Piecewise((15120*(1 - q)**7 - 8640*(2 - q)**7 + 3240*(3 - q)**7 - 720*(4 - q)**7 + 72*(5 - q)**7, q < 1), (-8640*(2 - q)**7 + 3240*(3 - q)**7 - 720*(4 - q)**7 + 72*(5 - q)**7, q < 2), (3240*(3 - q)**7 - 720*(4 - q)**7 + 72*(5 - q)**7, q < 3), (-720*(4 - q)**7 + 72*(5 - q)**7, q < 4), (72*(5 - q)**7, q < 5), (0, True))
phi_tilde_3d(q) = Piecewise((2*pi*(-63*q**11 + 378*q**10 - 3850*q**8 + 37620*q**6 - 291060*q**4 + 1718090*q**2 - 8766690)/33, q < 1), (2*pi*(42*q**12 - 756*q**11 + 5544*q**10 - 20020*q**9 + 31185*q**8 - 3960*q**7 + 38808*q**6 - 316008*q**5 + 10395*q**4 + 1715780*q**3 - 8766564*q - 21)/(33*q), q < 2), (2*pi*(-18*q**12 + 540*q**11 - 7128*q**10 + 53900*q**9 - 253935*q**8 + 756360*q**7 - 1380456*q**6 + 1508760*q**5 - 1510245*q**4 + 2391620*q**3 - 8914020*q + 49131)/(33*q), q < 3), (pi*(9*q**12 - 378*q**11 + 7128*q**10 - 79310*q**9 + 574695*q**8 - 2817540*q**7 + 9363816*q**6 - 20365884*q**5 + 26208765*q**4 - 14702930*q**3 - 8262102*q - 4684707)/(33*q), q < 4), (pi*(-q**12 + 54*q**11 - 1320*q**10 + 19250*q**9 - 185625*q**8 + 1237500*q**7 - 5775000*q**6 + 18562500*q**5 - 38671875*q**4 + 42968750*q**3 - 58593750*q + 28869725)/(33*q), q < 5), (-604800*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((2*pi*(-693*q**10 + 3780*q**9 - 30800*q**7 + 225720*q**5 - 1164240*q**3 + 3436180*q)/33, q < 1), (2*pi*(504*q**11 - 8316*q**10 + 55440*q**9 - 180180*q**8 + 249480*q**7 - 27720*q**6 + 232848*q**5 - 1580040*q**4 + 41580*q**3 + 5147340*q**2 - 8766564)/(33*q) - 2*pi*(42*q**12 - 756*q**11 + 5544*q**10 - 20020*q**9 + 31185*q**8 - 3960*q**7 + 38808*q**6 - 316008*q**5 + 10395*q**4 + 1715780*q**3 - 8766564*q - 21)/(33*q**2), q < 2), (2*pi*(-216*q**11 + 5940*q**10 - 71280*q**9 + 485100*q**8 - 2031480*q**7 + 5294520*q**6 - 8282736*q**5 + 7543800*q**4 - 6040980*q**3 + 7174860*q**2 - 8914020)/(33*q) - 2*pi*(-18*q**12 + 540*q**11 - 7128*q**10 + 53900*q**9 - 253935*q**8 + 756360*q**7 - 1380456*q**6 + 1508760*q**5 - 1510245*q**4 + 2391620*q**3 - 8914020*q + 49131)/(33*q**2), q < 3), (pi*(108*q**11 - 4158*q**10 + 71280*q**9 - 713790*q**8 + 4597560*q**7 - 19722780*q**6 + 56182896*q**5 - 101829420*q**4 + 104835060*q**3 - 44108790*q**2 - 8262102)/(33*q) - pi*(9*q**12 - 378*q**11 + 7128*q**10 - 79310*q**9 + 574695*q**8 - 2817540*q**7 + 9363816*q**6 - 20365884*q**5 + 26208765*q**4 - 14702930*q**3 - 8262102*q - 4684707)/(33*q**2), q < 4), (pi*(-12*q**11 + 594*q**10 - 13200*q**9 + 173250*q**8 - 1485000*q**7 + 8662500*q**6 - 34650000*q**5 + 92812500*q**4 - 154687500*q**3 + 128906250*q**2 - 58593750)/(33*q) - pi*(-q**12 + 54*q**11 - 1320*q**10 + 19250*q**9 - 185625*q**8 + 1237500*q**7 - 5775000*q**6 + 18562500*q**5 - 38671875*q**4 + 42968750*q**3 - 58593750*q + 28869725)/(33*q**2), q < 5), (604800*pi/q**2, True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M10
------------------------------------------
template<class Tscal>
class KernelDefM10 {
    public:
    inline static constexpr Tscal Rkern = 5; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 1.0/362880.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (11.0/2922230.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1.0/604800.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<9>(1 - q);
        Tscal t2 = sham::pow_constexpr<9>(2 - q);
        Tscal t3 = sham::pow_constexpr<9>(3 - q);
        Tscal t4 = sham::pow_constexpr<9>(4 - q);
        Tscal t5 = sham::pow_constexpr<9>(5 - q);
        if (q < 1) {
            return 210*t1 - 120*t2 + 45*t3 - 10*t4 + t5;
        }         else if (q < 2) {
            return -120*t2 + 45*t3 - 10*t4 + t5;
        }         else if (q < 3) {
            return 45*t3 - 10*t4 + t5;
        }         else if (q < 4) {
            return -10*t4 + t5;
        }         else if (q < 5) {
            return t5;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<8>(1 - q);
        Tscal t2 = sham::pow_constexpr<8>(2 - q);
        Tscal t3 = sham::pow_constexpr<8>(3 - q);
        Tscal t4 = sham::pow_constexpr<8>(4 - q);
        Tscal t5 = sham::pow_constexpr<8>(5 - q);
        if (q < 1) {
            return -1890*t1 + 1080*t2 - 405*t3 + 90*t4 - 9*t5;
        }         else if (q < 2) {
            return 1080*t2 - 405*t3 + 90*t4 - 9*t5;
        }         else if (q < 3) {
            return -405*t3 + 90*t4 - 9*t5;
        }         else if (q < 4) {
            return 90*t4 - 9*t5;
        }         else if (q < 5) {
            return -9*t5;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<7>(1 - q);
        Tscal t2 = sham::pow_constexpr<7>(2 - q);
        Tscal t3 = sham::pow_constexpr<7>(3 - q);
        Tscal t4 = sham::pow_constexpr<7>(4 - q);
        Tscal t5 = sham::pow_constexpr<7>(5 - q);
        if (q < 1) {
            return 15120*t1 - 8640*t2 + 3240*t3 - 720*t4 + 72*t5;
        }         else if (q < 2) {
            return -8640*t2 + 3240*t3 - 720*t4 + 72*t5;
        }         else if (q < 3) {
            return 3240*t3 - 720*t4 + 72*t5;
        }         else if (q < 4) {
            return -720*t4 + 72*t5;
        }         else if (q < 5) {
            return 72*t5;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<12>(q);
        Tscal t4 = sham::pow_constexpr<2>(q);
        Tscal t5 = sham::pow_constexpr<3>(q);
        Tscal t6 = sham::pow_constexpr<4>(q);
        Tscal t7 = sham::pow_constexpr<5>(q);
        Tscal t8 = sham::pow_constexpr<6>(q);
        Tscal t9 = sham::pow_constexpr<7>(q);
        Tscal t10 = sham::pow_constexpr<8>(q);
        Tscal t11 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (2.0/33.0)*shambase::constants::pi<Tscal>*(-63*t2 + 378*t1 - 3850*t10 + 37620*t8 - 291060*t6 + 1718090*t4 - 8766690);
        }         else if (q < 2) {
            return (2.0/33.0)*shambase::constants::pi<Tscal>*(42*t3 - 756*t2 + 5544*t1 - 20020*t11 + 31185*t10 - 3960*t9 + 38808*t8 - 316008*t7 + 10395*t6 + 1715780*t5 - 8766564*q - 21)/q;
        }         else if (q < 3) {
            return (2.0/33.0)*shambase::constants::pi<Tscal>*(-18*t3 + 540*t2 - 7128*t1 + 53900*t11 - 253935*t10 + 756360*t9 - 1380456*t8 + 1508760*t7 - 1510245*t6 + 2391620*t5 - 8914020*q + 49131)/q;
        }         else if (q < 4) {
            return (1.0/33.0)*shambase::constants::pi<Tscal>*(9*t3 - 378*t2 + 7128*t1 - 79310*t11 + 574695*t10 - 2817540*t9 + 9363816*t8 - 20365884*t7 + 26208765*t6 - 14702930*t5 - 8262102*q - 4684707)/q;
        }         else if (q < 5) {
            return (1.0/33.0)*shambase::constants::pi<Tscal>*(-t3 + 54*t2 - 1320*t1 + 19250*t11 - 185625*t10 + 1237500*t9 - 5775000*t8 + 18562500*t7 - 38671875*t6 + 42968750*t5 - 58593750*q + 28869725)/q;
        }
        else
            return -604800*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<12>(q);
        Tscal t4 = sham::pow_constexpr<2>(q);
        Tscal t5 = sham::pow_constexpr<3>(q);
        Tscal t6 = sham::pow_constexpr<4>(q);
        Tscal t7 = sham::pow_constexpr<5>(q);
        Tscal t8 = sham::pow_constexpr<6>(q);
        Tscal t9 = sham::pow_constexpr<7>(q);
        Tscal t10 = sham::pow_constexpr<8>(q);
        Tscal t11 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (2.0/33.0)*shambase::constants::pi<Tscal>*(-693*t1 + 3780*t11 - 30800*t9 + 225720*t7 - 1164240*t5 + 3436180*q);
        }         else if (q < 2) {
            return (2.0/33.0)*shambase::constants::pi<Tscal>*(504*t2 - 8316*t1 + 55440*t11 - 180180*t10 + 249480*t9 - 27720*t8 + 232848*t7 - 1580040*t6 + 41580*t5 + 5147340*t4 - 8766564)/q - (2.0/33.0)*shambase::constants::pi<Tscal>*(42*t3 - 756*t2 + 5544*t1 - 20020*t11 + 31185*t10 - 3960*t9 + 38808*t8 - 316008*t7 + 10395*t6 + 1715780*t5 - 8766564*q - 21)/t4;
        }         else if (q < 3) {
            return (2.0/33.0)*shambase::constants::pi<Tscal>*(-216*t2 + 5940*t1 - 71280*t11 + 485100*t10 - 2031480*t9 + 5294520*t8 - 8282736*t7 + 7543800*t6 - 6040980*t5 + 7174860*t4 - 8914020)/q - (2.0/33.0)*shambase::constants::pi<Tscal>*(-18*t3 + 540*t2 - 7128*t1 + 53900*t11 - 253935*t10 + 756360*t9 - 1380456*t8 + 1508760*t7 - 1510245*t6 + 2391620*t5 - 8914020*q + 49131)/t4;
        }         else if (q < 4) {
            return (1.0/33.0)*shambase::constants::pi<Tscal>*(108*t2 - 4158*t1 + 71280*t11 - 713790*t10 + 4597560*t9 - 19722780*t8 + 56182896*t7 - 101829420*t6 + 104835060*t5 - 44108790*t4 - 8262102)/q - (1.0/33.0)*shambase::constants::pi<Tscal>*(9*t3 - 378*t2 + 7128*t1 - 79310*t11 + 574695*t10 - 2817540*t9 + 9363816*t8 - 20365884*t7 + 26208765*t6 - 14702930*t5 - 8262102*q - 4684707)/t4;
        }         else if (q < 5) {
            return (1.0/33.0)*shambase::constants::pi<Tscal>*(-12*t2 + 594*t1 - 13200*t11 + 173250*t10 - 1485000*t9 + 8662500*t8 - 34650000*t7 + 92812500*t6 - 154687500*t5 + 128906250*t4 - 58593750)/q - (1.0/33.0)*shambase::constants::pi<Tscal>*(-t3 + 54*t2 - 1320*t1 + 19250*t11 - 185625*t10 + 1237500*t9 - 5775000*t8 + 18562500*t7 - 38671875*t6 + 42968750*t5 - 58593750*q + 28869725)/t4;
        }
        else
            return 604800*shambase::constants::pi<Tscal>/t4;
    }

};


------------------------------------------

Test the kernel

test_result_m10 = test_kernel(ret, tolerance=1e-8)

------------------------------------------
Testing kernel M10 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.7557319223985893e-06, expected=0.00000275573192239859, delta=0
Shamrock norm_2d = 1.1981975231318881e-06, expected=0.00000119819752313189, delta=0
Shamrock norm_3d = 5.263060287430401e-07, expected=5.26306028743040E-7, delta=1.05879118406788E-22

Testing kernel radius:
Shamrock Rkern = 5.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M10_f of PyCapsule object at 0x7f5b01781d70>
Shamrock df(q) = <built-in method M10_df of PyCapsule object at 0x7f5b01781e00>
Shamrock ddf(q) = <built-in method M10_ddf of PyCapsule object at 0x7f5b01781e90>
Shamrock phi_tilde_3d(q) = <built-in method M10_phi_tilde_3d of PyCapsule object at 0x7f5b01781ef0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M10_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01781f20>
Absolute error f(q) = 1.3969838619232178e-09
Absolute error df(q) = 2.928572939708829e-09
Absolute error ddf(q) = 3.725290298461914e-09
Absolute error phi_tilde_3d(q) = 0.00000052670637671644808483816909891633
Absolute error phi_tilde_3d_prime(q) = 0.00000075282548568718158583115382224247
------------------------------------------

C2 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(c2)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel C2
------------------------------------------
f(q)   = Piecewise(((1 - q/2)**4*(2*q + 1), q < 2), (0, True))
df(q)  = Piecewise((2*(1 - q/2)**4 - 2*(1 - q/2)**3*(2*q + 1), q < 2), (0, True))
ddf(q) = Piecewise((-8*(1 - q/2)**3 + 3*(1 - q/2)**2*(2*q + 1), q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(3*q**5 - 30*q**4 + 112*q**3 - 168*q**2 + 224) - 384)/336, q < 2), (-16*pi/(21*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(15*q**4 - 120*q**3 + 336*q**2 - 336*q) + 2*q*(3*q**5 - 30*q**4 + 112*q**3 - 168*q**2 + 224))/336, q < 2), (16*pi/(21*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel C2
------------------------------------------
template<class Tscal>
class KernelDefC2 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 3.0/4.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (7.0/4.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (21.0/16.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<4>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return t1*(2*q + 1);
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1 - (1.0/2.0)*q);
        Tscal t2 = sham::pow_constexpr<4>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return 2*t2 - 2*t1*(2*q + 1);
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - (1.0/2.0)*q);
        Tscal t2 = sham::pow_constexpr<3>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return -8*t2 + 3*t1*(2*q + 1);
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        if (q < 2) {
            return (1.0/336.0)*shambase::constants::pi<Tscal>*(t1*(3*t4 - 30*t3 + 112*t2 - 168*t1 + 224) - 384);
        }
        else
            return -(16.0/21.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        if (q < 2) {
            return (1.0/336.0)*shambase::constants::pi<Tscal>*(t1*(15*t3 - 120*t2 + 336*t1 - 336*q) + 2*q*(3*t4 - 30*t3 + 112*t2 - 168*t1 + 224));
        }
        else
            return (16.0/21.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_c2 = test_kernel(ret)

------------------------------------------
Testing kernel C2 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.75, expected=0.750000000000000, delta=0
Shamrock norm_2d = 0.5570423008216338, expected=0.557042300821634, delta=1.11022302462516E-16
Shamrock norm_3d = 0.4177817256162253, expected=0.417781725616225, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method C2_f of PyCapsule object at 0x7f5b017822b0>
Shamrock df(q) = <built-in method C2_df of PyCapsule object at 0x7f5b01782340>
Shamrock ddf(q) = <built-in method C2_ddf of PyCapsule object at 0x7f5b017823d0>
Shamrock phi_tilde_3d(q) = <built-in method C2_phi_tilde_3d of PyCapsule object at 0x7f5b01782430>
Shamrock phi_tilde_3d_prime(q) = <built-in method C2_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01782460>
Absolute error f(q) = 4.211561406673559480149857592443e-16
Absolute error df(q) = 4.0217953441708541414143680276179e-16
Absolute error ddf(q) = 1.371769129953306271714661459359e-15
Absolute error phi_tilde_3d(q) = 4.8594086645948423336187829411137e-15
Absolute error phi_tilde_3d_prime(q) = 6.7553673539479592524938178204882e-15
------------------------------------------

C4 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(c4)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel C4
------------------------------------------
f(q)   = Piecewise(((1 - q/2)**6*(35*q**2/12 + 3*q + 1), q < 2), (0, True))
df(q)  = Piecewise(((1 - q/2)**6*(35*q/6 + 3) - 3*(1 - q/2)**5*(35*q**2/12 + 3*q + 1), q < 2), (0, True))
ddf(q) = Piecewise((35*(1 - q/2)**6/6 - 6*(1 - q/2)**5*(35*q/6 + 3) + 15*(1 - q/2)**4*(35*q**2/12 + 3*q + 1)/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(105*q**8 - 1408*q**7 + 7700*q**6 - 21120*q**5 + 26400*q**4 - 29568*q**2 + 42240) - 56320)/63360, q < 2), (-256*pi/(495*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(840*q**7 - 9856*q**6 + 46200*q**5 - 105600*q**4 + 105600*q**3 - 59136*q) + 2*q*(105*q**8 - 1408*q**7 + 7700*q**6 - 21120*q**5 + 26400*q**4 - 29568*q**2 + 42240))/63360, q < 2), (256*pi/(495*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel C4
------------------------------------------
template<class Tscal>
class KernelDefC4 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 27.0/32.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (9.0/4.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (495.0/256.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return t2*((35.0/12.0)*t1 + 3*q + 1);
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<5>(1 - (1.0/2.0)*q);
        Tscal t3 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return t3*((35.0/6.0)*q + 3) - 3*t2*((35.0/12.0)*t1 + 3*q + 1);
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<4>(1 - (1.0/2.0)*q);
        Tscal t3 = sham::pow_constexpr<5>(1 - (1.0/2.0)*q);
        Tscal t4 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return (35.0/6.0)*t4 - 6*t3*((35.0/6.0)*q + 3) + (15.0/2.0)*t2*((35.0/12.0)*t1 + 3*q + 1);
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<4>(q);
        Tscal t3 = sham::pow_constexpr<5>(q);
        Tscal t4 = sham::pow_constexpr<6>(q);
        Tscal t5 = sham::pow_constexpr<7>(q);
        Tscal t6 = sham::pow_constexpr<8>(q);
        if (q < 2) {
            return (1.0/63360.0)*shambase::constants::pi<Tscal>*(t1*(105*t6 - 1408*t5 + 7700*t4 - 21120*t3 + 26400*t2 - 29568*t1 + 42240) - 56320);
        }
        else
            return -(256.0/495.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        if (q < 2) {
            return (1.0/63360.0)*shambase::constants::pi<Tscal>*(t1*(840*t6 - 9856*t5 + 46200*t4 - 105600*t3 + 105600*t2 - 59136*q) + 2*q*(105*t7 - 1408*t6 + 7700*t5 - 21120*t4 + 26400*t3 - 29568*t1 + 42240));
        }
        else
            return (256.0/495.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_c4 = test_kernel(ret)

------------------------------------------
Testing kernel C4 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.84375, expected=0.843750000000000, delta=0
Shamrock norm_2d = 0.716197243913529, expected=0.716197243913529, delta=0
Shamrock norm_3d = 0.6154820064881891, expected=0.615482006488189, delta=1.11022302462516E-16

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method C4_f of PyCapsule object at 0x7f5b01782820>
Shamrock df(q) = <built-in method C4_df of PyCapsule object at 0x7f5b017828b0>
Shamrock ddf(q) = <built-in method C4_ddf of PyCapsule object at 0x7f5b01782940>
Shamrock phi_tilde_3d(q) = <built-in method C4_phi_tilde_3d of PyCapsule object at 0x7f5b017829a0>
Shamrock phi_tilde_3d_prime(q) = <built-in method C4_phi_tilde_3d_prime of PyCapsule object at 0x7f5b017829d0>
Absolute error f(q) = 6.77057790166676168827606387144e-16
Absolute error df(q) = 7.3058119683808204154226284388856e-16
Absolute error ddf(q) = 3.9201297321255701516431654464333e-15
Absolute error phi_tilde_3d(q) = 3.647079135469228542501400203625e-14
Absolute error phi_tilde_3d_prime(q) = 6.0839977038123511945173638615324e-14
------------------------------------------

C6 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(c6)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel C6
------------------------------------------
f(q)   = Piecewise(((1 - q/2)**8*(4*q**3 + 25*q**2/4 + 4*q + 1), q < 2), (0, True))
df(q)  = Piecewise(((1 - q/2)**8*(12*q**2 + 25*q/2 + 4) - 4*(1 - q/2)**7*(4*q**3 + 25*q**2/4 + 4*q + 1), q < 2), (0, True))
ddf(q) = Piecewise(((1 - q/2)**8*(24*q + 25/2) - 8*(1 - q/2)**7*(12*q**2 + 25*q/2 + 4) + 14*(1 - q/2)**6*(4*q**3 + 25*q**2/4 + 4*q + 1), q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(480*q**11 - 8085*q**10 + 58240*q**9 - 229320*q**8 + 512512*q**7 - 560560*q**6 + 549120*q**4 - 768768*q**2 + 931840) - 1003520)/1397760, q < 2), (-512*pi/(1365*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(5280*q**10 - 80850*q**9 + 524160*q**8 - 1834560*q**7 + 3587584*q**6 - 3363360*q**5 + 2196480*q**3 - 1537536*q) + 2*q*(480*q**11 - 8085*q**10 + 58240*q**9 - 229320*q**8 + 512512*q**7 - 560560*q**6 + 549120*q**4 - 768768*q**2 + 931840))/1397760, q < 2), (512*pi/(1365*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel C6
------------------------------------------
template<class Tscal>
class KernelDefC6 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 15.0/16.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (39.0/14.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (1365.0/512.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<8>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return t3*(4*t2 + (25.0/4.0)*t1 + 4*q + 1);
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<7>(1 - (1.0/2.0)*q);
        Tscal t4 = sham::pow_constexpr<8>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return t4*(12*t1 + (25.0/2.0)*q + 4) - 4*t3*(4*t2 + (25.0/4.0)*t1 + 4*q + 1);
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
        Tscal t4 = sham::pow_constexpr<7>(1 - (1.0/2.0)*q);
        Tscal t5 = sham::pow_constexpr<8>(1 - (1.0/2.0)*q);
        if (q < 2) {
            return t5*(24*q + 25.0/2.0) - 8*t4*(12*t1 + (25.0/2.0)*q + 4) + 14*t3*(4*t2 + (25.0/4.0)*t1 + 4*q + 1);
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<4>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        Tscal t8 = sham::pow_constexpr<9>(q);
        if (q < 2) {
            return (1.0/1397760.0)*shambase::constants::pi<Tscal>*(t3*(480*t2 - 8085*t1 + 58240*t8 - 229320*t7 + 512512*t6 - 560560*t5 + 549120*t4 - 768768*t3 + 931840) - 1003520);
        }
        else
            return -(512.0/1365.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(q);
        Tscal t5 = sham::pow_constexpr<4>(q);
        Tscal t6 = sham::pow_constexpr<5>(q);
        Tscal t7 = sham::pow_constexpr<6>(q);
        Tscal t8 = sham::pow_constexpr<7>(q);
        Tscal t9 = sham::pow_constexpr<8>(q);
        Tscal t10 = sham::pow_constexpr<9>(q);
        if (q < 2) {
            return (1.0/1397760.0)*shambase::constants::pi<Tscal>*(t3*(5280*t1 - 80850*t10 + 524160*t9 - 1834560*t8 + 3587584*t7 - 3363360*t6 + 2196480*t4 - 1537536*q) + 2*q*(480*t2 - 8085*t1 + 58240*t10 - 229320*t9 + 512512*t8 - 560560*t7 + 549120*t5 - 768768*t3 + 931840));
        }
        else
            return (512.0/1365.0)*shambase::constants::pi<Tscal>/t3;
    }

};


------------------------------------------

Test the kernel

test_result_c6 = test_kernel(ret)

------------------------------------------
Testing kernel C6 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.9375, expected=0.937500000000000, delta=0
Shamrock norm_2d = 0.886720397226274, expected=0.886720397226274, delta=0
Shamrock norm_3d = 0.8486191301579575, expected=0.848619130157958, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method C6_f of PyCapsule object at 0x7f5b01782dc0>
Shamrock df(q) = <built-in method C6_df of PyCapsule object at 0x7f5b01782e50>
Shamrock ddf(q) = <built-in method C6_ddf of PyCapsule object at 0x7f5b01782ee0>
Shamrock phi_tilde_3d(q) = <built-in method C6_phi_tilde_3d of PyCapsule object at 0x7f5b01782f10>
Shamrock phi_tilde_3d_prime(q) = <built-in method C6_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01782f40>
Absolute error f(q) = 8.0028218067564718410675074309679e-16
Absolute error df(q) = 1.2403417665354433235366686954591e-15
Absolute error ddf(q) = 8.6752247783722774983167486346621e-15
Absolute error phi_tilde_3d(q) = 1.6291150200008408805291652636617e-13
Absolute error phi_tilde_3d_prime(q) = 3.4079679145226729610704815366427e-13
------------------------------------------

M4DH Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4dh)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4DH
------------------------------------------
f(q)   = Piecewise((q**2*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**2*(2 - q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((q**2*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 2*q*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**2*(2 - q)**2/4 + q*(2 - q)**3/2, q < 2), (0, True))
ddf(q) = Piecewise((q**2*(9*q/2 - 3) + 4*q*(3*(1 - q)**2 - 3*(2 - q)**2/4) - 2*(1 - q)**3 + (2 - q)**3/2, q < 1), (-3*q**2*(2*q - 4)/4 - 3*q*(2 - q)**2 + (2 - q)**3/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**4*(15*q**3 - 40*q**2 + 56) - 248)/280, q < 1), (pi*(-5*q**8 + 40*q**7 - 112*q**6 + 112*q**5 - 256*q + 4)/(280*q), q < 2), (-9*pi/(10*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**4*(45*q**2 - 80*q) + 4*q**3*(15*q**3 - 40*q**2 + 56))/280, q < 1), (pi*(-40*q**7 + 280*q**6 - 672*q**5 + 560*q**4 - 256)/(280*q) - pi*(-5*q**8 + 40*q**7 - 112*q**6 + 112*q**5 - 256*q + 4)/(280*q**2), q < 2), (9*pi/(10*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4DH
------------------------------------------
template<class Tscal>
class KernelDefM4DH {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 2;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (70.0/31.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (10.0/9.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(1 - q);
        Tscal t3 = sham::pow_constexpr<3>(2 - q);
        if (q < 1) {
            return t1*(-t2 + (1.0/4.0)*t3);
        }         else if (q < 2) {
            return (1.0/4.0)*t1*t3;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(1 - q);
        Tscal t5 = sham::pow_constexpr<3>(2 - q);
        if (q < 1) {
            return t3*(3*t1 - (3.0/4.0)*t2) + 2*q*(-t4 + (1.0/4.0)*t5);
        }         else if (q < 2) {
            return -(3.0/4.0)*t3*t2 + (1.0/2.0)*q*t5;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(1 - q);
        Tscal t5 = sham::pow_constexpr<3>(2 - q);
        if (q < 1) {
            return t3*((9.0/2.0)*q - 3) + 4*q*(3*t1 - (3.0/4.0)*t2) - 2*t4 + (1.0/2.0)*t5;
        }         else if (q < 2) {
            return -(3.0/4.0)*t3*(2*q - 4) - 3*q*t2 + (1.0/2.0)*t5;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        if (q < 1) {
            return (1.0/280.0)*shambase::constants::pi<Tscal>*(t3*(15*t2 - 40*t1 + 56) - 248);
        }         else if (q < 2) {
            return (1.0/280.0)*shambase::constants::pi<Tscal>*(-5*t7 + 40*t6 - 112*t5 + 112*t4 - 256*q + 4)/q;
        }
        else
            return -(9.0/10.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        if (q < 1) {
            return (1.0/280.0)*shambase::constants::pi<Tscal>*(t3*(45*t1 - 80*q) + 4*t2*(15*t2 - 40*t1 + 56));
        }         else if (q < 2) {
            return (1.0/280.0)*shambase::constants::pi<Tscal>*(-40*t6 + 280*t5 - 672*t4 + 560*t3 - 256)/q - (1.0/280.0)*shambase::constants::pi<Tscal>*(-5*t7 + 40*t6 - 112*t5 + 112*t4 - 256*q + 4)/t1;
        }
        else
            return (9.0/10.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4dh = test_kernel(ret)

------------------------------------------
Testing kernel M4DH matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.0, expected=2.00000000000000, delta=0
Shamrock norm_2d = 0.7187642591246886, expected=0.718764259124689, delta=1.11022302462516E-16
Shamrock norm_3d = 0.35367765131532297, expected=0.353677651315323, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4DH_f of PyCapsule object at 0x7f5b017832d0>
Shamrock df(q) = <built-in method M4DH_df of PyCapsule object at 0x7f5b01783360>
Shamrock ddf(q) = <built-in method M4DH_ddf of PyCapsule object at 0x7f5b01783390>
Shamrock phi_tilde_3d(q) = <built-in method M4DH_phi_tilde_3d of PyCapsule object at 0x7f5b017833c0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH_phi_tilde_3d_prime of PyCapsule object at 0x7f5b017833f0>
Absolute error f(q) = 1.2402696483592525566706857128555e-16
Absolute error df(q) = 5.1251967243981971244483580440983e-16
Absolute error ddf(q) = 8.1537932702989126670072438081336e-16
Absolute error phi_tilde_3d(q) = 7.8656660600690683965281955112701e-15
Absolute error phi_tilde_3d_prime(q) = 3.3207246830111135096394910341455e-14
------------------------------------------

M4DH3 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4dh3)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4DH3
------------------------------------------
f(q)   = Piecewise((q**3*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**3*(2 - q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((q**3*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 3*q**2*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**3*(2 - q)**2/4 + 3*q**2*(2 - q)**3/4, q < 2), (0, True))
ddf(q) = Piecewise((q**3*(9*q/2 - 3) + 6*q**2*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 6*q*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**3*(2*q - 4)/4 - 9*q**2*(2 - q)**2/2 + 3*q*(2 - q)**3/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**5*(35*q**3 - 90*q**2 + 112) - 756)/840, q < 1), (pi*(-35*q**9 + 270*q**8 - 720*q**7 + 672*q**6 - 2304*q + 20)/(2520*q), q < 2), (-127*pi/(126*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**5*(105*q**2 - 180*q) + 5*q**4*(35*q**3 - 90*q**2 + 112))/840, q < 1), (pi*(-315*q**8 + 2160*q**7 - 5040*q**6 + 4032*q**5 - 2304)/(2520*q) - pi*(-35*q**9 + 270*q**8 - 720*q**7 + 672*q**6 - 2304*q + 20)/(2520*q**2), q < 2), (127*pi/(126*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4DH3
------------------------------------------
template<class Tscal>
class KernelDefM4DH3 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 70.0/31.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (20.0/9.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (126.0/127.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1 - q);
        Tscal t2 = sham::pow_constexpr<3>(2 - q);
        Tscal t3 = sham::pow_constexpr<3>(q);
        if (q < 1) {
            return t3*(-t1 + (1.0/4.0)*t2);
        }         else if (q < 2) {
            return (1.0/4.0)*t3*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(1 - q);
        Tscal t5 = sham::pow_constexpr<3>(2 - q);
        Tscal t6 = sham::pow_constexpr<3>(q);
        if (q < 1) {
            return t6*(3*t1 - (3.0/4.0)*t2) + 3*t3*(-t4 + (1.0/4.0)*t5);
        }         else if (q < 2) {
            return -(3.0/4.0)*t6*t2 + (3.0/4.0)*t3*t5;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(1 - q);
        Tscal t5 = sham::pow_constexpr<3>(2 - q);
        Tscal t6 = sham::pow_constexpr<3>(q);
        if (q < 1) {
            return t6*((9.0/2.0)*q - 3) + 6*t3*(3*t1 - (3.0/4.0)*t2) + 6*q*(-t4 + (1.0/4.0)*t5);
        }         else if (q < 2) {
            return -(3.0/4.0)*t6*(2*q - 4) - (9.0/2.0)*t3*t2 + (3.0/2.0)*q*t5;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<5>(q);
        Tscal t4 = sham::pow_constexpr<6>(q);
        Tscal t5 = sham::pow_constexpr<7>(q);
        Tscal t6 = sham::pow_constexpr<8>(q);
        Tscal t7 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (1.0/840.0)*shambase::constants::pi<Tscal>*(t3*(35*t2 - 90*t1 + 112) - 756);
        }         else if (q < 2) {
            return (1.0/2520.0)*shambase::constants::pi<Tscal>*(-35*t7 + 270*t6 - 720*t5 + 672*t4 - 2304*q + 20)/q;
        }
        else
            return -(127.0/126.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        Tscal t8 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (1.0/840.0)*shambase::constants::pi<Tscal>*(t4*(105*t1 - 180*q) + 5*t3*(35*t2 - 90*t1 + 112));
        }         else if (q < 2) {
            return (1.0/2520.0)*shambase::constants::pi<Tscal>*(-315*t7 + 2160*t6 - 5040*t5 + 4032*t4 - 2304)/q - (1.0/2520.0)*shambase::constants::pi<Tscal>*(-35*t8 + 270*t7 - 720*t6 + 672*t5 - 2304*q + 20)/t1;
        }
        else
            return (127.0/126.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4dh3 = test_kernel(ret)

------------------------------------------
Testing kernel M4DH3 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.2580645161290325, expected=2.25806451612903, delta=4.44089209850063E-16
Shamrock norm_2d = 0.7073553026306459, expected=0.707355302630646, delta=0
Shamrock norm_3d = 0.31580350912722543, expected=0.315803509127225, delta=5.55111512312578E-17

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4DH3_f of PyCapsule object at 0x7f5b01783660>
Shamrock df(q) = <built-in method M4DH3_df of PyCapsule object at 0x7f5b01783690>
Shamrock ddf(q) = <built-in method M4DH3_ddf of PyCapsule object at 0x7f5b017836c0>
Shamrock phi_tilde_3d(q) = <built-in method M4DH3_phi_tilde_3d of PyCapsule object at 0x7f5b017836f0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH3_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01783720>
Absolute error f(q) = 1.2218718934561889186295748586557e-16
Absolute error df(q) = 8.8551048198837049750116589405827e-16
Absolute error ddf(q) = 1.0260752842383410777313944331012e-15
Absolute error phi_tilde_3d(q) = 1.5383844323012062258140079552524e-14
Absolute error phi_tilde_3d_prime(q) = 4.5790021856924746178274336277744e-14
------------------------------------------

M4DH5 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4dh5)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4DH5
------------------------------------------
f(q)   = Piecewise((q**5*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**5*(2 - q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((q**5*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 5*q**4*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**5*(2 - q)**2/4 + 5*q**4*(2 - q)**3/4, q < 2), (0, True))
ddf(q) = Piecewise((q**5*(9*q/2 - 3) + 10*q**4*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 20*q**3*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**5*(2*q - 4)/4 - 15*q**4*(2 - q)**2/2 + 5*q**3*(2 - q)**3, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**7*(63*q**3 - 154*q**2 + 165) - 2805)/2310, q < 1), (pi*(-21*q**11 + 154*q**10 - 385*q**9 + 330*q**8 - 2816*q + 7)/(2310*q), q < 2), (-511*pi/(330*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**7*(189*q**2 - 308*q) + 7*q**6*(63*q**3 - 154*q**2 + 165))/2310, q < 1), (pi*(-231*q**10 + 1540*q**9 - 3465*q**8 + 2640*q**7 - 2816)/(2310*q) - pi*(-21*q**11 + 154*q**10 - 385*q**9 + 330*q**8 - 2816*q + 7)/(2310*q**2), q < 2), (511*pi/(330*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4DH5
------------------------------------------
template<class Tscal>
class KernelDefM4DH5 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 252.0/127.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (28.0/17.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (330.0/511.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1 - q);
        Tscal t2 = sham::pow_constexpr<3>(2 - q);
        Tscal t3 = sham::pow_constexpr<5>(q);
        if (q < 1) {
            return t3*(-t1 + (1.0/4.0)*t2);
        }         else if (q < 2) {
            return (1.0/4.0)*t3*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<3>(1 - q);
        Tscal t4 = sham::pow_constexpr<3>(2 - q);
        Tscal t5 = sham::pow_constexpr<4>(q);
        Tscal t6 = sham::pow_constexpr<5>(q);
        if (q < 1) {
            return t6*(3*t1 - (3.0/4.0)*t2) + 5*t5*(-t3 + (1.0/4.0)*t4);
        }         else if (q < 2) {
            return -(3.0/4.0)*t6*t2 + (5.0/4.0)*t5*t4;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<3>(1 - q);
        Tscal t4 = sham::pow_constexpr<3>(2 - q);
        Tscal t5 = sham::pow_constexpr<3>(q);
        Tscal t6 = sham::pow_constexpr<4>(q);
        Tscal t7 = sham::pow_constexpr<5>(q);
        if (q < 1) {
            return t7*((9.0/2.0)*q - 3) + 10*t6*(3*t1 - (3.0/4.0)*t2) + 20*t5*(-t3 + (1.0/4.0)*t4);
        }         else if (q < 2) {
            return -(3.0/4.0)*t7*(2*q - 4) - (15.0/2.0)*t6*t2 + 5*t5*t4;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(q);
        Tscal t5 = sham::pow_constexpr<7>(q);
        Tscal t6 = sham::pow_constexpr<8>(q);
        Tscal t7 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (1.0/2310.0)*shambase::constants::pi<Tscal>*(t5*(63*t4 - 154*t3 + 165) - 2805);
        }         else if (q < 2) {
            return (1.0/2310.0)*shambase::constants::pi<Tscal>*(-21*t2 + 154*t1 - 385*t7 + 330*t6 - 2816*q + 7)/q;
        }
        else
            return -(511.0/330.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<2>(q);
        Tscal t4 = sham::pow_constexpr<3>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        Tscal t8 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (1.0/2310.0)*shambase::constants::pi<Tscal>*(t6*(189*t3 - 308*q) + 7*t5*(63*t4 - 154*t3 + 165));
        }         else if (q < 2) {
            return (1.0/2310.0)*shambase::constants::pi<Tscal>*(-231*t1 + 1540*t8 - 3465*t7 + 2640*t6 - 2816)/q - (1.0/2310.0)*shambase::constants::pi<Tscal>*(-21*t2 + 154*t1 - 385*t8 + 330*t7 - 2816*q + 7)/t3;
        }
        else
            return (511.0/330.0)*shambase::constants::pi<Tscal>/t3;
    }

};


------------------------------------------

Test the kernel

test_result_m4dh5 = test_kernel(ret)

------------------------------------------
Testing kernel M4DH5 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 1.984251968503937, expected=1.98425196850394, delta=0
Shamrock norm_2d = 0.5242751066556552, expected=0.524275106655655, delta=0
Shamrock norm_3d = 0.20556215741810357, expected=0.205562157418104, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4DH5_f of PyCapsule object at 0x7f5b01783990>
Shamrock df(q) = <built-in method M4DH5_df of PyCapsule object at 0x7f5b017839f0>
Shamrock ddf(q) = <built-in method M4DH5_ddf of PyCapsule object at 0x7f5b01783a20>
Shamrock phi_tilde_3d(q) = <built-in method M4DH5_phi_tilde_3d of PyCapsule object at 0x7f5b01783a50>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH5_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01783a80>
Absolute error f(q) = 1.1325334002180934793019450988886e-16
Absolute error df(q) = 3.5314759513233869425368604924612e-15
Absolute error ddf(q) = 3.6315357589065079464569273461409e-15
Absolute error phi_tilde_3d(q) = 3.0290833162156024982358814506602e-14
Absolute error phi_tilde_3d_prime(q) = 2.1751989573652618482335489019732e-13
------------------------------------------

M4DH7 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4dh7)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4DH7
------------------------------------------
f(q)   = Piecewise((q**7*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**7*(2 - q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((q**7*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 7*q**6*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**7*(2 - q)**2/4 + 7*q**6*(2 - q)**3/4, q < 2), (0, True))
ddf(q) = Piecewise((q**7*(9*q/2 - 3) + 14*q**6*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 42*q**5*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**7*(2*q - 4)/4 - 21*q**6*(2 - q)**2/2 + 21*q**5*(2 - q)**3/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**9*(495*q**3 - 1170*q**2 + 1144) - 53196)/25740, q < 1), (pi*(-165*q**13 + 1170*q**12 - 2808*q**11 + 2288*q**10 - 53248*q + 36)/(25740*q), q < 2), (-2047*pi/(715*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**9*(1485*q**2 - 2340*q) + 9*q**8*(495*q**3 - 1170*q**2 + 1144))/25740, q < 1), (pi*(-2145*q**12 + 14040*q**11 - 30888*q**10 + 22880*q**9 - 53248)/(25740*q) - pi*(-165*q**13 + 1170*q**12 - 2808*q**11 + 2288*q**10 - 53248*q + 36)/(25740*q**2), q < 2), (2047*pi/(715*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4DH7
------------------------------------------
template<class Tscal>
class KernelDefM4DH7 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 660.0/511.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (30.0/31.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (715.0/2047.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(1 - q);
        Tscal t2 = sham::pow_constexpr<3>(2 - q);
        Tscal t3 = sham::pow_constexpr<7>(q);
        if (q < 1) {
            return t3*(-t1 + (1.0/4.0)*t2);
        }         else if (q < 2) {
            return (1.0/4.0)*t3*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<3>(1 - q);
        Tscal t4 = sham::pow_constexpr<3>(2 - q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        Tscal t6 = sham::pow_constexpr<7>(q);
        if (q < 1) {
            return t6*(3*t1 - (3.0/4.0)*t2) + 7*t5*(-t3 + (1.0/4.0)*t4);
        }         else if (q < 2) {
            return -(3.0/4.0)*t6*t2 + (7.0/4.0)*t5*t4;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(1 - q);
        Tscal t2 = sham::pow_constexpr<2>(2 - q);
        Tscal t3 = sham::pow_constexpr<3>(1 - q);
        Tscal t4 = sham::pow_constexpr<3>(2 - q);
        Tscal t5 = sham::pow_constexpr<5>(q);
        Tscal t6 = sham::pow_constexpr<6>(q);
        Tscal t7 = sham::pow_constexpr<7>(q);
        if (q < 1) {
            return t7*((9.0/2.0)*q - 3) + 14*t6*(3*t1 - (3.0/4.0)*t2) + 42*t5*(-t3 + (1.0/4.0)*t4);
        }         else if (q < 2) {
            return -(3.0/4.0)*t7*(2*q - 4) - (21.0/2.0)*t6*t2 + (21.0/2.0)*t5*t4;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<12>(q);
        Tscal t4 = sham::pow_constexpr<13>(q);
        Tscal t5 = sham::pow_constexpr<2>(q);
        Tscal t6 = sham::pow_constexpr<3>(q);
        Tscal t7 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (1.0/25740.0)*shambase::constants::pi<Tscal>*(t7*(495*t6 - 1170*t5 + 1144) - 53196);
        }         else if (q < 2) {
            return (1.0/25740.0)*shambase::constants::pi<Tscal>*(-165*t4 + 1170*t3 - 2808*t2 + 2288*t1 - 53248*q + 36)/q;
        }
        else
            return -(2047.0/715.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<10>(q);
        Tscal t2 = sham::pow_constexpr<11>(q);
        Tscal t3 = sham::pow_constexpr<12>(q);
        Tscal t4 = sham::pow_constexpr<13>(q);
        Tscal t5 = sham::pow_constexpr<2>(q);
        Tscal t6 = sham::pow_constexpr<3>(q);
        Tscal t7 = sham::pow_constexpr<8>(q);
        Tscal t8 = sham::pow_constexpr<9>(q);
        if (q < 1) {
            return (1.0/25740.0)*shambase::constants::pi<Tscal>*(t8*(1485*t5 - 2340*q) + 9*t7*(495*t6 - 1170*t5 + 1144));
        }         else if (q < 2) {
            return (1.0/25740.0)*shambase::constants::pi<Tscal>*(-2145*t3 + 14040*t2 - 30888*t1 + 22880*t8 - 53248)/q - (1.0/25740.0)*shambase::constants::pi<Tscal>*(-165*t4 + 1170*t3 - 2808*t2 + 2288*t1 - 53248*q + 36)/t5;
        }
        else
            return (2047.0/715.0)*shambase::constants::pi<Tscal>/t5;
    }

};


------------------------------------------

Test the kernel

test_result_m4dh7 = test_kernel(ret)

------------------------------------------
Testing kernel M4DH7 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 1.2915851272015655, expected=1.29158512720157, delta=0
Shamrock norm_2d = 0.3080418253391523, expected=0.308041825339152, delta=5.55111512312578E-17
Shamrock norm_3d = 0.11118298418241834, expected=0.111182984182418, delta=1.38777878078145E-17

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4DH7_f of PyCapsule object at 0x7f5b01783cf0>
Shamrock df(q) = <built-in method M4DH7_df of PyCapsule object at 0x7f5b01783d50>
Shamrock ddf(q) = <built-in method M4DH7_ddf of PyCapsule object at 0x7f5b01783d80>
Shamrock phi_tilde_3d(q) = <built-in method M4DH7_phi_tilde_3d of PyCapsule object at 0x7f5b01783db0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH7_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01783de0>
Absolute error f(q) = 2.7745488273644668238686811902901e-16
Absolute error df(q) = 1.4083465578598287695160080384023e-14
Absolute error ddf(q) = 9.2957788638054756203862267414385e-15
Absolute error phi_tilde_3d(q) = 9.2324613461043945026233153008238e-14
Absolute error phi_tilde_3d_prime(q) = 4.9411274580246639734784111264843e-13
------------------------------------------

M4Shift2 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4shift2)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4Shift2
------------------------------------------
f(q)   = Piecewise((1, q < 1), (-(3 - 2*q)**3 + (4 - 2*q)**3/4, q < 3/2), ((4 - 2*q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((0, q < 1), (6*(3 - 2*q)**2 - 3*(4 - 2*q)**2/2, q < 3/2), (-3*(4 - 2*q)**2/2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 1), (36*q - 48, q < 3/2), (24 - 12*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(40*q**2 - 231)/60, q < 1), (pi*(48*q**6 - 288*q**5 + 600*q**4 - 440*q**3 - 39*q - 72)/(60*q), q < 3/2), (pi*(-32*q**6 + 288*q**5 - 960*q**4 + 1280*q**3 - 1536*q + 585)/(120*q), q < 2), (-439*pi/(120*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 1), (pi*(288*q**5 - 1440*q**4 + 2400*q**3 - 1320*q**2 - 39)/(60*q) - pi*(48*q**6 - 288*q**5 + 600*q**4 - 440*q**3 - 39*q - 72)/(60*q**2), q < 3/2), (pi*(-192*q**5 + 1440*q**4 - 3840*q**3 + 3840*q**2 - 1536)/(120*q) - pi*(-32*q**6 + 288*q**5 - 960*q**4 + 1280*q**3 - 1536*q + 585)/(120*q**2), q < 2), (439*pi/(120*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4Shift2
------------------------------------------
template<class Tscal>
class KernelDefM4Shift2 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 4.0/11.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (40.0/77.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (120.0/439.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(3 - 2*q);
        Tscal t2 = sham::pow_constexpr<3>(4 - 2*q);
        if (q < 1) {
            return 1;
        }         else if (q < 3.0/2.0) {
            return -t1 + (1.0/4.0)*t2;
        }         else if (q < 2) {
            return (1.0/4.0)*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(3 - 2*q);
        Tscal t2 = sham::pow_constexpr<2>(4 - 2*q);
        if (q < 1) {
            return 0;
        }         else if (q < 3.0/2.0) {
            return 6*t1 - (3.0/2.0)*t2;
        }         else if (q < 2) {
            return -(3.0/2.0)*t2;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        if (q < 1) {
            return 0;
        }         else if (q < 3.0/2.0) {
            return 36*q - 48;
        }         else if (q < 2) {
            return 24 - 12*q;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 1) {
            return (1.0/60.0)*shambase::constants::pi<Tscal>*(40*t1 - 231);
        }         else if (q < 3.0/2.0) {
            return (1.0/60.0)*shambase::constants::pi<Tscal>*(48*t5 - 288*t4 + 600*t3 - 440*t2 - 39*q - 72)/q;
        }         else if (q < 2) {
            return (1.0/120.0)*shambase::constants::pi<Tscal>*(-32*t5 + 288*t4 - 960*t3 + 1280*t2 - 1536*q + 585)/q;
        }
        else
            return -(439.0/120.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 1) {
            return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
        }         else if (q < 3.0/2.0) {
            return (1.0/60.0)*shambase::constants::pi<Tscal>*(288*t4 - 1440*t3 + 2400*t2 - 1320*t1 - 39)/q - (1.0/60.0)*shambase::constants::pi<Tscal>*(48*t5 - 288*t4 + 600*t3 - 440*t2 - 39*q - 72)/t1;
        }         else if (q < 2) {
            return (1.0/120.0)*shambase::constants::pi<Tscal>*(-192*t4 + 1440*t3 - 3840*t2 + 3840*t1 - 1536)/q - (1.0/120.0)*shambase::constants::pi<Tscal>*(-32*t5 + 288*t4 - 960*t3 + 1280*t2 - 1536*q + 585)/t1;
        }
        else
            return (439.0/120.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4shift2 = test_kernel(ret)

------------------------------------------
Testing kernel M4Shift2 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.36363636363636365, expected=0.363636363636364, delta=0
Shamrock norm_2d = 0.1653557850305406, expected=0.165355785030541, delta=0
Shamrock norm_3d = 0.08700953608668538, expected=0.0870095360866854, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift2_f of PyCapsule object at 0x7f5b01788060>
Shamrock df(q) = <built-in method M4Shift2_df of PyCapsule object at 0x7f5b01788090>
Shamrock ddf(q) = <built-in method M4Shift2_ddf of PyCapsule object at 0x7f5b017880c0>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift2_phi_tilde_3d of PyCapsule object at 0x7f5b017880f0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift2_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01788120>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 9.7144514654701197287067770957947e-16
Absolute error ddf(q) = 1.7763568394002505e-15
Absolute error phi_tilde_3d(q) = 4.2657861392141001571651397171902e-14
Absolute error phi_tilde_3d_prime(q) = 8.5828449702952621285761020609829e-14
------------------------------------------

M4Shift4 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4shift4)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4Shift4
------------------------------------------
f(q)   = Piecewise((1, q < 3/2), (-(7 - 4*q)**3 + (8 - 4*q)**3/4, q < 7/4), ((8 - 4*q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((0, q < 3/2), (12*(7 - 4*q)**2 - 3*(8 - 4*q)**2, q < 7/4), (-3*(8 - 4*q)**2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 3/2), (288*q - 480, q < 7/4), (192 - 96*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(160*q**2 - 1371)/240, q < 3/2), (pi*(1536*q**6 - 11520*q**5 + 31680*q**4 - 34400*q**3 + 25845*q - 14580)/(240*q), q < 7/4), (pi*(-2048*q**6 + 18432*q**5 - 61440*q**4 + 81920*q**3 - 98304*q + 59329)/(960*q), q < 2), (-2069*pi/(320*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 3/2), (pi*(9216*q**5 - 57600*q**4 + 126720*q**3 - 103200*q**2 + 25845)/(240*q) - pi*(1536*q**6 - 11520*q**5 + 31680*q**4 - 34400*q**3 + 25845*q - 14580)/(240*q**2), q < 7/4), (pi*(-12288*q**5 + 92160*q**4 - 245760*q**3 + 245760*q**2 - 98304)/(960*q) - pi*(-2048*q**6 + 18432*q**5 - 61440*q**4 + 81920*q**3 - 98304*q + 59329)/(960*q**2), q < 2), (2069*pi/(320*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4Shift4
------------------------------------------
template<class Tscal>
class KernelDefM4Shift4 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 8.0/27.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (160.0/457.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (320.0/2069.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(7 - 4*q);
        Tscal t2 = sham::pow_constexpr<3>(8 - 4*q);
        if (q < 3.0/2.0) {
            return 1;
        }         else if (q < 7.0/4.0) {
            return -t1 + (1.0/4.0)*t2;
        }         else if (q < 2) {
            return (1.0/4.0)*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(7 - 4*q);
        Tscal t2 = sham::pow_constexpr<2>(8 - 4*q);
        if (q < 3.0/2.0) {
            return 0;
        }         else if (q < 7.0/4.0) {
            return 12*t1 - 3*t2;
        }         else if (q < 2) {
            return -3*t2;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        if (q < 3.0/2.0) {
            return 0;
        }         else if (q < 7.0/4.0) {
            return 288*q - 480;
        }         else if (q < 2) {
            return 192 - 96*q;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 3.0/2.0) {
            return (1.0/240.0)*shambase::constants::pi<Tscal>*(160*t1 - 1371);
        }         else if (q < 7.0/4.0) {
            return (1.0/240.0)*shambase::constants::pi<Tscal>*(1536*t5 - 11520*t4 + 31680*t3 - 34400*t2 + 25845*q - 14580)/q;
        }         else if (q < 2) {
            return (1.0/960.0)*shambase::constants::pi<Tscal>*(-2048*t5 + 18432*t4 - 61440*t3 + 81920*t2 - 98304*q + 59329)/q;
        }
        else
            return -(2069.0/320.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 3.0/2.0) {
            return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
        }         else if (q < 7.0/4.0) {
            return (1.0/240.0)*shambase::constants::pi<Tscal>*(9216*t4 - 57600*t3 + 126720*t2 - 103200*t1 + 25845)/q - (1.0/240.0)*shambase::constants::pi<Tscal>*(1536*t5 - 11520*t4 + 31680*t3 - 34400*t2 + 25845*q - 14580)/t1;
        }         else if (q < 2) {
            return (1.0/960.0)*shambase::constants::pi<Tscal>*(-12288*t4 + 92160*t3 - 245760*t2 + 245760*t1 - 98304)/q - (1.0/960.0)*shambase::constants::pi<Tscal>*(-2048*t5 + 18432*t4 - 61440*t3 + 81920*t2 - 98304*q + 59329)/t1;
        }
        else
            return (2069.0/320.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4shift4 = test_kernel(ret)

------------------------------------------
Testing kernel M4Shift4 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.2962962962962963, expected=0.296296296296296, delta=0
Shamrock norm_2d = 0.11144328619126151, expected=0.111443286191262, delta=0
Shamrock norm_3d = 0.049231108544617215, expected=0.0492311085446172, delta=0

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift4_f of PyCapsule object at 0x7f5b01788360>
Shamrock df(q) = <built-in method M4Shift4_df of PyCapsule object at 0x7f5b01788390>
Shamrock ddf(q) = <built-in method M4Shift4_ddf of PyCapsule object at 0x7f5b017883c0>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift4_phi_tilde_3d of PyCapsule object at 0x7f5b017883f0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift4_phi_tilde_3d_prime of PyCapsule object at 0x7f5b01788420>
Absolute error f(q) = 2.7755575615628913510590791702271e-16
Absolute error df(q) = 1.8318679906315083e-15
Absolute error ddf(q) = 2.1316282072803006e-14
Absolute error phi_tilde_3d(q) = 6.2097746275599988316063341215081e-13
Absolute error phi_tilde_3d_prime(q) = 1.0539828786480343310056803196174e-12
------------------------------------------

M4Shift8 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4shift8)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4Shift8
------------------------------------------
f(q)   = Piecewise((1, q < 7/4), (-(15 - 8*q)**3 + (16 - 8*q)**3/4, q < 15/8), ((16 - 8*q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((0, q < 7/4), (24*(15 - 8*q)**2 - 6*(16 - 8*q)**2, q < 15/8), (-6*(16 - 8*q)**2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 7/4), (2304*q - 4224, q < 15/8), (1536 - 768*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(640*q**2 - 6531)/960, q < 7/4), (pi*(49152*q**6 - 405504*q**5 + 1236480*q**4 - 1504640*q**3 + 1491693*q - 907578)/(960*q), q < 15/8), (pi*(-131072*q**6 + 1179648*q**5 - 3932160*q**4 + 5242880*q**3 - 6291456*q + 4130001)/(7680*q), q < 2), (-64303*pi/(7680*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 7/4), (pi*(294912*q**5 - 2027520*q**4 + 4945920*q**3 - 4513920*q**2 + 1491693)/(960*q) - pi*(49152*q**6 - 405504*q**5 + 1236480*q**4 - 1504640*q**3 + 1491693*q - 907578)/(960*q**2), q < 15/8), (pi*(-786432*q**5 + 5898240*q**4 - 15728640*q**3 + 15728640*q**2 - 6291456)/(7680*q) - pi*(-131072*q**6 + 1179648*q**5 - 3932160*q**4 + 5242880*q**3 - 6291456*q + 4130001)/(7680*q**2), q < 2), (64303*pi/(7680*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4Shift8
------------------------------------------
template<class Tscal>
class KernelDefM4Shift8 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 16.0/59.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (640.0/2177.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (7680.0/64303.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(15 - 8*q);
        Tscal t2 = sham::pow_constexpr<3>(16 - 8*q);
        if (q < 7.0/4.0) {
            return 1;
        }         else if (q < 15.0/8.0) {
            return -t1 + (1.0/4.0)*t2;
        }         else if (q < 2) {
            return (1.0/4.0)*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(15 - 8*q);
        Tscal t2 = sham::pow_constexpr<2>(16 - 8*q);
        if (q < 7.0/4.0) {
            return 0;
        }         else if (q < 15.0/8.0) {
            return 24*t1 - 6*t2;
        }         else if (q < 2) {
            return -6*t2;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        if (q < 7.0/4.0) {
            return 0;
        }         else if (q < 15.0/8.0) {
            return 2304*q - 4224;
        }         else if (q < 2) {
            return 1536 - 768*q;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 7.0/4.0) {
            return (1.0/960.0)*shambase::constants::pi<Tscal>*(640*t1 - 6531);
        }         else if (q < 15.0/8.0) {
            return (1.0/960.0)*shambase::constants::pi<Tscal>*(49152*t5 - 405504*t4 + 1236480*t3 - 1504640*t2 + 1491693*q - 907578)/q;
        }         else if (q < 2) {
            return (1.0/7680.0)*shambase::constants::pi<Tscal>*(-131072*t5 + 1179648*t4 - 3932160*t3 + 5242880*t2 - 6291456*q + 4130001)/q;
        }
        else
            return -(64303.0/7680.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 7.0/4.0) {
            return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
        }         else if (q < 15.0/8.0) {
            return (1.0/960.0)*shambase::constants::pi<Tscal>*(294912*t4 - 2027520*t3 + 4945920*t2 - 4513920*t1 + 1491693)/q - (1.0/960.0)*shambase::constants::pi<Tscal>*(49152*t5 - 405504*t4 + 1236480*t3 - 1504640*t2 + 1491693*q - 907578)/t1;
        }         else if (q < 2) {
            return (1.0/7680.0)*shambase::constants::pi<Tscal>*(-786432*t4 + 5898240*t3 - 15728640*t2 + 15728640*t1 - 6291456)/q - (1.0/7680.0)*shambase::constants::pi<Tscal>*(-131072*t5 + 1179648*t4 - 3932160*t3 + 5242880*t2 - 6291456*q + 4130001)/t1;
        }
        else
            return (64303.0/7680.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4shift8 = test_kernel(ret)

------------------------------------------
Testing kernel M4Shift8 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.2711864406779661, expected=0.271186440677966, delta=0
Shamrock norm_2d = 0.09357755037098117, expected=0.0935775503709812, delta=1.38777878078145E-17
Shamrock norm_3d = 0.03801719866711526, expected=0.0380171986671153, delta=6.93889390390723E-18

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift8_f of PyCapsule object at 0x7f5b01788630>
Shamrock df(q) = <built-in method M4Shift8_df of PyCapsule object at 0x7f5b01788660>
Shamrock ddf(q) = <built-in method M4Shift8_ddf of PyCapsule object at 0x7f5b01788690>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift8_phi_tilde_3d of PyCapsule object at 0x7f5b017886c0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift8_phi_tilde_3d_prime of PyCapsule object at 0x7f5b017886f0>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 5.773159728050814e-15
Absolute error ddf(q) = 2.8421709430404007e-13
Absolute error phi_tilde_3d(q) = 5.6625300231436590457808797992127e-12
Absolute error phi_tilde_3d_prime(q) = 8.557371206771937550851606776476e-12
------------------------------------------

M4Shift16 Kernel

Generate c++ code for the kernel

ret = kernel_to_shamrock(m4shift16)

print_kernel_info(ret)

------------------------------------------
Sympy expression for kernel M4Shift16
------------------------------------------
f(q)   = Piecewise((1, q < 15/8), (-(31 - 16*q)**3 + (32 - 16*q)**3/4, q < 31/16), ((32 - 16*q)**3/4, q < 2), (0, True))
df(q)  = Piecewise((0, q < 15/8), (48*(31 - 16*q)**2 - 12*(32 - 16*q)**2, q < 31/16), (-12*(32 - 16*q)**2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 15/8), (18432*q - 35328, q < 31/16), (12288 - 6144*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(2560*q**2 - 28371)/3840, q < 15/8), (pi*(1572864*q**6 - 13565952*q**5 + 43315200*q**4 - 55293440*q**3 + 60721629*q - 38728125)/(3840*q), q < 31/16), (pi*(-8388608*q**6 + 75497472*q**5 - 251658240*q**4 + 335544320*q**3 - 402653184*q + 267853681)/(61440*q), q < 2), (-38785*pi/(4096*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 15/8), (pi*(9437184*q**5 - 67829760*q**4 + 173260800*q**3 - 165880320*q**2 + 60721629)/(3840*q) - pi*(1572864*q**6 - 13565952*q**5 + 43315200*q**4 - 55293440*q**3 + 60721629*q - 38728125)/(3840*q**2), q < 31/16), (pi*(-50331648*q**5 + 377487360*q**4 - 1006632960*q**3 + 1006632960*q**2 - 402653184)/(61440*q) - pi*(-8388608*q**6 + 75497472*q**5 - 251658240*q**4 + 335544320*q**3 - 402653184*q + 267853681)/(61440*q**2), q < 2), (38785*pi/(4096*q**2), True))
------------------------------------------

print_kernel_cpp_code(ret)

------------------------------------------
C++ generated code for kernel M4Shift16
------------------------------------------
template<class Tscal>
class KernelDefM4Shift16 {
    public:
    inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
    /// default hfact to be used for this kernel
    inline static constexpr Tscal hfactd = 1.2;

    /// 1D norm of the kernel
    inline static constexpr Tscal norm_1d = 32.0/123.0;
    /// 2D norm of the kernel
    inline static constexpr Tscal norm_2d = (2560.0/9457.0)/shambase::constants::pi<Tscal>;
    /// 3D norm of the kernel
    inline static constexpr Tscal norm_3d = (4096.0/38785.0)/shambase::constants::pi<Tscal>;

    inline static Tscal f(Tscal q) {
        Tscal t1 = sham::pow_constexpr<3>(31 - 16*q);
        Tscal t2 = sham::pow_constexpr<3>(32 - 16*q);
        if (q < 15.0/8.0) {
            return 1;
        }         else if (q < 31.0/16.0) {
            return -t1 + (1.0/4.0)*t2;
        }         else if (q < 2) {
            return (1.0/4.0)*t2;
        }
        else
            return 0;
    }


    inline static Tscal df(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(31 - 16*q);
        Tscal t2 = sham::pow_constexpr<2>(32 - 16*q);
        if (q < 15.0/8.0) {
            return 0;
        }         else if (q < 31.0/16.0) {
            return 48*t1 - 12*t2;
        }         else if (q < 2) {
            return -12*t2;
        }
        else
            return 0;
    }

    inline static Tscal ddf(Tscal q) {
        if (q < 15.0/8.0) {
            return 0;
        }         else if (q < 31.0/16.0) {
            return 18432*q - 35328;
        }         else if (q < 2) {
            return 12288 - 6144*q;
        }
        else
            return 0;
    }

    inline static Tscal phi_tilde_3d(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 15.0/8.0) {
            return (1.0/3840.0)*shambase::constants::pi<Tscal>*(2560*t1 - 28371);
        }         else if (q < 31.0/16.0) {
            return (1.0/3840.0)*shambase::constants::pi<Tscal>*(1572864*t5 - 13565952*t4 + 43315200*t3 - 55293440*t2 + 60721629*q - 38728125)/q;
        }         else if (q < 2) {
            return (1.0/61440.0)*shambase::constants::pi<Tscal>*(-8388608*t5 + 75497472*t4 - 251658240*t3 + 335544320*t2 - 402653184*q + 267853681)/q;
        }
        else
            return -(38785.0/4096.0)*shambase::constants::pi<Tscal>/q;
    }

    inline static Tscal phi_tilde_3d_prime(Tscal q) {
        Tscal t1 = sham::pow_constexpr<2>(q);
        Tscal t2 = sham::pow_constexpr<3>(q);
        Tscal t3 = sham::pow_constexpr<4>(q);
        Tscal t4 = sham::pow_constexpr<5>(q);
        Tscal t5 = sham::pow_constexpr<6>(q);
        if (q < 15.0/8.0) {
            return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
        }         else if (q < 31.0/16.0) {
            return (1.0/3840.0)*shambase::constants::pi<Tscal>*(9437184*t4 - 67829760*t3 + 173260800*t2 - 165880320*t1 + 60721629)/q - (1.0/3840.0)*shambase::constants::pi<Tscal>*(1572864*t5 - 13565952*t4 + 43315200*t3 - 55293440*t2 + 60721629*q - 38728125)/t1;
        }         else if (q < 2) {
            return (1.0/61440.0)*shambase::constants::pi<Tscal>*(-50331648*t4 + 377487360*t3 - 1006632960*t2 + 1006632960*t1 - 402653184)/q - (1.0/61440.0)*shambase::constants::pi<Tscal>*(-8388608*t5 + 75497472*t4 - 251658240*t3 + 335544320*t2 - 402653184*q + 267853681)/t1;
        }
        else
            return (38785.0/4096.0)*shambase::constants::pi<Tscal>/t1;
    }

};


------------------------------------------

Test the kernel

test_result_m4shift16 = test_kernel(ret)

------------------------------------------
Testing kernel M4Shift16 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.2601626016260163, expected=0.260162601626016, delta=5.55111512312578E-17
Shamrock norm_2d = 0.08616615296928246, expected=0.0861661529692825, delta=0
Shamrock norm_3d = 0.03361601891991251, expected=0.0336160189199125, delta=6.93889390390723E-18

Testing kernel radius:
Shamrock Rkern = 2.0, delta=0

Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift16_f of PyCapsule object at 0x7f5b01788930>
Shamrock df(q) = <built-in method M4Shift16_df of PyCapsule object at 0x7f5b01788960>
Shamrock ddf(q) = <built-in method M4Shift16_ddf of PyCapsule object at 0x7f5b01788990>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift16_phi_tilde_3d of PyCapsule object at 0x7f5b017889c0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift16_phi_tilde_3d_prime of PyCapsule object at 0x7f5b017889f0>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 5.329070518200751e-15
Absolute error ddf(q) = 2.2737367544323206e-12
Absolute error phi_tilde_3d(q) = 4.1121221586034293016703365082286e-11
Absolute error phi_tilde_3d_prime(q) = 6.4350934828587342081606789080707e-11
------------------------------------------

Plot the kernels

fig_sz = (6.4, 12)

Plotting helper functions

def create_kernel_plot_figure(fig_sz):
    """Create a figure with 5 subplots for kernel plotting"""
    plt.figure(figsize=fig_sz)
    ax_f = plt.subplot(5, 1, 1)
    ax_df = plt.subplot(5, 1, 2)
    ax_ddf = plt.subplot(5, 1, 3)
    ax_phi_tilde_3d = plt.subplot(5, 1, 4)
    ax_phi_tilde_3d_prime = plt.subplot(5, 1, 5)
    return ax_f, ax_df, ax_ddf, ax_phi_tilde_3d, ax_phi_tilde_3d_prime


def plot_kernel_result(axes, test_result, kernel_label):
    """Plot a single kernel result on the given axes"""
    ax_f, ax_df, ax_ddf, ax_phi_tilde_3d, ax_phi_tilde_3d_prime = axes
    q_arr = test_result["q_arr"]

    ax_f.plot(q_arr, test_result["shamrock_Cf"], label=f"C_3d f_{kernel_label}(q)")
    ax_df.plot(q_arr, test_result["shamrock_Cdf"], label=f"C_3d df_{kernel_label}(q)")
    ax_ddf.plot(q_arr, test_result["shamrock_Cddf"], label=f"C_3d ddf_{kernel_label}(q)")
    ax_phi_tilde_3d.plot(
        q_arr, test_result["shamrock_Cphi_tilde_3d"], label=f"C_3d phi_tilde_3d_{kernel_label}(q)"
    )
    ax_phi_tilde_3d_prime.plot(
        q_arr,
        test_result["shamrock_Cphi_tilde_3d_prime"],
        label=f"C_3d phi_tilde_3d_prime_{kernel_label}(q)",
    )


def finalize_kernel_plot(axes):
    """Add titles, labels, and legends to the kernel plot"""
    ax_f, ax_df, ax_ddf, ax_phi_tilde_3d, ax_phi_tilde_3d_prime = axes

    # Get current axis limits before adding reference line
    xlim = ax_phi_tilde_3d.get_xlim()
    ylim = ax_phi_tilde_3d.get_ylim()

    ylim = (1.2 * ylim[0], ylim[1])

    # Add -1/r reference line (only within current x range)
    q = np.linspace(max(1e-6, xlim[0]), xlim[1], 1000)
    one_over_r = -1 / q
    ax_phi_tilde_3d.plot(q, one_over_r, "--", color="grey", label="-1/r")

    # Restore original limits (ignore reference line for autoscaling)
    ax_phi_tilde_3d.set_xlim(xlim)
    ax_phi_tilde_3d.set_ylim(ylim)

    # Get current axis limits before adding reference line
    xlim = ax_phi_tilde_3d_prime.get_xlim()
    ylim = ax_phi_tilde_3d_prime.get_ylim()

    ylim = (ylim[0], 1.5 * ylim[1])

    # Add 1/r^2 reference line (only within current x range)
    one_over_r_squared = 1 / q**2
    ax_phi_tilde_3d_prime.plot(q, one_over_r_squared, "--", color="grey", label="1/r^2")

    # Restore original limits (ignore reference line for autoscaling)
    ax_phi_tilde_3d_prime.set_xlim(xlim)
    ax_phi_tilde_3d_prime.set_ylim(ylim)

    ax_f.set_title("C_3d f(q)", fontsize=10)
    ax_df.set_title("C_3d df(q)", fontsize=10)
    ax_ddf.set_title("C_3d ddf(q)", fontsize=10)
    ax_phi_tilde_3d.set_title("C_3d phi_tilde_3d(q)", fontsize=10)
    ax_phi_tilde_3d_prime.set_title("C_3d phi_tilde_3d_prime(q)", fontsize=10)

    ax_f.set_xlabel("q")
    ax_df.set_xlabel("q")
    ax_ddf.set_xlabel("q")
    ax_phi_tilde_3d.set_xlabel("q")
    ax_phi_tilde_3d_prime.set_xlabel("q")

    ax_f.legend(loc="right", fontsize=8)
    ax_df.legend(loc="right", fontsize=8)
    ax_ddf.legend(loc="right", fontsize=8)
    ax_phi_tilde_3d.legend(loc="right", fontsize=8)
    ax_phi_tilde_3d_prime.legend(loc="right", fontsize=8)

    plt.tight_layout()
    plt.show()

Cubic splines kernels (M-series)

axes = create_kernel_plot_figure(fig_sz)
plot_kernel_result(axes, test_result_m4, "m4")
plot_kernel_result(axes, test_result_m5, "m5")
plot_kernel_result(axes, test_result_m6, "m6")
plot_kernel_result(axes, test_result_m7, "m7")
plot_kernel_result(axes, test_result_m8, "m8")
plot_kernel_result(axes, test_result_m9, "m9")
plot_kernel_result(axes, test_result_m10, "m10")
finalize_kernel_plot(axes)
plt.show()

Wendland kernels (C-series)

axes = create_kernel_plot_figure(fig_sz)
plot_kernel_result(axes, test_result_c2, "c2")
plot_kernel_result(axes, test_result_c4, "c4")
plot_kernel_result(axes, test_result_c6, "c6")
finalize_kernel_plot(axes)
plt.show()

Double hump kernels

axes = create_kernel_plot_figure(fig_sz)
plot_kernel_result(axes, test_result_m4dh, "m4dh")
plot_kernel_result(axes, test_result_m4dh3, "m4dh3")
plot_kernel_result(axes, test_result_m4dh5, "m4dh5")
plot_kernel_result(axes, test_result_m4dh7, "m4dh7")
finalize_kernel_plot(axes)
plt.show()

M4Shift kernels

axes = create_kernel_plot_figure(fig_sz)
plot_kernel_result(axes, test_result_m4shift2, "m4shift2")
plot_kernel_result(axes, test_result_m4shift4, "m4shift4")
plot_kernel_result(axes, test_result_m4shift8, "m4shift8")
plot_kernel_result(axes, test_result_m4shift16, "m4shift16")
finalize_kernel_plot(axes)
plt.show()

Estimated memory usage: 210 MB