Pull to refresh

Comments 188

> 0x5f3759df (это название статьи)
Содержание — супер.
… А в плане SEO название статьи — предельно неудачное: ) Когда клоны наедятся — стоит переименовать. Может, так и было задумано?
Это название оригинальной статьи. Ради какого-то там SEO портить авторскую задумку — несправедливо.
Должен отметить вы сломали TMFeed image
P.S. Привет из SEnC :)
Но SEO в данном случае я видел «в благородном смысле» — ради пользования. Читатель ни за что не запомнит название статьи и не введёт в поиск. Значит, возможность SEO для него пропала. Если бы статья называлась "Быстрый алгоритм квадратного корня и 0x5f3759df" — не пострадали бы ни поиск, ни задумка автора.

(Или, хуже вариант — «Зачем 0x5f3759df для квадратного корня»)
Читатель ни за что не запомнит название статьи и не введёт в поиск.

Поисковый запрос вроде «магическая константа обратный корень» вполне себе нормально гуглится.
Со всеми этими выкладками мы приблизились к пониманию того, что в коде в начале статьи нет никакой “магии”, “фокусов”, “интуиции”, “подбора” и прочих грязных хаков, а есть лишь чистая математика, во всей своей первозданной красоте. Для меня главным выводом из этой истории стала новость о том, что преобразование float в int и наоборот путем переинтерпретации одного и того же набора бит – это не ошибка программиста и не “хак”, а вполне себе иногда разумная операция.
Ошибаетесь. В этом коде есть чудовищный и катастрофический хак, который просто обязательно надо было упомянуть. А именно, нарушение правила strict aliasing-а при разадресации указателя. Все беды в коде возникают оттого, что кто-то проникся «элегантным решением», не разобравшись во всех деталях, а потом начинает применять его к месту и не к месту.

Эту и другие темы я разбирал в своем докладе на конференции C++ Russia. Кому интересно, можете заглянуть.
а что поделать, тут пользуются аппаратным ускорением логарифмирования
Вопрос не в том, пользуются или не пользуются. Любое трюкачество должно быть уместно. А если его таки применяют, должен быть большой жирный комментарий, почему оно здесь и какие условия нужны для того чтобы все работало, включая фазы луны и прочее.

Код в статье будет работать, если компилятору передали ключ -fno-strict-aliasing, либо если эта опция умолчательная для компилятора.

Если через n лет проект начнут портировать на другую архитектуру (или переедут на другой компилятор) и забудут про этот маленький трюк, последствия могут оказаться самыми печальными.
Чтобы печальных последствий не было, достаточно 1-2 юнит-тестов. Если компилятор вдруг решит поменять размер int или стрюкачить с алиасингом — это произойдёт на этапе компиляции и легко отловится тестом.
К сожалению, этого не достаточно. Программа, содержащая неопределенное поведение — некорректна. Не существует никаких критериев, позволяющих доказать корректность программы, кроме формального анализа и ручного разбора ассемблера в данном конкретном случае. Даже компилятор сам не знает, есть ли там UB или нет.

Поменяв одно выражение в любой части кода вы можете повлиять на инлайнинг и как следствие на порядок выполнения проходов компилятора. Вам не доводилось видеть ситуацию, когда в программе выполняются обе ветви условия одновременно? Это один из возможных вариантов реакции компилятора на UB.

Выполняя оптимизации, компилятор исходит из предположения, что неопределенного поведения в программе нет. Соответственно абсолютно все действия, выполняемые компилятором, заботятся только об обеспечении наблюдаемого поведения. Код из статьи таковым не является.

Компилятор может изменять структуру программы так, чтобы сохранить наблюдаемое поведение и упростить себе жизнь. Под нож могут попадать выражения, условия, и даже целые функции. Представьте, что одним из таких условий был assert в вашем юнит тесте.

Резюмируя, юнит тест в таком случае даст не больше, чем комментарий «мамой клянусь, оно работает!». Даже хуже, проход тестов может вселить только ложную уверенность в том, что все хорошо.

assert в тесте как раз защитить можно — надо только модуль с тестами компилировать без оптимизации, и линковать с остальной программой динамически (или статически, но в два этапа — второй этап линковки уже без оптимизации).


Хуже другое — UB может никак не проявляться в тестах, но проявиться уже в реальной программе. Вот тут и правда никак не закрыться.

Идя таким путем, вы вступаете на очень зыбкую почву. Неизвестно еще, что будет хуже.

Если возникает объективная потребность в небезопасных манипуляциях, их нужно изолировать, вплоть до ручного написания функции на ассемблере.

Другое дело, что вся эта машинерия может негативно сказаться на производительности и потерять в итоге смысл.

Современные компиляторы достаточно умные, чтобы использовать весь набор доступных инструкций, включая уже названную rsqrtps. Главное просто не мешать и не считать себя умнее компилятора.

Так я потому и пишу: хуже другое — UB может никак не проявляться в тестах, но проявиться уже в реальной программе. Вот тут и правда никак не закрыться.

Я видел. Это были мысли вслух, а не наезд :)
а что поделать, тут пользуются аппаратным ускорением логарифмирования

можно написать так:


float FastInvSqrt(float x) {
  float xhalf = 0.5f * x;
  int32_t i;
  memcpy(&i, &x, sizeof(i));  // представим биты float в виде целого числа
  i = 0x5f3759df - (i >> 1);  // какого черта здесь происходит ?
  memcpy(&x, &i, sizeof(x));
  x = x*(1.5f-(xhalf*x*x));
  return x;
}

современный компилятор заинлайнит memcpy с фиксированным размером,
и в результате будет тот же код, что и задумывался,
зато ни UB, ни прочих непрятных эффектов

Вы нашли правильное место, но потеряли весь контекст на том сайте.

Проблема в том, что вот это вот (код из примера с вашего сайта):
auto p = reinterpret_cast<unsigned char*>(&x);
это — отлично, правильно, великолепно.

А вот это вот (минимальная модификация, напрашивающаяся для кода в статье):
auto p = reinterpret_cast<unsigned int*>(&x);
это — кровь, кишки, расчленёнка, все горько плачут над отформатированным винчестером.

Ибо для char (или unsigned char) в разделе 6.5.7 An object shall have its stored value accessed only by an lvalue expression that has one of the following types сделано специальное исключение, но ни для intа, ни для floatа — исключений нет.

А дальше, если мы хотим что-то полезное сделать — то нам нужно либо использовать memcpy, либо читать по одному байту и собирать в int (надеясь на компилятор), либо ещё как-нибудь…
UFO just landed and posted this here
что конкретно не так в коде функции в статье?
Код в статье нарушает стандарт С++ и провоцирует UB. Дальше можно уже не разбираться, ибо на этом месте заканчиваются все гарантии.

Типы примитивные, массивов нет, компилятор поместит i и x в регистры или на стеке… и что не так?
Что именно компилятор будет делать со значениями зависит от сотни-другой параметров, включая оценку весов при инлайнинге и register pressure в данной точке. Сама попытка объяснить, что сделает компилятор из общих соображений и «здравого смысла» уже обречена на провал.

Ради интереса скомпилировал тест с -O3 и -O0 с включенным и выключенным алиасингом на g++ и clang++ (std=c++14) — результат выполнения одинаков.
То что он одинаков здесь и сейчас не значит, что он будет одинаков всегда. В этом и засада. Пример того, что может случиться, озвучивается в докладе в районе 28 и 29 слайдов.

К сожалению, самое страшное, что может сделать программист в непонятной ситуации, — «это проверить на практике и убедиться, что все хорошо».
UFO just landed and posted this here

Код основан на трюке с обращением к float через указатель на int: int i = *(int*)&x;, x = *(float*)&i;. Вот это и есть "не так".


Подвох в том, что эта операция — по стандарту UB (нарушение strict aliasing rule). Некоторые компиляторы имеют ключи для отключения этого правила или вообще не учитывают его при оптимизациях, потому код и работает. Но такой код является непереносимым не только на другие платформы, но даже на новую версию компилятора.

Это всё так, но вернёмся к тому же Quake. Это коммерческий продукт. Игра, в которую должно быть можно играть. В нём можно было бы использовать «правильный код без UB» — но потерять на этом пару fps (а может и пару десятков). В результате код был бы с точки зрения компилятора правильным, но бизнес-задачу не решал бы. Люди бы в эту игру не играли, деньги бы за неё не платили. Ну или можно было бы пойти на UB, рискнуть непереносимостью — и создать тот Quake, которым он был создан.

Так и для кого же мы, программисты, пишем код — для компилятора, или всё-таки для других людей?

Зачем вы читаете у меня между строк то, чего я не писал? Я отвечал на конкретный вопрос:


Материалы по ссылкам прочёл, и даже понимание пришло в целом, но не могу понять, что конкретно не так в коде функции в статье?
Я это понял и именно на это и отвечал. Вы описали своё понимание «что не так в коде» с точки зрения того, что код не соответствует стандарту языка и, по Вашему мнению, именно это «не так». Я же попытался объяснить своё виденье вопроса: соответствие стандарту не есть необходимым условиям для кода. Необходимым условием есть практическая польза этого кода.

От того, что код решает свою бизнес-задачу, он не уходит из категории "хитрых хаков", нуждающихся в дополнительной поддержке.

От того, что код остаётся в категории «хитрых хаков» и нуждается в дополнительной поддержке он не перестаёт решать свою бизнес-задачу.
:)
Стоит добавить, что решает на конкретном компиляторе, при конкретных его ключах и на конкретной архитектуре. Если вдруг забудут про этот хак (или про другой аналогичный) и переедут на другой компилятор/платформу он уже задачу решать не будет.
Стоит отметить, процессоры тоже на месте не стоят, и сейчас очень не факт, что он решает свою бизнес-задачу.
Стоит отметить, что «хитрый хак» использован в игре — т.е. в коде, который поддерживаться будет как максимум пару лет.
> Ну или можно было бы пойти на UB, рискнуть непереносимостью — и создать тот Quake, которым он был создан.

Авторы Quake, я предполагаю, вообще не озадачивались подобным выбором, зная свою целевую обстановку. Агрессивные оптимизации, которые могли бы испортить такой код, тогда не применялись ни одним компилятором, доступным для Windows; тогда только-только начинала появляться теория таких оптимизаций. Например, книга Мучника — одно из классических изданий такого рода — это 1997-й, но реальное появление чего-то такого в доступных компиляторах это уже 2000-е. Формализация проблемы это тоже стандарт C99, не раньше(?)
Зная, что тут проблемы нет, авторы Quake могли применять подобные фишки в полный рост, не опасаясь последствий.
Кроме того, сам факт бинарной поставки мог позволить доточить по вкусу, где нужно, даже после компилятора :)

