```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 <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAStream.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <vector>

constexpr int THREADS = 256;

__device__ __forceinline__ float block_sum(float value, float* shared) {
    const int tid = threadIdx.x;
    shared[tid] = value;
    __syncthreads();

    for (int offset = THREADS / 2; offset > 0; offset >>= 1) {
        if (tid < offset)
            shared[tid] += shared[tid + offset];
        __syncthreads();
    }

    const float result = shared[0];
    __syncthreads();
    return result;
}

template<bool Training>
__global__ void spatial_stats_kernel(
    const float* __restrict__ x,
    float* __restrict__ means,
    float* __restrict__ variances,
    int64_t spatial)
{
    __shared__ float shared[THREADS];

    const int row = blockIdx.x;
    const float* input = x + static_cast<int64_t>(row) * spatial;

    // Centering the sum improves accuracy when activations have a large mean.
    const float anchor = input[0];
    float sum = 0.0f;
    for (int64_t i = threadIdx.x; i < spatial; i += THREADS)
        sum += input[i] - anchor;

    const float mean = anchor + block_sum(sum, shared) / float(spatial);

    if (threadIdx.x == 0)
        means[row] = mean;

    if (Training) {
        float squared_sum = 0.0f;
        for (int64_t i = threadIdx.x; i < spatial; i += THREADS) {
            const float delta = input[i] - mean;
            squared_sum = fmaf(delta, delta, squared_sum);
        }

        const float variance = block_sum(squared_sum, shared) / float(spatial);
        if (threadIdx.x == 0)
            variances[row] = variance;
    }
}

template<bool Training>
__global__ void normalize_pool_kernel(
    const float* __restrict__ means,
    const float* __restrict__ variances,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ output,
    float* __restrict__ saved_stats,
    int batch,
    int channels,
    int64_t spatial,
    float scale,
    float eps,
    float momentum)
{
    __shared__ float shared[THREADS];
    const int channel = blockIdx.x;

    float scaled_mean;
    float invstd;

    if (Training) {
        const float anchor = means[channel];
        float mean_sum = 0.0f;

        for (int n = threadIdx.x; n < batch; n += THREADS)
            mean_sum += means[n * channels + channel] - anchor;

        const float mean = anchor + block_sum(mean_sum, shared) / float(batch);

        float variance_sum = 0.0f;
        for (int n = threadIdx.x; n < batch; n += THREADS) {
            const int index = n * channels + channel;
            const float delta = means[index] - mean;
            variance_sum += variances[index] + delta * delta;
        }

        const float variance =
            block_sum(variance_sum, shared) / float(batch);
        const float scaled_variance = variance * scale * scale;

        scaled_mean = mean * scale;
        invstd = rsqrtf(scaled_variance + eps);

        if (threadIdx.x == 0) {
            const double count = double(batch) * double(spatial);
            const float unbiased_variance =
                scaled_variance * float(count / (count - 1.0));

            running_mean[channel] =
                (1.0f - momentum) * running_mean[channel]
                + momentum * scaled_mean;
            running_var[channel] =
                (1.0f - momentum) * running_var[channel]
                + momentum * unbiased_variance;
        }
    } else {
        scaled_mean = running_mean[channel];
        invstd = rsqrtf(running_var[channel] + eps);
    }

    if (threadIdx.x == 0) {
        saved_stats[channel] = scaled_mean;
        saved_stats[channels + channel] = invstd;
    }

    const float gamma = weight[channel];
    const float beta = bias[channel];

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        const int index = n * channels + channel;
        output[index] =
            ((scale * means[index] - scaled_mean) * invstd) * gamma + beta;
    }
}

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.scalar_type() == at::kFloat,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");

    const c10::cuda::CUDAGuard 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);

    TORCH_CHECK(batch > 0 && spatial > 0, "Input must be nonempty");
    TORCH_CHECK(!training || int64_t(batch) * spatial > 1,
                "Expected more than one value per channel when training");

    auto means = torch::empty({batch, channels}, x.options());
    auto variances = torch::empty({batch, channels}, x.options());
    auto stats = torch::empty({2, channels}, x.options());
    auto output = torch::empty({batch, channels, 1, 1, 1}, x.options());

    const auto stream = c10::cuda::getCurrentCUDAStream();

    if (training) {
        spatial_stats_kernel<true><<<batch * channels, THREADS, 0, stream>>>(
            x.data_ptr<float>(), means.data_ptr<float>(),
            variances.data_ptr<float>(), spatial);

        normalize_pool_kernel<true><<<channels, THREADS, 0, stream>>>(
            means.data_ptr<float>(), variances.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>(), stats.data_ptr<float>(),
            batch, channels, spatial,
            float(scale), float(eps), float(momentum));
    } else {
        spatial_stats_kernel<false><<<batch * channels, THREADS, 0, stream>>>(
            x.data_ptr<float>(), means.data_ptr<float>(),
            variances.data_ptr<float>(), spatial);

        normalize_pool_kernel<false><<<channels, THREADS, 0, stream>>>(
            means.data_ptr<float>(), variances.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>(), stats.data_ptr<float>(),
            batch, channels, spatial,
            float(scale), float(eps), float(momentum));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return {output, means, stats};
}
"""

_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 _FusedBNPool(torch.autograd.Function):
    @staticmethod
    def forward(
        ctx, x, weight, bias, running_mean, running_var,
        scale, eps, momentum, training
    ):
        output, means, stats = _fused_ops.fused_bn_pool_cuda(
            x, weight, bias, running_mean, running_var,
            scale, eps, momentum, training
        )
        ctx.save_for_backward(x, weight, means, stats)
        ctx.scale = scale
        ctx.training = training
        ctx.spatial = x.shape[2] * x.shape[3] * x.shape[4]
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, means, stats = ctx.saved_tensors
        grad = grad_output.reshape_as(means)
        mean, invstd = stats[0], stats[1]
        centered_means = ctx.scale * means - mean.unsqueeze(0)

        grad_x = grad_weight = grad_bias = None

        if ctx.needs_input_grad[0]:
            factor = ctx.scale * weight * invstd / ctx.spatial

            if ctx.training:
                centered_grad = grad - grad.mean(dim=0, keepdim=True)
                covariance = (grad * centered_means).mean(dim=0)
                correction = covariance * invstd.square()

                grad_x = (
                    centered_grad[:, :, None, None, None]
                    - (ctx.scale * x - mean[None, :, None, None, None])
                    * correction[None, :, None, None, None]
                ) * factor[None, :, None, None, None]
            else:
                grad_x = (
                    grad * factor.unsqueeze(0)
                )[:, :, None, None, None].expand_as(x)

        if ctx.needs_input_grad[1]:
            grad_weight = (
                grad * centered_means * invstd.unsqueeze(0)
            ).sum(dim=0)

        if ctx.needs_input_grad[2]:
            grad_bias = grad.sum(dim=0)

        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 bn.weight is None
            or bn.weight.dtype != torch.float32
            or bn.running_mean is None
            or bn.running_var is None
            or x.numel() == 0
        ):
            return self.global_avg_pool(bn(x * self.scale_factor))

        momentum = bn.momentum if bn.momentum is not None else 0.0

        if bn.training:
            if x.numel() // x.shape[1] <= 1:
                return self.global_avg_pool(bn(x * self.scale_factor))

            with torch.no_grad():
                bn.num_batches_tracked.add_(1)

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

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