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


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

torch::Tensor transpose_conv64_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    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.h>
#include <cuda_runtime.h>
#include <mma.h>

namespace {

constexpr int TILE_M = 128;
constexpr int TILE_N = 64;
constexpr int TILE_K = 32;

union SharedStorage {
    struct {
        float ah[TILE_M * TILE_K];
        float al[TILE_M * TILE_K];
        float bh[TILE_K * TILE_N];
        float bl[TILE_K * TILE_N];
    } operands;
    float result[TILE_M * TILE_N];
};

__global__ __launch_bounds__(256)
void transpose_conv64_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int height,
    int width,
    int output_height,
    int output_width)
{
#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 800
    using namespace nvcuda;

    __shared__ SharedStorage shared;

    const int tid = threadIdx.x;
    const int warp = tid / 32;
    const int batch = blockIdx.y;
    const int input_area = height * width;
    const int output_area = output_height * output_width;
    const int spatial_start = blockIdx.x * TILE_M;

    wmma::fragment<
        wmma::accumulator, 16, 16, 8, float
    > accum[4];

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        wmma::fill_fragment(accum[j], 0.0f);
    }

    const int spatial_row = tid % TILE_M;
    const int spatial_index = spatial_start + spatial_row;
    const int oh = spatial_index / output_width;
    const int ow = spatial_index % output_width;

    for (int kh = 0; kh < 3; ++kh) {
        const int ih = oh - kh;

        for (int kw = 0; kw < 3; ++kw) {
            const int iw = ow - kw;
            const bool valid =
                spatial_index < output_area &&
                static_cast<unsigned>(ih) <
                    static_cast<unsigned>(height) &&
                static_cast<unsigned>(iw) <
                    static_cast<unsigned>(width);

            for (int channel_base = 0; channel_base < 64;
                 channel_base += TILE_K) {

                // Neighboring threads read neighboring spatial positions.
                #pragma unroll
                for (int k = tid / TILE_M; k < TILE_K; k += 2) {
                    float value = 0.0f;
                    if (valid) {
                        const int channel = channel_base + k;
                        const int64_t offset =
                            (static_cast<int64_t>(batch) * 64 + channel) *
                            input_area + ih * width + iw;
                        value = x[offset];
                    }

                    const float high = __float_to_tf32(value);
                    const float low = __float_to_tf32(value - high);
                    const int index = spatial_row * TILE_K + k;

                    shared.operands.ah[index] = high;
                    shared.operands.al[index] = low;
                }

                #pragma unroll
                for (int index = tid; index < TILE_K * TILE_N;
                     index += 256) {
                    const int k = index / TILE_N;
                    const int output_channel = index % TILE_N;
                    const int input_channel = channel_base + k;
                    const float value =
                        weight[(input_channel * 64 + output_channel) * 9
                               + kh * 3 + kw];

                    const float high = __float_to_tf32(value);
                    shared.operands.bh[index] = high;
                    shared.operands.bl[index] =
                        __float_to_tf32(value - high);
                }

                __syncthreads();

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

                    wmma::load_matrix_sync(
                        ah,
                        shared.operands.ah + warp * 16 * TILE_K + k,
                        TILE_K);
                    wmma::load_matrix_sync(
                        al,
                        shared.operands.al + warp * 16 * TILE_K + k,
                        TILE_K);

                    #pragma unroll
                    for (int j = 0; j < 4; ++j) {
                        wmma::fragment<
                            wmma::matrix_b, 16, 16, 8,
                            wmma::precision::tf32, wmma::row_major
                        > bh, bl;

                        wmma::load_matrix_sync(
                            bh,
                            shared.operands.bh + k * TILE_N + j * 16,
                            TILE_N);
                        wmma::load_matrix_sync(
                            bl,
                            shared.operands.bl + k * TILE_N + j * 16,
                            TILE_N);

                        // Three products recover substantially more precision
                        // than a single TF32 multiplication.
                        wmma::mma_sync(accum[j], ah, bl, accum[j]);
                        wmma::mma_sync(accum[j], al, bh, accum[j]);
                        wmma::mma_sync(accum[j], ah, bh, accum[j]);
                    }
                }

                __syncthreads();
            }
        }
    }

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        wmma::store_matrix_sync(
            shared.result + warp * 16 * TILE_N + j * 16,
            accum[j],
            TILE_N,
            wmma::mem_row_major);
    }

    __syncthreads();

    // Transpose the tile into coalesced NCHW output stores.
    #pragma unroll
    for (int channel = tid / TILE_M; channel < TILE_N; channel += 2) {
        if (spatial_index < output_area) {
            float value =
                shared.result[spatial_row * TILE_N + channel];
            if (bias != nullptr) {
                value += bias[channel];
            }

            const int64_t offset =
                (static_cast<int64_t>(batch) * 64 + channel) *
                output_area + spatial_index;
            output[offset] = value;
        }
    }
