```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 triplet_margin_cuda(
    torch::Tensor anchor,
    torch::Tensor positive,
    torch::Tensor negative,
    double margin);
"""

_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 <cstdint>

__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 triplet_rows(
    const float* __restrict__ anchor,
    const float* __restrict__ positive,
    const float* __restrict__ negative,
    float* __restrict__ losses,
    int64_t width,
    float margin)
{
    const int tid = threadIdx.x;
    const int64_t base = static_cast<int64_t>(blockIdx.x) * width;
    float p0 = 0.f, p1 = 0.f, p2 = 0.f, p3 = 0.f;
    float n0 = 0.f, n1 = 0.f, n2 = 0.f, n3 = 0.f;

    if (Vectorized) {
        const float4* a = reinterpret_cast<const float4*>(anchor + base);
        const float4* p = reinterpret_cast<const float4*>(positive + base);
        const float4* n = reinterpret_cast<const float4*>(negative + base);

        for (int64_t i = tid; i < width / 4; i += blockDim.x) {
            const float4 av = a[i];
            const float4 pv = p[i];
            const float4 nv = n[i];

            const float dp0 = (av.x - pv.x) + 1.e-6f;
            const float dp1 = (av.y - pv.y) + 1.e-6f;
            const float dp2 = (av.z - pv.z) + 1.e-6f;
            const float dp3 = (av.w - pv.w) + 1.e-6f;
            const float dn0 = (av.x - nv.x) + 1.e-6f;
            const float dn1 = (av.y - nv.y) + 1.e-6f;
            const float dn2 = (av.z - nv.z) + 1.e-6f;
            const float dn3 = (av.w - nv.w) + 1.e-6f;

            p0 += dp0 * dp0;
            p1 += dp1 * dp1;
            p2 += dp2 * dp2;
            p3 += dp3 * dp3;
            n0 += dn0 * dn0;
            n1 += dn1 * dn1;
            n2 += dn2 * dn2;
            n3 += dn3 * dn3;
        }
    } else {
        for (int64_t i = tid; i < width; i += blockDim.x) {
            const float av = anchor[base + i];
            const float dp = (av - positive[base + i]) + 1.e-6f;
            const float dn = (av - negative[base + i]) + 1.e-6f;
            p0 += dp * dp;
            n0 += dn * dn;
        }
    }

    float ps = warp_sum((p0 + p1) + (p2 + p3));
    float ns = warp_sum((n0 + n1) + (n2 + n3));

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

    if (lane == 0) {
        partial_p[warp] = ps;
        partial_n[warp] = ns;
    }
    __syncthreads();

    if (warp == 0) {
        ps = warp_sum(lane < 8 ? partial_p[lane] : 0.f);
        ns = warp_sum(lane < 8 ? partial_n[lane] : 0.f);
        if (lane == 0) {
            const float loss = (margin + sqrtf(ps)) - sqrtf(ns);
            losses[blockIdx.x] = loss < 0.f ? 0.f : loss;
        }
    }
}

__global__ void mean_loss(
    const float* __restrict__ losses,
    float* __restrict__ result,
    int64_t rows)
{
    float sum = 0.f;
    for (int64_t i = threadIdx.x; i < rows; i += blockDim.x)
        sum += losses[i];

    sum = warp_sum(sum);
    __shared__ float partial[8];

    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;
    if (lane == 0)
        partial[warp] = sum;
    __syncthreads();

    if (warp == 0) {
        sum = warp_sum(lane < 8 ? partial[lane] : 0.f);
        if (lane == 0)
            result[0] = sum / static_cast<float>(rows);
    }
}

torch::Tensor triplet_margin_cuda(
    torch::Tensor anchor,
    torch::Tensor positive,
    torch::Tensor negative,
    double margin)
{
    TORCH_CHECK(anchor.is_cuda() && positive.is_cuda() && negative.is_cuda(),
                "Inputs must be CUDA tensors");
    TORCH_CHECK(anchor.device() == positive.device() &&
                anchor.device() == negative.device(),
                "Inputs must be on the same device");
    TORCH_CHECK(anchor.scalar_type() == at::kFloat &&
                positive.scalar_type() == at::kFloat &&
                negative.scalar_type() == at::kFloat,
                "Inputs must have dtype float32");
    TORCH_CHECK(anchor.dim() == 2 &&
                anchor.sizes() == positive.sizes() &&
                anchor.sizes() == negative.sizes(),
                "Inputs must have identical [batch, features] shapes");

    const c10::cuda::CUDAGuard guard(anchor.device());
    anchor = anchor.contiguous();
    positive = positive.contiguous();
    negative = negative.contiguous();

    const int64_t rows = anchor.size(0);
    const int64_t width = anchor.size(1);
    TORCH_CHECK(rows <= 2147483647LL, "Batch size is too large");

    auto losses = torch::empty({rows}, anchor.options());
    auto result = torch::empty({}, anchor.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    if (rows > 0) {
        if (width % 4 == 0) {
            triplet_rows<true><<<static_cast<unsigned int>(rows), 256, 0, stream>>>(
                anchor.data_ptr<float>(), positive.data_ptr<float>(),
                negative.data_ptr<float>(), losses.data_ptr<float>(),
                width, static_cast<float>(margin));
        } else {
            triplet_rows<false><<<static_cast<unsigned int>(rows), 256, 0, stream>>>(
                anchor.data_ptr<float>(), positive.data_ptr<float>(),
                negative.data_ptr<float>(), losses.data_ptr<float>(),
                width, static_cast<float>(margin));
        }
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    mean_loss<<<1, 256, 0, stream>>>(
        losses.data_ptr<float>(), result.data_ptr<float>(), rows);
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return result;
}
"""

_triplet_ext = load_inline(
    name="fused_triplet_margin_cuda",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["triplet_margin_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self, margin=1.0):
        super().__init__()
        if margin <= 0:
            raise ValueError("margin must be greater than zero")
        self.margin = float(margin)

    def forward(self, anchor, positive, negative):
        return _triplet_ext.triplet_margin_cuda(
            anchor, positive, negative, self.margin
        )
```