А вот для современного применения, с этим образцовым кодом уже проблемы. И если то, что коллега Halt взвился на этот код по нынешним меркам, было просто реакцией на персональную разновидность красной тряпки (у кого что болит...), и за пределами основной темы статьи, то Ваше предложение проверить на юнит-тестах, увы, совсем некорректно.
Добавлю также, что раньше игры довольно плохо поддерживались после выпуска, т.к. это было довольно бессмысленно. Так что все баги в записанном на болванку бинарнике внезапно превращались в фичи. После если повезёт, м.б. два — три багафикс патча выпустят, которые никто не заметит.
Также в те времена в игровой индустрии госпотствовал Windows, по сему речи о портировании в основном не возникало. Стоит отметить, и сейчас довольно редко возникает.
Как следствие: неизменность компилятора вплоть до версии, библиотек, платформ и т.д.
Также в те времена в игровой индустрии госпотствовал Windows, по сему речи о портировании в основном не возникало. Стоит отметить, и сейчас довольно редко возникает.

Неправда: игровые консоли уже тогда были популярны не меньше, чем игры на PC. Ну и в частности — Quake III была выпущена, кроме Windows, под Dreamcast, PlayStation и Xbox.
Подвох в том, что по стандарту это UB. Несмотря на то что указатель скастили вот только что — на указателях этих есть пометка, что они гарантированно не пересекаются (когда они гарантированно пересекаются), нарушается правило strict aliasing стандарта. Кстати в статье не указан еще один метод: пометить указатели как могущие пересекаться __attribute__((__may_alias__)), но это уже расширение.

Как минимум, автор данного метода заложился на то, что представление плавающих чисел будет именно таким:


Для начала очевидно, что этот трюк работает лишь на платформах, где и целые, и плавающие числа имеют размер = 32 бита; это уже сказано в комментариях. Однако, тут надо добавить, что ещё и порядок битов д.б. именно таким, а не другим: например, на произвольном бстрактном процессоре экспонента м.б. справа, а и мантисса — слева. Или экспонента может занимать не восемь бит, а как-то иначе.


По хорошему этот код надо предаварить проверкой, выполняемой при компиляции — например, что битовое представление, изображённое на картинке "0,01111100,01000000000000000000000", даёт именно это плавающее число "0.15625". И для разных других чисел — тоже.


Да и представление целых чисел бывает разное: little-endian и big-endian.

> Как минимум, автор данного метода заложился на то, что представление плавающих чисел будет именно таким:

Да, но IEEE754 сейчас чуть более, чем везде. А где есть другие варианты (IBM zSeries, VAX) — IEEE754 тоже есть.

> Да и представление целых чисел бывает разное: little-endian и big-endian.

На платформах с целым LE плавающие тоже в LE, и наоборот.
Исключений пока не видно.

Вы правы в том, что этот трюк платформенно-зависим. Изначально он предполагался исключительно для x86. На другие платформы его надо переносить с проверкой этих условий. Но сейчас Вы вряд ли найдёте платформу общего назначения (не какой-то специальный embedded), где бы это не работало.

Это точно так же, как с дополнительным кодом для отрицательных целых (в английском — жаргонное twoʼs complement) — C ещё допускает иные варианты, а вот Java, C#, Go — нет, они реализуют только такой вариант.

> По хорошему этот код надо предаварить проверкой, выполняемой при компиляции

Думаю, в кроссплатформенных версиях так и делают. Я бы делал. Но шансы, что сработает защита, считаю нулевыми :)
В ряде источников (например, тут) сказано, что обходить проблему strict aliasing можно надёжно через unionʼы — правда, сказано это для C, а не C++. У Вас утверждается, что union только условно безопасен. Кому верить и отчего это зависит?
На том же Stackoverflow есть хороший вопрос после которого понятно, что ничего непонятно.

Очевидно, что верить нужно стандарту, вот только пользуемся мы продуктами интерпретации оного, так что правильного ответа я не знаю.

Но судя по изменениям, движемся мы в сторону разрешения type punning через union.
Кому верить и отчего это зависит?
От компилятора. GCC разрешает подобные трюки, правда с оговорками, clang — нет.

Непонятно, впрочем, зачем всё это, если есть переносимая реализация:

Смотрите сами
$ cat cast.cc
#include <string.h>

int float_to_int(float x) {
  int result;
  memcpy(&result, &x, sizeof(float));
  return result;
}
khim@khim:~/work/android/android-standalone/prebuilts/clang/host/linux-x86/clang-4053586/bin$ ./clang++ -O3 -S cast.cc -o- | c++filt
        .text
        .file   "cast.cc"
        .globl  float_to_int(float)
        .p2align        4, 0x90
        .type   float_to_int(float),@function
float_to_int(float):                      # @float_to_int(float)
        .cfi_startproc
# BB#0:
        movd    %xmm0, %eax
        retq
.Lfunc_end0:
        .size   float_to_int(float), .Lfunc_end0-float_to_int(float)
        .cfi_endproc

        .ident  "Android clang version 5.0.300080  (based on LLVM 5.0.300080)"
        .section        ".note.GNU-stack","",@progbits


Как видите никакого вызова memcpy нигде не остаётся и всё «честно».
А самое главное, любой читающий код поймет, что тут происходит, даже без пояснений. Главное не подпускать к коду тех, кто вопит про неэффективность и упущенную прибыль возможность оптимизации :)

Хотя, конечно, по-прежнему остается проблема разных размеров. Кэп подсказывает, что надо добавить ассерт, плюс использовать типы с фиксированным размером, вроде uint32_t.
Главное не подпускать к коду тех, кто вопит про неэффективность и упущенную прибыль возможность оптимизации :)
Нет никакой «упущенной прибыли». Даже самый последний MSVC научился с memcpy работать несколько лет назад. Покажите им сгененированный ассеблерный код — и они умолкнут.

Кэп подсказывает, что надо добавить ассерт, плюс использовать типы с фиксированным размером, вроде uint32_t.
Не бывает флоатов фиксированного размера в стандарте, увы. А вот static_assert — это да.
UFO just landed and posted this here
IEEE754 binary32 AKA single — это не фиксированный размер в 32 бита?
Фиксированный — только не в C. ISO/IEC TS 18661-3:2015, увы, мало какие компиляторы поддерживают. Хотя GCC и Clang вроде умеют __float32
Стандарты C, C++ описывают некий общий вариант, в котором жёстко не прописаны ни IEEE754, ни двоичная иерархия размеров, ни даже представление чисел в дополнительном коде. Вот в C99:

For signed integer types, the bits of the object representation shall be divided into three groups: value bits, padding bits, and the sign bit. There need not be any padding bits; there shall be exactly one sign bit. Each bit that is a value bit shall have the same value as the same bit in the object representation of the corresponding unsigned type (if there are
M value bits in the signed type and N in the unsigned type, then M ≤ N ). If the sign bit is zero, it shall not affect the resulting value. If the sign bit is one, the value shall be modified in one of the following ways:
— the corresponding value with sign bit 0 is negated (sign and magnitude);
— the sign bit has the value −(2^N) (two’s complement);
— the sign bit has the value −(2^(N−1)) (one’s complement).


Читая это, мне реально интересно, где они находят машины, у которых для данного описания M <= N, или где есть что-то кроме дополнительного кода для представления отрицательных. Возможно, Вы, если описываете себя как embedded специалист, назовёте пару имён?

Видимо, для тех, кто устал от этого всего, в Java, C#, Go сделано (цитирую по Go)

> The value of an n-bit integer is n bits wide and represented using two's complement arithmetic.
(никаких тебе «знаковый может быть у́же беззнакового»)
> In a function call, the function value and arguments are evaluated in the usual order.
(тут есть отсылка на вычисление слева направо; никаких «в любом порядке, как нам удобно», как в C)
> Shifts behave as if the left operand is shifted n times by 1 for a shift count of n.
(никаких «если сдвиг на величину большую либо равную ширине сдвигаемого, результат не определён»)
> when evaluating the operands of an expression, assignment, or return statement, all function calls, method calls, and communication operations are evaluated in lexical left-to-right order.
(слева направо, считаем, в большинстве случаев)
> For unsigned integer values, the operations +, -, *, and << are computed modulo 2^n, where n is the bit width of the unsigned integer's type.
(тут так же, как в C)
> For signed integers, the operations +, -, *, and << may legally overflow and the resulting value exists and is deterministically defined by the signed integer representation, the operation, and its operands. No exception is raised as a result of overflow. A compiler may not optimize code under the assumption that overflow does not occur. For instance, it may not assume that x < x + 1 is always true.
(принципиальное отличие от C; как если бы gcc вызван с -fwrapv)

возможно, результирующий код и не так эффективен, но зато надёжен в плане отсутствия эффектов, неожиданно вышибающих почву из-под ног.

И struct aliasing rule там тоже нет, насколько я знаю.

