```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 split_gemm_cuda(torch::Tensor A, torch::Tensor B);
"""

_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 <cuda_fp16.h>
#include <cublas_v2.h>
#include <climits>

namespace {

constexpr float RESIDUAL_SCALE = 4096.0f;

__global__ void split_inputs(
    const float* __restrict__ a,
    const float* __restrict__ b,
    __half* __restrict__ ah,
    __half* __restrict__ al,
    __half* __restrict__ bh,
    __half* __restrict__ bl,
    int64_t na,
    int64_t nb)
{
    const int64_t total = na + nb;
    for (int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
                     + threadIdx.x;
         i < total;
         i += static_cast<int64_t>(blockDim.x) * gridDim.x) {
        const bool first = i < na;
        const int64_t j = first ? i : i - na;
        const float x = first ? a[j] : b[j];

        const __half hi = __float2half_rn(x);
        const __half lo = __float2half_rn(
            (x - __half2float(hi)) * RESIDUAL_SCALE);

        if (first) {
            ah[j] = hi;
            al[j] = lo;
        } else {
            bh[j] = hi;
            bl[j] = lo;
        }
    }
}

void check_blas(cublasStatus_t status) {
    TORCH_CHECK(status == CUBLAS_STATUS_SUCCESS,
                "cuBLAS failed with status ", static_cast<int>(status));
}

class MathModeGuard {
public:
    explicit MathModeGuard(cublasHandle_t handle) : handle_(handle) {
        check_blas(cublasGetMathMode(handle_, &previous_));
        check_blas(cublasSetMathMode(
            handle_,
            static_cast<cublasMath_t>(
                CUBLAS_DEFAULT_MATH |
                CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION)));
    }

    ~MathModeGuard() {
        cublasSetMathMode(handle_, previous_);
    }

private:
    cublasHandle_t handle_;
    cublasMath_t previous_;
};

} // namespace

torch::Tensor split_gemm_cuda(torch::Tensor A, torch::Tensor B) {
    TORCH_CHECK(A.is_cuda() && B.is_cuda(), "Inputs must be CUDA tensors");
    TORCH_CHECK(A.device() == B.device(), "Inputs must share a device");
    TORCH_CHECK(A.dim() == 2 && B.dim() == 2, "Inputs must be matrices");
    TORCH_CHECK(A.size(1) == B.size(1), "Incompatible inner dimensions");

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

    if (A.scalar_type() != at::kFloat ||
        B.scalar_type() != at::kFloat) {
        return at::matmul(A, B.transpose(0, 1));
    }

    const int64_t m64 = A.size(0);
    const int64_t n64 = B.size(0);
    const int64_t k64 = A.size(1);

    TORCH_CHECK(m64 <= INT_MAX && n64 <= INT_MAX && k64 <= INT_MAX,
                "Matrix dimensions exceed cuBLAS integer limits");

    const int M = static_cast<int>(m64);
    const int N = static_cast<int>(n64);
    const int K = static_cast<int>(k64);

    if (M == 0 || N == 0) {
        return torch::empty({M, N}, A.options());
    }
    if (K == 0) {
        return torch::zeros({M, N}, A.options());
    }

    auto a = A.contiguous();
    auto b = B.contiguous();

    auto half_options = A.options().dtype(at::kHalf);
    auto ah = torch::empty({M, K}, half_options);
    auto al = torch::empty({M, K}, half_options);
    auto bh = torch::empty({N, K}, half_options);
    auto bl = torch::empty({N, K}, half_options);
    auto out = torch::empty({M, N}, A.options());

    const int64_t total = a.numel() + b.numel();
    const int blocks = static_cast<int>(
        std::min<int64_t>((total + 255) / 256, 65535));

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

    split_inputs<<<blocks, 256, 0, stream>>>(
        a.data_ptr<float>(),
        b.data_ptr<float>(),
        reinterpret_cast<__half*>(ah.data_ptr<at::Half>()),
        reinterpret_cast<__half*>(al.data_ptr<at::Half>()),
        reinterpret_cast<__half*>(bh.data_ptr<at::Half>()),
        reinterpret_cast<__half*>(bl.data_ptr<at::Half>()),
        a.numel(),
        b.numel());
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    cublasHandle_t handle = at::cuda::getCurrentCUDABlasHandle();
    MathModeGuard math_guard(handle);

    const float zero = 0.0f;
    const float one = 1.0f;
    const float residual_alpha = 1.0f / RESIDUAL_SCALE;

    // Row-major C = A B^T is column-major C^T = B A^T.
    auto gemm = [&](const torch::Tensor& left,
                    const torch::Tensor& right,
                    float alpha,
                    float beta) {
        check_blas(cublasGemmEx(
            handle,
            CUBLAS_OP_T,
            CUBLAS_OP_N,
            N, M, K,
            &alpha,
            left.data_ptr<at::Half>(), CUDA_R_16F, K,
            right.data_ptr<at::Half>(), CUDA_R_16F, K,
            &beta,
            out.data_ptr<float>(), CUDA_R_32F, N,
            CUBLAS_COMPUTE_32F,
            CUBLAS_GEMM_DEFAULT_TENSOR_OP));
    };

    // Preserve FP32 accuracy using high components and both cross terms.
    gemm(bl, ah, residual_alpha, zero);
    gemm(bh, al, residual_alpha, one);
    gemm(bh, ah, one, one);

    return out;
}
"""

_extension = load_inline(
    name="split_tensorcore_matmul",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["split_gemm_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    extra_ldflags=["-lcublas"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self):
        super().__init__()

    def forward(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
        if not A.is_cuda or not B.is_cuda:
            return torch.matmul(A, B.T)
        return _extension.split_gemm_cuda(A, B)
```