Hello!
I always used this code based from Intel code to compute the Screen Space AABB from the Bounding Sphere.
It's based from this sample code:
https://www.intel.com/content/www/us/en/developer/articles/technical/deferred-rendering-for-current-and-future-rendering-pipelines.html

void UpdateClipRegionRoot(in float nc, in float lc, in float lz, in float Radius, in float CameraScale, inout float ClipMin, inout float ClipMax)
{
    float nz = (Radius - nc * lc) / lz;
    float pz = (lc * lc + lz * lz - Radius * Radius) / (lz - (nz / nc) * lc);
    if (pz > 0.0f)
    {
        float c = -nz * CameraScale / nc;
        if (nc > 0.0f)
            ClipMin = max(ClipMin, c);
        else
            ClipMax = min(ClipMax, c);
    }
}

void UpdateClipRegion(in float lc, in float lz, in float Radius, in float CameraScale, inout float ClipMin, inout float ClipMax)
{
    float rSq = Radius * Radius;
    float lcSqPluslzSq = lc * lc + lz * lz;
    float d = rSq * lc * lc - lcSqPluslzSq * (rSq - lz * lz);
    if (d > 0.0f)
    {
        float a = Radius * lc;
        float b = sqrt(d);
        float nx0 = (a + b) / lcSqPluslzSq;
        float nx1 = (a - b) / lcSqPluslzSq;
        UpdateClipRegionRoot(nx0, lc, lz, Radius, CameraScale, ClipMin, ClipMax);
        UpdateClipRegionRoot(nx1, lc, lz, Radius, CameraScale, ClipMin, ClipMax);
    }
}

void ComputeClipRegion(in float3 Center, in float Radius, out float4 ClipRegion)
{
    ClipRegion = float4(1.0f, 1.0f, 0.0f, 0.0f);
    if ((Center.z + Radius) >= CameraNear)
    {
        float2 ClipMin = float2(-1.0f, -1.0f);
        float2 ClipMax = float2(+1.0f, +1.0f);
        UpdateClipRegion(Center.x, Center.z, Radius, Projection[0].x, ClipMin.x, ClipMax.x);
        UpdateClipRegion(Center.y, Center.z, Radius, Projection[1].y, ClipMin.y, ClipMax.y);
        ClipRegion = float4(ClipMin, ClipMax);
    }
}

void ComputeBoundingBox(in float3 Center, in float Radius, out float4 Bounds)
{
    ComputeClipRegion(Center, Radius, Bounds);
    Bounds = 0.5f * float4(Bounds.x, -Bounds.w, Bounds.z, -Bounds.y) + 0.5f;
}

It works very good but I wonder if there is a more modern and efficient way to do it?
Something to note is you have Point Light, Spot Light, Tube Light and Rect Light in lighting calculations.
I actually compute the bounding sphere in all these types so it's not the most possible efficient way.
But for point light in all cases you would need a Sphere to Screen Space AABB.

Thanks a lot for all the help and resources about this topic!
If you can recommend another way for the sphere to screen space AABB, it's very welcome too!

// returns the screen-space (normalized device coordinates) bounds of a projected sphere
// center = view-space center of the sphere
// radius = world or view space radius of the sphere
// boxMin = bottom left of projected bounds
// boxMax = top right of projected bounds

bool ProjectedSphereAABB(const vec3 &center, const float radius, vec3 &boxMin, vec3 &boxMax)
{
    /// sphere frustum culling 
    if (CullSphere(center,radius))
        return false;

    float d2 = dot(center,center);

    float a = sqrt(d2 - radius * radius);

    /// view-aligned "right" vector (right angle to the view plane from the center of the sphere. Since  "up" is always (0,n,0), replaced cross product with vec3(-c.z, 0, c.x)
    vec3 right = (radius / a) * vec3(-center.z, 0, center.x);
    vec3 up = vec3(0,radius,0);

    vec4 projectedRight  = Projection * vec4(right,0);
    vec4 projectedUp     = Projection * vec4(up,0);

    vec4 projectedCenter = Projection * vec4(center,1);

    vec4 north  = projectedCenter + projectedUp;
    vec4 east   = projectedCenter + projectedRight;
    vec4 south  = projectedCenter - projectedUp;
    vec4 west   = projectedCenter - projectedRight;

    north /= north.w ;
    east  /= east.w  ;
    west  /= west.w  ;
    south /= south.w ;

    boxMin = min(min(min(east,west),north),south).xyz;
    boxMax = max(max(max(east,west),north),south).xyz;


    return true;
}