Если бы можно было гарантированно и по стандарту получать подобные ограничения в C, этим бы пользовалось >90% его пользователей, и им бы было несущественно, что от этого они теряют 10-20-30% скорости. Но сейчас C развивается немножечко так в другую сторону :(
Уточнение про ЯВУ — не во всех все такие правила (например, для сдвигов нет гарантий «как N раз по 1 биту»; это специфика Go). Речь про общее направление и режим по умолчанию. Но оно достаточно показательно. Обычно мы ожидаем меньше «умничания» от компилятора для языка более низкого уровня, и больше — для высокого (вплоть до изменения алгоритма), здесь же наоборот.
UFO just landed and posted this here
> Так что если у меня написано i = *(uint32 *)&f то оно проверено и это никогда уже не поменяется.

Проблема в том, что тут стандарты «плывут» и компиляторы вслед за ними. В сверхбольшом проекте каждая смена версии компилятора выявляет пачки проблемных мест. Вот в этом случае GCC 4.8 внёс по сравнению с 4.7 сильно больше агрессивности в том, что до того в принципе уже считалось существующим в нём и дало неожиданность на ровном месте.
И если ваше «никогда не поменяется» включает в себя фиксацию версии компилятора — ok, но тогда вы иначе ограничены в развитии. Если нет — будьте готовы к сюрпризам.

> Но вообще интересно, как «по-правильному» передавать информацию в кодограммах между двумя машинами с заранее неизвестными компиляторами.

Если речь о взаимодействии протоколами, то тут есть несколько универсальных, но весьма громоздких подходов.
Например, стандарт обещает отсутствие проблем с алиасингом при доступе через char*. Если вам надо прочитать uint32_t, то гарантированно остаётся только что-то вроде

unsigned char *p;
return ((uint32_t)p[3]<<24) + ((uint32_t)p[2]<<16) +
  ((uint32_t)p[1]<<8) + (uint32_t)p[0];

(заметим, что конверсию в uint32_t тут нужно делать явно! по умолчанию оно стремится расширить в signed int, если влезают все значения, а не в unsigned int — ещё одна засада...)

Сможет компилятор понять, что это можно на x86/ARM/MIPS/etc. превратить в mov в регистр — не знаю заранее. В моих тестах сейчас ни GCC, ни Clang этого не сумели. Но можно хотя бы зависимость от платформы вписать в #ifdef.

> А ещё интереснее на архитектурах, где на одной стороне байт 32-битный.

А данные из сети поступают на такие машины как? Я думаю, если интерфейс соответствует BSD sockets, то в один машинный 32-битный байт укладывается 1 сетевой октет (потеря 75% места, а что поделаешь). Тогда код взятия 2 абзацами выше — будет работать, и аналогичная обратная раскладка — тоже. А вот прямая запись в unsigned-поле в предположении, что в нём 4 байта — уже не сработает.
UFO just landed and posted this here
В моих тестах сейчас ни GCC, ни Clang этого не сумели.
Плохо тестировали. Clang так умеет, а для GCC нужно чуть-чуть по другому написать. Пипхольные оптимизации — они такие, да. Шаг влево, шаг вправо — и всё развалилось…

А ещё интереснее на архитектурах, где на одной стороне байт 32-битный.
Для таких архитектур придётся отдельный код писать — а куда деваться? Напишите «if constexpr (CHAR_BITS == 32)» и засунете туда свой вариант. Или, на старых версиях — шаблонный тип от CHAR_BITS. Или, уж совсем по старинке, ifdef или memcpy…
> Плохо тестировали. Clang так умеет, а для GCC нужно чуть-чуть по другому написать.
У меня основной GCC пока 4.8, а Clang — района 3.8. Даже 4.0 эту свёртку ещё не умеет. Но за идею про "|" спасибо, тут было неочевидно, что это ещё может на него повлиять :)
Но за идею про "|" спасибо, тут было неочевидно, что это ещё может на него повлиять :)
Типичная ситуация когда вы думаете о компиляторе как о живом, почти разумном, существе. Который пытается понять смысл вашей программы и сделать её лучше. Чего и в помине нету.

Битовые манипуляции — по определению работают с полями, кусками машинных слов. И кроме простейшего варианта, который тут случился могут быть много других — и там есть где развернуться оптимизации. Которая, в вырожденном случае, превратится в ничто. А когда вы используете плюс то из-за переноса могут быть разные странные эффекты. Оптимизировать такое — гораздо сложнее, потому разработчики GCC и не пытались.

P.S. Не удивлюсь если clang только вот эту вот связку распознаёт и больше — почти ничего. «В пользу бедных», так сказать…
> Типичная ситуация когда вы думаете о компиляторе как о живом, почти разумном, существе. Который пытается понять смысл вашей программы и сделать её лучше. Чего и в помине нету.

Мысленное одушевление столь сложной сущности это вообще-то неизбежно для человека. И смысл он вполне может понять, в простых случаях. Но я не это имел в виду.

> А когда вы используете плюс то из-за переноса могут быть разные странные эффекты.

Вот именно.
UFO just landed and posted this here
Ну вот я считаю, по своему прочтению доступных данных (стандарты и final drafts), что решает, хоть и чуть дороже (за счёт обязательного размещения в оперативной памяти, регистры не годятся). А коллеги Halt и khim этот метод не отвергают, но просто игнорируют; зато настаивают, что memcpy в современном мире сам по себе достаточен, и даёт даже работу через регистры, где компилятор это сумел понять. Я понимаю — там, где вариант с memcpy работает, он эффективнее. Но прямого ответа нет…
UFO just landed and posted this here
А коллеги Halt и khim этот метод не отвергают, но просто игнорируют
Лечить головную боль с помощью гильотины? Увольте. Средство надёжное, спору нет, просто… несколько радикальное.

Я понимаю — там, где вариант с memcpy работает, он эффективнее.
Он, собственно, работает во всех рапространённых компиляторах и является единственным переносимым способом сделать то, что вы хотите (можно, конечно, реализовать memcpy руками черещ указатели на char, но это, скорее всего, будет менее эффективно).

Вот тут есть дискуссия, тут — ещё одна, здесь — формальный proposal, тут — шаблончик.

Это уже обсуждалось 100 раз.

Но прямого ответа нет…
Прямого ответа не будет пока bit_cast в C++ не добавят…
> Он, собственно, работает во всех рапространённых компиляторах и является единственным переносимым способом сделать то, что вы хотите

И ещё раз. Тут Вы пишете, что это _единственный_ (сами выделили) переносимый. А выше пишете

> Средство надёжное, спору нет, просто… несколько радикальное.

Если надёжное — то тоже работает или переносимое? Или нет? Я очень прошу дать прямой ответ, без всяких отклонений на bit_cast и тому подобное.
Если надёжное — то тоже работает или переносимое? Или нет? Я очень прошу дать прямой ответ, без всяких отклонений на bit_cast и тому подобное.
Ответ простой: конечно же volatile работает. С чего б ему не работать? Однако в качестве побочного эффекта — возможно замедление программы. Не на проценты — в разы. В 10-20 раз если сильно не повезёт.

Так как подобные трюки обычно проделываются с целью ускорения программы, то volatile там не подходит от слова «совсем».
> Ответ простой: конечно же volatile работает.

Спасибо за прямой ответ. Жаль, что его буквально пришлось вытягивать.

> Однако в качестве побочного эффекта — возможно замедление программы. Не на проценты — в разы.

Какое будет замедление в конкретном случае — позвольте судить именно по этому случаю, а не предсказывать, не зная, что где и как работает.
Меня интересовало, подходит ли этот метод на пусть медленную, но гарантированную базу. А уже имея её в запасе — можно выбирать и более эффективные методы.
Меня интересовало, подходит ли этот метод на пусть медленную, но гарантированную базу.
Не подходит. Volatile — это ошибка в дизайне C. Его невозможно реализовать корректно и эффективно и обсуждать скорость работы программы, использующей volatile ни один разработчик компилятора с вами не будет. Ибо бессмысленно.

А уже имея её в запасе — можно выбирать и более эффективные методы.
Volatile — это не база. Volatile — это ошибка.

Даже в ядре, при общении с железом (а это то, ради чего придуман volatile) его рекомендуют избегать, а уж вне ядра — тем более.

И уж конечно никакому нормальному человеку не придёт в голову его использовать там, где важна скорость ибо единственная цель volatile — это запрет на оптимизации.

У нас в кодовой базе на несколько миллионов строк volatile встречается примерно в тысяче мест и, разумеется, только и исключительно в многопоточных структурах. Нигде больше volatile никому и никогда даже в голову использовать не пришло.
> Volatile — это ошибка в дизайне C. Его невозможно реализовать корректно и эффективно

Его возможно реализовать корректно: данное в памяти надо читать, не обращая внимание на все предыдущие идеи о том, что там должно быть, и на последующие о том, что туда записано.
Его возможно реализовать эффективно: надо всего лишь прочитать и записать. Ничего больше.
«Ошибки» здесь нет.

> Даже в ядре, при общении с железом (а это то, ради чего придуман volatile) его рекомендуют избегать

Я в курсе этого сообщения. Там возражения против него рассматриваются исключительно в свете синхронизации общения с другими параллельно исполняющимися задачами. Для этого он действительно не нужен.
Но там не рассматривается тема проблем неожиданного алиасинга.

> там, где важна скорость

Вы зациклились на понятии скорости и ничего другого не желаете слышать. Я же всё время говорю об обеспечении _корректности_ операций такого «сомнительного» типа. И только после обеспечения корректности, считаю, можно начинать рассматривать тему скорости.
Давайте, пока Вы не согласитесь с последним абзацем, не будете мне отвечать:) потому что иначе мы впадём в вечный цикл. Молчание будет достаточным индикатором несогласия.
Его возможно реализовать эффективно: надо всего лишь прочитать и записать.
Это, к сожалению, не так: нужно также отказаться от возможности переносить обращения к другим переменным через точку в программе, где вы вставили обращение к volatile.

Но там не рассматривается тема проблем неожиданного алиасинга.
Угу — там указывается на то, что volatile плохо решает даже те задачи, для которых он изначально предназначен. Задачи же, для которых он не был предназначен он решает ещё хуже.

Я же всё время говорю об обеспечении _корректности_ операций такого «сомнительного» типа.
Для этого существует char* и memcpy.

Вы зациклились на понятии скорости и ничего другого не желаете слышать.
Ok. Принято. Рассмотрим отличный переносимый вариант:
int x;
float y;
static_assert(sizeof(x) == sizeof(y));
...
FILE* f=tmpfile();
fwrite(&x, sizeof(x), 1, f);
fseek(f, 0, SEEK_SET);
fread(&y, sizeof(y), 1, f);
fclose(f);
Чем этот вариант хуже, чем volatile? В его корректности, как бы, тоже никто не сомневается. Почему вы не хотите его в качестве базового рассмотреть?
> отказаться от возможности переносить обращения к другим переменным через точку в программе, где вы вставили обращение к volatile.

Прошу обоснования цитатами из стандарта.

> Чем этот вариант хуже, чем volatile?

Завязкой на сущности, которые заметно выходят за пределы собственно обстановки рантайма C.
Прошу обоснования цитатами из стандарта.
Ась?
5.1.2.3 Program execution

In the abstract machine, all expressions are evaluated as specified by the semantics. An actual implementation need not evaluate part of an expression if it can deduce that its value is not used and that no needed side effects are produced (including any caused by calling a function or accessing a volatile object).


Volatile разрушает оптимизации так же, как вызов функции без исходников (за исключением известных компилятору типа memcpy) — а вы его во внутренний цикл пихать хотите…

Завязкой на сущности, которые заметно выходят за пределы собственно обстановки рантайма C.
Извините, но все эти сущности являются частью языка C. Такой же неотьемлемой, как volatile или memcpy.

Если вы работаете с чем-то, что не реализуется язык C, как он описан в стандарте — то у вас всё, что угодно может быть, конечно…
> 5.1.2.3 Program execution
Ну, не этот пункт, но я нашёл нужное совсем рядом. Спасибо. Стандарты такого рода, увы, очень неудобны для чтения.

