É possível escrever a função InvSqrt () rápida do Quake no Rust?

101

Isso é apenas para satisfazer minha própria curiosidade.

Existe uma implementação disso:

float InvSqrt (float x)
{
   float xhalf = 0.5f*x;
   int i = *(int*)&x;
   i = 0x5f3759df - (i>>1);
   x = *(float*)&i;
   x = x*(1.5f - xhalf*x*x);
   return x;
}

em Rust? Se existir, publique o código.

Eu tentei e falhei. Não sei como codificar o número da bóia usando o formato inteiro. Aqui está a minha tentativa:

fn main() {
    println!("Hello, world!");
    println!("sqrt1: {}, ",sqrt2(100f64));
}

fn sqrt1(x: f64) -> f64 {
    x.sqrt()
}

fn sqrt2(x: f64) -> f64 {
    let mut x = x;
    let xhalf = 0.5*x;
    let mut i = x as i64;
    println!("sqrt1: {}, ", i);

    i = 0x5f375a86 as i64 - (i>>1);

    x = i as f64;
    x = x*(1.5f64 - xhalf*x*x);
    1.0/x
}

Referência:
1. Origem do Rápido InvSqrt do Quake3 () - Página 1
2. Entendendo a Raiz Quadrada Inversa Rápida do Quake3
3. RÁPIDO QUADRADO INVERSO RÁPIDO.pdf
4. código fonte: q_math.c # L552-L572

Flyq
fonte
4
Pelo que entendi, esse código é UB em C devido à violação da regra estrita de alias . A maneira abençoada padrão de executar esse tipo de punção de tipo é com a union.
trentcl
4
@trentcl: Eu também não acho que unionfunciona. memcpydefinitivamente funciona, embora seja detalhado.
Matthieu M.
14
@MatthieuM. Punição de tipo com uniões é C perfeitamente válido , mas não C ++ válido.
Moira
4
Suponho que essa questão seja boa do ponto de vista da pura curiosidade, mas entenda que os tempos mudaram. No x86, as instruções rsqrtsse rsqrtps, introduzidas no Pentium III em 1999, são mais rápidas e precisas que esse código. BRAÇO NEON tem o vrsqrteque é semelhante. E qualquer que seja o cálculo usado pelo Quake III, provavelmente seria feito na GPU hoje em dia.
benrg

Respostas:

87

Não sei como codificar o número da bóia usando o formato inteiro.

Existe uma função para isso: f32::to_bitsque retorna um u32. Há também a função para a outra direção: f32::from_bitsque aceita um u32argumento as. Essas funções são preferidas, mem::transmutepois a última é unsafedifícil de usar.

Com isso, aqui está a implementação de InvSqrt:

fn inv_sqrt(x: f32) -> f32 {
    let i = x.to_bits();
    let i = 0x5f3759df - (i >> 1);
    let y = f32::from_bits(i);

    y * (1.5 - 0.5 * x * y * y)
}

( Parque infantil )


Essa função é compilada no seguinte assembly em x86-64:

.LCPI0_0:
        .long   3204448256        ; f32 -0.5
.LCPI0_1:
        .long   1069547520        ; f32  1.5
example::inv_sqrt:
        movd    eax, xmm0
        shr     eax                   ; i << 1
        mov     ecx, 1597463007       ; 0x5f3759df
        sub     ecx, eax              ; 0x5f3759df - ...
        movd    xmm1, ecx
        mulss   xmm0, dword ptr [rip + .LCPI0_0]    ; x *= 0.5
        mulss   xmm0, xmm1                          ; x *= y
        mulss   xmm0, xmm1                          ; x *= y
        addss   xmm0, dword ptr [rip + .LCPI0_1]    ; x += 1.5
        mulss   xmm0, xmm1                          ; x *= y
        ret

Não encontrei nenhum conjunto de referência (se tiver, por favor me diga!), Mas me parece bastante bom. Só não sei por que o flutuador foi movido eaxapenas para fazer a subtração de deslocamento e número inteiro. Talvez os registros SSE não suportem essas operações?

O clang 9.0 with -O3compila o código C basicamente para o mesmo assembly . Então esse é um bom sinal.


Vale ressaltar que, se você realmente deseja usar isso na prática: por favor, não. Como o benrg apontou nos comentários , as modernas CPUs x86 têm uma instrução especializada para essa função, que é mais rápida e precisa do que esse hack. Infelizmente, 1.0 / x.sqrt() parece não otimizar essa instrução . Portanto, se você realmente precisa da velocidade, usar os _mm_rsqrt_psintrínsecos provavelmente é o caminho a percorrer. Isso, no entanto, exige novamente unsafecódigo. Não vou entrar em muitos detalhes nesta resposta, pois uma minoria de programadores realmente precisará dela.