I guess the CullSphere is simply culling with the near distance:

if ((Center.z + Radius) < CameraNear) {
return false;
}

The code you posted here is a lot more small and direct.
Is it as safe as the Intel version?
I wonder the performance difference between the two methods.

I understand the end of the algorithm which is to project on screen the right and up direction + project the center position on screen to compute the 4 points.
Actually, I think we could just project 2 points, the bottom left and upper right.

I have difficulty to understand the entire calculation of the up and right direction vector before to project them.
An explanation would be very appreciated.
Thanks!

If you project only the first two point you will get a perfect circle, a sphere in 3d is a round sphere but when projected into a 2d screen it may be distorted into an ellipsoid.

Cullsphere is a simple culling , correct on this, note that projectedcenter need your projection matrix, the rest is simple vector addition

I tested the two algorithms and I don't get the same result at all, did I do something wrong?

#include <iostream>
#include <cmath>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>

const float CameraNear = 0.1f;
const glm::mat4 Projection = glm::perspective(glm::radians(45.0f), 1280.0f / 720.0f, CameraNear, 100.0f);

void UpdateClipRegionRoot(float nc, float lc, float lz, float Radius, float CameraScale, float& ClipMin, float& ClipMax)
{
    float nz = (Radius - nc * lc) / lz;
    float pz = (lc * lc + lz * lz - Radius * Radius) / (lz - (nz / nc) * lc);
    if (pz > 0.0f)
    {
        float c = -nz * CameraScale / nc;
        if (nc > 0.0f)
            ClipMin = std::max(ClipMin, c);
        else
            ClipMax = std::min(ClipMax, c);
    }
}

void UpdateClipRegion(float lc, float lz, float Radius, float CameraScale, float& ClipMin, float& ClipMax)
{
    float rSq = Radius * Radius;
    float lcSqPluslzSq = lc * lc + lz * lz;
    float d = rSq * lc * lc - lcSqPluslzSq * (rSq - lz * lz);
    if (d > 0.0f)
    {
        float a = Radius * lc;
        float b = sqrt(d);
        float nx0 = (a + b) / lcSqPluslzSq;
        float nx1 = (a - b) / lcSqPluslzSq;
        UpdateClipRegionRoot(nx0, lc, lz, Radius, CameraScale, ClipMin, ClipMax);
        UpdateClipRegionRoot(nx1, lc, lz, Radius, CameraScale, ClipMin, ClipMax);
    }
}

void ComputeClipRegion(glm::vec3 Center, float Radius, glm::vec4& ClipRegion)
{
    ClipRegion = glm::vec4(1.0f, 1.0f, 0.0f, 0.0f);
    if ((Center.z + Radius) >= CameraNear)
    {
        glm::vec2 ClipMin = glm::vec2(-1.0f, -1.0f);
        glm::vec2 ClipMax = glm::vec2(+1.0f, +1.0f);
        UpdateClipRegion(Center.x, Center.z, Radius, Projection[0].x, ClipMin.x, ClipMax.x);
        UpdateClipRegion(Center.y, Center.z, Radius, Projection[1].y, ClipMin.y, ClipMax.y);
        ClipRegion = glm::vec4(ClipMin, ClipMax);
    }
}

