```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 bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    double eps);
"""

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

constexpr int THREADS = 256;
constexpr int CHUNK = 8192;

struct Welford {
    float mean;
    float m2;
    int count;
};

__device__ __forceinline__ Welford combine(Welford a, Welford b) {
    if (b.count == 0) return a;
    if (a.count == 0) return b;

    const int count = a.count + b.count;
    const float delta = b.mean - a.mean;
    const float fraction = float(b.count) / float(count);

    a.mean += delta * fraction;
    a.m2 += b.m2 + delta * delta * float(a.count) * fraction;
    a.count = count;
    return a;
}

__device__ __forceinline__ Welford warp_reduce(Welford a) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Welford b;
        b.mean = __shfl_down_sync(0xffffffff, a.mean, offset);
        b.m2 = __shfl_down_sync(0xffffffff, a.m2, offset);
        b.count = __shfl_down_sync(0xffffffff, a.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            a = combine(a, b);
    }
    return a;
}

__device__ __forceinline__ Welford block_reduce(Welford a) {
    __shared__ float means[8];
    __shared__ float m2s[8];
    __shared__ int counts[8];

    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    a = warp_reduce(a);
    if (lane == 0) {
        means[warp] = a.mean;
        m2s[warp] = a.m2;
        counts[warp] = a.count;
    }
    __syncthreads();

    Welford result = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8) {
            result.mean = means[lane];
            result.m2 = m2s[lane];
            result.count = counts[lane];
        }
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_statistics(
    const float* __restrict__ x,
    float2* __restrict__ partials,
    int channels,
    int spatial,
    int chunks,
    int batch)
{
    const int part = blockIdx.x;
    const int channel = blockIdx.y;
    const int n = part / chunks;
    const int chunk = part - n * chunks;

    const int begin = chunk * CHUNK;
    const int end = min(begin + CHUNK, spatial);
    const int64_t base = (int64_t(n) * channels + channel) * spatial;

    Welford acc = {0.0f, 0.0f, 0};
    for (int i = begin + threadIdx.x; i < end; i += THREADS) {
        const float value = x[base + i];
        ++acc.count;
        const float delta = value - acc.mean;
        acc.mean += delta / float(acc.count);
        acc.m2 += delta * (value - acc.mean);
    }

    acc = block_reduce(acc);
    if (threadIdx.x == 0)
        partials[int64_t(channel) * batch * chunks + part] =
            make_float2(acc.mean, acc.m2);
}

__global__ void finalize_statistics(
    const float2* __restrict__ partials,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float2* __restrict__ affine,
    int batch,
    int spatial,
    int chunks,
    bool training,
    float momentum,
    float eps)
{
    const int channel = blockIdx.x;

    if (!training) {
        if (threadIdx.x == 0) {
            const float scale =
                weight[channel] * rsqrtf(running_var[channel] + eps);
            affine[channel] =
                make_float2(running_mean[channel], scale);
        }
        return;
    }

    const int parts = batch * chunks;
    Welford acc = {0.0f, 0.0f, 0};

    for (int part = threadIdx.x; part < parts; part += THREADS) {
        const int chunk = part % chunks;
        const float2 value = partials[int64_t(channel) * parts + part];
        Welford next;
        next.mean = value.x;
        next.m2 = value.y;
        next.count = min(CHUNK, spatial - chunk * CHUNK);
        acc = combine(acc, next);
    }

    acc = block_reduce(acc);

    if (threadIdx.x == 0) {
        const float variance = fmaxf(acc.m2 / float(acc.count), 0.0f);
        const float unbiased = fmaxf(
            acc.m2 / float(acc.count - 1), 0.0f);

        running_mean[channel] =
            (1.0f - momentum) * running_mean[channel] +
            momentum * acc.mean;
        running_var[channel] =
            (1.0f - momentum) * running_var[channel] +
            momentum * unbiased;

        affine[channel] = make_float2(
            acc.mean, weight[channel] * rsqrtf(variance + eps));
    }
}

__global__ void normalized_pool4(
    const float* __restrict__ x,
    const float2* __restrict__ affine,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int channels,
    int depth,
    int height,
    int width,
    int out_depth,
    int out_height,
    int out_width,
    int64_t elements)
{
    const int64_t out_spatial =
        int64_t(out_depth) * out_height * out_width;
    const int64_t in_spatial = int64_t(depth) * height * width;

    for (int64_t index = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         index < elements;
         index += int64_t(blockDim.x) * gridDim.x) {
        const int64_t nc = index / out_spatial;
        const int channel = nc % channels;
        int position = int(index - nc * out_spatial);

        const int ow = position % out_width;
        position /= out_width;
        const int oh = position % out_height;
        const int od = position / out_height;

        const int64_t base =
            nc * in_spatial +
            int64_t(od * 4) * height * width +
            int64_t(oh * 4) * width + ow * 4;

        float sum = 0.0f;
        #pragma unroll
        for (int dz = 0; dz < 4; ++dz) {
            #pragma unroll
            for (int dy = 0; dy < 4; ++dy) {
                const int64_t p =
                    base + int64_t(dz) * height * width + dy * width;
                sum += x[p] + x[p + 1] + x[p + 2] + x[p + 3];
            }
        }

        const float2 stats = affine[channel];
        output[index] =
            (sum * (1.0f / 64.0f) - stats.x) * stats.y + bias[channel];
    }
}

torch::Tensor bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    double eps)
{
    c10::cuda::CUDAGuard guard(x.device());
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == torch::kFloat32,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");

    const int batch = x.size(0);
    const int channels = x.size(1);
    const int depth = x.size(2);
    const int height = x.size(3);
    const int width = x.size(4);
    const int spatial = depth * height * width;
    const int chunks = (spatial + CHUNK - 1) / CHUNK;
    const int od = depth / 4;
    const int oh = height / 4;
    const int ow = width / 4;

    TORCH_CHECK(od > 0 && oh > 0 && ow > 0,
                "Input is too small for two pooling layers");
    TORCH_CHECK(!training || int64_t(batch) * spatial > 1,
                "Batch normalization requires more than one value per channel");

    auto output = torch::empty({batch, channels, od, oh, ow}, x.options());
    auto affine = torch::empty({channels, 2}, x.options());
    auto stream = at::cuda::getCurrentCUDAStream();

    torch::Tensor partials;
    const float2* partial_ptr = nullptr;

    if (training) {
        partials = torch::empty(
            {channels, batch * chunks, 2}, x.options());
        partial_ptr =
            reinterpret_cast<const float2*>(partials.data_ptr<float>());

        partial_statistics<<<dim3(batch * chunks, channels), THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partials.data_ptr<float>()),
            channels, spatial, chunks, batch);
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    finalize_statistics<<<channels, THREADS, 0, stream>>>(
        partial_ptr,
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        reinterpret_cast<float2*>(affine.data_ptr<float>()),
        batch, spatial, chunks, training, float(momentum), float(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const int64_t elements = output.numel();
    const int blocks = int(std::min<int64_t>(
        (elements + THREADS - 1) / THREADS, 65535));

    if (blocks > 0) {
        normalized_pool4<<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<const float2*>(affine.data_ptr<float>()),
            bias.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, depth, height, width, od, oh, ow, elements);
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    return output;
}
"""

_bn_pool_extension = load_inline(
    name="fused_bn_double_avg_pool3d_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["bn_pool_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_channels, out_channels, kernel_size, stride, padding, bias_shape
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
        )
        self.batch_norm = nn.BatchNorm3d(out_channels)
        self.avg_pool1 = nn.AvgPool3d(kernel_size=2)
        self.avg_pool2 = nn.AvgPool3d(kernel_size=2)

    def forward(self, x):
        x = self.conv_transpose(x)
        bn = self.batch_norm

        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean.dtype != torch.float32
        ):
            return self.avg_pool2(self.avg_pool1(bn(x)))

        momentum = bn.momentum
        if bn.training:
            bn.num_batches_tracked.add_(1)
            if momentum is None:
                momentum = 1.0 / float(bn.num_batches_tracked.item())

        return _bn_pool_extension.bn_pool_cuda(
            x.contiguous(),
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            bn.training,
            0.0 if momentum is None else momentum,
            bn.eps,
        )
```