Обрезать числа с плавающей запятой и удвоения после определенных пользователем точек в FPU X87 и SSE

Как отметили несколько человек, более раннее округление не помогает окончательному результату быть более точным. Если вы хотите узнать больше о сравнениях с плавающей запятой и странностях / подводных камнях, я настоятельно рекомендую серию статей Брюса Доусона о плавающей запятой. Вот цитата из с индексом < / а>

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

Математика [с плавающей точкой] сложна.

Вы просто не поверите, насколько это невероятно сложно. Я имею в виду, вы можете подумать, что сложно рассчитать, когда поезда из Чикаго и Лос-Анджелеса столкнутся, но это пустяки в математике с плавающей запятой.

(Бонусные баллы, если вы узнаете этот последний абзац как пересказ известной строчки о космосе.)

Как на самом деле можно было реализовать свою плохую идею:

Нет никаких машинных инструкций или функций стандартной библиотеки C, которые можно было бы усечь или округлить до чего-либо, кроме целого числа.

Обратите внимание, что есть машинные инструкции (и функции C), которые округляют double до ближайшего (представимого) целого числа без преобразования его в intmax_t или что-то еще, просто _3 _-> _ 4_. Таким образом, нет двустороннего обхода через целое число с дополнением 2 фиксированной ширины.

Таким образом, чтобы использовать их, вы можете масштабировать плавающий объект в несколько раз, округлить до ближайшего целого числа, а затем снова уменьшить его. (например, функция chux на основе round(), но я бы рекомендовал C99 double rint(double) вместо round(). round имеет странную семантику округления, которая не соответствует ни одному из доступных режимов округления на x86, поэтому он компилируется в худший код.

В инструкциях asm для x86, которые вы все время упоминаете, нет ничего особенного, и они не делают ничего такого, о чем вы не можете попросить компилятор сделать с чистым C.

FISTP (Float Integer STore (и Pop the x87 stack) - один из способов компилятор или asm-программист для реализации long lrint(double) или (int)nearbyint(double). Некоторые компиляторы улучшают код для того или другого. Он округляется с использованием текущего режим округления x87 (по умолчанию: округление до ближайшего), который имеет ту же семантику, что и ISO Стандартные функции C.

часть SSE3, даже если он работает со стеком x87. Это позволяет компиляторам избегать установки режима округления на усечение (round -towards-zero), чтобы реализовать семантику усечения C (long)x, а затем восстановить старый режим округления.

В этом и заключается суть «не изменять контрольное слово». Ни одна из инструкций этого не делает, но для реализации (int)x без FISTTP компилятор должен использовать другие инструкции для изменения и восстановления режима округления вокруг инструкции FIST. Или просто используйте SSE2 CVTTSD2SI для преобразования двойного в регистре xmm с усечением вместо Значение FP в устаревшем стеке x87.

Поскольку FISTTP доступен только с SSE3, вы должны использовать его только для long double или в 32-битном коде, который в любом случае имеет значения FP в регистрах x87 из-за устаревшего 32-битного ABI, который возвращает значения FP в стеке x87.

PS. если вы не узнали HHGTG ссылку Брюса, оригинал:

Пространство большое. Действительно большой. Вы просто не поверите, насколько он невероятно велик. Я имею в виду, вы можете подумать, что до аптеки еще далеко, но в космос это пустяки.

как мне это сделать в X87 (скалярный FPU) и SSE (векторный FPU)?

Следующее не использует ни X87, ни SSE. Я включил его как ссылку сообщества на код общего назначения. Во всяком случае, его можно использовать для тестирования решения X87.

1. Любое "измельчение" результата g() по крайней мере незначительно увеличит ошибку, что, надеюсь, допустимо, как сказал OP: "Чтобы избежать распространения ошибки на другие вычисления ..."

2. Неясно, хочет ли OP «точные результаты с точностью до 5 знаков после запятой» для отражения абсолютной точности (+/- 0,000005) или относительной точности (+/- 0,000005 * результат). Предположим «абсолютную точность».

3. Поскольку float, double очень часто являются двоичными числами с плавающей запятой, любое «чоп» будет отражать число FP, ближайшее к кратному 0,00001.

4. Текстовый метод:

   //       - x    xxx...xxx    . xxxxx \0
   char buf[1+1+ DBL_MAX_10_EXP+1  +5   +1];
   sprintf(buf, "%.5f", x);
   x = atof(buf);
   

5. round() rint() метод:

   #define SCALE 100000.0
   if (fabs(x) < DBL_MAX/SCALE) {
     x = x*SCALE;
     x = rint(x)/SCALE;
   }
   

6. Прямая битовая обработка x. Просто обнуляйте биты выбора в мантиссе.

   TBD code.