```python
import torch
import torch.nn as nn
from torch.autograd.function import once_differentiable
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>

__device__ __forceinline__ float warp_sum(float v) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        v += __shfl_down_sync(0xffffffff, v, offset);
    return v;
}

template<bool Training>
__global__ void spatial_moments(
    const float* __restrict__ x,
    float* __restrict__ moments,
    int spatial,
    int nc,
    float scale)
{
    __shared__ float shared_sum[8];
    __shared__ float shared_square[8];

    const int row = blockIdx.x;
    const int tid = threadIdx.x;
    const int lane = tid & 31;
    const int warp = tid >> 5;
    const float* src = x + static_cast<int64_t>(row) * spatial;

    float sum = 0.0f;
    float square = 0.0f;

    for (int i = tid; i < spatial; i += blockDim.x) {
        const float v = src[i] * scale;
        sum += v;
        if (Training)
            square = fmaf(v, v, square);
    }

    sum = warp_sum(sum);
    if (Training)
        square = warp_sum(square);

    if (lane == 0) {
        shared_sum[warp] = sum;
        if (Training)
            shared_square[warp] = square;
    }
    __syncthreads();

    if (warp == 0) {
        sum = lane < 8 ? shared_sum[lane] : 0.0f;
        sum = warp_sum(sum);
        if (Training) {
            square = lane < 8 ? shared_square[lane] : 0.0f;
            square = warp_sum(square);
        }
        if (lane == 0) {
            moments[row] = sum / static_cast<float>(spatial);
            if (Training)
                moments[nc + row] = square / static_cast<float>(spatial);
        }
    }
}

template<bool Training>
__global__ void normalize_and_pool(
    const float* __restrict__ moments,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ stats,
    float* __restrict__ out,
    int batch,
    int channels,
    int spatial,
    float eps,
    float momentum)
{
    const int c = blockIdx.x;
    const int lane = threadIdx.x;
    const int nc = batch * channels;

    float mean;
    float variance;

    if (Training) {
        float sum = 0.0f;
        float square = 0.0f;
        for (int n = lane; n < batch; n += 32) {
            sum += moments[n * channels + c];
            square += moments[nc + n * channels + c];
        }

        sum = warp_sum(sum);
        square = warp_sum(square);
        mean = __shfl_sync(0xffffffff, sum, 0) / batch;
        const float second =
            __shfl_sync(0xffffffff, square, 0) / batch;
        variance = fmaxf(fmaf(-mean, mean, second), 0.0f);

        if (lane == 0) {
            const float count = static_cast<float>(
                static_cast<int64_t>(batch) * spatial);
            const float unbiased = variance * (count / (count - 1.0f));
            running_mean[c] =
                (1.0f - momentum) * running_mean[c] + momentum * mean;
            running_var[c] =
                (1.0f - momentum) * running_var[c] + momentum * unbiased;
        }
    } else {
        mean = running_mean[c];
        variance = running_var[c];
    }

    const float invstd = rsqrtf(variance + eps);
    const float gamma = weight[c];
    const float beta = bias[c];

    if (lane == 0) {
        stats[c] = mean;
        stats[channels + c] = invstd;
    }

    for (int n = lane; n < batch; n += 32) {
        const float sample_mean = moments[n * channels + c];
        out[n * channels + c] =
            ((sample_mean - 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() == torch::kFloat32,
                "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 = x.size(0);
    const int channels = x.size(1);
    const int spatial = x.size(2) * x.size(3) * x.size(4);
    const int nc = batch * channels;

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

    auto moments = torch::empty({2, batch, channels}, x.options());
    auto stats = torch::empty({2, channels}, x.options());
    auto out = torch::empty({batch, channels, 1, 1, 1}, x.options());
    const cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        spatial_moments<true><<<nc, 256, 0, stream>>>(
            x.data_ptr<float>(), moments.data_ptr<float>(),
            spatial, nc, static_cast<float>(scale));
        normalize_and_pool<true><<<channels, 32, 0, stream>>>(
            moments.data_ptr<float>(),
            weight.data_ptr<float>(), bias.data_ptr<float>(),
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            stats.data_ptr<float>(), out.data_ptr<float>(),
            batch, channels, spatial,
            static_cast<float>(eps), static_cast<float>(momentum));
    } else {
        spatial_moments<false><<<nc, 256, 0, stream>>>(
            x.data_ptr<float>(), moments.data_ptr<float>(),
            spatial, nc, static_cast<float>(scale));
        normalize_and_pool<false><<<channels, 32, 0, stream>>>(
            moments.data_ptr<float>(),
            weight.data_ptr<float>(), bias.data_ptr<float>(),
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            stats.data_ptr<float>(), out.data_ptr<float>(),
            batch, channels, spatial,
            static_cast<float>(eps), static_cast<float>(momentum));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return {out, moments, stats};
}
"""

_fused_ops = load_inline(
    name="fused_scale_bn_global_pool_3d",
    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
    ):
        out, moments, stats = _fused_ops.fused_bn_pool_cuda(
            x, weight, bias, running_mean, running_var,
            scale, eps, momentum, training
        )
        ctx.save_for_backward(x, weight, moments, stats)
        ctx.scale = scale
        ctx.training = training
        return out

    @staticmethod
    @once_differentiable
    def backward(ctx, grad_output):
        x, weight, moments, stats = ctx.saved_tensors
        n, c = x.shape[:2]
        spatial = x.shape[2] * x.shape[3] * x.shape[4]
        g = grad_output.reshape(n, c)
        mean, invstd = stats[0], stats[1]
        pooled_normalized = (moments[0] - mean) * invstd

        grad_weight = (
            (g * pooled_normalized).sum(dim=0)
            if ctx.needs_input_grad[1] else None
        )
        grad_bias = g.sum(dim=0) if ctx.needs_input_grad[2] else None
        grad_x = None

        if ctx.needs_input_grad[0]:
            factor = (
                ctx.scale * weight * invstd / spatial
            ).view(1, c, 1, 1, 1)

            if ctx.training:
                normalized = (
                    x * ctx.scale - mean.view(1, c, 1, 1, 1)
                ) * invstd.view(1, c, 1, 1, 1)
                correction = (g * pooled_normalized).mean(dim=0)
                grad_x = (
                    (g - g.mean(dim=0)).view(n, c, 1, 1, 1)
                    - normalized * correction.view(1, c, 1, 1, 1)
                ) * factor
            else:
                grad_x = (
                    g.view(n, c, 1, 1, 1) * factor
                ).expand_as(x).contiguous()

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

        if bn.training and x.numel() // x.shape[1] <= 1:
            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:
            bn.num_batches_tracked.add_(1)
            if bn.momentum is None:
                momentum = 1.0 / bn.num_batches_tracked.item()

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