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

torch.backends.cudnn.benchmark = True
torch.backends.cudnn.allow_tf32 = True

# Fused "bias fill + strided scatter" kernels.
# A stride-1 reduced conv (cuDNN, TF32 tensor cores) produces tmp in a compact
# coordinate space; these kernels expand it into the full strided
# transposed-conv output in a single bandwidth-bound pass: strided positions
# get conv value + bias, everything else gets bias/zero.
cuda_src = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda_runtime.h>

// Exact unsigned division by runtime constant via precomputed magic numbers:
// q = umulhi(pos, M) >> sh  (exact for 0 <= pos < 2^31, d >= 2; d==1 handled).
__device__ __forceinline__ void dvm(unsigned int pos, int d, unsigned int M, int sh,
                                    int& q, int& r) {
    if (d == 1) { q = (int)pos; r = 0; }
    else {
        q = (int)(__umulhi(pos, M) >> sh);
        r = (int)(pos - (unsigned int)(q * d));
    }
}

// Specialized kernel for sd == sh == sw == 2: cheap parity checks instead of
// full divisions. Each thread produces one aligned float4 (flat mapping, so
// stores are always 16B aligned); groups crossing a row boundary take a short
// scalar path. Requires total < 2^31.
__global__ void __launch_bounds__(256)
tconv_scatter_s2_kernel(
    const float* __restrict__ tmp, const float* __restrict__ bias, float* __restrict__ out,
    int Co, int Od, int Oh, int Ow,
    int MDc, int MHc, int MWc,
    int md0, int mh0, int mw0,
    int mdn, int mhn, int mwn,
    int pd, int ph, int pw,
    unsigned int MOw, int sOw,
    unsigned int MOh, int sOh,
    unsigned int MOd, int sOd,
    unsigned int MCo, int sCo,
    int ngroups, int total)
{
    const int gid = blockIdx.x * blockDim.x + threadIdx.x;
    if (gid >= ngroups) return;
    const int n0 = gid << 2;

    int t, ow0;
    dvm((unsigned int)n0, Ow, MOw, sOw, t, ow0);
    int t2, oh, z, od;
    dvm((unsigned int)t, Oh, MOh, sOh, t2, oh);
    dvm((unsigned int)t2, Od, MOd, sOd, z, od);

    float biasv = 0.0f;
    if (bias != nullptr) {
        int bb, co;
        dvm((unsigned int)z, Co, MCo, sCo, bb, co);
        biasv = __ldg(bias + co);
    }

    // Row activeness: valid iff (od+pd) even, (oh+ph) even, and in range.
    bool act = false;
    long long tbase = 0;
    if (((od ^ pd) & 1) == 0) {
        const int md = ((od + pd) >> 1) - md0;
        if ((unsigned int)md < (unsigned int)mdn) {
            if (((oh ^ ph) & 1) == 0) {
                const int mh = ((oh + ph) >> 1) - mh0;
                if ((unsigned int)mh < (unsigned int)mhn) {
                    act = true;
                    tbase = (((long long)z * MDc + md) * MHc + mh) * (long long)MWc;
                }
            }
        }
    }

    if (ow0 + 3 < Ow && n0 + 3 < total) {
        // Fast path: whole float4 inside one row.
        float4 v = make_float4(biasv, biasv, biasv, biasv);
        if (act) {
            const int q = (pw ^ ow0) & 1;   // valid lane parity within the group
            const int mw = ((ow0 + q + pw) >> 1) - mw0;
            if (q == 0) {
                if ((unsigned int)mw < (unsigned int)mwn)       v.x = biasv + __ldg(tmp + tbase + mw);
                if ((unsigned int)(mw + 1) < (unsigned int)mwn) v.z = biasv + __ldg(tmp + tbase + mw + 1);
            } else {
                if ((unsigned int)mw < (unsigned int)mwn)       v.y = biasv + __ldg(tmp + tbase + mw);
                if ((unsigned int)(mw + 1) < (unsigned int)mwn) v.w = biasv + __ldg(tmp + tbase + mw + 1);
            }
        }
        __stcs(reinterpret_cast<float4*>(out + n0), v);
    } else {
        // Slow path: group crosses a row boundary or the end of the tensor.
        #pragma unroll
        for (int l = 0; l < 4; ++l) {
            const int n = n0 + l;
            if (n >= total) break;
            const int ow = (l == 0) ? ow0 : ow0 + l;   // may wrap; recompute below if so
            int owr = ow0 + l;
            int ohr = oh, zr = z, odr = od;
            if (owr >= Ow) {                            // wrapped into next row(s)
                owr -= Ow;
                ohr = oh + 1;
                if (ohr >= Oh) { ohr -= Oh; odr = od + 1; if (odr >= Od) { odr -= Od; zr = z + 1; } }
            }
            float v = biasv;
            if (((odr ^ pd) & 1) == 0 && ((ohr ^ ph) & 1) == 0 && ((owr ^ pw) & 1) == 0) {
                const int md = ((odr + pd) >> 1) - md0;
                const int mh = ((ohr + ph) >> 1) - mh0;
                const int mw = ((owr + pw) >> 1) - mw0;
                if ((unsigned int)md < (unsigned int)mdn &&
                    (unsigned int)mh < (unsigned int)mhn &&
                    (unsigned int)mw < (unsigned int)mwn) {
                    v += __ldg(tmp + (((long long)zr * MDc + md) * MHc + mh) * (long long)MWc + mw);
                }
            }
            __stcs(out + n, v);
        }
    }
}

