Можно ли написать быструю функцию InvSqrt () Quake в Rust?

101

Это просто для удовлетворения моего собственного любопытства.

Есть ли реализация этого:

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;
}

в ржавчине? Если он существует, отправьте код.

Я попробовал это и потерпел неудачу. Я не знаю, как кодировать число с плавающей точкой, используя целочисленный формат. Вот моя попытка:

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
}

Ссылка:
1. Происхождение Быстрого InvSqrt () Quake3 - Страница 1
2. Понимание быстрого обратного квадратного корня Quake
3. БЫСТРАЯ ОБРАТНАЯ КВАДРАТНАЯ КОРНА.pdf
4. Исходный код: q_math.c # L552-L572

Flyq
источник
4
Насколько я понимаю, этот код UB в C из-за нарушения строгого правила псевдонимов . Стандартный способ выполнения этого вида наказания - с union.
trentcl
4
@trentcl: я тоже не думаю, что unionработает. memcpyопределенно работает, хотя это многословно.
Матье М.
14
@MatthieuM. Тип Punning с помощью союзов является совершенно допустимым C , но не допустимым C ++.
Мойра
4
Я полагаю, что этот вопрос хорош с точки зрения чистого любопытства, но, пожалуйста, поймите, что времена изменились. На x86, то rsqrtssи rsqrtpsинструкция, введенная с Pentium III в 1999 году, быстрее и точнее , чем этот код. ARM NEON имеет vrsqrteчто похоже. И какие бы вычисления Quake III ни использовал для этого, в любом случае, вероятно, в наши дни будет сделано на GPU.
Бенг

Ответы:

87

Я не знаю, как кодировать число с плавающей точкой, используя целочисленный формат.

Для этого есть функция, f32::to_bitsкоторая возвращает u32. Существует также функция для другого направления: она f32::from_bitsпринимает u32аргумент as. Эти функции предпочтительнее, чем mem::transmuteпоследние, unsafeи их сложно использовать.

С этим, вот реализация 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)
}

( Детская площадка )


Эта функция компилируется в следующую сборку на 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

Я не нашел ни одной справочной сборки (если есть, скажите, пожалуйста!), Но она мне кажется довольно хорошей. Я просто не уверен, почему поплавок был перемещен eaxтолько для того, чтобы сделать сдвиг и целочисленное вычитание. Может быть, регистры SSE не поддерживают эти операции?

clang 9.0 с -O3компилирует код C в основном в ту же сборку . Так что это хороший знак.


Стоит отметить, что если вы действительно хотите использовать это на практике: пожалуйста, не делайте этого. Как отметил Бенгр в комментариях , современные процессоры x86 имеют специальную инструкцию для этой функции, которая быстрее и точнее, чем этот хак. К сожалению, 1.0 / x.sqrt() похоже, не оптимизировать под эту инструкцию . Так что, если вам действительно нужна скорость, возможно,_mm_rsqrt_ps стоит использовать встроенные функции. Это, однако, снова требует unsafeкода. Я не буду вдаваться в подробности в этом ответе, так как на самом деле это нужно меньшинству программистов.

