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

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

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

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

Разложение синуса в ряд Тейлора.

Функция синуса раскладывается в бесконечный ряд Тейлора.

$\$\$display\$\$\backslash sin(x)=x-\{\backslash frac\; \{x^\{3\}\}\{3!\}\}+\{\backslash frac\; \{x^\{5\}\}\{5!\}\}-\{\backslash frac\; \{x^\{7\}\}\{7!\}\}+\{\backslash frac\; \{x^\{9\}\}\{9!\}\}-\backslash cdots\; \$\$display\$\$$

Понятно, что бесконечный ряд мы посчитать не можем, кроме случаев, когда есть аналитическая формула бесконечной суммы. Но это не наш случай))) Предположим, что мы хотим посчитать синус в интервале $\$inline\$[0,\; \backslash frac\{\backslash pi\}\{2\}]\$inline\$$. Более подробно работу с интервалами обсудим в части 3. Зная, что $\$inline\$\backslash sin(\backslash frac\{\backslash pi\}\{2\})=1\$inline\$$ оценим найдём первый член который можно отбросить исходя из условия, что , где $\$inline\$e\$inline\$$ это разница между числом 1 и наименьшем числом, которое больше 1. Грубо говоря это последний бит мантиссы (wiki). Решить данное уравнение проще перебором. Для $\$inline\$e\; \backslash approx\; 2.22\backslash times10^\{-16\}\$inline\$$. У меня получилось $\$inline\$n=23\$inline\$$ уже можно отбросить. Правильный выбор количества слагаемых будет обсуждено в одной из следующей частей, поэтому на сегодня «перестрахуемся» и возьмём слагаемых до $\$inline\$n=25\$inline\$$ включительно.
Последнее слагаемое приблизительно в 10000 раз меньше, чем $\$inline\$e\$inline\$$.

Простейшее решение

Руки уже чешутся, пишем:

#include <iostream> #include <iomanip> #include <cmath> #include <array> #include <bitset> #include <quadmath.h> // Полный путь к файлу для clang //#include "/usr/lib/gcc/x86_64-linux-gnu/10/include/quadmath.h" #include <numeric> #include <limits> #include <vector>  #include <boost/timer/timer.hpp> #include <boost/math/special_functions/factorials.hpp>  namespace bm = boost::math;  using namespace std;  typedef union { uint32_t i[2]; double x; } mynumber;  array<double, 25> fc;  double sin_e1(double x) {   double result = 0;   int sign = 1;   for(int i = 1; i < 25; i += 2) {     result += sign * pow(x, i) / bm::unchecked_factorial<double>(i);     sign = -sign;   }   return result; }  double sin_e2(double x) {   double result = 0;   int sign = 1;   double xx = x * x;   double pw = x;   double fti = 1.0;   for(int i = 1; i < 25; i += 2) {     fti /= i;     result += sign * pw * fti;     fti /= ( i + 1 );     sign = -sign;     pw  *= xx;   }   return result; }  double sin_e3(double x) {   double result = 0;   for(int i = 25; i >= 1; i -= 2) {     result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);   }   return result; }  double sin_e4(double x) {   double xx = x * x;   double res = fc[25];   for(int i = 23; i >= 1; i -= 2) {     res = fc[i] + xx * res;   }   return x * res; }  double sin_e5(double x) {   double xx = x * x;   double res = fc[25];   for(int i = 23; i >= 3; i -= 2) {     res = fc[i] + xx * res;   }   return x + x * xx * res; }  #define SIN(a) sin_e5(a) // ^^ Изменить функцию для вычисления здесь. ^^  int main() {    __uint128_t ft = 1;   fc[1] = 1.0; //3 * 5;   for(int i = 2; i < fc.size(); i++) {     ft *= i;     // factorial with sign for Taylor series     fc[i] = (1.0 / ft) * (( (i - 2) % 4 < 2) ? -1 : 1);   }   vector<double> xv;   xv.resize(8 * 2000000);   // Линейное заполнение массива значениями от 0 до M_PI/2   for (int i = 0; i < xv.size(); i++) {     // Максимальное значение в массиве изменять здесь.     xv[i] = (M_PI / 2) * i / double(xv.size());   }    {     boost::timer::auto_cpu_timer at;     double res = 0;     for(int i = 0; i < xv.size(); i++) {       res += SIN(xv[i]);     }   }    int co = 0, cn = 0;   // Используем числа четверной точности как эталон.   __float128 avg = 0.0, div = 0.0;   for(int i = 0; i < xv.size(); i++) {     mynumber dold, dnew;     dold.x = sin(xv[i]);     dnew.x = SIN(xv[i]);     __float128 q = sinq(xv[i]); // <= sinq считаем эталоном.     __float128 dd = __float128(dnew.x) - q;     // Вычисляем среднее и стандартное отклонение.     div += dd * dd;     avg += dd;     // Сравниваем побитово, что значения синуса от встроенной функции и от нашей.     // Если они различаются, то выясняем какая из функций даёт более правильный результат.     if( dold.i[0] != dnew.i[0] || dold.i[1] != dnew.i[1] ) {       if( fabsq(q - dold.x) <= fabsq(q - dnew.x) )         co++;       else         cn++;     }   }   avg /= xv.size();   div /= xv.size();    // Количество случаев, когда внутренняя функция дала более правильный результат к общему количеству вычислений.   cout << "Better libm: " <<  co << " / " << xv.size() << "(" << 100.0 * co / xv.size() << "%)" << endl;    // Количество случаев, когда "наша" функция дала более правильный результат к общему количеству вычислений.   cout << "Better new: " <<  cn << " / " << xv.size() << "(" << 100.0 * cn / xv.size() << "%)" << endl;    // Среднее отклонения и отклонение отклонения нашей функции от эталона.   cout << "  Avg / std new: " << double(avg) << " / " << double(sqrtq( div - avg * avg )) << endl;   return 0; } 