// Fast path: grid = (ceil(Oh*Ow/256), Od, B*Co). One output element per
// thread, coalesced streaming stores, magic-multiply division only.
__global__ void tconv_scatter_fast_kernel(
    const float* __restrict__ tmp, const float* __restrict__ bias, float* __restrict__ out,
    int Co, int Od, int Oh, int Ow,
    int MDc, int MHc, int MWc,
    int md0, int mh0, int mw0,
    int mdn, int mhn, int mwn,
    int sd, int sh_, int sw,
    int pd, int ph, int pw,
    unsigned int MCo, int sCo,
    unsigned int MOw, int sOw,
    unsigned int Msd, int ssd,
    unsigned int Msh, int ssh,
    unsigned int Msw, int ssw)
{
    const int z = blockIdx.z;               // b*Co + co
    const int od = blockIdx.y;
    const int p = blockIdx.x * blockDim.x + threadIdx.x;
    const int plane = Oh * Ow;
    if (p >= plane) return;

    int oh, ow;
    dvm((unsigned int)p, Ow, MOw, sOw, oh, ow);

    int bidx, co;
    dvm((unsigned int)z, Co, MCo, sCo, bidx, co);

    float v = (bias != nullptr) ? __ldg(bias + co) : 0.0f;

    int qd, rd;
    dvm((unsigned int)(od + pd), sd, Msd, ssd, qd, rd);
    if (rd == 0) {
        int qh, rh;
        dvm((unsigned int)(oh + ph), sh_, Msh, ssh, qh, rh);
        if (rh == 0) {
            const int md = qd - md0;
            const int mh = qh - mh0;
            if ((unsigned int)md < (unsigned int)mdn &&
                (unsigned int)mh < (unsigned int)mhn) {
                int qw, rw;
                dvm((unsigned int)(ow + pw), sw, Msw, ssw, qw, rw);
                if (rw == 0) {
                    const int mw = qw - mw0;
                    if ((unsigned int)mw < (unsigned int)mwn) {
                        const long long tb =
                            (((long long)z * MDc + md) * MHc + mh) * (long long)MWc;
                        v += __ldg(tmp + tb + mw);
                    }
                }
            }
        }
    }
    const long long ob = ((long long)z * Od + od) * (long long)plane + p;
    __stwt(out + ob, v);
}