Лукас Калбертодт
источник
4
Согласно Intel Intrinsics Guide, нет операции целочисленного сдвига, которая только сдвигает самый младший 32-разрядный из 128-разрядного регистра в аналог addssили mulss. Но если другие 96 бит xmm0 можно игнорировать, тогда можно использовать psrldинструкцию. То же самое касается целочисленного вычитания.
fsasm
Я признаю, что почти ничего не знал о ржавчине, но разве «небезопасно» не является основным свойством fast_inv_sqrt? С его полным неуважением к типам данных и тому подобное.
Gloweye
12
@Gloweye Это другой тип "небезопасных", о которых мы говорим. Быстрое приближение, которое получает плохую ценность слишком далеко от сладкого места, по сравнению с чем-то быстрым и свободным с неопределенным поведением.
дедупликатор
8
@Gloweye: математически, последняя часть этого fast_inv_sqrt- всего лишь один шаг итерации Ньютона-Рафсона, чтобы найти лучшее приближение inv_sqrt. В этой части нет ничего опасного. Обман в первой части, которая находит хорошее приближение. Это работает, потому что это делает целочисленное деление на 2 в экспоненциальной части числа с плавающей точкой, и действительноsqrt(pow(0.5,x))=pow(0.5,x/2)
MSalters
1
@fsasm: это правильно; movdв EAX и обратно - это пропущенная оптимизация текущих компиляторов. (И да, соглашения о вызовах передают / возвращают скаляр floatв элементе low XMM и допускают, что старшие биты могут быть мусором. Но обратите внимание, что, если он был расширен с нуля, он может легко остаться таким: сдвиг вправо не приводит к ноль элементов и ни одно из них не вычитает из _mm_set_epi32(0,0,0,0x5f3759df), т. е. movdзагружает. Вам нужно было movdqa xmm1,xmm0бы скопировать регистр раньше psrld. Обойти задержку от FP, переадресация команды на целое число и наоборот скрыта за mulssзадержкой
Питер Кордес
37

Этот реализован с менее известным unionв 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)
    }
}

Сделал некоторые микро тесты, используя criterioncrate на x86-64 Linux box. Удивительно, но у Руста sqrt().recip()самое быстрое. Но, конечно, любой результат микро-теста должен быть взят с зерном соли.

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
источник
22
Я нисколько не удивлен sqrt().inv(), быстрее всех. И sqrt, и inv - это отдельные инструкции, и они идут довольно быстро. Doom был написан в те дни, когда было небезопасно предполагать, что аппаратная плавающая точка вообще существует, а трансцендентные функции, такие как sqrt, определенно были бы программными. +1 за отметки.
Мартин Боннер поддерживает Монику
4
Что меня удивляет , что по- transmuteвидимому , отличаются от to_и from_bits- я бы ожидать тех научения эквивалента даже до оптимизации.
trentcl
2
@MartinBonner (Кроме того, это не имеет значения, но sqrt не является трансцендентной функцией .)
benrg
4
@MartinBonner: любой аппаратный FPU, который поддерживает деление, обычно также поддерживает sqrt. IEEE «базовые» операции (+ - * / sqrt) необходимы для получения правильно округленного результата; Вот почему SSE предоставляет все эти операции, но не exp, sin или что-то еще. На самом деле, split и sqrt обычно работают на одном и том же исполнительном модуле, сконструированном аналогичным образом. См. HW div / sqrt unit details . В любом случае, они все еще не быстрые по сравнению с умножением, особенно в латентном режиме.
Питер Кордес
1
В любом случае, Skylake имеет значительно лучшую конвейеризацию для div / sqrt, чем предыдущие Uarches. Посмотрите деление с плавающей точкой против умножения с плавающей точкой для некоторых выдержек из таблицы Агнер Фог. Если вы не выполняете много другой работы в цикле, так что sqrt + div является узким местом, вы можете использовать HW быстрый взаимный sqrt (вместо взлома Quake) + итерация Ньютона. Особенно с FMA, это хорошо для пропускной способности, если не задержки. Быстрая векторизация rsqrt и обратная с SSE / AVX в зависимости от точности
Питер Кордес
10

Вы можете использовать, std::mem::transmuteчтобы сделать необходимое преобразование:

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
}

Вы можете посмотреть живой пример здесь: здесь

Настоящий Свежий
источник
4
В небезопасном нет ничего плохого, но есть способ сделать это без явного небезопасного блока, поэтому я бы предложил переписать этот ответ, используя f32::to_bitsи f32::from_bits. Он также несет в себе цель, явно отличную от трансмутации, которую большинство людей, вероятно, считают «магией».
Sahsahae
5
@Sahsahae Я только что опубликовал ответ, используя две функции, которые вы упомянули :) И я согласен, здесь unsafeследует избегать, так как в этом нет необходимости.
Лукас Калбертодт