Можно ли реализовать быстрый обратный квадратный корень из Quake на Rust?

Безопасная и эффективная реализация InvSqrt на Rust

Да, реализовать InvSqrt() можно, и Rust предлагает несколько способов. Начнём с главного: unsafe здесь оправдан, потому что мы работаем на уровне представления чисел в памяти — это низкоуровневая оптимизация, которую компилятор не может проверить статически. Однако есть способы сделать код более идиоматичным.

1. Решение через f32::to_bits() и f32::from_bits() (рекомендуется)

В стандартной библиотеке есть безопасные методы для reinterpret_cast между f32 и u32 (не i32!):

fn inv_sqrt(x: f32) -> f32 {
    let x2 = x * 0.5;
    let mut i = x.to_bits();  // f32 -> u32, безопасно
    i = 0x5f3759df - (i >> 1);
    let mut y = f32::from_bits(i); // u32 -> f32, безопасно
    y = y * (1.5 - x2 * y * y);
    y
}

Это полностью безопасный Rust без unsafe. Методы to_bits() и from_bits() определены в стандартной библиотеке и гарантированно корректны. Единственный нюанс: мы используем u32, а не i32, но для побитовых операций это неважно.

2. Почему transmute — не лучший выбор?

• transmute требует unsafe и может привести к UB, если типы не совпадают по размеру.
• to_bits()/from_bits() делают то же самое, но безопасно и с проверками на этапе компиляции.
• Производительность идентична — оба метода компилируются в одинаковый машинный код.

3. Полная версия с const и тестами

const MAGIC: u32 = 0x5f3759df;
const THREE_HALVES: f32 = 1.5;

fn inv_sqrt(x: f32) -> f32 {
    let x2 = x * 0.5;
    let i = MAGIC - (x.to_bits() >> 1);
    let y = f32::from_bits(i);
    y * (THREE_HALVES - x2 * y * y) // одна итерация Ньютона
}

#[cfg(test)]
mod tests {
    use super::*;
    
    #[test]
    fn test_accuracy() {
        let x = 4.0;
        let result = inv_sqrt(x);
        let expected = 1.0 / x.sqrt();
        let diff = (result - expected).abs();
        assert!(diff < 1e-6, "Разница слишком велика: {}", diff);
    }
    
    #[test]
    fn test_positive() {
        for i in 1..1000 {
            let x = i as f32;
            let result = inv_sqrt(x);
            assert!(result.is_finite() && result > 0.0);
        }
    }
}

4. Производительность

На современных процессорах аппаратная инструкция rsqrtss (доступна через _mm_rsqrt_ss из std::arch) может быть быстрее. Но Quake-версия остаётся актуальной для embedded-систем без FPU или для консольных процессоров.

Вывод: Используйте to_bits() и from_bits() — это безопасно, идиоматично для Rust и даёт ту же производительность, что и оригинальная C-версия.

#[cfg(target_arch = "x86_64")]
fn inv_sqrt_simd(x: f32) -> f32 {
    unsafe {
        let a = std::arch::x86_64::_mm_set_ss(x);
        let r = std::arch::x86_64::_mm_rsqrt_ss(a);
        std::arch::x86_64::_mm_cvtss_f32(r)
    }
}

let approx = inv_sqrt_simd(x);
approx * (1.5 - 0.5 * x * approx * approx)