> Volatile разрушает оптимизации так же, как вызов функции без исходников (за исключением известных компилятору типа memcpy) — а вы его во внутренний цикл пихать хотите…

Простите, про внутренний цикл и т.п. — это Ваши необоснованные домыслы. Я нигде не предлагал использовать его в критических по скорости местах; на практике я бы там скорее применил ассемблерные вставки (типа movd на x86), не надеясь на оптимизации всяких memcpy(). В остальном я всё время вёл речь про то, как сделать относительно быстро при надёжной гарантии (и без закапывания в особенности каждого компилятора). Ответы я в итоге получил, нового уже, скорее всего, тут не будет совсем.

> Извините, но все эти сущности являются частью языка C. Такой же неотьемлемой, как volatile или memcpy.

Не во freestanding случае. Я его учитываю как важный.

В общем, спасибо, можно закрывать ветку.
Я нигде не предлагал использовать его в критических по скорости местах; на практике я бы там скорее применил ассемблерные вставки (типа movd на x86), не надеясь на оптимизации всяких memcpy().
И получили бы замедление раз в два (хорошего если в 10) через несколько лет при смене архитектуры процессора. Или если вы считаете что на каждый критичный участок нужно посадить программиста, который будет следить за изменениями в мире процессоров и постоянно подкручивать ваш код?
> И получили бы замедление раз в два (хорошего если в 10) через несколько лет при смене архитектуры процессора.

На movd с аналогом? Это что же должно случиться, чтобы тут возникло замедление в 10 раз?

> вы считаете что на каждый критичный участок нужно посадить программиста, который будет следить за изменениями в мире процессоров и постоянно подкручивать ваш код?

Я не «считаю», я вижу, что оно сейчас именно так и происходит. Каждый год от смены чего-то «утекают» как минимум несколько процентов от результата. Иногда — в разы (как с приходом GPU).
На movd с аналогом? Это что же должно случиться, чтобы тут возникло замедление в 10 раз?
AVX512, например. На RSQRTPS компилятор этот код заменить не сможет, но если функция выполняется в цикле — вполне может её векторизовать. А прямое использование movd, как и volatile, приведут к тому, что у него этой возможности не будет.

Я не «считаю», я вижу, что оно сейчас именно так и происходит.
Вы видите, я — нет. Может быть мы просто над разными вещами работаем, впрочем.

По моим наблюдениями всякие eigenы и TensorFlow (заточенные под то, что их компилятор «правильно» скомпилирует и очень-очень редко применяющие ассемюлерные трюки) используются на несколько порядков чаще, чем вещи типа кодеков (а эти, последние, часто перескакивают стадию «а напишем-ка мы это на ассеблере, чтобы иметь работу до гроба по переписыванию» и реализуются сразу «в железе»).

Иногда — в разы (как с приходом GPU).
GPU — это хороший пример. Ни один GPU не даёт вам возможности «штатно» программировать его на ассемблере. Ни один. А языки, на которых их-таки программируют чем дальше, тем больше становятся похожи на старый добрый C.

Программирование на ассемблере — это временный костыль, но так как развитие идёт по спирали, то эти подходы снова «всплывают» и «умирают» самых разных, зачастую весьма неожиданных местах.

Обсуждаемый пример — в этом смысле черезвычайно показателен: уверяю вас — Кармак отлично знает ассемблер (мы говори о человеке, который начинал свою карьеру с программирования регистров EGA на низком уровне, напоминаю!), однако в данном конкретном случае он-таки предпочёл остаться в рамках C.
> но если функция выполняется в цикле — вполне может её векторизовать

Случай векторизации я тут не рассматривал. Вообще у нас дискуссия как-то странно идёт — одни и те же аргументы применяются несколькими сторонами то в контексте конкретного примера, то в контексте высокоскоростного кода, который надо писать, чтобы компилятор его оптимизировал и векторизовал… половина путаницы из-за этого.

> используются на несколько порядков чаще, чем вещи типа кодеков

Ну вот моя специфика тут пары последних лет это где-то «кодеки». Параллельности немного (ещё и навязанные средства мешают), а вот путаницы в форматах разборе данных ой-ой. Пока что по результатам этой дискуссии пришлось добавить -fno-strict-aliasing на пару исходников, причём я подозреваю, что если я туда повлепливаю memcpy, то меня не поймут кое-какие коллеги.

> GPU — это хороший пример. Ни один GPU не даёт вам возможности «штатно» программировать его на ассемблере.

И снова — я имел в виду вообще переписывать код, а Вы вдруг ассемблер вспомнили :) почему?

> Кармак отлично знает ассемблер
> однако в данном конкретном случае он-таки предпочёл остаться в рамках C.

А вот тут уже очень интересный момент. Если он «предпочёл остаться в рамках C» (за это Вы его хвалите), но применил вместо разрешённого подхода — трюк, про который говорите «никогда-никогда!» (и который не был допустим даже по C89, насколько я нагуглил), о чём это говорит? Он был прав или нет?
Вообще у нас дискуссия как-то странно идёт — одни и те же аргументы применяются несколькими сторонами то в контексте конкретного примера, то в контексте высокоскоростного кода, который надо писать, чтобы компилятор его оптимизировал и векторизовал… половина путаницы из-за этого.
Вы саму статью-то читали? Конкретно — раздел «зачем»?

И снова — я имел в виду вообще переписывать код, а Вы вдруг ассемблер вспомнили :) почему?
Потому что если бы вы писали на C, то код мог бы быть перенесён на GPU «штатно», а вот код на ассемблере пришлось бы понять и переписать на совсем другой язык.

А вот тут уже очень интересный момент. Если он «предпочёл остаться в рамках C» (за это Вы его хвалите), но применил вместо разрешённого подхода — трюк, про который говорите «никогда-никогда!» (и который не был допустим даже по C89, насколько я нагуглил), о чём это говорит? Он был прав или нет?
Он был прав в рамках того мира, в котором он этот код писал. memcpy тогда не был «бесплатным», но и проблемы strict-aliasing'а — это GCC 4.0, 2005й год! Хотя стандарт — да, стандарт об этом говорил и в 1989м, но если вы не использовали компиляторы для суперкомпьютеров типа Крей, то в 1997м у вас проблемы алиасинга не возникали…

я подозреваю, что если я туда повлепливаю memcpy, то меня не поймут кое-какие коллеги.
Скопировать из Хрома bit_cast вместе с комментарием — нет?

То есть в его случае — оставалось только комментарий там написать. И жаль, что этого не было сделано.

Но у нас-то на дворе 2017й!

Это не ошибка языка, а средство для описания работы с отображенными в память устройствами.

Для работы с отображаемыми в память устройствами volatile тоже с большим трудом подходит, так как во многих случаях нужно более точно контролировать, что происходит. На каком-нибудь 8088 при записи простого intа было важно указать — какой из байтов пойдёт в память первым. А в современных системах — нужно всё равно барьеры явно указывать и с кешами с помощью вне-языковых средств работать, так что тольку от volatile — с гулькин нос, а зато соблазну использовать его не по назначению — вагон.

Там, где порядок байт был важен — надо было использовать volatile char, только и всего.

Ну так и копирование из int во float тоже можно побайтно организовать — и это будет работать на всех компиляторах — но люди этого мало почему-то.

А я что, с этим спорю? :-)

volatile решает проблему "неожиданных" оптимизаций, но при этом замедляет программу: без volatile толковый компилятор мог всю работу в регистрах провернуть, а с volatile компилятор вынужден использовать память. В итоге при нуедачном стечении обстоятельств "быстрый" обратный квадратный корень, который использует volatile, может оказаться даже медленнее чем простой 1/fsqrt(x).

UFO just landed and posted this here
> Мне у ADSP нравится: один и тот же регистр может быть float, а может быть int в зависимости от имени, по которому к нему обращаешься. :)

Вот народ из RISC-V разработки пишет:

>> We considered a unified register file for both integer and floating-point values as this simplifies software register allocation and calling conventions, and reduces total user state. However, a split organization increases the total number of registers accessible with a given instruction width, simplifies provision of enough regfile ports for wide superscalar issue, supports decoupled floating-point-unit architectures, and simplifies use of internal floating-point encoding techniques. Compiler support and calling conventions for split register file architectures are well understood, and using dirty bits on floating-point register file state can reduce context-switch overhead.

Аргументы, как для не-крошечных систем, вполне разумны.

> ЗЗЫ: Интел — отстой ;)

Я согласен во многом, но не в этом случае разделения целочисленных и {x,y,z}mm регистров.
UFO just landed and posted this here
> В итоге при нуедачном стечении обстоятельств

А как может возникнуть такое стечение обстоятельств, если содержимое не пойдёт дальше L1?
Затраты на вытеснение чего-то другого в DRAM? Ну так их лучше считать размазанными по всей программе, а локальные переменные всё равно возле вершины стека и наверняка уже давно кэшированы.
Другого источника с ходу не вижу.

Если бы с L1 все было так просто — регистры общего назначения у процессора были бы не нужны.

Регистры нужны не потому, что они не L1, а потому, что они не память — их мало (проще адресовать), и они не ограничены требованием представления в памяти (например, можно переименовывать внутри процессора, или когда память ещё не настроена). К обсуждаемому это не относится.
К обсуждаемому это не относится.
Как это не относится? Вряд ли вы будете так уж сильно оптимизировать вычисление одного корня — скорее всего их у вас будут тучи. А какой-нибудь Cannonlake может обрабатывать 16 int или floatов одновременно. Заставив процессор положить всё на память и запретив, соответсвенно, все оптимизации вы закрыли для компилятора это путь. Раз — и навсегда.

Оставьте вы эту затею с volatile — он не для этого предназанчен.
А как может возникнуть такое стечение обстоятельств, если содержимое не пойдёт дальше L1?
Да легко. Сделав переменную volatile вы, фактически, заставляете компилятор работать в режиме -O0.

Простейший пример:
int foo(int* p, volatile long x) {
    x *= 2;
    return p[x];
}
int bar(int* p, long x) {
    x *= 2;
    return p[x];
}
Сравните код. Вместо одной инструкции вы получаете 4 — неужели вы думаете, что это не отразится на скорости работы программы? И это ещё мы до векторизации не добрались, которую volatile вырубает сразу — и «на корню».

Ну честное слово, volatile — это ну вот совсем не то, что вы хотите увидеть в кусочке программы, который вы всеми правдами и неправдами пытаетесь ускорить!
> Вместо одной инструкции вы получаете 4 — неужели вы думаете, что это не отразится на скорости работы программы?