double sin_e1(double x) {   double result = 0;   int sign = 1;   for(int i = 1; i < 25; i += 2) {     result += sign * pow(x, i) / bm::factorial<double>(i);     sign = -sign;   }   return result; } 

Как ускорить программу я думаю, что многие сообразили сразу. Как вы думаете, во сколько раз ваши изменения могут ускорить программу? Оптимизированная версия и ответ на вопрос под спойлером.

double sin_e2(double x) {   double result = 0;   int sign = 1;   double xx = x * x;   double pw = x;   double fti = 1.0;   for(int i = 1; i < 25; i += 2) {     fti /= i;     result += sign * pw * fti;     fti /= ( i + 1 );     sign = -sign;     pw  *= xx;   }   return result; } 

Улучшение точности

Методика

Точность вычисления функции будем определять 2-мя стандартными параметрами.

Среднеквадратичное отклонение от истинного значения sin(float128) и среднее данного отклонения. Последний параметр может дать важную информацию о том, как ведёт себя наша функция. Она может систематически занижать или завышать результат.

В дополнение к данным параметрам ввёдём еще два. Вместе с нашей функции мы вызываем ещё встроенную в библиотеку функцию sin(double). Если результаты двух функций: нашей и встроенной не совпадают (побитово), то добавляем в статистику, какая из двух функций дальше от истинного значения.

Порядок суммирования

Вернёмся снова к исходному примеру. Как можно увеличить его точность «по-быстренькому»? Те, кто внимательно читал статью Можно ли сложить N чисел типа double наиболее точно? скорее всего дадут ответ сразу. Надо крутить цикл в обратную сторону. Чтобы складывать от наименьших по-модулю, к наибольшим.

double sin_e3(double x) {   double result = 0;   for(int i = 25; i >= 1; i -= 2) {     result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);   }   return result; } 

Результаты приведены в табличке.

Может показаться, что использование алгоритмов «умного» суммирования уберёт ошибку практически до 0, но это не так. Конечно эти алгоритмы дадут увеличение точности, но для полного избавления от ошибок требуются ещё и алгоритмы умного умножения. Они существуют, но очень накладны: очень много лишних операций. Применение их здесь не оправдано. Впрочем позднее мы к ним вернёмся в другом контексте.

Осталось совсем немного. Объединить быстрый и точный алгоритмы. Для этого снова вернёмся к ряду Тейлора. Ограничем его для примера 4-мя членами и сделаем следующее преобразование.

$\$\$display\$\$\backslash sin(x)\backslash approx\; x(1+x^2(-1/3!+x^2(1/5!+x^2(-1/7!+x^2\backslash cdot1/9!))))\$\$display\$\$$

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

double sin_e4(double x) {   double xx = x * x;   double res = fc[25];   for(int i = 23; i >= 1; i -= 2) {     res = fc[i] + xx * res;   }   return x * res; } 

Работает быстро, но потеряли точность, по сравнению с e3. Опять же проблема в округлении. Давайте рассмотрим последний шаг цикла и немного преобразуем исходное выражение

$\$\$display\$\$\backslash sin(x)\backslash approx\; x+x\; \backslash cdot\; x^2(-1/3!+\; \backslash cdots))\$\$display\$\$$

И соответствующий код.

double sin_e5(double x) {   double xx = x * x;   double res = fc[25];   for(int i = 23; i >= 3; i -= 2) {     res = fc[i] + xx * res;   }   return x + x * xx * res; } 

Точность в сравнении с libm увеличилась в 2 раза. Если догадываетесь почему точность увеличилась, пишите в комментариях. К тому же есть ещё одна, гораздо более неприятная вещь у sin_e4, которая отсутствует у sin_e5, связанная с точностью. Попробуйте догадаться в чём проблема. В следующей части я обязательно о ней расскажу подробно.

Если статья Вам понравится, то в следующей я расскажу, как в GNU libc считается синус с максимальным ULP в 0.548.