Loading Document Achar o vetor normal é uma tarefa fácil e computacionalmente simples (basta realizar o produto vetorial entre dois vetores que definem o plano), entretanto, o vetor normal deve ser normalizado para que haja consistência no resultado em diferentes superfícies. Para isso se faz necessário saber a recíproca da raiz quadrada de uma variável, o que é computacionalmente custuso. Aqui veremos um algoritmo capaz de realizar essa operação de forma significativamente mais rápida, sem perda relevante de precisão em relação às ferramentas nativas do Python. Usaremos a biblioteca numpy para ter acesso a algumas funcionalidades.
Existe uma forma padronizada de representar valores flutuantes na computação, baseada na notação científica em base 2, já que toda informação é representada em bits. Essa notação é composta por três partes: o sinal, o expoente e a mantissa (ou fração).
Neste projeto, utilizaremos números de ponto flutuante de 32 bits, conforme o padrão IEEE 754. Segundo esse padrão:
Como lidaremos apenas com números positivos, sabemos de antemão que o primeiro bit de todos os números com que trabalharmos será sempre 0.
Como o expoente tem 8 bits, sabemos que ele pode assumir valores de 0 a 255. No entanto, também queremos representar expoentes negativos. Para isso, utilizamos um viés (bias) de 127, ou seja, subtraímos 127 do valor armazenado para obter o valor real do expoente.
Isso nos dá um intervalo efetivo de aproximadamente −126 a +127, já que os valores 0 e 255 são reservados para representações especiais (como zero, infinitos e NaNs).
Já no caso da mantissa, como o primeiro dígito da notação científica em base 2 é sempre 1, podemos representá-la armazenando apenas a parte fracionária. Esse primeiro dígito é chamado de bit implícito e não precisa ser guardado, o que economiza um bit e permite maior precisão na parte decimal.
Assim, tendo uma mantissa M e um expoente E, sabemos que o valor representado pelo número flutuante é:
V=(1+223M)⋅2E−127.
Sabemos também que se os mesmo bits fossem colocados na mesma ordem, mas representassem um número inteiro, o seu valor seria:
R=223E+M.
Esses dois resultados sao muito importantes pois eles tem um relação imprevisível que vai nos ajudar a resolver o problema inicial.
Queremos achar Vy=√Vx1. Tirando o logaritmo na base 2 dos dois lados, podemos deixar de nos preocupar com a raiz.
Vy=√Vx1⟹log2(Vx)=−21log2(Vx).
Em azul: f(x)=log2(x+1); Em verde: f(x)=x+μ
Com isso, temos tudo o que precisamos para resolver o problema inical. Vamos passo a passo para englobar todas as ideias.
Vy=√Vx1⟹log2(Vy)=−21log2(Vx)⟹223Ry+μ−127≈−21(223Rx+μ−127)⟹Ry≈3⋅222(127−μ)−2Rx⟹Ry≈C−2Rx
Isso no levou a uma relação da disposição dos bits do argumento com a disposição dos bits da resposta, já que C é uma constante. Usando o numpy , podemos alternar como lemos os bits (como flutuante ou inteiro) e, portanto, resolver se forma extremamente rápida o problema.
A última coisa que precisamos é achar essa constante C. Vale ressaltar que, apesar de ser possível encontrar matematicamente o melhor C a diminuir o erro, dado a precisão da máquina em determinadas distribuições de entrada, faz-se mais preciso usar um valor empírico diferente para essa constante. Entretanto, o valor matemático de C ainda é excelente na prática.
Vamos passar o passo a passo de como encontrar matematicamente o valor de C. Como C depende de μ, precisamos primeiro definir μ. Essa constante deve ser a que melhor aproxima as funções f(x)=log2(x+1) e g(x)=x+μ no intervalo [0,1]. Ou seja, queremos μ tal que a f(x)−g(x) seja 0 no intervalo [0,1]. Dessa forma, temos:
∫01(f(x)−g(x))dx=0⟹∫01(log2(x+1)−x−μ)dx=0⟹μ=ln21∫01ln(x+1)dx−21⟹μ=ln21(2ln(2)−1)−21⟹μ=23−ln21
Como C=3⋅222(127−μ), temos C=3⋅222(127−23+ln21)⟹C=3⋅221(251+ln22). Agora temos tudo que precisamos para encontrar uma ótima aproximação da recíproca da raiz quadrada de qualquer valor.
Apesar de já termos um ótimo método para resolver nosso problema, ainda podemos melhorá-lo usando o método Newton-Raphson. Pense da seguinte forma: ao aplicarmos nosso método, temos Vy tal que:
Vy≈√Vx1⟹Vy21−Vx≈0
dVdf(Vy)≈Vy−V0Vy21−Vx−0⟹−Vy32≈Vy−V0Vy21−Vx⟹V0−Vy≈2Vy(1−VxVy2)⟹V0≈Vy(23−2VxVy2)
Assim, com a nossa aproximação, nos conseguimos melhorar a precisão do resultado. Note que podemos reutilizar esse método quantas vezes quisermos, entretanto normalmente utilizar-lo uma vez já nos dá um resultado com precisão excepcional.
Focando na implementação, temos:
import numpy as np
def fastInvSqrtRoot(n: float) -> float:
# aproximação da constante C passado em hexadecimal
# valor matemático: 0x5f34ff58
# valor empírico: 0x5f3759df
C = np.int32(0x5f3759df)
# convertento n em np.float32 para acessar a leitura de bits
n = np.float32(n)
# interpretando os bits como inteiro
i = n.view(int32)
# encontrando os bits do resultado
i = C - (i >> 1)
# interpretando os bits do resultado como flutuante
y = i.view(float32)
# método Newton-Raphson
y *= (1.5 - 0.5 * n * y * y)
return float(y)