```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 transposed_conv_cuda(
    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>
#include <mma.h>

namespace wmma = nvcuda::wmma;

constexpr int CHANNELS = 64;
constexpr int REDUCTION = 64 * 3 * 3;
constexpr int TILE_M = 32;
constexpr int TILE_N = 64;
constexpr int TILE_K = 32;

// Pack the original [input_channel, output_channel, kh, kw] weights
// into a matrix, retaining a TF32 high part and a TF32 residual.
__global__ void pack_weights(
    const float* __restrict__ weight,
    float* __restrict__ packed)
{
#if __CUDA_ARCH__ >= 800
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index < CHANNELS * REDUCTION) {
        int co = index / REDUCTION;
        int k = index % REDUCTION;
        int ci = k % CHANNELS;
        int tap = k / CHANNELS;
        float value = weight[(ci * CHANNELS + co) * 9 + tap];
        float high = wmma::__float_to_tf32(value);
        packed[index] = high;
        packed[CHANNELS * REDUCTION + index] =
            wmma::__float_to_tf32(value - high);
    }
#endif
}

__global__ void transposed_conv_kernel(
    const float* __restrict__ input,
    const float* __restrict__ weights,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int height,
    int width)
{
#if __CUDA_ARCH__ >= 800
    __shared__ __align__(32) float a_high[TILE_M * TILE_K];
    __shared__ __align__(32) float a_low[TILE_M * TILE_K];
    __shared__ __align__(32) float b_high[TILE_K * TILE_N];
    __shared__ __align__(32) float b_low[TILE_K * TILE_N];
    __shared__ __align__(32) float result[TILE_M * TILE_N];

    const int tid = threadIdx.x;
    const int warp = tid / 32;
    const int warp_m = warp / 4;
    const int warp_n = warp % 4;

    const int channel_base = blockIdx.y * TILE_M;
    const int64_t position_base = int64_t(blockIdx.x) * TILE_N;
    const int batch = blockIdx.z;
    const int output_height = height + 2;
    const int output_width = width + 2;
    const int64_t input_area = int64_t(height) * width;
    const int64_t output_area =
        int64_t(output_height) * output_width;

    wmma::fragment<wmma::accumulator, 16, 16, 8, float> accumulator;
    wmma::fill_fragment(accumulator, 0.0f);

    for (int k_base = 0; k_base < REDUCTION; k_base += TILE_K) {
        for (int index = tid; index < TILE_M * TILE_K; index += 256) {
            int row = index / TILE_K;
            int k = index % TILE_K;
            int offset = (channel_base + row) * REDUCTION + k_base + k;
            a_high[index] = weights[offset];
            a_low[index] = weights[CHANNELS * REDUCTION + offset];
        }

        for (int index = tid; index < TILE_K * TILE_N; index += 256) {
            int k = k_base + index / TILE_N;
            int64_t position = position_base + index % TILE_N;
            int ci = k % CHANNELS;
            int tap = k / CHANNELS;

            float value = 0.0f;
            if (position < output_area) {
                int oy = int(position / output_width);
                int ox = int(position % output_width);
                int iy = oy - tap / 3;
                int ix = ox - tap % 3;

                if (iy >= 0 && iy < height && ix >= 0 && ix < width) {
                    int64_t offset =
                        (int64_t(batch) * CHANNELS + ci) * input_area
                        + int64_t(iy) * width + ix;
                    value = input[offset];
                }
            }

            float high = wmma::__float_to_tf32(value);
            b_high[index] = high;
            b_low[index] = wmma::__float_to_tf32(value - high);
        }
        __syncthreads();

        #pragma unroll
        for (int step = 0; step < TILE_K; step += 8) {
            wmma::fragment<
                wmma::matrix_a, 16, 16, 8,
                wmma::precision::tf32, wmma::row_major> ah, al;
            wmma::fragment<
                wmma::matrix_b, 16, 16, 8,
                wmma::precision::tf32, wmma::row_major> bh, bl;

            int a_offset = warp_m * 16 * TILE_K + step;
            int b_offset = step * TILE_N + warp_n * 16;

            wmma::load_matrix_sync(ah, a_high + a_offset, TILE_K);
            wmma::load_matrix_sync(al, a_low + a_offset, TILE_K);
            wmma::load_matrix_sync(bh, b_high + b_offset, TILE_N);
            wmma::load_matrix_sync(bl, b_low + b_offset, TILE_N);

            // Residual products recover substantially more precision
            // than a single TF32 matrix multiplication.
            wmma::mma_sync(accumulator, al, bh, accumulator);
            wmma::mma_sync(accumulator, ah, bl, accumulator);
            wmma::mma_sync(accumulator, ah, bh, accumulator);
        }
        __syncthreads();
    }

    wmma::store_matrix_sync(
        result + warp_m * 16 * TILE_N + warp_n * 16,
        accumulator,
        TILE_N,
        wmma::mem_row_major);
    __syncthreads();

    for (int index = tid; index < TILE_M * TILE_N; index += 256) {
        int co = channel_base + index / TILE_N;
        int64_t position = position_base + index % TILE_N;
        if (position < output_area) {
            float value = result[index];
            if (bias != nullptr) {
                value += bias[co];
            }
            output[(int64_t(batch) * CHANNELS + co) * output_area
                   + position] = value;
        }
    }
#endif
}

torch::Tensor transposed_conv_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    c10::optional<torch::Tensor> bias)
{
    TORCH_CHECK(x.is_cuda() && weight.is_cuda(), "Expected CUDA tensors");
    TORCH_CHECK(
        x.scalar_type() == at::kFloat &&
        weight.scalar_type() == at::kFloat,
        "Expected float32 tensors");
    TORCH_CHECK(x.is_contiguous() && weight.is_contiguous(),
                "Expected contiguous tensors");
    TORCH_CHECK(x.dim() == 4 && x.size(1) == CHANNELS,
                "Expected NCHW input with 64 channels");
    TORCH_CHECK(
        weight.dim() == 4 &&
        weight.size(0) == CHANNELS &&
        weight.size(1) == CHANNELS &&
        weight.size(2) == 3 &&
        weight.size(3) == 3,
        "Expected weight shape [64, 64, 3, 3]");
    TORCH_CHECK(x.device() == weight.device(), "Device mismatch");

    c10::cuda::CUDAGuard guard(x.device());
    TORCH_CHECK(at::cuda::getCurrentDeviceProperties()->major >= 8,
                "This kernel requires compute capability 8.0 or newer");

    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.numel() == CHANNELS,
            "Invalid bias");
        bias_ptr = b.data_ptr<float>();
    }

    const int batch = int(x.size(0));
    const int height = int(x.size(2));
    const int width = int(x.size(3));
    auto output = torch::empty(
        {x.size(0), CHANNELS, x.size(2) + 2, x.size(3) + 2},
        x.options());

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

    auto packed = torch::empty({2, CHANNELS, REDUCTION}, weight.options());
    cudaStream_t stream = c10::cuda::getCurrentCUDAStream();

    pack_weights<<<(CHANNELS * REDUCTION + 255) / 256, 256, 0, stream>>>(
        weight.data_ptr<float>(),
        packed.data_ptr<float>());
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    int64_t output_area = int64_t(height + 2) * (width + 2);
    dim3 grid(
        unsigned((output_area + TILE_N - 1) / TILE_N),
        CHANNELS / TILE_M,
        unsigned(batch));

    transposed_conv_kernel<<<grid, 256, 0, stream>>>(
        x.data_ptr<float>(),
        packed.data_ptr<float>(),
        bias_ptr,
        output.data_ptr<float>(),
        height,
        width);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

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


class ModelNew(nn.Module):
    def __init__(
        self,
        in_channels: int,
        out_channels: int,
        kernel_size: int,
        stride: int = 1,
        padding: int = 0,
        output_padding: int = 0,
        groups: int = 1,
        bias: bool = False,
    ):
        super().__init__()
        self.conv_transpose2d = nn.ConvTranspose2d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
            output_padding=output_padding,
            groups=groups,
            bias=bias,
        )
        conv = self.conv_transpose2d
        self._specialized = (
            in_channels == 64
            and out_channels == 64
            and conv.kernel_size == (3, 3)
            and conv.stride == (1, 1)
            and conv.padding == (0, 0)
            and conv.output_padding == (0, 0)
            and groups == 1
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        conv = self.conv_transpose2d
        needs_grad = torch.is_grad_enabled() and (
            x.requires_grad
            or conv.weight.requires_grad
            or (conv.bias is not None and conv.bias.requires_grad)
        )
        if (
            self._specialized
            and not needs_grad
            and x.is_cuda
            and x.dtype == torch.float32
            and conv.weight.dtype == torch.float32
            and x.is_contiguous()
            and conv.weight.is_contiguous()
            and x.ndim == 4
            and x.shape[1] == 64
            and x.shape[0] <= 65535
            and x.shape[2] > 0
            and x.shape[3] > 0
            and not torch.is_autocast_enabled()
            and torch.cuda.get_device_capability(x.device)[0] >= 8
        ):
            return _extension.transposed_conv_cuda(x, conv.weight, conv.bias)
        return conv(x)
```