function I = get_integral(cache_filename, var_name, tau_range, M)
    syms x f;

    if exist(cache_filename, "file")
        data = load(cache_filename, var_name);
        I = data.(var_name);
    else
        I = zeros(length(tau_range), 0);

        for tau_idx = 1:length(tau_range)
            tau = tau_range(tau_idx);
            f = (x - tau) ^ 2 * exp(-x ^ 2/2) * x ^ (M - 1) / (2 ^ (M / 2 - 1) * gamma(M / 2));
            I(tau_idx) = double(int(f, [tau, +inf]));
        end

        S.(var_name) = I;
        save(cache_filename, '-struct', 'S', var_name);
    end

end