#endif
}

} // namespace

torch::Tensor transpose_conv64_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    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) == 64,
                "Expected NCHW input with 64 channels");
    TORCH_CHECK(weight.dim() == 4 &&
                weight.size(0) == 64 && weight.size(1) == 64 &&
                weight.size(2) == 3 && weight.size(3) == 3,
                "Expected weight shape (64, 64, 3, 3)");
    TORCH_CHECK(weight.device() == x.device() &&
                bias.device() == x.device(),
                "All tensors must be on the same device");
    TORCH_CHECK(bias.numel() == 0 ||
                (bias.numel() == 64 &&
                 bias.scalar_type() == at::kFloat &&
                 bias.is_contiguous()),
                "Invalid bias");

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

    const int height = static_cast<int>(x.size(2));
    const int width = static_cast<int>(x.size(3));
    const int output_height = height + 2;
    const int output_width = width + 2;
    const int64_t output_area =
        static_cast<int64_t>(output_height) * output_width;

    TORCH_CHECK(output_area <= 2147483647LL,
                "Spatial dimensions exceed kernel indexing limits");

    auto output = torch::empty(
        {x.size(0), 64, output_height, output_width}, x.options());

    if (x.size(0) == 0) {
        return output;
    }

    const dim3 grid(
        static_cast<unsigned>((output_area + TILE_M - 1) / TILE_M),
        static_cast<unsigned>(x.size(0)));

    transpose_conv64_kernel<<<
        grid, 256, 0, at::cuda::getCurrentCUDAStream()
    >>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.numel() ? bias.data_ptr<float>() : nullptr,
        output.data_ptr<float>(),
        height,
        width,
        output_height,
        output_width);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = None


def _get_extension():
    global _extension
    if _extension is None:
        _extension = load_inline(
            name="transpose_conv64_tf32_compensated_v1",
            cpp_sources=_cpp_source,
            cuda_sources=_cuda_source,
            functions=["transpose_conv64_cuda"],
            extra_cflags=["-O3"],
            extra_cuda_cflags=["-O3"],
            verbose=False,
        )
    return _extension


class _TransposeConv64(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias):
        if any(ctx.needs_input_grad):
            ctx.save_for_backward(x, weight, bias)
        return _get_extension().transpose_conv64_cuda(x, weight, bias)

    @staticmethod
    @once_differentiable
    def backward(ctx, grad_output):
        x, weight, bias = ctx.saved_tensors
        needed = ctx.needs_input_grad

        with torch.enable_grad():
            x_copy = x.detach().requires_grad_(needed[0])
            weight_copy = weight.detach().requires_grad_(needed[1])
            bias_copy = bias.detach().requires_grad_(needed[2])

            output = F.conv_transpose2d(
                x_copy,
                weight_copy,
                bias_copy if bias.numel() else None,
            )

            candidates = (x_copy, weight_copy, bias_copy)
            indices = [
                i for i in range(3)
                if needed[i] and (i != 2 or bias.numel() != 0)
            ]
            gradients = torch.autograd.grad(
                output,
                tuple(candidates[i] for i in indices),
                grad_output,
            )

        result = [None, None, None]
        for index, gradient in zip(indices, gradients):
            result[index] = gradient
        return tuple(result)


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,
        )
        self.register_buffer("_empty_bias", torch.empty(0), persistent=False)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        conv = self.conv_transpose2d

        supported = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 4
            and x.is_contiguous()
            and conv.weight.dtype == torch.float32
            and conv.weight.is_contiguous()
            and conv.in_channels == 64
            and conv.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 conv.dilation == (1, 1)
            and conv.groups == 1
            and x.shape[0] <= 65535
            and x.shape[2] > 0
            and x.shape[3] > 0
            and not torch.is_autocast_enabled()
        )

        if supported and torch.cuda.get_device_capability(x.device)[0] >= 8:
            bias = conv.bias if conv.bias is not None else self._empty_bias
            return _TransposeConv64.apply(x, conv.weight, bias)

        return conv(x)
```