inline float DifferenceOfProducts(float a, float b, float c, float d) {
float cd = c * d;
float err = std::fma(-c, d, cd);
float dop = std::fma(a, b, -cd);
return dop + err;
}
DifferenceOfProducts()
вычисляет таким образом, что избегает катастрофического сокращения. Впервые эта методика была описана легендарным Уильямом Кэхэном в статье On the Cost of Floating-PointComputation Without Extra-Precise Arithmetic. Стоит заметить, что работы Кэхэна интересно читать в целом, в них есть множество комментариев о текущем состоянии мира плавающих запятых, а также математических и технических рассуждений. Вот один из его выводов:Те из нас, кто сражался с превратностями арифметики с плавающей запятой и плохо продуманными «оптимизациями» компиляторов, могут справедливо испытывать гордость за победу в этой битве. Но если мы передадим продолжение этой битвы последующим поколениям, это будет противоречить всей сути цивилизации. Наш опыт говорит, что языки программирования и системы разработки являются источниками слишком большого количества хаоса, с которым нам приходится бороться. Слишком от многих ошибок вполне можно обойтись, как и от некоторых привлекательных «оптимизаций», безопасных для целых чисел, но иногда оказывающихся фатальными для чисел с плавающей запятой.
DifferenceOfProducts()
: в основе искусности этой функции лежит использование инструкций совмещённого умножения-сложения (fused multiply-add, FMA)3. С математической точки зрения FMA(a,b,c)
является , поэтому поначалу кажется, что эта операция полезна только в качестве микрооптимизаци: одна инструкция вместо двух. Однако FMADifferenceOfProducts()
. Снова покажу первые две её строки: float cd = c * d;
float err = std::fma(-c, d, cd);
err
всегда будет равна нулю. Но при работе с FMA вторая строка на самом деле извлекает величину погрешности округления в вычисленном значении и сохраняет её в err
. После этого на выходе всё получается очень просто: float dop = std::fma(a, b, -cd);
return dop + err;
return
«патчит» эту проблему, прибавляя погрешность, выделенную во второй строке. В статье Jeannenrod et al. показано, что результат верен до 1,5 ulps (мер единичной точности), что великолепно: FMA и простые операции с плавающей запятой точны до 0,5 ulps, поэтому алгоритм практически идеален.DifferenceOfProducts()
, она оказывается на удивление полезной. Вычисляете дискриминанта квадратного уравнения? Вызывайте DifferenceOfProducts(b, b, 4 * a, c)
4. Вычисляете определитель матрицы 2x2? Алгоритм решит эту задачу. В реализации следующей версии pbrt я нашёл ему около 80 применений. Из всех них самой любимой является функция векторного произведения. Она всегда была источником проблем, из-за которого приходилось поднимать руки вверх и использовать в реализации числа double, чтобы избежать катастрофического сокращения:inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) {
double v1x = v1.x, v1y = v1.y, v1z = v1.z;
double v2x = v2.x, v2y = v2.y, v2z = v2.z;
return Vector3f(v1y * v2z - v1z * v2y,
v1z * v2x - v1x * v2z,
v1x * v2y - v1y * v2x);
}
DifferenceOfProducts()
.inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) {
return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y),
DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z),
DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x));
}
DifferenceOfProducts()
и значениями float? . Значения и соответствуют двойной точности, а слегка смещён — отсюда и взялась лишняя ulp.DifferenceOfProducts()
выполняет две FMA, а также умножение и сложение. Наивный алгоритм можно реализовать с одной FMA и одним умножением, что, казалось бы, должно занять в два раза меньше времени. Но на практике оказывается, что после получения значений из регистров это неважно: в синтетическом бенчмарке, проведённом на моём ноутбуке, DifferenceOfProducts()
всего в 1,09 раза затратнее наивного алгоритма. Работа с двойной точностью была медленнее в 2,98 раза.DifferenceOfProducts()
является хорошим лекарством против большей их части. Она проста в использовании, и у нас нет особых причин не применять её.SumOfProducts()
, обеспечивающую защиту от катастрофического сокращения. Если вы хотите усложнить задачу, то объясните, почему в DifferenceOfProducts()
, dop + err == dop
, ведь знаки a*b
и c*d
различаются.std::fma()
; в gcc и clang работает -march=native
.4 * a
не вызывает никакой ошибки округления, если только не происходит переполнения.К сожалению, не доступен сервер mySQL