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


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

std::vector<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_runtime.h>
#include <vector>

constexpr int THREADS = 256;
constexpr int TILE = 4096;

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

__device__ __forceinline__ Moments merge_moments(Moments a, Moments b) {
    if (b.count == 0) return a;
    if (a.count == 0) return b;

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

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

__device__ __forceinline__ Moments warp_reduce(Moments value) {
    int lane = threadIdx.x & 31;
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Moments other;
        other.mean = __shfl_down_sync(0xffffffff, value.mean, offset);
        other.m2 = __shfl_down_sync(0xffffffff, value.m2, offset);
        other.count = __shfl_down_sync(0xffffffff, value.count, offset);
        if (lane + offset < 32)
            value = merge_moments(value, other);
    }
    return value;
}

// Only thread zero receives the complete result.
__device__ __forceinline__ Moments block_reduce(Moments value) {
    __shared__ float means[8];
    __shared__ float m2s[8];
    __shared__ int counts[8];

    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    value = warp_reduce(value);

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

    Moments 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_moments_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_mean,
    float* __restrict__ partial_m2,
    int spatial,
    int tiles,
    float scale)
{
    int tile = blockIdx.x % tiles;
    int nc = blockIdx.x / tiles;
    int begin = tile * TILE;
    int end = min(begin + TILE, spatial);
    int64_t base = int64_t(nc) * spatial;

    Moments value = {0.0f, 0.0f, 0};

    for (int i = begin + threadIdx.x; i < end; i += THREADS) {
        float v = x[base + i] * scale;
        ++value.count;
        float delta = v - value.mean;
        value.mean += delta / float(value.count);
        value.m2 += delta * (v - value.mean);
    }

    value = block_reduce(value);
    if (threadIdx.x == 0) {
        partial_mean[blockIdx.x] = value.mean;
        partial_m2[blockIdx.x] = value.m2;
    }
}

__global__ void finish_bn_pool_kernel(
    const float* __restrict__ partial_mean,
    const float* __restrict__ partial_m2,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ output,
    float* __restrict__ saved_mean,
    float* __restrict__ saved_invstd,
    int batch,
    int channels,
    int spatial,
    int tiles,
    float eps,
    float momentum,
    bool training)
{
    int c = blockIdx.x;
    __shared__ float channel_mean;
    __shared__ float channel_invstd;

    if (training) {
        Moments value = {0.0f, 0.0f, 0};
        for (int j = threadIdx.x; j < batch * tiles; j += THREADS) {
            int n = j / tiles;
            int tile = j % tiles;
            int index = (n * channels + c) * tiles + tile;
            Moments part;
            part.mean = partial_mean[index];
            part.m2 = partial_m2[index];
            part.count = min(TILE, spatial - tile * TILE);
            value = merge_moments(value, part);
        }

        value = block_reduce(value);
        if (threadIdx.x == 0) {
            channel_mean = value.mean;
            channel_invstd = rsqrtf(value.m2 / float(value.count) + eps);

            running_mean[c] =
                (1.0f - momentum) * running_mean[c] + momentum * value.mean;
            running_var[c] =
                (1.0f - momentum) * running_var[c]
                + momentum * (value.m2 / float(value.count - 1));
        }
    } else if (threadIdx.x == 0) {
        channel_mean = running_mean[c];
        channel_invstd = rsqrtf(running_var[c] + eps);
    }

    __syncthreads();

    if (threadIdx.x == 0) {
        saved_mean[c] = channel_mean;
        saved_invstd[c] = channel_invstd;
    }

    // Pooling the normalized tensor only requires each sample's mean.
    for (int n = threadIdx.x; n < batch; n += THREADS) {
        float centered_sum = 0.0f;
        int base = (n * channels + c) * tiles;

        for (int tile = 0; tile < tiles; ++tile) {
            int count = min(TILE, spatial - tile * TILE);
            centered_sum +=
                (partial_mean[base + tile] - channel_mean) * float(count);
        }

        float normalized_mean =
            (centered_sum / float(spatial)) * channel_invstd;
        output[n * channels + c] =
            normalized_mean * weight[c] + bias[c];
    }
}

std::vector<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.is_contiguous(),
                "Expected contiguous CUDA input");
    TORCH_CHECK(x.scalar_type() == torch::kFloat32 && x.dim() == 5,
                "Expected five-dimensional float32 input");

    c10::cuda::CUDAGuard guard(x.device());

    int batch = x.size(0);
    int channels = x.size(1);
    int spatial = x.size(2) * x.size(3) * x.size(4);
    int tiles = (spatial + TILE - 1) / TILE;

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

    auto partial_mean = torch::empty({batch, channels, tiles}, x.options());
    auto partial_m2 = torch::empty_like(partial_mean);
    auto output = torch::empty({batch, channels, 1, 1, 1}, x.options());
    auto mean = torch::empty({channels}, x.options());
    auto invstd = torch::empty_like(mean);

    auto stream = at::cuda::getCurrentCUDAStream();

    partial_moments_kernel<<<batch * channels * tiles, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        partial_mean.data_ptr<float>(),
        partial_m2.data_ptr<float>(),
        spatial, tiles, float(scale));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_bn_pool_kernel<<<channels, THREADS, 0, stream>>>(
        partial_mean.data_ptr<float>(),
        partial_m2.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        output.data_ptr<float>(),
        mean.data_ptr<float>(),
        invstd.data_ptr<float>(),
        batch, channels, spatial, tiles,
        float(eps), float(momentum), training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return {output, mean, invstd};
}
"""

_fused_ops = load_inline(
    name="fused_scaled_bn3d_global_pool_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_bn_pool_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class _FusedBatchNormPool(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, running_mean, running_var,
                scale, eps, momentum, training):
        output, mean, invstd = _fused_ops.fused_bn_pool_cuda(
            x, weight, bias, running_mean, running_var,
            scale, eps, momentum, training,
        )
        ctx.save_for_backward(x, weight, mean, invstd)
        ctx.scale = scale
        ctx.training = training
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, mean, invstd = ctx.saved_tensors
        shape = (1, -1, 1, 1, 1)
        spatial = x.shape[2] * x.shape[3] * x.shape[4]
        batch = x.shape[0]

        invstd = invstd.view(shape)
        normalized = (x * ctx.scale - mean.view(shape)) * invstd
        sample_mean = normalized.mean(dim=(2, 3, 4), keepdim=True)

        grad_weight = (grad_output * sample_mean).sum(
            dim=(0, 2, 3, 4)
        )
        grad_bias = grad_output.sum(dim=(0, 2, 3, 4))

        if ctx.training:
            centered_grad = (
                grad_output - grad_output.mean(dim=0, keepdim=True)
            ) / spatial
            grad_x = (
                centered_grad
                - normalized * grad_weight.view(shape) / (batch * spatial)
            )
        else:
            grad_x = (grad_output / spatial).expand_as(x)

        grad_x = grad_x * (ctx.scale * weight.view(shape) * invstd)
        return (
            grad_x, grad_weight, grad_bias,
            None, None, None, None, None, None,
        )


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 x.numel() == 0
            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.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()
        if momentum is None:
            momentum = 0.0

        return _FusedBatchNormPool.apply(
            x.contiguous(),
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            float(self.scale_factor),
            float(bn.eps),
            float(momentum),
            bn.training,
        )
```