Я в курсе, спасибо.

> Ну честное слово, volatile — это ну вот совсем не то, что вы хотите увидеть в кусочке программы, который вы всеми правдами и неправдами пытаетесь ускорить!

Позвольте уж каждому на месте решать, насколько и что ему надо ускорять, и является ли вообще скорость тут проблемой. «Tools, not policy.»
А в реальности — скорее всего будут применены результаты сравнения нескольких методов для конкретного таргета и компилятора, даже если лучшим окажется то, что казалось худшим.
Позвольте уж каждому на месте решать, насколько и что ему надо ускорять, и является ли вообще скорость тут проблемой.
Но… зачем? Компиляторы делятся на четыре большие группы:
1. Современные компиляторы — знают о memcpy и могут вызывать проблемы с алиасингом. Volatile не нужен, нужно использовать memcpy.
2. Чуть более старые компиляторы — знают о memcpy, но не используют информацию о типах для алиасинг оптимизации. Volatile не нужен, можно использовать memcpy.
3. Совсем старые компиляторы — вообще имеют мало оптимизаций и проблемы алиасинга от них бесконечно далеко. Volatile не нужен, хотя и memcpy тоже использовать нельзя (хотя можно использовать каламбур типизации в чистом виде).
4. Ещё более древние компиляторы, не соблюдающие стандарт и требующие вставки специальных инструкций при необходимости обращения к памяти то как к целому числу, то как к числу с плавающей точкой. Volatile не нужен — ибо не работает, нужны специальные, зависящие от компилятора инструкции (на старых компиляторах для DOS'а нужно было ручками вставлять wait в нужные участки кода).

А в реальности — скорее всего будут применены результаты сравнения нескольких методов для конкретного таргета и компилятора, даже если лучшим окажется то, что казалось худшим.
Вот найдёте хоть один компилятор где volatile нужен, спасает и работает быстрее, чем memcpy — будет о чём поговорить.

Какой смысл нам, живущим на Шарообразной Земле обсуждать «краевую радугу»? Жили бы мы на диске — другое дело…
Ну, во-первых,
> Современные компиляторы — знают о memcpy и могут вызывать проблемы с алиасингом. Volatile не нужен, нужно использовать memcpy.

это знание отрубается нафиг простым включением freestanding :) после чего надо ещё искать, как ему обратно объяснить, что он должен знать, что такое memcpy.
Но это несколько побочная ветвь, я её упоминаю только для того, что Ваша линейная картина далеко не такая линейная.

> Чуть более старые компиляторы — знают о memcpy, но не используют информацию о типах для алиасинг оптимизации. Volatile не нужен, можно использовать memcpy.

Вот GCC 4.2 — он ещё идёт в текущей FreeBSD как альтернатива clangʼу. Он уже умеет -fstrict-aliasing, но при -O0 он на Вашем варианте исходной функции явно зовёт библиотечную memcpy(), а не встроенный вариант. Значит, уже реальный пример картины наоборот — алиасинг есть, «хорошего» memcpy нет.

Кстати, gcc всех версий у меня показал в варианте с memcpy, что float->int он копирует через регистр, а финальный int->float через стек, например вот для 6.4.0 и -O3:

FastInvSqrt:
        movd    %xmm0, %edx
        movl    $1597463007, %eax
        sarl    %edx
        subl    %edx, %eax
        movl    %eax, -4(%rsp)
        movss   -4(%rsp), %xmm0
        ret


и запись результата такая же, как если бы была через volatile, через union…

> Вот найдёте хоть один компилятор где volatile нужен, спасает и работает быстрее, чем memcpy — будет о чём поговорить.

Вот, считайте, уже нашёл — GCC 4.2 без -O. Если Вы не считаете этот вариант значимым именно из-за неоптимизации — я не могу с этим согласиться, скорость при отладке тоже имеет значение, и получать заведомые тормоза на ровном месте как-то не хочется. И на «не медленнее» тоже примеры уже показал.

Так что как минимум с одной стороны наш мир плоский :)
это знание отрубается нафиг простым включением freestanding :)
С чего вдруг?

после чего надо ещё искать, как ему обратно объяснить, что он должен знать, что такое memcpy
Не нужно. -fbuiltin-memcpy нужно задавать только если кто-то зачем-то сказал -fno-builtin.

Он уже умеет -fstrict-aliasing, но при -O0 он на Вашем варианте исходной функции явно зовёт библиотечную memcpy(), а не встроенный вариант.
Если вы задали -O0, то это значит что скорость работы программы вас не волнует, уж извините.

запись результата такая же, как если бы была через volatile, через union…
Что это доказывает? Что проблемы с оптимизацией существуют в GCC — а не в memcpy. Вот если бы вариант с volatile оказался не таким же, а быстрее — было бы о чём говорить.

Если Вы не считаете этот вариант значимым именно из-за неоптимизации — я не могу с этим согласиться, скорость при отладке тоже имеет значение, и получать заведомые тормоза на ровном месте как-то не хочется.
Если вам нужна скорость — используйте хотя бы -Og. -O0 — это принципиальный отказ от всех и всяческии оптимизаций.

Так что как минимум с одной стороны наш мир плоский :)
Со стороны людей, которые явно выключают оптимизации, а потом жалуются на то, что программы работают межденно. Извините, но я с этим миром стараюсь не пересекаться.
> Не нужно. -fbuiltin-memcpy нужно задавать только если кто-то зачем-то сказал -fno-builtin.

-ffreestanding включает -fno-builtin.

> Если вы задали -O0, то это значит что скорость работы программы вас не волнует, уж извините.

Всё равно волнует, не извиню. Когда что-то просто не оптимизируется — это одно, а когда явно включается противодействие — это уже другое.

> Если вам нужна скорость — используйте хотя бы -Og.

-Og появился только в 4.8 и отсутствует вообще в Clang.

> Вот если бы вариант с volatile оказался не таким же, а быстрее — было бы о чём говорить.

При 4.2 и -O0 он быстрее — за счёт отсутствия тяжёлого вызова библиотечной функции.

> Со стороны людей, которые явно выключают оптимизации, а потом жалуются на то, что программы работают межденно. Извините, но я с этим миром стараюсь не пересекаться.

Когда появится уровень оптимизации, который никак не вредит отладке (действия не переходят границу исходной строки) — я займу такую же позицию. Пока что даже -Og такого не даёт, хоть и заметно приближается.
UFO just landed and posted this here
UFO just landed and posted this here
Вот везде написано «undefined behaviour».
Угу.

А какие есть варианты этого behaviour?
Да какие угодно, блин. Вот прекрасный пример. Модификация в коде сохранилась, а результат… результат возвращать — не нужно. Всё равно он «undefined», какая разница что там будет?

Чтобы behaviour было undefined
… нужно, чтобы вы сделали что-то, что объявлено как undefined behavior в стандарте.

Чтобы behaviour было undefined, требуется наличие более одного варианта, а я что-то не вижу.
Чего вы не видите? Что обращаются к одному и тому же участку памяти то как к intу, то как к floatу?

Стандарт совершенно явно говорит о том, что если к одному и тому же участку памяти обращаются то как к floatу, то как к intу, то результат может быть каким угодно и может не иметь вообще никакого отношения к тем данным, что вы передали на вход. Я как-то статью писал про это и там даже было обьяснено «откуда ноги» у этого запрета растут. Но всё не впрок. «Embedded программисты» по прежнему задают вопрос: «а я что-то не вижу» какие тут могут быть проблемы. Если вы не видите — то это не значит что и компилятор не видит их тоже! И даже если сегодня не видит — завтра увидит! Вы вызвали «undefined behavior» и ваша программа самоуничтожилась? Ну эт нормально, бывает, дело-то житейское… Хорошо винчестер не отформатировала, а только образ в памяти… Никто чинить это не будет.
UFO just landed and posted this here

Там этих ссылок в других ветках накидали уже до черта. Можно было там читать, а не спрашивать уже спрошенное.

UFO just landed and posted this here
> Со всеми этими выкладками мы приблизились к пониманию того, что в коде в начале статьи нет никакой “магии”

Скорее, наоборот, я ещё больше укрепился в том, что математика — это магия.
Третий закон Кларка, однако.

Any sufficiently advanced technology is indistinguishable from magic
нам не хватает лишь уверенности в том, что вычисленное таким образом приближенное значение может быть столь же эффективно улучшено алгоритмом Ньютона

Я тут на прикинул на бумажке. Получается, что для возведения x в степень p, где p не ноль, шаг по Ньютону будет image, где z это аргумент функции (который используется в оригинальном коде для вычисления xhalf). Похоже, что удобно будет только когда -1/p — положительное целое.

В блоке "интересные публикации" название статьи отображается переведённым в десятичную систему счисления. Exploitable?

Вот это да. Ну, заодно и потестили Хабр.

В TMFeed аналогично — «1597463007».

Возможно там при рендеринге есть что-то типа:


if (is_numeric($value)) echo (0 + $value);

Работает на PHP 5.6

Администрация Хабра переименовала статью :)
Ну ок.

Интересно какие ксплойты повылазят, если опубликовать перевод к "0x5f3759df(appendix)"...


ReferenceError: appendix is not defined?

Ради какого-то там SEO портить авторскую задумку — несправедливо.

+1
Про заворот слепых кишок у SEO-шников опосля второй статьи я помолчу лучше...

Хорошо, что на данной платформе sizeof(int) == sizeof(float).
Это надо в начале алгоритма поставить как assert, иначе долго будут искать причину гуру программирования когда неправильно считать будет.
С современными компиляторами этот код вообще использовать нельзя, нарушения стандарта — они такие, да. Компилятор вполне может «соптимизировать» весь этот код в «return 0».
Что то меня терзают сомнения на тему таких умных компиляторов. Может конечно алгоритмы которые в шахматы выигрывают и скомпилируют в return 0, но те что в шашки — точно так не смогут. Да и отладчик с #pragma останется всегда. На крайняк и на ассемблере можно наколбасить. На нём уж точно замучаются оптимизировать.

Good news, everyone.


И хотя более новые версии GCC так (пока) не делают, такое поведение формально не противоречит стандарту, потому что в результате undefined behaviour может получиться все, что угодго (по факту оптимизатор решил, что (unsigned short *)&a и &a не алиасятся, и решил переставить инструкции местами, потому что почему бы и нет).

такое чудо вылазит, начиная с О2
А Clang сразу пишет 4, забив на 5.

Я ненастоящий сварщик, но — volatile проставить и все, не?

