Симметричный Lerp & amp; оптимизация компилятора

Если повсеместно используется один двоичный формат с плавающей точкой IEEE-754 (например, базовый 32-разрядный двоичный файл, этот формат обычно используется в C ++). float) со всеми операторами C ++, сопоставленными с операциями IEEE-754 прямым и простым способом, затем lerp_symmetric(alpha, x0, x1) (далее упоминается как A) равно lerp_symmetric(1-alpha, x1, x0) (B)

Доказательство:

• Если alpha, который мы предполагаем, находится в [0, 1], больше или равно ½, то 1-alpha точно по лемме Стербенса. (Под «точным» мы подразумеваем, что вычисленный результат с плавающей запятой равен математическому результату; ошибки округления нет.) Затем при вычислении A, w0 точный, так как это 1-alpha, а также w1 является точным, поскольку его математический результат alpha, так что это точно представимо. И в вычислительной технике B, w0 является точным, поскольку его математический результат alpha, а также w1 точно, так как это снова 1-alpha,
• Если alpha меньше ½, то 1-alpha может иметь некоторую ошибку округления. Пусть результат будет beta, Затем в A, w0 является beta, Сейчас ½ ≤ betaТаким образом, лемма Стербенса относится к оценке w1 = 1.0f - w0, так w1 является точным (и равен математическому результату 1-beta). И в B, w0 точно, опять же по лемме Стербенса, и равняется w1 из A, а также w1 (из B) является точным, поскольку его математический результат beta, что точно представимо.

Теперь мы можем видеть, что w0 в A равняется w1 в B а также w1 в A равняется w0 в B, Позволить beta быть 1-alpha в любом из вышеуказанных случаев, A а также B поэтому вернитесь (1-beta)*x0 + beta*x1 а также beta*x1 + (1-beta)*x0соответственно. Добавление IEEE-754 является коммутативным (за исключением полезных нагрузок NaN), поэтому A а также B вернуть идентичные результаты.

Отвечая на вопросы:

1. Я бы сказал, что это разумный подход. Я бы не утверждал, что нет улучшений, которые можно было бы сделать без дальнейших размышлений.

2. Нет, вы не можете доверять вашему компилятору:

   • C ++ позволяет реализациям использовать избыточную точность при оценке арифметики с плавающей точкой. таким образом w0*x0 + w1*x1 можно оценить с помощью double, long doubleили другая точность, даже если все операнды float,
   • C ++ допускает сокращения, если они отключены, поэтому w0*x0 + w1*x1 может оцениваться как fmaf(w0, x0, w1*x1), таким образом, используя точную арифметику для одного из умножений, но не другого.

Вы можете частично обойти это, используя:

float w0 = 1.0f - alpha;
float w1 = 1.0f - w0;
float t0 = w0*x0;
float t1 = w1*x1;
return t0+t1;

Стандарт C ++ требует, чтобы избыточная точность отбрасывалась в присваиваниях и приведениях. Это распространяется на возврат функций. (Я сообщаю об этой и других спецификациях C ++ из памяти; стандарт должен быть проверен.) Таким образом, каждый из приведенных выше округляет свой результат до float даже если изначально была использована дополнительная точность. Это предотвратит сокращение.

(Следует также иметь возможность отключить сокращение путем включения <cmath> и вставив директиву препроцессора #pragma STDC FP_CONTRACT off, Некоторые компиляторы могут не поддерживать это.)

Одна из проблем, связанных с вышеупомянутым обходным решением, заключается в том, что значения сначала округляются до точности оценки, а затем округляются до float, Существуют математические значения, для которых, для такого значения Икс, округление Икс первым double (или другой точности), а затем float дает другой результат, чем округление Икс прямо к float, Диссертация Строгая структура для полной поддержки стандарта IEEE для арифметики с плавающей точкой в ​​языках программирования высокого уровня Самуэль А. Фигероа дель Сид устанавливает, что при оценке одной операции умножения или сложения в базовой 64-битной IEEE-754 с плавающей запятой (обычно используется для double), а затем округление до 32-разрядного формата никогда не приводит к ошибке двойного округления (поскольку эти операции с учетом входных данных, которые являются элементами 32-разрядного формата, никогда не могут вызвать одну из проблемных Икс значения описаны выше).¹

Если я прав в отношении спецификаций C ++, сообщаемых из памяти, то описанный выше обходной путь должен быть завершен до тех пор, пока реализация C ++ оценивает выражения с плавающей запятой либо с номинальным форматом, либо с форматом, достаточно широким для удовлетворения требований, предъявляемых Figueroa del Cid. ,

¹ За Фигероа дель Сид, если x а также y иметь п-немного значащих, и x+y или же x*y вычисляется точно, а затем округляется до Q мест, второе округление до п места будут иметь тот же ответ, как если бы результат был непосредственно округлен до п места, если п ≤ (Q — 1) / 2. Это выполняется для базового 32-разрядного двоичного числа с плавающей запятой IEEE-754 (п = 24) и 64-разрядная (Q = 53) Эти форматы обычно используются для float а также doubleи обходного пути выше должно быть достаточно в реализации C ++, которая их использует. Если реализация C ++ оценивается float используя точность, которая не удовлетворяет условию, которое дает Фигероа дель Сид, тогда могут возникнуть ошибки двойного округления.