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

_CPP = 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 = 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>
#include <cstdint>

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

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

__device__ __forceinline__ void block_sum_pair(
    float& sum, float& square_sum, float* shared) {
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    sum = warp_sum(sum);
    square_sum = warp_sum(square_sum);

    if (lane == 0) {
        shared[warp] = sum;
        shared[8 + warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = lane < 8 ? shared[lane] : 0.0f;
        square_sum = lane < 8 ? shared[8 + lane] : 0.0f;
        sum = warp_sum(sum);
        square_sum = warp_sum(square_sum);
    }
}

template<bool TRAINING>
__global__ void partial_moments(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int64_t spatial,
    int tiles,
    int64_t jobs,
    float scale) {
    __shared__ float shared[16];

    for (int64_t job = blockIdx.x; job < jobs; job += gridDim.x) {
        const int tile = job % tiles;
        const int64_t row = job / tiles;
        const int64_t begin = static_cast<int64_t>(tile) * TILE;
        const int64_t end = min(begin + TILE, spatial);

        float sum = 0.0f;
        float square_sum = 0.0f;

        for (int64_t i = begin + threadIdx.x; i < end; i += THREADS) {
            const float value = x[row * spatial + i] * scale;
            sum += value;
            if (TRAINING)
                square_sum = fmaf(value, value, square_sum);
        }

        block_sum_pair(sum, square_sum, shared);
        if (threadIdx.x == 0)
            partial[job] = make_float2(sum, square_sum);
        __syncthreads();
    }
}

template<bool TRAINING>
__global__ void finish_bn_pool(
    const float2* __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,
    int64_t spatial,
    int tiles,
    float eps,
    float momentum) {
    const int channel = blockIdx.x;
    __shared__ float shared[16];
    __shared__ float mean;
    __shared__ float inverse_std;

    if (TRAINING) {
        float sum = 0.0f;
        float square_sum = 0.0f;
        const int64_t count = static_cast<int64_t>(batch) * tiles;

        for (int64_t i = threadIdx.x; i < count; i += THREADS) {
            const int n = i / tiles;
            const int tile = i % tiles;
            const float2 p =
                partial[(static_cast<int64_t>(n) * channels + channel)
                        * tiles + tile];
            sum += p.x;
            square_sum += p.y;
        }

        block_sum_pair(sum, square_sum, shared);

        if (threadIdx.x == 0) {
            const float elements =
                static_cast<float>(static_cast<int64_t>(batch) * spatial);
            const float m = sum / elements;
            const float variance =
                fmaxf(square_sum / elements - m * m, 0.0f);

            mean = m;
            inverse_std = rsqrtf(variance + eps);

            running_mean[channel] =
                (1.0f - momentum) * running_mean[channel] + momentum * m;
            running_var[channel] =
                (1.0f - momentum) * running_var[channel]
                + momentum * variance * (elements / (elements - 1.0f));
        }
    } else if (threadIdx.x == 0) {
        mean = running_mean[channel];
        inverse_std = rsqrtf(running_var[channel] + eps);
    }
    __syncthreads();

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

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        float sample_sum = 0.0f;
        const int64_t base =
            (static_cast<int64_t>(n) * channels + channel) * tiles;

        for (int tile = 0; tile < tiles; ++tile)
            sample_sum += partial[base + tile].x;

        const float sample_mean = sample_sum / static_cast<float>(spatial);
        output[static_cast<int64_t>(n) * channels + channel] =
            ((sample_mean - mean) * inverse_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) {
    const c10::cuda::CUDAGuard guard(x.device());
    x = x.contiguous();

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

    auto partial = torch::empty({jobs, 2}, x.options());
    auto output = torch::empty({batch, channels, 1, 1, 1}, x.options());

    const int blocks = static_cast<int>(std::min<int64_t>(jobs, 65535));
    auto stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        partial_moments<true><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, tiles, jobs, static_cast<float>(scale));

        finish_bn_pool<true><<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const float2*>(partial.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>(), batch, channels, spatial, tiles,
            static_cast<float>(eps), static_cast<float>(momentum));
    } else {
        partial_moments<false><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, tiles, jobs, static_cast<float>(scale));

        finish_bn_pool<false><<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const float2*>(partial.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>(), batch, channels, spatial, tiles,
            static_cast<float>(eps), static_cast<float>(momentum));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_EXT = load_inline(
    name="fused_scale_bn_pool3d_ext",
    cpp_sources=_CPP,
    cuda_sources=_CUDA,
    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

        # Preserve autograd and support other devices and data types.
        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or bn.weight.dtype != torch.float32
            or x.numel() == 0
            or (torch.is_grad_enabled() and (
                x.requires_grad
                or bn.weight.requires_grad
                or bn.bias.requires_grad
            ))
        ):
            return self.global_avg_pool(bn(x * self.scale_factor))

        if bn.training and x.numel() // x.size(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 _EXT.fused_bn_pool_cuda(
            x,
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            float(self.scale_factor),
            bn.eps,
            momentum,
            bn.training,
        )
```