// Generic flat 64-bit fallback for exotic shapes (huge B*Co or Od).
__global__ void tconv_scatter_kernel(
    const float* __restrict__ tmp, const float* __restrict__ bias, float* __restrict__ out,
    int Co, int Od, int Oh, int Ow,
    int MDc, int MHc, int MWc,
    int md0, int mh0, int mw0,
    int mdn, int mhn, int mwn,
    int sd, int sh_, int sw,
    int pd, int ph, int pw,
    long long total)
{
    const long long idx = (long long)blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= total) return;

    const int ow = (int)(idx % Ow);
    const long long t1 = idx / Ow;
    const int oh = (int)(t1 % Oh);
    const long long t2 = t1 / Oh;
    const int od = (int)(t2 % Od);
    const long long t3 = t2 / Od;
    const int co = (int)(t3 % Co);

    float v = (bias != nullptr) ? __ldg(bias + co) : 0.0f;
    const int td = od + pd, th = oh + ph, tw = ow + pw;
    if ((td % sd == 0) && (th % sh_ == 0) && (tw % sw == 0)) {
        const int md = td / sd - md0;
        const int mh = th / sh_ - mh0;
        const int mw = tw / sw - mw0;
        if ((unsigned int)md < (unsigned int)mdn &&
            (unsigned int)mh < (unsigned int)mhn &&
            (unsigned int)mw < (unsigned int)mwn) {
            v += tmp[(((long long)t3 * MDc + md) * MHc + mh) * (long long)MWc + mw];
        }
    }
    out[idx] = v;
}

static inline void magic_u32(int d, unsigned int& M, int& sh) {
    if (d <= 1) { M = 0; sh = 0; return; }
    int l = 0;
    while ((1LL << l) < d) ++l;                 // ceil(log2(d))
    unsigned long long T = 31 + l;              // 2^T/d >= 2^31 -> exact for pos < 2^31
    M = (unsigned int)(((1ULL << T) + (unsigned long long)d - 1) / (unsigned long long)d);
    sh = l - 1;
}