Не уверен, станет ли валиднее код, но медленнее — точно.
Доступ к памяти всё равно пройдёт через кэш, потеря будет в пару тактов на операции — это малосущественно.

А ещё — для реализации под конкретный процессор можно вставить операцию местного ассемблера, и volatile не потребуется :)

volatile нужен platform-agnostic, multi-threaded ONLY.
Грубо говоря не даст прочитать некоторому потоку (соптимизированное компилятором) промежуточное значение (явно не указанный программистом v1 = v2) или не даст вырезать (оптимизировать до 0 инструкций) неочевидную для компилятора установку такой переменной.
Грубо говоря, он думает что переменная тут (или ниже) не используется (этим потоком), зато может вполне себе — другим.
Т.е. конкретно тут — ни о чем...


khim
С современными компиляторами этот код вообще использовать нельзя
Компилятор вполне может «соптимизировать» весь этот код в «return 0».

С первым высказыванием согласиться могу, условно. Только не "нельзя" а не комильфо… И assert-ами обвязать желательно…
Со вторым не согласен в корне — с какого это перепуга он его вырежет?
В худшем случае UB (да и UB на касте int/float мне очень сомнительно)...

> Т.е. конкретно тут — ни о чем…

Конкретно тут — о том, что компилятор обязан записать в память то, что пишут по указателю volatile, сразу, не задерживаясь; и прочитать по другому volatile указателю, не делая никаких предположений, что по этому указателю и почему.
Этого достаточно, чтобы убрать проблему strict aliasing в компиляторе.

> volatile нужен platform-agnostic, multi-threaded ONLY.

К тредовости оно вообще-то отношения в C/C++ не имеет от слова «никак», пока пользуетесь стандартными механизмами вроде mutex lock/unlock.
(Это отличается, кстати, от Java — где на volatile нагрузили синхронизацию доступа между тредами.)
Компилятор не имеет права кэшировать чтение или запись переменной даже без volatile через вызов любой неизвестной компилятору функции с побочным эффектом изменения чего-то в памяти. Операции с мьютексами декларируются со свойством такого изменения, даже если они инлайнятся.

> В худшем случае UB (да и UB на касте int/float мне очень сомнительно)…

Таки может :(
К тредовости оно вообще-то отношения в C/C++ не имеет от слова «никак»

Да ну? А вы про атомарные операции (конструкты, инструкции) что-нибудь слышали?
Т.е. простейшие вещи (ака одна переменная типа int) вы правда мютексами защищаете?


Мютекс кстати тут вообще никаким местом и именно что "от слова никак", ибо если явно не указан volatile, то мютекс вокруг вам в случае не очевидной (нежелательной) оптимизации не поможет "от слова ничем"…
Например если компилятор выпилит установку переменной, как неиспользуемой далее.


Про остальное же я промолчу, да…
Например:


Компилятор не имеет права кэшировать чтение или запись переменной даже без volatile через вызов любой неизвестной компилятору функции

А то что эту функцию (неявно для компилятора) исполняет прямо сейчас другой поток, это куда?
Я вам еще раз повторю volatile однопоточно не интересен (ну если не нативно какой либо контроллер с мапой переменная<->чего-нибудь железное, что сравнимо с многопоточностью, ибо есть еще один актер)...


И да здесь не про Java...

> Да ну? А вы про атомарные операции (конструкты, инструкции) что-нибудь слышали?

«Мы» слышали. В том числе то, что атомарные переменные вообще-то не объявляют volatile, у них подобная функциональность заложена в штатные средства доступа.

> Т.е. простейшие вещи (ака одна переменная типа int) вы правда мютексами защищаете?

Давайте Вы внимательно прочтёте, на что отвечаете. А то такие Ваши домыслы, мягко говоря, удивляют.

> Например если компилятор выпилит установку переменной, как неиспользуемой далее.

Он не выпилит установку переменной, которая не локальная, потому что не следит за тем, кто её дальше будет использовать. Локальные для функции — они в таком обычно всё равно не участвуют (а если участвуют, это вообще-то не лучший стиль).

> А то что эту функцию (неявно для компилятора) исполняет прямо сейчас другой поток, это куда?

Какое отношение вообще имеет _исполнение_ функции в другом потоке? Или в Вашей обстановке у них при этом общие локальные данные?

> Я вам еще раз повторю volatile однопоточно не интересен

И я повторю, что он решает проблему со strict aliasing.

> Про остальное же я промолчу, да…

Настоятельно прошу вместо подобных «многозначительных» «умолчаний» вести обсуждение по сути.
В том числе то, что атомарные переменные вообще-то не объявляют volatile
Давайте Вы внимательно прочтёте, на что отвечаете.

Вы про что вообще? Здесь вообще-то речь про int i (который атомарный) и про функцию FastInvSqrt (в которой ни в каком её месте не нужен volatile от того самого слова).


Настоятельно прошу… вести обсуждение по сути.

По сути: расскажите где в этом куске из трех строчек нужен volatile и почему и где "компилятор обязан записать в память сразу"...


float FastInvSqrt(float x) {
  ...
  int i = *(int*)&x;
  i = 0x5f3759df - (i >> 1);
  x = *(float*)&i;
  ...
}

Даже если компилятор соптимизирует это совсем без использования памяти (т.е. будет использовать например регистр EAX для преобразования — он обязан следовать логике преобразования для каста как написано (а не как вы думаете).


Т.е. для того-же x86:


; int i = *(int*)&x;
mov         eax,dword ptr [x]
; i = 0x5f3759df - (i >> 1);
sar         eax,1
mov         ecx,5F3759DFh
sub         ecx,eax
; x = *(float*)&i;
mov         dword ptr [x],ecx

ВСЕ!


Причем здесь вообще "strict aliasing"? Тип int (как и float) не являются структурами...


Если вы про это, почему mov dword ptr [x],ecx вместо например:


; x = *(float*)&i;
mov         dword ptr [i],ecx
fld         dword ptr [i]
fstp        dword ptr [x]

То и это совсем из другой оперы и не имеет ни к volatile ни к "strict aliasing", ни даже к выравниванию никакого отношения.


Под UB же я имел ввиду как раз например если типы для платформы есть не равной длинны, и/или для этой платформы биты перевернуты LSB/MSB или little-endian vs. big-endian или подобное (однако и тут при желании оно лечится двумя строчками обернутыми в #if… #endif).

> По сути: расскажите где в этом куске из трех строчек нужен volatile

Вы сузили контекст до конкретного куска кода. Этого не было в предыдущем обсуждении.

Если говорить об этом коде, то именно за счёт того, что x безусловно перетирается, проблемы нет.
Если же он был бы написан чуть-чуть иначе — например, было бы так

  int i = *(int*)&x;
  i = 0x5f3759df - (i >> 1);
  *(int*)&x = i;

(одна только последняя строчка поменялась)

то проблема уже, очень вероятно, вылезла бы в полный рост, — компилятор мог бы просто не заметить, что x изменено. (Мало ли кто пишет чего по постороннему адресу...)

И GCC начинает об этом предупреждать — мне он на неё дословно сказал

fsq3.cc:3:19: warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]
   int i = *(int*)&x;  // представим биты float в виде цел
                   ^
fsq3.cc:5:11: warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]
   *(int*)&x = i;
           ^


К вопросу о том, структуры это или нет, это не имеет ни малейшего отношения. Коллега Halt в соседних комментариях приводил пример (в видео), как это реально влияло на каком-то компиляторе — аналогом для данной функции было бы, что код свернулся до { return x; }

К этому, все Ваши реплики как
> он обязан следовать логике преобразования для каста как написано (а не как вы думаете).
> Тип int (как и float) не являются структурами…

относятся к совсем другим аспектам, но не к этому.

А вот volatile помогло бы именно за счёт того, что запись в *(int*)&x переписало бы независимо от того, что компилятор думает о том, как относятся x и i друг к другу — он бы просто обязан был держать обоих в памяти и не делать никаких предположений о том, что x не изменилось, раз он не видит данного изменения. И это даже для локальных переменных. В выходном коде видно, что они явно пишут в память (в стек) и тут же из неё читают. Именно об этом я говорил как о «паре лишних тактов».

К счастью, пока что (?) встреченные Clang и GCC во всех тестах сумели как-то опознать, что это одна и та же сущность в памяти. Но гарантировать это я бы не пытался.

Если хотите возражать — прошу это делать со ссылками на стандарты (на публичные final drafts, если иного источника нет). Любая отсылка вида «у меня всё работает» тут, увы, совершенно непригодна.
Вы сузили контекст до конкретного куска кода. Этого не было в предыдущем обсуждении.

Чего это? Велось обсуждение конкретного куска кода из статьи — ничего я не "сужал"!
Мы что обсуждаем-то и где?
Во вторых там выше черным по белому стоит "конкретно тут — ни о чем…".


Вы же влезли в эту ветку с доказательствами "О чем, еще как о чем!" А потом выясняется, что для какого-то отстраненного теоретического случая (который здесь ни к месту и ни про то вовсе).


К этому, все Ваши реплики как относятся к совсем другим аспектам, но не к этому.

(Подавившись) вы очень нахальны, молодой человек...


*(int*)&x = i;
то проблема уже, очень вероятно, вылезла бы в полный рост, — компилятор мог бы просто не заметить, что x изменено.

Да ну? А ничего что в данном случае его должно интересовать что i (а не x) изменено, и "не заметить" этого компилятор ну никак не сможет (сдается мне сударь, вы не совсем понимаете о чем здесь пишете)…
Действительно есть такие компиляторы, которые при использования подобного "тухлого" каста (ударение на подобного) имеют UB (который описан у них в доку и как минимум выкинет warning при агрессивной оптимизации), но...


Но во вторых и у ув. мистера Уолша и у мистера Кармака хватило ума не использовать такой "гнилой" каст, и они таки написали x = *(float*)&i; вместо вот этого вот *(int*)&x = i;.


Если хотите возражать — прошу это делать со ссылками на стандарты.

Я вам вовсе ничего не собираюсь не доказывать не возражать...


Для начала, что бы не быть голословным быть может вы мне приведете пример ("стандарт, публичные final drafts") где описан такой UB и на каком основании компилятор не заметит изменения "x" после даже такого "гнилого" каста и использует старое (не измененное) значение.


*(int*)&x = i;
x = x*(1.5f-(xhalf*x*x));
Для начала, что бы не быть голословным быть может вы мне приведете пример («стандарт, публичные final drafts») где описан такой UB и на каком основании компилятор не заметит изменения «x» после даже такого «гнилого» каста и использует старое (не измененное) значение.


Да пожалуйста:
C99 – 6.5 Expressions

An object shall have its stored value accessed only by an lvalue expression that has one of the following types:
  • a type compatible with the effective type of the object,
  • a qualified version of a type compatible with the effective type of the object,
  • a type that is the signed or unsigned type corresponding to the effective type of the object,
  • a a type that is the signed or unsigned type corresponding to a qualified version of the effective type of the object,
  • an aggregate or union type that includes one of the aforementioned types among its members (including, recursively, a member of a subaggregate or contained union), or
  • a character type.


И сноска: The intent of this list is to specify those circumstances in which an object may or may not be aliased.


То есть стандарт явно разрешает считать что значение x и значение *(int*)&x никак, совсем никак не связаны между собой. Равно как и значения i и *(int*)&i. Тем самым вторая строка в этом примере — не зависит от первой и от параметров функции. А четвёртая строка — не зависит от третьей. То есть пару из второй и третьей строки можно вынести куда угодно. C учётом того, что значение i, вычисленное в строке 4 ни на что не влияет вообще — можно эту пару строк и удалить…

Да кстати о птичках...


GCC… warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]

