```python
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.cpp_extension import load_inline


_cpp_source = """
#include <torch/extension.h>

torch::Tensor fused_epilogue_cuda(
    torch::Tensor output,
    torch::Tensor bias,
    double multiplier,
    double negative_slope);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/AccumulateType.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>

__device__ __forceinline__ float activate(
    float value, float bias, float multiplier, float slope) {
    float z = (value + bias) * multiplier;
    return z > 0.0f ? z : z * slope;
}

__global__ void epilogue_float4(
    float* __restrict__ output,
    const float* __restrict__ bias,
    int64_t vector_count,
    int64_t width4,
    float multiplier,
    float slope) {
    auto* output4 = reinterpret_cast<float4*>(output);
    const auto* bias4 = reinterpret_cast<const float4*>(bias);

    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < vector_count;
         i += int64_t(blockDim.x) * gridDim.x) {
        float4 v = output4[i];
        float4 b = bias4[i % width4];

        v.x = activate(v.x, b.x, multiplier, slope);
        v.y = activate(v.y, b.y, multiplier, slope);
        v.z = activate(v.z, b.z, multiplier, slope);
        v.w = activate(v.w, b.w, multiplier, slope);

        output4[i] = v;
    }
}

template <typename scalar_t>
__global__ void epilogue_scalar(
    scalar_t* __restrict__ output,
    const scalar_t* __restrict__ bias,
    int64_t count,
    int64_t width,
    double multiplier,
    double slope) {
    using acc_t = at::acc_type<scalar_t, true>;

    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < count;
         i += int64_t(blockDim.x) * gridDim.x) {
        // Preserve the intermediate rounding of separate tensor operators.
        scalar_t sum = static_cast<scalar_t>(
            static_cast<acc_t>(output[i]) +
            static_cast<acc_t>(bias[i % width]));
        scalar_t scaled = static_cast<scalar_t>(
            static_cast<acc_t>(sum) * static_cast<acc_t>(multiplier));
        acc_t z = static_cast<acc_t>(scaled);
        output[i] = static_cast<scalar_t>(
            z > acc_t(0) ? z : z * static_cast<acc_t>(slope));
    }
}

torch::Tensor fused_epilogue_cuda(
    torch::Tensor output,
    torch::Tensor bias,
    double multiplier,
    double negative_slope) {
    TORCH_CHECK(output.is_cuda() && bias.is_cuda(), "Expected CUDA tensors");
    TORCH_CHECK(output.device() == bias.device(), "Device mismatch");
    TORCH_CHECK(output.scalar_type() == bias.scalar_type(), "Dtype mismatch");
    TORCH_CHECK(output.is_contiguous() && bias.is_contiguous(),
                "Expected contiguous tensors");
    TORCH_CHECK(output.dim() == 2 && bias.dim() == 1 &&
                output.size(1) == bias.numel(), "Invalid tensor shapes");

    c10::cuda::CUDAGuard guard(output.device());
    const int64_t count = output.numel();
    if (count == 0) {
        return output;
    }

    const int64_t width = bias.numel();
    constexpr int threads = 256;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (output.scalar_type() == at::kFloat && width % 4 == 0) {
        const int64_t vector_count = count / 4;
        const int blocks = static_cast<int>(
            std::min<int64_t>((vector_count + threads - 1) / threads, 4096));

        epilogue_float4<<<blocks, threads, 0, stream>>>(
            output.data_ptr<float>(),
            bias.data_ptr<float>(),
            vector_count,
            width / 4,
            static_cast<float>(multiplier),
            static_cast<float>(negative_slope));
    } else {
        const int blocks = static_cast<int>(
            std::min<int64_t>((count + threads - 1) / threads, 4096));

        AT_DISPATCH_FLOATING_TYPES_AND2(
            at::ScalarType::Half,
            at::ScalarType::BFloat16,
            output.scalar_type(),
            "fused_linear_epilogue",
            [&] {
                epilogue_scalar<scalar_t><<<blocks, threads, 0, stream>>>(
                    output.data_ptr<scalar_t>(),
                    bias.data_ptr<scalar_t>(),
                    count,
                    width,
                    multiplier,
                    negative_slope);
            });
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="linear_bias_mul_leaky_epilogue",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_epilogue_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3", "--fmad=false"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self, in_features, out_features, multiplier, negative_slope):
        super().__init__()
        self.gemm = nn.Linear(in_features, out_features)
        self.multiplier = multiplier
        self.leaky_relu = nn.LeakyReLU(negative_slope)

    def forward(self, x):
        if (
            not x.is_cuda
            or torch.is_autocast_enabled()
            or (
                torch.is_grad_enabled()
                and (
                    x.requires_grad
                    or self.gemm.weight.requires_grad
                    or self.gemm.bias.requires_grad
                )
            )
        ):
            return F.leaky_relu(
                self.gemm(x) * self.multiplier,
                negative_slope=self.leaky_relu.negative_slope,
            )

        original_shape = x.shape[:-1]
        flat_x = x.reshape(-1, self.gemm.in_features)
        output = torch.mm(flat_x, self.gemm.weight.t())
        output = _fused_ops.fused_epilogue_cuda(
            output,
            self.gemm.bias,
            self.multiplier,
            self.leaky_relu.negative_slope,
        )
        return output.view(*original_shape, self.gemm.out_features)
```