torch::Tensor tconv_scatter(
    torch::Tensor tmp, torch::Tensor bias,
    int64_t B, int64_t Co,
    int64_t Od, int64_t Oh, int64_t Ow,
    int64_t md0, int64_t mh0, int64_t mw0,
    int64_t mdn, int64_t mhn, int64_t mwn,
    int64_t sd, int64_t sh, int64_t sw,
    int64_t pd, int64_t ph, int64_t pw)
{
    auto out = torch::empty({B, Co, Od, Oh, Ow}, tmp.options());
    const long long total = (long long)out.numel();
    if (total == 0) return out;

    const bool has_bias = bias.defined() && bias.numel() > 0;
    const float* bias_ptr = has_bias ? bias.data_ptr<float>() : nullptr;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    const int MDc = (int)tmp.size(2), MHc = (int)tmp.size(3), MWc = (int)tmp.size(4);
    const long long plane = Oh * Ow;

    if (sd == 2 && sh == 2 && sw == 2 && total < (1LL << 31) - 8) {
        // Specialized parity-based float4 kernel.
        unsigned int MOw, MOh, MOd, MCo;
        int sOw, sOh, sOd, sCo;
        magic_u32((int)Ow, MOw, sOw);
        magic_u32((int)Oh, MOh, sOh);
        magic_u32((int)Od, MOd, sOd);
        magic_u32((int)Co, MCo, sCo);

        const int ngroups = (int)((total + 3) >> 2);
        const int threads = 256;
        const int blocks = (ngroups + threads - 1) / threads;
        tconv_scatter_s2_kernel<<<blocks, threads, 0, stream>>>(
            tmp.data_ptr<float>(), bias_ptr, out.data_ptr<float>(),
            (int)Co, (int)Od, (int)Oh, (int)Ow,
            MDc, MHc, MWc,
            (int)md0, (int)mh0, (int)mw0,
            (int)mdn, (int)mhn, (int)mwn,
            (int)pd, (int)ph, (int)pw,
            MOw, sOw, MOh, sOh, MOd, sOd, MCo, sCo,
            ngroups, (int)total);
    } else if (B * Co <= 65535 && Od <= 65535 && plane > 0 && plane < 0x7fffffffLL) {
        unsigned int MCo, MOw, Msd, Msh, Msw;
        int sCo, sOw, ssd, ssh, ssw;
        magic_u32((int)Co, MCo, sCo);
        magic_u32((int)Ow, MOw, sOw);
        magic_u32((int)sd, Msd, ssd);
        magic_u32((int)sh, Msh, ssh);
        magic_u32((int)sw, Msw, ssw);

        const int threads = 256;
        dim3 grid((unsigned int)((plane + threads - 1) / threads),
                  (unsigned int)Od, (unsigned int)(B * Co));
        tconv_scatter_fast_kernel<<<grid, threads, 0, stream>>>(
            tmp.data_ptr<float>(), bias_ptr, out.data_ptr<float>(),
            (int)Co, (int)Od, (int)Oh, (int)Ow,
            MDc, MHc, MWc,
            (int)md0, (int)mh0, (int)mw0,
            (int)mdn, (int)mhn, (int)mwn,
            (int)sd, (int)sh, (int)sw,
            (int)pd, (int)ph, (int)pw,
            MCo, sCo, MOw, sOw, Msd, ssd, Msh, ssh, Msw, ssw);
    } else {
        const int threads = 256;
        const long long blocks = (total + threads - 1) / threads;
        tconv_scatter_kernel<<<(unsigned int)blocks, threads, 0, stream>>>(
            tmp.data_ptr<float>(), bias_ptr, out.data_ptr<float>(),
            (int)Co, (int)Od, (int)Oh, (int)Ow,
            MDc, MHc, MWc,
            (int)md0, (int)mh0, (int)mw0,
            (int)mdn, (int)mhn, (int)mwn,
            (int)sd, (int)sh, (int)sw,
            (int)pd, (int)ph, (int)pw,
            total);
    }
    return out;
}
"""

cpp_src = r"""
torch::Tensor tconv_scatter(
    torch::Tensor tmp, torch::Tensor bias,
    int64_t B, int64_t Co,
    int64_t Od, int64_t Oh, int64_t Ow,
    int64_t md0, int64_t mh0, int64_t mw0,
    int64_t mdn, int64_t mhn, int64_t mwn,
    int64_t sd, int64_t sh, int64_t sw,
    int64_t pd, int64_t ph, int64_t pw);