$ cd /tmp/
$ echo '#include <stdio.h>

float FastInvSqrt(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;
}

int main()
{
  printf("fastinvsqrt of 9.0: %10.5f\n", FastInvSqrt(9.0));
  return 0;
}
' > test.c
$ gcc -Wstrict-aliasing test.c -o test
$ # и даже совсем усе-усе:
$ gcc -Wall -Wextra test.c -o test
$ ./test
fastinvsqrt of 9.0:    0.33295

Чисто для самообразования… На которой версии gcc (и/или платформе) вы видите такое?
А то я перепробовал от 3.4 до 6.2 (linux, win, и от отчаянья даже mac и solaris), нигде не увидел...

Сорри за спам, забылось -O2 (вот блин что autoconf с людьми делает)…
Т.е. воспроизводится везде, начиная с 3.4…
Еще раз звиняюсь...

g++ -Wall -Wextra -O2 test.c -o test
test.c: In function ‘float FastInvSqrt(float)’:
test.c:5:19: warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]
   int i = *(int*)&x;
                   ^
test.c:7:17: warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]
   x = *(float*)&i;

Правда, скомпилировать код так, чтобы он сломался, я не смог.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609

Покреативив вокруг да около, сделал чуть более "покладистый" вариант, вдруг кому буде нужно:


Example
$ echo '
#include <stdio.h>
#include <assert.h>

#if defined(static_assert) || (__STDC_VERSION__ >= 201100L)
static_assert(sizeof(int) == sizeof(float), "Cannot cast types of different sizes (int/float)");
#else
#warning Unknown sizes of int/float to the compile time: assumed are equal (4/4).
#endif

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

int main()
{
  printf("fastinvsqrt of 9.0: %10.5f\n", FastInvSqrt(9.0));
  return 0;
}
' > test.c; gcc -Ofast -Wall -Wextra test.c -o test; ./test

Как видно собирается даже с fast оптимизацией без всяких...


По хорошему к assert по размерности, сюда еще нужен бы assert, проверяющий endianness (если вдруг оне различные для int и float)…
Ну да это сильно муторно компиляторо-платформозависимо… потому, пусть будет так.

Причем здесь вообще «strict aliasing»? Тип int (как и float) не являются структурами...
А какая разница? Strict aliasing гласит, что записать по адресу "(int*)&x" никак не может повлять на значение x!

он обязан следовать логике преобразования для каста как написано (а не как вы думаете).
С кастом — всё хорошо. С обращением в память — плохо.

В следующей программе:
  float x;
  float *q=&x
  int *p=(int*)&x;
  ...
Компилятор имеет право считать, что запись через "*p" и чтение через "*q" (и наоборот) никак друг с другом не связаны. То же самое — в исходной программе. В строке x = *(float*)&i; x получает некоторое значение — но, с точки зрения компилятора, точно не имеющее никакого отношения к тому, что находится в i. Стало быть манипуляции с i — можно перенести ниже. А поскольку строка int i = *(int*)&x; инициализирует i каким-то значением, не имеющим отношения к тому, что находится в i — то это значение тоже можно вынести «вниз». А другая фаза того же компилятора может заметить, после всего этого, что мы пишем в x значение, которое ещё ничем не инициализировано вообще — ну так чего бы туда 0 не записать?
Со вторым не согласен в корне — с какого это перепуга он его вырежет?
С целью оптимизации кода. Код, который нифига не делает — точно короче того, который что-то делает.

В худшем случае UB (да и UB на касте int/float мне очень сомнительно)...
Тут не каст int/float. Тут Type punning. В результате компилятор может, ну, например, переставить работу с целыми числами до работы с плавучкой или после.

То есть будет что-то подобное:
float FastInvSqrt(float x) {
  int i;
  float xhalf = 0.5f * x;
  x = *(float*)&i;
  x = x*(1.5f-(xhalf*x*x));
  i = *(int*)&x;  // представим биты float в виде целого числа
  i = 0x5f3759df - (i >> 1);  // какого черта здесь происходит ?
  return x;
}
А после можно и «заметить» что мы обращаемся к неинициализированной переменной. Ну и раз мы всё равно можем вернуть что угодно — вернуть 0.
У меня дежавю, кажется что-то похожее я видел в хабе «Разработка игр».
UFO just landed and posted this here
Есть какие то замеры производительности по сравнению с традиционным способом? С 2005 года и уж с тем более с 80х прошло много лет, может быть сегодня овчинка уже не стоит выделки?
Что-то с 2005 года поменялось в теории сложности или архитектуре процессоров, что операции корня и деления стали работать быстрее суммирования и умножения? Вряд ли.
Что-то с 2005 года поменялось в теории сложности или архитектуре процессоров
Таки да
Вот это да. И всё же вопрос был о сравнении с «традиционными способами», а вшитая в проц аппроксимация — это не оно.
Мнэээ…

1. Таки нет, чисто формально, ибо RSQRTPS — команда SSE1, а это 1999, а не 2005 :)

2. Вопрос был про «корня и деления» против «суммирования и умножения». Эти фокусы позволяют, например, ускорять деление через обратный корень, но с потерей точности — вот Intel про это пишет в 248966:

> In microarchitectures that provide DIVPS/SQRTPS with high latency and low throughput, it is possible to speed up single-precision divide and square root calculations using the (V)RSQRTPS and (V)RCPPS instructions. For example, with 128-bit RCPPS/RSQRTPS at 5-cycle latency and 1-cycle throughput or with 256-bit implementation of these instructions at 7-cycle latency and 2-cycle throughput, a single Newton-Raphson iteration or Taylor approximation can achieve almost the same precision as the (V)DIVPS and (V)SQRTPS instructions.

но всё равно это сильно дороже умножения…

В этом плане лучше поставить на оптимизацию Goldschmidt method — AMD клянётся, что получает на делении до 16 циклов на single, 20 на double, 24 на extended — это заведомо быстрее типичного SRT. И наверняка там ещё можно что-то выжать. (А для квадратного корня ещё быстрее получается.)
Обратная сторона — всё в плавучке — прибавлять 10 тактов на конверсию в обе стороны для целых это уже сравнимо с интеловским SRT…
1. Таки нет, чисто формально, ибо RSQRTPS — команда SSE1, а это 1999, а не 2005 :)
AMD K6-III выпускались точно до 2000го года, а компьютеры, я думаю, продавались и позже.

но всё равно это сильно дороже умножения…
Не сильно. На каком-нибудь K7 RSQRTPS вообще быстрее умножения. Можете в таблички посмотреть.

На современных процессорах на оптимизацию подобных инструкций «забили», потому что кто же считает подобные вещи на CPU? Для этого GPU есть…
> На каком-нибудь K7 RSQRTPS вообще быстрее умножения.

Смешно — но показывает странное направление работы AMD.

> На современных процессорах на оптимизацию подобных инструкций «забили», потому что кто же считает подобные вещи на CPU? Для этого GPU есть…

Ну, похоже, не совсем. Например, после Nehalem повырезали много идиотизмов работы с денормализованными (там, где наследие K-C-S тратило тысячи тактов). Как минимум в SSE оба вендора не игнорируют плавучку. Но это происходит медленнее, чем хотелось бы…
Есть же еще математические сопроцессоры (FPU). Например, в STM32 вычисление квадратного корня занимает 14 тактов, что не так много. Деление — еще 14. Так что если нужен обратный квадратный корень — лучше взять этот алгоритм, если прямой — лучше вычислить на FPU.
Вот как раз в шейдерах-то таких хаков не надо, они только вредят. GPU обычно существенно быстрее работают с плавающей точкой, они на это заточены.

Шикарно!


Поделюсь трюком, который я вычитал в исходниках Elite для zx spectrum: быстрое умножение двух однобайтных чисел.
ab = (a^2 + b^2 - (a-b)^2)/2
Эта формула элементарно выводится из школьного разложения квадрата разности. Т.е. можно быстро и точно умножать однобайтные числа, используя таблицу квадратов (512 байт)

Это было реально круто на процессорах Z80 или 6502, на которых отсутствовала операция умножения. Современные же процессоры перемножают числа за 2-3 такта, так что необходимость в подобных трюках, в общем и целом, отпала…

Современные процессоры реализуют эту оптимизацию в железе, потому и так мало тактов требуется.

C современными возможностями по количеству вентилей вообще умножение строят на сетках и затем параллельных сумматорах… так что обычно не «эта» оптимизация. Хотя результат ещё лучше :)
Однако блоков умножения у современных процессоров обычно меньше чем ALU, и при out-of-order исполнении давление на эти блоки может оказаться высоким, а ALU блоки простаивать. Так что такие оптимизации все еще могут иметь смысл.
Не имеют. Архитектура процессоров меняется настолько часто, что вы не успеете программу пересобирать под 100 разных архитектур…
Меняются то они часто, но вот принципиально блоков умножения больше чем блоков ALU не становится, т.к. эти блоки дорогие по площади кристалла. Вот например Skylake посмотрите https://en.wikichip.org/wiki/intel/microarchitectures/skylake_(client)#Execution_engine
ALU по прежнему больше и тому объективные причины, а не потому что не хочется.
Три дереференса + две суммы + одно вычитание + 1 побитовый сдвиг вместо одного умножения? оО
Там всегда a не меньше b?
Иначе ещё нужна проверка. (Или я чего-то не понимаю?)

Программная реализация умножения 8-битных чисел выполняется за пару сотен тактов на z80. умножение с использованием таблиц квадратов — несколько десятков тактов. При частоте процессора в 3,5 MHz ещё как ощутимо

