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


_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 <math_constants.h>

constexpr int THREADS = 256;
constexpr int CHUNK = 4096;

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

template<bool VECTORIZED>
__global__ void partial_moments(
    const float* __restrict__ x,
    const float* __restrict__ conv_bias,
    float2* __restrict__ partial,
    int spatial,
    int groups,
    int group_size,
    int chunks
) {
    const int row = blockIdx.y;
    const int tid = threadIdx.x;
    const int start = blockIdx.x * CHUNK;
    const int first_channel = (row % groups) * (64 / groups);
    const float* base = x + static_cast<int64_t>(row) * group_size;

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

    if constexpr (VECTORIZED) {
        #pragma unroll
        for (int j = 0; j < CHUNK / (4 * THREADS); ++j) {
            const int index = start + 4 * (tid + j * THREADS);
            if (index < group_size) {
                float4 v = *reinterpret_cast<const float4*>(base + index);
                const float b = conv_bias[first_channel + index / spatial];
                v.x += b;
                v.y += b;
                v.z += b;
                v.w += b;
                sum += v.x;
                sum += v.y;
                sum += v.z;
                sum += v.w;
                square_sum = fmaf(v.x, v.x, square_sum);
                square_sum = fmaf(v.y, v.y, square_sum);
                square_sum = fmaf(v.z, v.z, square_sum);
                square_sum = fmaf(v.w, v.w, square_sum);
            }
        }
    } else {
        #pragma unroll
        for (int j = 0; j < CHUNK / THREADS; ++j) {
            const int index = start + tid + j * THREADS;
            if (index < group_size) {
                const float v =
                    base[index] + conv_bias[first_channel + index / spatial];
                sum += v;
                square_sum = fmaf(v, v, square_sum);
            }
        }
    }

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

    __shared__ float sums[8];
    __shared__ float squares[8];
    const int lane = tid & 31;
    const int warp = tid >> 5;

    if (lane == 0) {
        sums[warp] = sum;
        squares[warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = warp_sum(lane < 8 ? sums[lane] : 0.0f);
        square_sum = warp_sum(lane < 8 ? squares[lane] : 0.0f);
        if (lane == 0) {
            partial[static_cast<int64_t>(row) * chunks + blockIdx.x] =
                make_float2(sum, square_sum);
        }
    }
}

__global__ void finish_moments(
    const float2* __restrict__ partial,
    float2* __restrict__ moments,
    int group_size,
    int chunks,
    float eps
) {
    const int row = blockIdx.x;
    const int lane = threadIdx.x;
    float sum = 0.0f;
    float square_sum = 0.0f;

    for (int i = lane; i < chunks; i += 32) {
        const float2 v = partial[static_cast<int64_t>(row) * chunks + i];
        sum += v.x;
        square_sum += v.y;
    }

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

    if (lane == 0) {
        const float mean = sum / static_cast<float>(group_size);
        const float variance = fmaxf(
            fmaf(-mean, mean, square_sum / static_cast<float>(group_size)),
            0.0f
        );
        moments[row] = make_float2(mean, rsqrtf(variance + eps));
    }
}

template<int FIXED_GROUPS>
__global__ void fused_output(
    const float* __restrict__ x,
    const float* __restrict__ conv_bias,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float2* __restrict__ moments,
    float* __restrict__ output,
    int spatial,
    int runtime_groups
) {
    const int position = blockIdx.x * blockDim.x + threadIdx.x;
    const int batch = blockIdx.y;
    if (position >= spatial)
        return;

    const int groups = FIXED_GROUPS ? FIXED_GROUPS : runtime_groups;
    const int channels_per_group = 64 / groups;
    const float* input =
        x + static_cast<int64_t>(batch) * 64 * spatial + position;

    float values[64];
    float maximum = -CUDART_INF_F;

    #pragma unroll
    for (int channel = 0; channel < 64; ++channel) {
        const float original =
            input[static_cast<int64_t>(channel) * spatial] + conv_bias[channel];
        const float2 stats =
            moments[batch * groups + channel / channels_per_group];

        const float normalized = fmaf(
            (original - stats.x) * stats.y,
            weight[channel],
            bias[channel]
        );
        const float activated = tanhf(normalized);
        const float hard_swish =
            (activated * (activated + 3.0f)) * (1.0f / 6.0f);
        const float residual = original + hard_swish;

        values[channel] = residual;
        maximum = fmaxf(maximum, residual);
    }

    float exponential_sum = 0.0f;
    #pragma unroll
    for (int channel = 0; channel < 64; ++channel)
        exponential_sum += __expf(values[channel] - maximum);

    output[static_cast<int64_t>(batch) * spatial + position] =
        maximum + __logf(exponential_sum);
}

torch::Tensor fused_post_conv_cuda(
    torch::Tensor x,
    torch::Tensor conv_bias,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps
) {
    TORCH_CHECK(
        x.is_cuda() && x.scalar_type() == torch::kFloat32 &&
        x.dim() == 4 && x.size(1) == 64 && x.is_contiguous(),
        "Expected contiguous CUDA float32 NCHW input with 64 channels"
    );
    TORCH_CHECK(groups > 0 && 64 % groups == 0, "Invalid group count");
    for (const auto& parameter : {conv_bias, weight, bias}) {
        TORCH_CHECK(
            parameter.device() == x.device() &&
            parameter.scalar_type() == torch::kFloat32 &&
            parameter.is_contiguous() && parameter.numel() == 64,
            "Expected 64 contiguous float32 parameters on the input device"
        );
    }

    const c10::cuda::CUDAGuard guard(x.device());
    const auto stream =
        c10::cuda::getCurrentCUDAStream(x.get_device()).stream();

    const int batch = static_cast<int>(x.size(0));
    const int spatial = static_cast<int>(x.size(2) * x.size(3));
    const int group_count = static_cast<int>(groups);
    auto output = torch::empty(
        {x.size(0), 1, x.size(2), x.size(3)}, x.options()
    );
    if (batch == 0 || spatial == 0)
        return output;

    const int group_size = (64 / group_count) * spatial;
    const int chunks = (group_size + CHUNK - 1) / CHUNK;
    auto partial = torch::empty(
        {static_cast<int64_t>(batch) * group_count, chunks, 2}, x.options()
    );
    auto moments = torch::empty(
        {static_cast<int64_t>(batch) * group_count, 2}, x.options()
    );

    const dim3 stats_grid(chunks, batch * group_count);
    if (spatial % 4 == 0) {
        partial_moments<true><<<stats_grid, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            conv_bias.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, group_count, group_size, chunks
        );
    } else {
        partial_moments<false><<<stats_grid, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            conv_bias.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, group_count, group_size, chunks
        );
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_moments<<<batch * group_count, 32, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        reinterpret_cast<float2*>(moments.data_ptr<float>()),
        group_size, chunks, static_cast<float>(eps)
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const dim3 output_grid((spatial + 127) / 128, batch);
    if (group_count == 16) {
        fused_output<16><<<output_grid, 128, 0, stream>>>(
            x.data_ptr<float>(),
            conv_bias.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            reinterpret_cast<const float2*>(moments.data_ptr<float>()),
            output.data_ptr<float>(),
            spatial, group_count
        );
    } else {
        fused_output<0><<<output_grid, 128, 0, stream>>>(
            x.data_ptr<float>(),
            conv_bias.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            reinterpret_cast<const float2*>(moments.data_ptr<float>()),
            output.data_ptr<float>(),
            spatial, group_count
        );
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

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

torch::Tensor fused_post_conv_cuda(
    torch::Tensor x,
    torch::Tensor conv_bias,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps
);
"""

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


class ModelNew(nn.Module):
    def __init__(
        self, in_channels, out_channels, kernel_size, groups, eps=1e-5
    ):
        super().__init__()
        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size)
        self.group_norm = nn.GroupNorm(groups, out_channels, eps=eps)
        self.tanh = nn.Tanh()
        self.hard_swish = nn.Hardswish()

    def forward(self, x):
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and self.conv.out_channels == 64
            and self.conv.weight.dtype == torch.float32
            and self.conv.weight.is_contiguous()
            and self.conv.bias is not None
            and self.conv.bias.dtype == torch.float32
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
            and not torch.is_autocast_enabled()
        ):
            x_conv = F.conv2d(
                x,
                self.conv.weight,
                None,
                self.conv.stride,
                self.conv.padding,
                self.conv.dilation,
                self.conv.groups,
            )
            return _fused_ops.fused_post_conv_cuda(
                x_conv,
                self.conv.bias,
                self.group_norm.weight,
                self.group_norm.bias,
                self.group_norm.num_groups,
                self.group_norm.eps,
            )

        x_conv = self.conv(x)
        x_norm = self.group_norm(x_conv)
        x_activated = self.hard_swish(self.tanh(x_norm))
        return torch.logsumexp(x_conv + x_activated, dim=1, keepdim=True)
```