Como esta aproximação de raiz quadrada flutuante funciona?

Encontrei uma aproximação de raiz quadrada bastante estranha, mas funcional, parafloats; Eu realmente não entendo. Alguém pode me explicar por que esse código funciona?

float sqrt(float f)
{
    const int result = 0x1fbb4000 + (*(int*)&f >> 1);
    return *(float*)&result;   
}

Eu testei um pouco eit gera valores fora destd::sqrt() em cerca de 1 a 3%. Eu sei do @ do Quake Iaiz quadrada inversa rápida e acho que é algo semelhante aqui (sem a iteração newton), mas eu realmente aprecio uma explicação sobrecomo funcion.

(nota: eu marquei ambos osc ec ++ desde que é válido-ish (ver comentários) código C e C ++)