Даже больше — умножение порядка 800 тактов, Z-80 весьма долго (в плане количества тактов) выполняет команды — минимум 4, максимум больше 20. Решил вспомнить старое и накидать умножение по таблице — вышло около 170 тактов (зависит от переходов).
Заголовок спойлера
	; HL = L*C ~ 171 clocks
    
	LD	H, high (SQRTABLE); 7

	; DE = L^2
	LD	E, (HL)		; 7
	INC	H		; 4
	LD	D, (HL)		; 7

	; A = abs (L - C)
	LD	A, L		; 4
	SUB	C		; 4
	JR	NC, L1		; 7/12
	NEG			; 8
L1:
	LD	L, A		; 4
	LD	A, C		; 4

	; BC = (L - C)^2
	LD	B, (HL)		; 7
	DEC	H		; 4
	LD	C, (HL)		; 7

	; HL = L^2 - (L - C)^2
	EX	DE, HL		; 4
	AND	A		; 4
	SBC	HL, BC		; 15
	EX	DE, HL		; 4

	; BC = C^2
	LD	L, A		; 4
	LD	C, (HL)		; 7
	INC	H		; 4
	LD	B, (HL)		; 7

	; HL = L^2 + C^2 - (L - C)^2
	EX	DE, HL		; 4
	JR	NC, L2		; 7/12
	
	; If L^2 < (L - C)^2,
	; L*C cannot be more than 32767.
	ADD	HL, BC		; 11

	; /2
	SRL	H		; 8
	RR	L		; 8

	RET			; 10

L2:
	ADD	HL, BC		; 11

	; /2 using the carry flag as 17th bit
	RR	H		; 8
	RR	L		; 8

	RET			; 10

я вычитал в исходниках Elite для zx spectrum

Где где? Можно ссылку на исходники?

Ну это просто. Загружаешь MONS на 58500, запускаешь его. вбиваешь загрузчик кода Elite, он с адреса 5B00, длину блока указываешь так, чтобы MONS не перетереть. Потом переходишь на 7CB7, это точка входа в игру. И жмешь комбинацию клавиш для показа листинга. Вуаля — на экране исходник.


20 лет прошло, до сих пор помню :-)

Есть более известный вариант:

ab = (a+b)^2/4 - (a-b)^2/4


только два лукапа по таблице, для a+b и a-b; деление на 4 выполняется нацело, дробные взаимно компенсируются.

Тогда a+b может быть больше 255, т.е. нужно оперировать словами и таблицу квадратов удваивать по размеру.

Да. Тут уже надо компромисс искать.
Подозреваю, именно для Z80-based сокращение таблицы было важнее, но точными данными не владею.
Такой вариант на Z80 будет чуть быстрее (137) — к счастью, (a+b)^2/4 все еще влазит в 16 бит, поэтому таблица будет всего 1К. Даже для спектрума 48К это вполне допустимо для такой ключевой функции, как умножение.
Дополнительный 9-й бит результата (a + b) попадет во флаг переноса, можно сделать условный переход и, таким образом, продолжить решать задачу в 8-битной арифметике.
Заголовок спойлера
	; HL = B*C

	LD	A, B	; 4
	ADD	C	; 4
	JR	NC, L1	; 7/12
	INC	H	; 4
	INC	H	; 4
L1:
	; DE = (B + C)^2/4
	LD	L, A	; 4
	LD	H, high (SQRTABLE); 7
	LD	E, (HL)	; 7
	INC	H	; 4
	LD	D, (HL)	; 7

	LD	A, B	; 4
	SUB	C	; 4
	JR	NC, L2	; 7/12
	NEG		; 8
L2:
    ; BC = (B - C)^2/4
	LD	L, A	; 4
	LD	H, high (SQRTABLE); 7
	LD	C, (HL)	; 7
	INC	H	; 4
	LD	B, (HL)	; 7

	EX	DE, HL	; 4
	AND	A	; 4
	SBC	HL, BC	; 15
	RET		; 10

Кто бы еще смог понять, что здесь происходит —
Эвклид без Эвклида ?
Понятно, что здесь происходит вычисление обратного числа, Y = X^-1
Но почему именно так? Я хотел было призвать алгоритм Эвклида, но наткнулся на это!
И почему для 64х бит достаточно всего лишь одной лишней итерации?
И что за странный комментарий 3 bits correct? Это потому что odd?

// modular inverse of odd 32-bit value
uint32_t modinv32( uint32_t x ) {
        uint32_t y = x;      // 3 bits correct
        y *= 2 - x * y; // 6
        y *= 2 - x * y; // 12
        y *= 2 - x * y; // 24
        y *= 2 - x * y; // 32
        return y;
}

// modular inverse of odd 64-bit value
uint64_t modinv64( uint64_t x ) {
        uint64_t y = modinv32( (uint32_t) x );
        y *= 2 - x * y;
        return y;
}



Это не алгоритм Евклида, а существенно более оптимальный, применимый для такого частного случая — получается сложность O(log k) вместо O(k) для Евклида. Переход можно обосновать так: если y = x^-1 mod 2^k, то xy = a 2^k + 1; домножим, как в вашей функции, на 2 - xy: получим xy * (2 - xy) = (a 2^k + 1) * (2 - a 2^k - 1) = -a^2 2^(2k) + a 2^k - a 2^k + 1 = 1 mod 2^(2k), то есть из обратного к х значения по модулю 2^k получили обратное по модулю 2^(2k). Делая такие переходы как раз из обратного по модулю 2^3 получим по 2^6 и т.д. Также отсюда ясно, почему для перехода к 64 битам достаточно одного шага.

А, про первые 3 бита забыл написать. Тут можно очень просто убедиться, перебрав все нечётные числа от 1 до 7 (т.к. рассматриваем обратное по модулю 2^3 = 8) — как видно, каждое число является обратным самому себе. Ну и ещё — к чётным числам обратного ясное дело не существует, поэтому их вообще не рассматриваем.

А в С++ в те времена int уже был 32-битным? Или это не имеет значения для указателей?
В какие «те»? Quake III Arena, о котором идёт речь в статье, вышел в 1999 году, а int стал 32-битным (сначала это было даже опционально, в зависимости от модели памяти) с распространением 32-битных операционок, т.е. ещё на пару лет раньше.
Quake III Arena, о котором идёт речь в статье, вышел в 1999 году, а int стал 32-битным (сначала это было даже опционально, в зависимости от модели памяти) с распространением 32-битных операционок, т.е. ещё на пару лет раньше.
Какие, к терапевту, 32-битные операционки? 32-битный int — это Doom, 1993й год.

Только для денормализованных значений аргумента (близких к нулю) в результате, кажется, будет мусор, а не ожидаемое +inf.

да, формула только для x>1
А побитовое извлечение корня разве не быстрее? Там только целочисленные операции, без операций с плавающей точкой. Кто-нибудь сравнивал?

Как вы будете извлекать квадратный корень из числа с плавающей точкой без операций с плавающей точкой?

Любое число с плавающей точкой f представляется в компе в виде f = z * 2^k. Его всегда легко привести к виду f = z' * 2^(2*k'). Тогда sqrt(f) = sqrt(z') * 2^k'.

z, z', k, k' — целые. Таким образом, вычисление корня из f легко сводится к вычислению корня из целого числа z'.
Ну а из целого числа как вы корень извлечёте? И ещё, чтобы такое провернуть нужно отдельно проводить операции над мантиссой и экспонентой, а процессоры хуже работают с числами нестандартного размера, да ещё и не выровненными.
Из целого числа есть поразрядный алгоритм, наподобие умножения столбиком. В десятичной системе он сложный, после адаптации к двоичной превращается в двоичные побитовые сдвиги и сравнения. Погуглите «вычисление квадратного корня столбиком». Например, dxdy.ru/topic24572.html только нужно алгоритм из десятичной системы перевести в двоичную, он сильно упрощается.

Причём тут невыровненные числа? В отдельные переменные (регистры) выделяется мантисса и экспонента, затем с каждыми из них операции ведутся как с обычными машинными словами. После вычислений float/double собирается обратно побитовыми операциями.
Если квадратный корень — целое число, то он является суммой нечетных чисел от 1 до 2n-1 (добавляем 2 порядка, от полученного значения убираем 1 порядок). Замерял скорость на скриптовых языках. Не впечатлён. Для C/C++ алгоритмов много.
Для вычисления целочисленных корней на 8 бит PIC микроконтроллере давно использую Hardware algorithm [GLS] из hackersdelight. Очень компактный и шустрый.
Есть ещё целочисленный алгоритм, там только побитовые сдвиги и сравнения, нет умножений. Когда-то писал на ассемблере, но листочек затерялся, не могу найти. Но можно попробовать воспроизвести.
Огромное спасибо! Ради таких статей я и пришёл в свое время на этот ресурс.
Математика === Магия. Классная статья, спасибо!
Проверил на GCC 4.9 (pentium 2020)
Быстрый метод в 27 раз быстрее, чем cmath-pow.
Средняя накопительная погрешность 0,1%
Припоминается с древних времен:

«Во имя эффективности в программировании было совершено больше прегрешений (причем не всегда ее удавалось достичь), чем по какой-либо другой причине, включая непроходимую глупость»
В. Вульф

Если в вашей программе есть место, которым вы особенно гордитесь, выбросьте его.
Э. Йордан.

Немного опоздал к обсуждению, но вот подробный разбор 0x5f3759df c проверкой оптимальности этой константы для начального приближения к значению обратного корня.
Автор показывает, что 0x5f375a86 на самом деле даёт немного более точные результаты.

Интересная работа, но там (как и в этом хабратопике) точность рассматривается только для первой половины Q_rsqrt — до ньютоновской итерации.
Как ни парадоксально, повышение точности начального приближения не повышает точность окончательного результата (после ньютоновской итерации).
Интересно было бы разобраться, почему так :)

>Но и это ещё не всё!
>И даже более того!

К концу статьи стал подозревать, что это магазин на диване и вы собираетесь мне что-то продать.
Я возьму число 0.0450465, как дающее неплохое приближение и использованное в оригинальной реализации.

Почему именно это число? Как получить такое число, не зная о нем заранее?
как дающее неплохое приближение

Видимо, каким-нибудь линейным или градиентным спуском перебирались числа в близком диапазоне для поиска наилучшего приближения (представимого с помощью float).
Я выше кидал ссылку на подробный разбор.
Судя по тому, что это не наилучшее приближение — это число действительно подбиралось методом проб и ошибок.
Sign up to leave a comment.