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

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

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

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

в Rust? Если он существует, опубликуйте код.

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

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. Origin of Quake3's Fast InvSqrt () - страница 1 < br> 2. Что такое быстрый обратный квадратный корень Quake
3. БЫСТРЫЙ ОБРАТНЫЙ КВАДРАТНЫЙ КОРЕНЬ.pdf
4. исходный код: q_math.c # L552 -L572


rust floating-point sqrt type-punning
person Flyq    schedule 28.11.2019    source источник
comment
Версия C #: Можно ли написать Quake быстрая функция InvSqrt () в C #?   -  person Flyq    schedule 28.11.2019
comment
Насколько я понимаю, этот код является UB на C из-за нарушения правила строгого псевдонима . Стандартный способ использовать этот тип каламбура - использовать union.   -  person trentcl    schedule 28.11.2019
comment
@trentcl: Я тоже не думаю, что union работает. memcpy определенно работает, хотя и многословно.   -  person Matthieu M.    schedule 28.11.2019
comment
@MatthieuM. Типовой каламбур с объединениями - это совершенно допустимый C, но не действительный C ++.   -  person Moira    schedule 28.11.2019
comment
@Bergi Threre - логическая ошибка, надеюсь, она вас не побеспокоит: неправильный код   -  person Flyq    schedule 29.11.2019
comment
Я полагаю, что этот вопрос хорош с точки зрения чистого любопытства, но, пожалуйста, поймите, что времена изменились. На x86 инструкции rsqrtss и rsqrtps, представленные в Pentium III в 1999 году, работают быстрее и точнее, чем этот код. ARM NEON имеет vrsqrte, что похоже. И какие бы вычисления Quake III ни использовал для этого, вероятно, в наши дни все равно будет выполняться на GPU.   -  person benrg    schedule 30.11.2019
comment
Термин, который вы ищете, - это ввести каламбур с плавающей точкой в ​​целое число.   -  person Peter Cordes    schedule 30.11.2019

Ответы (3)


arrow_upward
95
arrow_downward

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

Для этого есть функция: _1 _ , который возвращает u32. Также есть функция для другого направления: f32::from_bits, который принимает u32 в качестве аргумента. Эти функции предпочтительнее, чем 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 в в основном ту же сборку. Так что это хороший знак.

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

person Lukas Kalbertodt    schedule 28.11.2019
comment
Согласно руководству Intel Intrinsics Guide не существует операции целочисленного сдвига, которая сдвигает только самые младшие 32-битные из 128-битных аналоговых регистров на addss или mulss. Но если остальные 96 бит xmm0 можно игнорировать, тогда можно использовать инструкцию psrld. То же самое и с целочисленным вычитанием. - person fsasm; 28.11.2019
comment
Признаюсь, я почти ничего не знаю о ржавчине, но разве небезопасно в основном свойство fast_inv_sqrt? С полным неуважением к типам данных и тому подобному. - person Gloweye; 28.11.2019
comment
@Gloweye Но мы говорим о другом типе небезопасности. Быстрое приближение, которое получает плохую ценность слишком далеко от сладкого места, по сравнению с чем-то быстрым и бесполезным с неопределенным поведением. - person Deduplicator; 28.11.2019
comment
@Gloweye: Математически последняя часть этого fast_inv_sqrt - это всего лишь один шаг итерации Ньютона-Рафсона, чтобы найти лучшее приближение inv_sqrt. В этой части нет ничего опасного. Уловка заключается в первой части, которая находит хорошее приближение. Это работает, потому что он делает целочисленное деление на 2 в экспоненциальной части числа с плавающей запятой, и действительно sqrt(pow(0.5,x))=pow(0.5,x/2) - person MSalters; 29.11.2019
comment
@fsasm: Верно; movd в EAX и обратно - это упущенная оптимизация текущими компиляторами. (И да, соглашения о вызовах передают / возвращают скаляр float в нижнем элементе XMM и позволяют старшим битам быть мусором. Но обратите внимание, что если он был расширен нулем, он может легко остаться таким: сдвиг вправо не вводит ненулевые элементы, как и вычитание из _mm_set_epi32(0,0,0,0x5f3759df), т.е. загрузка movd. Вам понадобится movdqa xmm1,xmm0, чтобы скопировать регистр перед psrld. Обход задержки при пересылке инструкций FP в целое число и наоборот скрывается задержкой mulss . - person Peter Cordes; 30.11.2019

arrow_upward
41
arrow_downward

Этот реализован с помощью менее известного 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)
    }
}