bool ProjectedSphereAABB(glm::vec3 center, float radius, glm::vec3& boxMin, glm::vec3& boxMax)
{
    if (center.z + radius < CameraNear)
        return false;

    float d2 = dot(center,center);

    float a = sqrt(d2 - radius * radius);

    glm::vec3 right = (radius / a) * glm::vec3(-center.z, 0, center.x);
    glm::vec3 up = glm::vec3(0,radius,0);

    glm::vec4 projectedRight  = Projection * glm::vec4(right,0);
    glm::vec4 projectedUp     = Projection * glm::vec4(up,0);

    glm::vec4 projectedCenter = Projection * glm::vec4(center,1);

    glm::vec4 north  = projectedCenter + projectedUp;
    glm::vec4 east   = projectedCenter + projectedRight;
    glm::vec4 south  = projectedCenter - projectedUp;
    glm::vec4 west   = projectedCenter - projectedRight;

    north /= north.w ;
    east  /= east.w  ;
    west  /= west.w  ;
    south /= south.w ;

    boxMin = glm::min(glm::min(glm::min(east,west),north),south);
    boxMax = glm::max(glm::max(glm::max(east,west),north),south);
    
    return true;
}

int main()
{
    glm::vec3 pos = glm::vec3(1.0f, 1.0f, 25.0f);

    glm::vec4 bounds;
    ComputeClipRegion(pos, 5.0f, bounds);
    std::printf("%f, %f, %f, %f\n", bounds.x, bounds.y, bounds.z, bounds.w);

    glm::vec3 min, max;
    ProjectedSphereAABB(pos, 5.0f, min, max);
    std::printf("%f, %f, %f, %f\n", min.x, min.y, max.x, max.y);
}

I get this result:

-0.220847, -0.392618, 0.334014, 0.593802
-0.333779, -0.579411, 0.220621, 0.386274

If we invert the min and max of Y and inverting the sign we are close but not the same:

-0.220847, -0.392618, 0.334014, 0.593802
-0.333779, -0.386274, 0.220621, 0.579411

I did a benchmark of the two methods and it looks like the Intel method is a lot faster:

int main()
{
    glm::vec3 pos = glm::vec3(1.0f, 1.0f, 25.0f);

    glm::vec4 bounds;
    uint64_t totalTime = 0;
    for (int i = 0; i < 10'000; ++i)
    {
        auto start = std::chrono::system_clock::now();

        ComputeClipRegion(pos, 5.0f, bounds);

        auto end = std::chrono::system_clock::now();
        auto duration = end - start;
        totalTime += std::chrono::duration_cast<std::chrono::nanoseconds>(duration).count();
    }
    totalTime /= 10'000;
    std::printf("%d\n", totalTime);

    glm::vec3 min, max;
    totalTime = 0;
    for (int i = 0; i < 10'000; ++i)
    {
        auto start = std::chrono::system_clock::now();

        ProjectedSphereAABB(pos, 5.0f, min, max);

        auto end = std::chrono::system_clock::now();
        auto duration = end - start;
        totalTime += std::chrono::duration_cast<std::chrono::nanoseconds>(duration).count();
    }
    totalTime /= 10'000;
    std::printf("%d\n", totalTime);
}

Result:

148
946

Are you sure that you have included in your computations also the sphere culling ? also glm is slow

Yes the sphere culling is here in the Intel method:

void ComputeClipRegion(glm::vec3 Center, float Radius, glm::vec4& ClipRegion)
{
    ClipRegion = glm::vec4(1.0f, 1.0f, 0.0f, 0.0f);
    if ((Center.z + Radius) >= CameraNear)
    {

In the method you presented here I added this for this purpose:

bool ProjectedSphereAABB(glm::vec3 center, float radius, glm::vec3& boxMin, glm::vec3& boxMax)
{
    if (center.z + radius < CameraNear)
        return false;

Maybe your method give more precise result for the cost of performance?
---
I understood finally the calculation of the square root, if you draw on paper you will see that a right triangle is formed, the square root allows you to compute the distance to the border of the circle.

In my method I have used the squared distance , probably something else is taking its toll, I will look into this immediately

My bottleneck was in the projection function, I have adopted the intel's method