```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_min_sub_cuda(
    torch::Tensor output,
    torch::Tensor bias,
    torch::Tensor constant);
"""

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

__device__ __forceinline__ float finish(
    float value, float bias, float constant)
{
    value = value + bias;
    float minimum = isnan(value) ? value :
                    (isnan(constant) ? constant : fminf(value, constant));
    return minimum - constant;
}

__global__ void fused_bias_min_sub_vectorized(
    float* __restrict__ output,
    const float* __restrict__ bias,
    const float* __restrict__ constant,
    int64_t vectors,
    int64_t width_vectors)
{
    const float c = constant[0];
    for (int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
                     + threadIdx.x;
         i < vectors;
         i += static_cast<int64_t>(blockDim.x) * gridDim.x) {
        float4 value = reinterpret_cast<const float4*>(output)[i];
        float4 b = reinterpret_cast<const float4*>(bias)[i % width_vectors];

        value.x = finish(value.x, b.x, c);
        value.y = finish(value.y, b.y, c);
        value.z = finish(value.z, b.z, c);
        value.w = finish(value.w, b.w, c);

        reinterpret_cast<float4*>(output)[i] = value;
    }
}

__global__ void fused_bias_min_sub_scalar(
    float* __restrict__ output,
    const float* __restrict__ bias,
    const float* __restrict__ constant,
    int64_t count,
    int64_t width)
{
    const float c = constant[0];
    for (int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
                     + threadIdx.x;
         i < count;
         i += static_cast<int64_t>(blockDim.x) * gridDim.x) {
        output[i] = finish(output[i], bias[i % width], c);
    }
}

torch::Tensor fused_bias_min_sub_cuda(
    torch::Tensor output,
    torch::Tensor bias,
    torch::Tensor constant)
{
    TORCH_CHECK(output.is_cuda() && bias.is_cuda() && constant.is_cuda(),
                "All tensors must be CUDA tensors");
    TORCH_CHECK(output.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat &&
                constant.scalar_type() == at::kFloat,
                "All tensors must be float32");
    TORCH_CHECK(output.device() == bias.device() &&
                output.device() == constant.device(),
                "All tensors must share a device");
    TORCH_CHECK(output.is_contiguous() && bias.is_contiguous() &&
                constant.is_contiguous(), "Tensors must be contiguous");
    TORCH_CHECK(output.dim() > 0 && bias.dim() == 1 &&
                output.size(-1) == bias.numel() && constant.numel() == 1,
                "Invalid tensor shapes");

    const 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();

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

    const int64_t work = vectorized ? count / 4 : count;
    const int64_t needed_blocks = (work + threads - 1) / threads;
    const int blocks = static_cast<int>(
        needed_blocks < 4096 ? needed_blocks : 4096);

    if (vectorized) {
        fused_bias_min_sub_vectorized<<<blocks, threads, 0, stream>>>(
            output.data_ptr<float>(),
            bias.data_ptr<float>(),
            constant.data_ptr<float>(),
            count / 4,
            width / 4);
    } else {
        fused_bias_min_sub_scalar<<<blocks, threads, 0, stream>>>(
            output.data_ptr<float>(),
            bias.data_ptr<float>(),
            constant.data_ptr<float>(),
            count,
            width);
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="linear_bias_min_sub_fused_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_bias_min_sub_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self, in_features, out_features, constant):
        super().__init__()
        self.linear = nn.Linear(in_features, out_features)
        self.constant = nn.Parameter(torch.tensor(constant))

    def forward(self, x):
        # Preserve autograd and mixed-precision behavior when needed.
        needs_grad = torch.is_grad_enabled() and (
            x.requires_grad
            or self.linear.weight.requires_grad
            or self.linear.bias.requires_grad
            or self.constant.requires_grad
        )
        if (
            needs_grad
            or not x.is_cuda
            or x.dtype != torch.float32
            or self.linear.weight.dtype != torch.float32
            or self.linear.bias.dtype != torch.float32
            or self.constant.dtype != torch.float32
            or torch.is_autocast_enabled()
        ):
            y = self.linear(x)
            return torch.minimum(y, self.constant) - self.constant

        # Keep the optimized GEMM and fuse its remaining epilogue in place.
        y = F.linear(x, self.linear.weight, None)
        return _fused_ops.fused_bias_min_sub_cuda(
            y, self.linear.bias, self.constant
        )
```