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


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

torch::Tensor fused_bias_scale_leaky_cuda(
    torch::Tensor input,
    torch::Tensor bias,
    double multiplier,
    double negative_slope);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAStream.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <algorithm>

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

__global__ void fused_vector_kernel(
    const float4* __restrict__ input,
    const float4* __restrict__ bias,
    float4* __restrict__ output,
    int64_t vectors,
    int64_t channel_vectors,
    float multiplier,
    float negative_slope) {

    int64_t index = static_cast<int64_t>(blockIdx.x) * blockDim.x
                    + threadIdx.x;
    int64_t stride = static_cast<int64_t>(blockDim.x) * gridDim.x;

    for (; index < vectors; index += stride) {
        float4 value = input[index];
        float4 offset = bias[index % channel_vectors];

        float4 result;
        result.x = activate(value.x, offset.x, multiplier, negative_slope);
        result.y = activate(value.y, offset.y, multiplier, negative_slope);
        result.z = activate(value.z, offset.z, multiplier, negative_slope);
        result.w = activate(value.w, offset.w, multiplier, negative_slope);
        output[index] = result;
    }
}

__global__ void fused_scalar_kernel(
    const float* __restrict__ input,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int64_t elements,
    int64_t channels,
    float multiplier,
    float negative_slope) {

    int64_t index = static_cast<int64_t>(blockIdx.x) * blockDim.x
                    + threadIdx.x;
    int64_t stride = static_cast<int64_t>(blockDim.x) * gridDim.x;

    for (; index < elements; index += stride) {
        output[index] = activate(
            input[index], bias[index % channels],
            multiplier, negative_slope);
    }
}

torch::Tensor fused_bias_scale_leaky_cuda(
    torch::Tensor input,
    torch::Tensor bias,
    double multiplier,
    double negative_slope) {

    TORCH_CHECK(input.is_cuda() && bias.is_cuda(),
                "Input and bias must be CUDA tensors");
    TORCH_CHECK(input.device() == bias.device(),
                "Input and bias must be on the same device");
    TORCH_CHECK(input.scalar_type() == torch::kFloat32 &&
                bias.scalar_type() == torch::kFloat32,
                "Input and bias must have float32 dtype");
    TORCH_CHECK(input.is_contiguous() && bias.is_contiguous(),
                "Input and bias must be contiguous");
    TORCH_CHECK(input.dim() >= 1 && bias.dim() == 1 &&
                input.size(-1) == bias.numel(),
                "Bias must match the last input dimension");

    const c10::cuda::CUDAGuard device_guard(input.device());
    auto output = torch::empty_like(input);

    const int64_t elements = input.numel();
    if (elements == 0) {
        return output;
    }

    const int64_t channels = bias.numel();
    constexpr int threads = 256;
    const auto stream =
        c10::cuda::getCurrentCUDAStream(input.get_device());

    const bool vectorized =
        channels % 4 == 0 &&
        reinterpret_cast<uintptr_t>(input.data_ptr<float>()) % 16 == 0 &&
        reinterpret_cast<uintptr_t>(bias.data_ptr<float>()) % 16 == 0;

    if (vectorized) {
        const int64_t vectors = elements / 4;
        const int blocks = static_cast<int>(
            std::min<int64_t>((vectors + threads - 1) / threads, 4096));

        fused_vector_kernel<<<blocks, threads, 0, stream.stream()>>>(
            reinterpret_cast<const float4*>(input.data_ptr<float>()),
            reinterpret_cast<const float4*>(bias.data_ptr<float>()),
            reinterpret_cast<float4*>(output.data_ptr<float>()),
            vectors, channels / 4,
            static_cast<float>(multiplier),
            static_cast<float>(negative_slope));
    } else {
        const int blocks = static_cast<int>(
            std::min<int64_t>((elements + threads - 1) / threads, 4096));

        fused_scalar_kernel<<<blocks, threads, 0, stream.stream()>>>(
            input.data_ptr<float>(),
            bias.data_ptr<float>(),
            output.data_ptr<float>(),
            elements, channels,
            static_cast<float>(multiplier),
            static_cast<float>(negative_slope));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="gemm_bias_scale_leaky_fused",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_bias_scale_leaky_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 x.dtype != torch.float32
            or self.gemm.weight.dtype != torch.float32
            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 self.leaky_relu(self.gemm(x) * self.multiplier)

        # Keep the GEMM in cuBLAS and fuse bias, scaling, and activation.
        product = F.linear(x, self.gemm.weight, bias=None)
        return _fused_ops.fused_bias_scale_leaky_cuda(
            product.contiguous(),
            self.gemm.bias,
            float(self.multiplier),
            float(self.leaky_relu.negative_slope),
        )
```