"""

tconv_scatter = load_inline(
    name="tconv_scatter_v6",
    cpp_sources=cpp_src,
    cuda_sources=cuda_src,
    functions=["tconv_scatter"],
    verbose=False,
    extra_cuda_cflags=["-O3"],
)


class ModelNew(nn.Module):
    """
    ConvTranspose3d implemented (when dilation is an integer multiple of
    stride per dimension) as:
      1) a regular stride-1 cuDNN conv3d on a cropped view of the input in a
         "reduced" coordinate space (kernel flipped, dilation = dil/stride),
         which runs on tensor cores with TF32, followed by
      2) a custom fused CUDA kernel that writes the full output in one
         bandwidth-bound pass: strided positions get conv value + bias,
         everything else gets bias/zero. A specialized parity-based float4
         kernel handles the common stride-2 case.
    Falls back to the PyTorch op for unsupported configurations.
    """

    def __init__(self, in_channels: int, out_channels: int, kernel_size: int,
                 stride: int = 1, padding: int = 0, dilation: int = 1,
                 bias: bool = False):
        super(ModelNew, self).__init__()
        self.conv_transpose3d = nn.ConvTranspose3d(
            in_channels, out_channels,
            kernel_size=(kernel_size, kernel_size, kernel_size),
            stride=stride, padding=padding, dilation=dilation, bias=bias)
        s = self.conv_transpose3d.stride
        p = self.conv_transpose3d.padding
        d = self.conv_transpose3d.dilation
        self.sd, self.sh, self.sw = int(s[0]), int(s[1]), int(s[2])
        self.pd, self.ph, self.pw = int(p[0]), int(p[1]), int(p[2])
        self.dd, self.dh, self.dw = int(d[0]), int(d[1]), int(d[2])
        self.kernel_size = int(kernel_size)

        self._fast = (self.dd % self.sd == 0 and
                      self.dh % self.sh == 0 and
                      self.dw % self.sw == 0)
        self.rd = self.dd // self.sd
        self.rh = self.dh // self.sh
        self.rw = self.dw // self.sw

        self._wf = None
        self._wf_ver = -1
        self._geom = {}
        self._empty_bias = None

    def _geometry(self, D, H, W):
        key = (D, H, W)
        g = self._geom.get(key)
        if g is not None:
            return g
        K = self.kernel_size

        def axis(n, s, p, dil, r):
            O = (n - 1) * s - 2 * p + dil * (K - 1) + 1
            m0 = (p + s - 1) // s
            mn = (O - 1 + p) // s - m0 + 1
            rK = r * (K - 1)
            i0 = m0 - rK
            i1 = m0 + mn - 1
            c0 = i0 if i0 > 0 else 0
            c1 = i1 + 1
            if c1 > n:
                c1 = n
            pp = rK + c0 - m0
            return O, m0, mn, c0, c1, pp

        g = (axis(D, self.sd, self.pd, self.dd, self.rd),
             axis(H, self.sh, self.ph, self.dh, self.rh),
             axis(W, self.sw, self.pw, self.dw, self.rw))
        self._geom[key] = g
        return g

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        if (not self._fast) or (not x.is_cuda) or x.dtype != torch.float32:
            return self.conv_transpose3d(x)

        w = self.conv_transpose3d.weight
        bias = self.conv_transpose3d.bias

        wf = self._wf
        if wf is None or self._wf_ver != w._version or wf.device != w.device:
            # [Ci, Co, Kd, Kh, Kw] -> flip kernel -> [Co, Ci, Kd, Kh, Kw]
            wf = w.detach().flip([2, 3, 4]).permute(1, 0, 2, 3, 4).contiguous()
            self._wf = wf
            self._wf_ver = w._version

        B = x.shape[0]
        D, H, W = int(x.shape[2]), int(x.shape[3]), int(x.shape[4])
        (Od, md0, mdn, dc0, dc1, pdp), \
            (Oh, mh0, mhn, hc0, hc1, php), \
            (Ow, mw0, mwn, wc0, wc1, pwp) = self._geometry(D, H, W)

        Co = self.conv_transpose3d.out_channels

        if mdn <= 0 or mhn <= 0 or mwn <= 0 or \
                dc1 <= dc0 or hc1 <= hc0 or wc1 <= wc0 or \
                Od <= 0 or Oh <= 0 or Ow <= 0:
            out = torch.empty((B, Co, Od, Oh, Ow), device=x.device, dtype=x.dtype)
            if out.numel() > 0:
                if bias is None:
                    out.zero_()
                else:
                    out.copy_(bias.view(1, -1, 1, 1, 1))
            return out

        xc = x[:, :, dc0:dc1, hc0:hc1, wc0:wc1]
        if not xc.is_contiguous():
            xc = xc.contiguous()

        # No bias here: bias is added by the scatter kernel (also at the
        # non-strided positions), avoiding a double-add.
        tmp = F.conv3d(xc, wf, None, stride=1,
                       padding=(pdp, php, pwp),
                       dilation=(self.rd, self.rh, self.rw))
        if not tmp.is_contiguous():
            tmp = tmp.contiguous()

        # Safety: never index outside the produced conv output.
        mdn = min(mdn, int(tmp.shape[2]))
        mhn = min(mhn, int(tmp.shape[3]))
        mwn = min(mwn, int(tmp.shape[4]))

        if bias is None:
            eb = self._empty_bias
            if eb is None or eb.device != x.device:
                eb = torch.empty(0, device=x.device, dtype=x.dtype)
                self._empty_bias = eb
            bias = eb

        return tconv_scatter.tconv_scatter(
            tmp, bias, B, Co, Od, Oh, Ow,
            md0, mh0, mw0, mdn, mhn, mwn,
            self.sd, self.sh, self.sw,
            self.pd, self.ph, self.pw)
```