Провел несколько микротестов с использованием criterion crate на Linux x86-64. Удивительно, но собственный 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]
person edwardw    schedule 28.11.2019
comment
Я нисколько не удивлен, что sqrt().inv() самый быстрый. В наши дни и sqrt, и inv представляют собой отдельные инструкции, которые выполняются довольно быстро. Doom был написан в те дни, когда было небезопасно предполагать, что вообще было аппаратное обеспечение с плавающей запятой, а такие трансцендентные функции, как sqrt, определенно были программными. +1 за тесты. - person Martin Bonner supports Monica; 28.11.2019
comment
Что меня удивляет, так это то, что transmute явно отличается от to_ и from_bits - я бы ожидал, что они будут эквивалентны инструкциям даже до оптимизации. - person trentcl; 29.11.2019
comment
@MartinBonnersupportsMonica Следует отметить, что это было верно для x86 в течение долгого времени, но на ARM нам действительно нужно проверять версии. Cortex-M4 имеет полную аппаратную поддержку только для чисел с плавающей запятой одинарной точности, тогда как Cortex-M7 обеспечивает одинарную и двойную точность. Так что современный планшет должен быть в порядке, а вот мелкие встроенные элементы - нет. Такого рода уловки полезно запомнить тем из нас, кому все еще нужно считать ОЗУ в байтах. :) - person Graham; 29.11.2019
comment
@Graham Отличный момент! (Тест почти наверняка будет проводиться на платформе с аппаратным FP, но, вероятно, существует больше установленных процессоров, которые не выполняют FP, чем они.) Прошло время, так как мне трудно заботиться о пространстве - а затем я считал (16-битными) словами. - person Martin Bonner supports Monica; 29.11.2019
comment
@Graham: эта конкретная реализация также ограничена числами с плавающей запятой одинарной точности, поэтому даже на M4 она имеет ограниченное значение. - person MSalters; 29.11.2019
comment
@MSalters, если не требуется диапазон чисел с плавающей запятой двойной точности, это алгоритм аппроксимации, поэтому я решил, что одинарная точность, вероятно, имеет больше смысла. - person edwardw; 29.11.2019
comment
@MartinBonner Все FPU x86, начиная с оригинального 8087, поддерживали аппаратное обеспечение fsqrt, а также трансцендентные функции вроде fsin и т. д. Однако они были микрокодированы и намного медленнее, чем InvSqrt. Сегодня вы должны использовать rsqrtss или rsqrtps, если вам нужен быстрый приблизительный обратный квадратный корень. - person benrg; 30.11.2019
comment
@MartinBonner (Кроме того, это не имеет значения, но sqrt не является трансцендентной функцией.) - person benrg; 30.11.2019
comment
@MartinBonner: Любой аппаратный FPU, поддерживающий разделение, обычно также поддерживает sqrt. Основные операции IEEE (+ - * / sqrt) необходимы для получения правильно округленного результата; вот почему SSE предоставляет все эти операции, но не exp, sin или что-то еще. Фактически, div и sqrt обычно выполняются в одном и том же исполнительном модуле, спроектированном аналогичным образом. См. HW сведения об модуле div / sqrt. В любом случае, они все еще не быстры по сравнению с умножением, особенно по задержке. - person Peter Cordes; 30.11.2019
comment
@benrg: fsqrt / sqrtps не микрокодированы; это операции IEEE Basic, напрямую поддерживаемые HW как единый uop, например rsqrtps. См. Мой предыдущий комментарий. - person Peter Cordes; 30.11.2019
comment
@edwardw был ли ваш микробенчмарк для измерения задержки или пропускной способности? Современные блоки div / sqrt несколько конвейерны, но не полностью. (например, задержка в Skylake ~ 12 циклов, один на каждые 3 цикла для независимых входов float или несколько хуже для double). Но multiple более конвейерный, например Задержка 4 цикла и обратная пропускная способность 0,5 цикла (то есть 8 в полете одновременно). Время 1,6 нс составляет всего ~ 6 тактовых циклов на частоте 4 ГГц, что неоправданно мало; вы, вероятно, измеряете среднюю пропускную способность sqrt + div по повторяющемуся циклу. Задержка для HW sqrt + div тоже может быть лучше. - person Peter Cordes; 30.11.2019
comment
В любом случае, Skylake имеет значительно лучшую конвейерную обработку для div / sqrt, чем предыдущие uarches. См. Деление с плавающей запятой против умножения с плавающей запятой для получения некоторых выдержек из таблицы Агнера Фога. Если вы не выполняете много другой работы в цикле, поэтому sqrt + div является узким местом, вы можете использовать HW-быстрый обратный sqrt (вместо quake hack) + итерацию Newton. Особенно с FMA, которая хороша для пропускной способности, если не для задержки. Быстрый векторизованный rsqrt и обратный с SSE / AVX в зависимости от точности - person Peter Cordes; 30.11.2019
comment
@PeterCordes - это действительно пропускная способность. Я буквально завернул вызов функции в замыкание и передал его итератору, как того требует criterion API, который, в свою очередь, позаботился о разогреве, измерении и статистическом анализе. Это мандалорский путь ^ H ^ H ^ H ^ H criterion. - person edwardw; 01.12.2019

arrow_upward
10
arrow_downward

Вы можете использовать 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
}

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

person Deedee Megadoodoo    schedule 28.11.2019
comment
В unsafe нет ничего плохого, но есть способ сделать это без явного небезопасного блока, поэтому я предлагаю переписать этот ответ, используя _ 1_ и _ 2_. Он также несет в себе цель, в отличие от трансмутации, которую большинство людей, вероятно, считают магией. - person Sahsahae; 28.11.2019
comment
@Sahsahae Я только что отправил ответ, используя две упомянутые вами функции :) И я согласен, что unsafe здесь следует избегать, поскольку в этом нет необходимости. - person Lukas Kalbertodt; 28.11.2019