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


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

torch::Tensor conv3d_specialized(
    torch::Tensor x,
    torch::Tensor weight,
    c10::optional<torch::Tensor> bias);
"""

_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>

namespace {

constexpr int TILE_N = 128;
constexpr int TILE_K = 16;

// Implicit GEMM: 64 output channels by 128 spatial positions.
// Each thread computes eight channels and four spatial positions.
__global__ __launch_bounds__(256)
void conv3d_3x3x3_c3_o64_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int depth,
    int height,
    int width,
    int out_height,
    int out_width,
    int64_t output_volume)
{
    __shared__ float shared_weight[TILE_K][65];
    __shared__ float shared_input[TILE_K][TILE_N];

    const int lane = threadIdx.x;
    const int channel_group = threadIdx.y;
    const int tid = channel_group * 32 + lane;
    const int batch = blockIdx.y;
    const int64_t spatial_base = int64_t(blockIdx.x) * TILE_N;

    float accum[8][4];

    #pragma unroll
    for (int m = 0; m < 8; ++m) {
        #pragma unroll
        for (int n = 0; n < 4; ++n) {
            accum[m][n] = 0.0f;
        }
    }

    // Spatial coordinates are invariant across reduction tiles.
    int iz[8];
    int iy[8];
    int ix[8];
    bool valid_position[8];

    #pragma unroll
    for (int t = 0; t < 8; ++t) {
        const int q = tid + t * 256;
        const int spatial_offset = q % TILE_N;
        const int64_t p = spatial_base + spatial_offset;

        valid_position[t] = p < output_volume;
        iz[t] = int(p / (int64_t(out_height) * out_width));
        iy[t] = int((p / out_width) % out_height);
        ix[t] = int(p % out_width);
    }

    #pragma unroll
    for (int reduction_base = 0; reduction_base < 81;
         reduction_base += TILE_K) {
        #pragma unroll
        for (int t = 0; t < 4; ++t) {
            const int q = tid + t * 256;
            const int oc = q / TILE_K;
            const int local_k = q % TILE_K;
            const int k = reduction_base + local_k;

            shared_weight[local_k][oc] =
                k < 81 ? weight[oc * 81 + k] : 0.0f;
        }

        #pragma unroll
        for (int t = 0; t < 8; ++t) {
            const int q = tid + t * 256;
            const int local_k = q / TILE_N;
            const int spatial_offset = q % TILE_N;
            const int k = reduction_base + local_k;

            float value = 0.0f;
            if (k < 81 && valid_position[t]) {
                const int ic = k / 27;
                const int kz = (k / 9) % 3;
                const int ky = (k / 3) % 3;
                const int kx = k % 3;

                const int64_t input_index =
                    (((int64_t(batch) * 3 + ic) * depth
                      + iz[t] + kz) * height
                     + iy[t] + ky) * width
                    + ix[t] + kx;

                value = x[input_index];
            }
            shared_input[local_k][spatial_offset] = value;
        }

        __syncthreads();

        #pragma unroll
        for (int k = 0; k < TILE_K; ++k) {
            float inputs[4];

            #pragma unroll
            for (int n = 0; n < 4; ++n) {
                inputs[n] = shared_input[k][lane + n * 32];
            }

            #pragma unroll
            for (int m = 0; m < 8; ++m) {
                const float w =
                    shared_weight[k][channel_group + m * 8];

                #pragma unroll
                for (int n = 0; n < 4; ++n) {
                    accum[m][n] = fmaf(w, inputs[n], accum[m][n]);
                }
            }
        }

        __syncthreads();
    }

    #pragma unroll
    for (int m = 0; m < 8; ++m) {
        const int oc = channel_group + m * 8;
        const float b = bias == nullptr ? 0.0f : bias[oc];
        const int64_t output_base =
            (int64_t(batch) * 64 + oc) * output_volume;

        #pragma unroll
        for (int n = 0; n < 4; ++n) {
            const int64_t p = spatial_base + lane + n * 32;
            if (p < output_volume) {
                output[output_base + p] = accum[m][n] + b;
            }
        }
    }
}

} // namespace

torch::Tensor conv3d_specialized(
    torch::Tensor x,
    torch::Tensor weight,
    c10::optional<torch::Tensor> bias)
{
    TORCH_CHECK(x.is_cuda(), "Expected CUDA input");
    TORCH_CHECK(x.scalar_type() == at::kFloat,
                "Expected float32 input");
    TORCH_CHECK(x.is_contiguous(), "Expected contiguous input");
    TORCH_CHECK(x.dim() == 5 && x.size(1) == 3,
                "Expected input shape [N, 3, D, H, W]");
    TORCH_CHECK(x.size(2) >= 3 && x.size(3) >= 3 && x.size(4) >= 3,
                "Input spatial dimensions must be at least 3");
    TORCH_CHECK(
        weight.device() == x.device() &&
        weight.scalar_type() == at::kFloat &&
        weight.is_contiguous() &&
        weight.dim() == 5 &&
        weight.size(0) == 64 &&
        weight.size(1) == 3 &&
        weight.size(2) == 3 &&
        weight.size(3) == 3 &&
        weight.size(4) == 3,
        "Expected contiguous float32 weight with shape [64, 3, 3, 3, 3]");

    const c10::cuda::CUDAGuard device_guard(x.device());

    const int64_t batch = x.size(0);
    const int64_t depth = x.size(2);
    const int64_t height = x.size(3);
    const int64_t width = x.size(4);
    const int64_t out_depth = depth - 2;
    const int64_t out_height = height - 2;
    const int64_t out_width = width - 2;
    const int64_t output_volume = out_depth * out_height * out_width;

    TORCH_CHECK(batch <= 65535, "Batch exceeds kernel grid limit");
    TORCH_CHECK(depth <= INT_MAX && height <= INT_MAX && width <= INT_MAX,
                "Spatial dimension exceeds kernel limit");
    TORCH_CHECK((output_volume + TILE_N - 1) / TILE_N <= INT_MAX,
                "Output exceeds kernel grid limit");

    const float* bias_ptr = nullptr;
    if (bias.has_value()) {
        const auto& b = bias.value();
        TORCH_CHECK(
            b.device() == x.device() &&
            b.scalar_type() == at::kFloat &&
            b.is_contiguous() &&
            b.dim() == 1 &&
            b.numel() == 64,
            "Expected contiguous float32 bias with shape [64]");
        bias_ptr = b.data_ptr<float>();
    }

    auto output = torch::empty(
        {batch, 64, out_depth, out_height, out_width}, x.options());

    if (batch == 0) {
        return output;
    }

    const dim3 block(32, 8);
    const dim3 grid(
        static_cast<unsigned int>((output_volume + TILE_N - 1) / TILE_N),
        static_cast<unsigned int>(batch));

    conv3d_3x3x3_c3_o64_kernel<<<
        grid, block, 0, at::cuda::getCurrentCUDAStream()
    >>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias_ptr,
        output.data_ptr<float>(),
        static_cast<int>(depth),
        static_cast<int>(height),
        static_cast<int>(width),
        static_cast<int>(out_height),
        static_cast<int>(out_width),
        output_volume);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = load_inline(
    name="conv3d_c3_o64_tiled_cuda",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["conv3d_specialized"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class _Conv3dFunction(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias):
        ctx.save_for_backward(x, weight)
        ctx.has_bias = bias is not None
        return _extension.conv3d_specialized(x, weight, bias)

    @staticmethod
    def backward(ctx, grad_output):
        x, weight = ctx.saved_tensors
        grad_output = grad_output.contiguous()
        grad_x = grad_weight = grad_bias = None

        if ctx.needs_input_grad[0]:
            grad_x = torch.nn.grad.conv3d_input(
                x.shape, weight, grad_output,
                stride=1, padding=0, dilation=1, groups=1,
            )

        if ctx.needs_input_grad[1]:
            grad_weight = torch.nn.grad.conv3d_weight(
                x, weight.shape, grad_output,
                stride=1, padding=0, dilation=1, groups=1,
            )

        if ctx.has_bias and ctx.needs_input_grad[2]:
            grad_bias = grad_output.sum(dim=(0, 2, 3, 4))

        return grad_x, grad_weight, grad_bias


class ModelNew(nn.Module):
    def __init__(
        self,
        in_channels: int,
        out_channels: int,
        kernel_size: int,
        stride: int = 1,
        padding: int = 0,
        dilation: int = 1,
        groups: int = 1,
        bias: bool = False,
    ):
        super().__init__()
        self.conv3d = nn.Conv3d(
            in_channels,
            out_channels,
            (kernel_size, kernel_size, kernel_size),
            stride=stride,
            padding=padding,
            dilation=dilation,
            groups=groups,
            bias=bias,
        )

        self._specialized = (
            in_channels == 3
            and out_channels == 64
            and self.conv3d.kernel_size == (3, 3, 3)
            and self.conv3d.stride == (1, 1, 1)
            and self.conv3d.padding == (0, 0, 0)
            and self.conv3d.dilation == (1, 1, 1)
            and groups == 1
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        weight = self.conv3d.weight
        bias = self.conv3d.bias

        if (
            self._specialized
            and x.is_cuda
            and x.dtype == torch.float32
            and x.dim() == 5
            and x.size(1) == 3
            and x.size(0) <= 65535
            and min(x.shape[2:]) >= 3
            and x.is_contiguous()
            and weight.dtype == torch.float32
            and weight.device == x.device
            and weight.is_contiguous()
            and (
                bias is None
                or (
                    bias.dtype == torch.float32
                    and bias.device == x.device
                    and bias.is_contiguous()
                )
            )
            and not torch.is_autocast_enabled()
        ):
            return _Conv3dFunction.apply(x, weight, bias)

        return self.conv3d(x)
```