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

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

constexpr int THREADS = 256;
constexpr int TILE = 2048;

__device__ __forceinline__ double2 warp_sum(double2 value) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        value.x += __shfl_down_sync(0xffffffff, value.x, offset);
        value.y += __shfl_down_sync(0xffffffff, value.y, offset);
    }
    return value;
}

__device__ __forceinline__ double2 block_sum(
    double2 value, double2* shared) {
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    value = warp_sum(value);
    if (lane == 0)
        shared[warp] = value;
    __syncthreads();

    value = threadIdx.x < 8
        ? shared[threadIdx.x]
        : make_double2(0.0, 0.0);

    if (warp == 0)
        value = warp_sum(value);

    return value;
}

template<bool TRAINING>
__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    double2* __restrict__ partial,
    int64_t spatial,
    int parts,
    float scale) {
    __shared__ double2 scratch[8];

    const int64_t block = blockIdx.x;
    const int64_t row = block / parts;
    const int part = block % parts;
    const int64_t start = static_cast<int64_t>(part) * TILE;
    const int64_t base = row * spatial;

    double2 moments = make_double2(0.0, 0.0);

    #pragma unroll
    for (int i = 0; i < TILE / THREADS; ++i) {
        const int64_t s = start + threadIdx.x + i * THREADS;
        if (s < spatial) {
            // Preserve the float32 rounding of the original scaling operation.
            const float scaled = x[base + s] * scale;
            const double value = static_cast<double>(scaled);
            moments.x += value;
            if (TRAINING)
                moments.y += value * value;
        }
    }

    moments = block_sum(moments, scratch);
    if (threadIdx.x == 0)
        partial[block] = moments;
}

template<bool TRAINING>
__global__ void normalize_pool_kernel(
    const double2* __restrict__ partial,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ output,
    int batch,
    int channels,
    int parts,
    int64_t spatial,
    double eps,
    double momentum) {
    __shared__ double2 scratch[8];
    __shared__ double mean;
    __shared__ double inv_std;

    const int channel = blockIdx.x;

    if (TRAINING) {
        double2 total = make_double2(0.0, 0.0);
        const int64_t count = static_cast<int64_t>(batch) * parts;

        for (int64_t i = threadIdx.x; i < count; i += THREADS) {
            const int64_t n = i / parts;
            const int p = i % parts;
            const double2 value =
                partial[(n * channels + channel) * parts + p];
            total.x += value.x;
            total.y += value.y;
        }

        total = block_sum(total, scratch);

        if (threadIdx.x == 0) {
            const double elements =
                static_cast<double>(batch) * static_cast<double>(spatial);
            mean = total.x / elements;
            const double variance =
                fmax(total.y / elements - mean * mean, 0.0);
            inv_std = 1.0 / sqrt(variance + eps);

            const double unbiased_variance =
                variance * elements / (elements - 1.0);

            running_mean[channel] = static_cast<float>(
                (1.0 - momentum) * running_mean[channel] +
                momentum * mean);
            running_var[channel] = static_cast<float>(
                (1.0 - momentum) * running_var[channel] +
                momentum * unbiased_variance);
        }
    } else {
        if (threadIdx.x == 0) {
            mean = static_cast<double>(running_mean[channel]);
            inv_std = 1.0 / sqrt(
                static_cast<double>(running_var[channel]) + eps);
        }
    }

    __syncthreads();

    const double gamma = static_cast<double>(weight[channel]);
    const double beta = static_cast<double>(bias[channel]);

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        double sample_sum = 0.0;
        const int64_t offset =
            (static_cast<int64_t>(n) * channels + channel) * parts;

        for (int p = 0; p < parts; ++p)
            sample_sum += partial[offset + p].x;

        const double sample_mean =
            sample_sum / static_cast<double>(spatial);

        output[static_cast<int64_t>(n) * channels + channel] =
            static_cast<float>(
                (sample_mean - mean) * inv_std * gamma + beta);
    }
}

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

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

    const int batch = static_cast<int>(x.size(0));
    const int channels = static_cast<int>(x.size(1));
    const int64_t spatial = x.size(2) * x.size(3) * x.size(4);
    const int parts = static_cast<int>((spatial + TILE - 1) / TILE);

    TORCH_CHECK(
        !training || static_cast<int64_t>(batch) * spatial > 1,
        "Batch normalization requires more than one value per channel");

    auto output = torch::empty(
        {batch, channels, 1, 1, 1}, x.options());
    if (batch == 0)
        return output;

    auto partial = torch::empty(
        {batch, channels, parts, 2},
        x.options().dtype(torch::kFloat64));

    auto stream = at::cuda::getCurrentCUDAStream();
    const int64_t blocks =
        static_cast<int64_t>(batch) * channels * parts;

    if (training) {
        partial_moments_kernel<true><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<double2*>(partial.data_ptr<double>()),
            spatial, parts, static_cast<float>(scale));

        normalize_pool_kernel<true><<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const double2*>(partial.data_ptr<double>()),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            output.data_ptr<float>(),
            batch, channels, parts, spatial, eps, momentum);
    } else {
        partial_moments_kernel<false><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<double2*>(partial.data_ptr<double>()),
            spatial, parts, static_cast<float>(scale));

        normalize_pool_kernel<false><<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const double2*>(partial.data_ptr<double>()),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            output.data_ptr<float>(),
            batch, channels, parts, spatial, eps, momentum);
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="convtranspose3d_fused_bn_global_pool",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_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,
        scale_factor,
        eps=1e-5,
        momentum=0.1,
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size
        )
        self.scale_factor = scale_factor
        self.batch_norm = nn.BatchNorm3d(
            out_channels, eps=eps, momentum=momentum
        )
        self.global_avg_pool = nn.AdaptiveAvgPool3d((1, 1, 1))

    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
        ):
            return self.global_avg_pool(bn(x * self.scale_factor))

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

        return _fused_ops.fused_bn_pool_cuda(
            x.contiguous(),
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            float(self.scale_factor),
            float(bn.eps),
            float(momentum) if momentum is not None else 0.0,
            bn.training,
        )
```