Lukas Kalbertodt
fonte
4
De acordo com o Intel Intrinsics Guide, não há operação de deslocamento inteiro que desloque apenas os 32 bits mais baixos do registrador analógico de 128 bits para addssou mulss. Mas se os outros 96 bits de xmm0 puderem ser ignorados, pode-se usar a psrldinstrução. O mesmo vale para subtração de número inteiro.
Fsasm 28/11/19
Admito que não sei quase nada sobre ferrugem, mas não é "inseguro" basicamente uma propriedade essencial do fast_inv_sqrt? Com seu total desrespeito por tipos de dados e coisas do gênero.
Gloweye
12
@ Gloweye É um tipo diferente de "inseguro" sobre o qual falamos. Uma aproximação rápida que obtém um valor ruim muito longe do ponto ideal, em comparação com algo jogando rápido e solto com um comportamento indefinido.
Deduplicator
8
@Gloweye: Matematicamente, a última parte disso fast_inv_sqrté apenas uma etapa da iteração de Newton-Raphson para encontrar uma melhor aproximação de inv_sqrt. Não há nada de inseguro nessa parte. O truque está na primeira parte, que encontra uma boa aproximação. Isso funciona porque ele está fazendo uma divisão inteira por 2 da parte expoente do flutuador, e de fatosqrt(pow(0.5,x))=pow(0.5,x/2)
MSalters
11
@ fsasm: Isso está correto; movdpara EAX e vice-versa é uma otimização perdida pelos compiladores atuais. (E sim, as convenções de chamada passam por escalar / retornar escalar floatno elemento baixo de um XMM e permitem que os bits altos sejam lixo. Mas observe que, se ele foi estendido para zero, pode facilmente ficar assim: a mudança à direita não introduz não- zero elementos e nem subtração _mm_set_epi32(0,0,0,0x5f3759df), ou seja, uma movdcarga Você iria precisar de um. movdqa xmm1,xmm0para copiar o reg antes psrldBypass latência encaminhe instrução FP para inteiro e vice-versa está oculta por. mulsslatência.
Peter Cordes
37

Este é implementado com menos conhecido unionem Rust:

union FI {
    f: f32,
    i: i32,
}

fn inv_sqrt(x: f32) -> f32 {
    let mut u = FI { f: x };
    unsafe {
        u.i = 0x5f3759df - (u.i >> 1);
        u.f * (1.5 - 0.5 * x * u.f * u.f)
    }
}

Fiz alguns micro benchmarks usando criterionengradado em uma caixa Linux x86-64. Surpreendentemente, o próprio Rust sqrt().recip()é o mais rápido. Mas é claro que qualquer resultado de micro benchmark deve ser obtido com um grão de sal.

inv sqrt with transmute time:   [1.6605 ns 1.6638 ns 1.6679 ns]
inv sqrt with union     time:   [1.6543 ns 1.6583 ns 1.6633 ns]
inv sqrt with to and from bits
                        time:   [1.7659 ns 1.7677 ns 1.7697 ns]
inv sqrt with powf      time:   [7.1037 ns 7.1125 ns 7.1223 ns]
inv sqrt with sqrt then recip
                        time:   [1.5466 ns 1.5488 ns 1.5513 ns]
edwardw
fonte
22
Não estou nem um pouco surpreso sqrt().inv()é o mais rápido. Atualmente, sqrt e inv são instruções únicas e são bastante rápidas. Doom foi escrito nos dias em que não era seguro supor que houvesse ponto flutuante de hardware, e funções transcendentais como o sqrt definitivamente seriam software. +1 para os benchmarks.
Martin Bonner apoia Monica
4
O que me surpreende é que transmuteaparentemente é diferente to_e from_bits- eu esperaria que eles fossem equivalentes a instruções antes mesmo da otimização.
trentcl
2
@MartinBonner (Também, não que isso importe, mas sqrt não é uma função transcendente .)
benrg
4
@ MartinBonner: Qualquer FPU de hardware que suporte a divisão normalmente também suporta o sqrt. As operações "básicas" do IEEE (+ - * / sqrt) são necessárias para produzir um resultado arredondado corretamente; é por isso que o SSE fornece todas essas operações, mas não exp, sin ou qualquer outra coisa. De fato, a divisão e o sqrt geralmente são executados na mesma unidade de execução, projetada de maneira semelhante. Consulte os detalhes da unidade div / sqrt de HW . De qualquer forma, eles ainda não são rápidos em comparação à multiplicação, especialmente em latência.
Peter Cordes
11
De qualquer forma, a Skylake possui um pipelining significativamente melhor para div / sqrt do que os uarches anteriores. Veja Divisão de ponto flutuante versus multiplicação de ponto flutuante para alguns extratos da tabela de Agner Fog. Se você não está fazendo muito outro trabalho em loop, então sqrt + div é um gargalo, convém usar o sqrt recíproco rápido do HW (em vez do hack do terremoto) + uma iteração de Newton. Especialmente com FMA, isso é bom para taxa de transferência, senão latência. Rápido vetorizado RSQRT e recíproca com SSE / AVX dependendo precisão
Peter Cordes
10

Você pode usar std::mem::transmutepara fazer a conversão necessária:

fn inv_sqrt(x: f32) -> f32 {
    let xhalf = 0.5f32 * x;
    let mut i: i32 = unsafe { std::mem::transmute(x) };
    i = 0x5f3759df - (i >> 1);
    let mut res: f32 = unsafe { std::mem::transmute(i) };
    res = res * (1.5f32 - xhalf * res * res);
    res
}

Você pode procurar um exemplo ao vivo aqui: aqui

Real Fresh
fonte
4
Não há nada errado com inseguro, mas há uma maneira de fazer isso sem bloqueio explícito inseguro, então sugiro reescrever esta resposta usando f32::to_bitse f32::from_bits. Ele também carrega a intenção claramente diferente de transmutar, que a maioria das pessoas provavelmente considera "mágica".
Sahsahae 28/11/19
5
@Sahsahae Acabei de postar uma resposta usando as duas funções que você mencionou :) E eu concordo, unsafedeve ser evitado aqui, pois não é necessário.
Lukas Kalbertodt 28/11/19