math::Sphere::OptimalEnclosingSphere() breaks down for torus with sufficient detail about mathgeolib HOT 6 OPEN

#include <stdio.h>
#include <stdlib.h>

#include "../src/MathGeoLib.h"
#include "../src/Geometry/Sphere.h"
#include "../src/Math/myassert.h"
#include "TestRunner.h"
#include "TestData.h"

MATH_IGNORE_UNUSED_VARS_WARNING

MATH_BEGIN_NAMESPACE

std::vector<math::vec> GenerateTorus(
    double major_radius,
    double minor_radius,
    int    major_axis_steps,
    int    minor_axis_steps
)
{
    std::vector<math::vec> points;
    constexpr double pi = 3.14159265358979323846264338327950288;
    const double R = major_radius;
    const double r = minor_radius;
    for (int major = 0; major < major_axis_steps; ++major)
    {
        const double rel_major = static_cast<double>(major) / static_cast<double>(major_axis_steps);
        for (int minor = 0; minor < minor_axis_steps; ++minor)
        {
            const double rel_minor = static_cast<double>(minor) / static_cast<double>(minor_axis_steps);
            const double theta     = (pi * 2.0 * rel_major);
            const double phi       = (pi * 2.0 * rel_minor);
            const double sin_theta = std::sin(theta);
            const double cos_theta = std::cos(theta);
            const double sin_phi   = std::sin(phi);
            const double cos_phi   = std::cos(phi);
            const double x         = (R + r * cos_phi) * cos_theta;
            const double z         = (R + r * cos_phi) * sin_theta;
            const double y         =      r * sin_phi;
            const float  fx        = static_cast<float>(x);
            const float  fy        = static_cast<float>(y);
            const float  fz        = static_cast<float>(z);
            points.push_back(math::vec{fx, fy, fz});
        }
    }
    return points;
}

UNIQUE_TEST(BoundingSphere)
{
    const auto points = GenerateTorus(0.5, 0.125, 40, 32);

    const math::Sphere sphere = math::Sphere::OptimalEnclosingSphere(
        points.data(),
        static_cast<int>(points.size())
    );
    assert(sphere.r < (0.5 + 0.125 + 0.001));
}

MATH_END_NAMESPACE