В качестве примера я выбрал тригонометрическую функцию синуса. Это широко распространённая функция, математика которой хорошо известна со школы и университета. В тоже время при её имплементации появятся много ярких примеров правильной работы с числами. В качестве числа с плавающей точкой я буду использовать double.
В данном цикле статей планируется много всего начиная от математики, заканчивая машинными кодами и опциями компилятора. Язык написания статьи С++, но без излишеств. В отличии от С работающие примеры будут более удобочитаемыми даже для людей не знакомым с этим языком и занимать меньше строк.
Статьи будут написаны методом погружения. Будут обсуждаться подзадачи, которые потом соберутся вместе в единое решение проблемы.
Разложение синуса в ряд Тейлора.
Функция синуса раскладывается в бесконечный ряд Тейлора.
$$display$$\sin(x)=x-{\frac {x^{3}}{3!}}+{\frac {x^{5}}{5!}}-{\frac {x^{7}}{7!}}+{\frac {x^{9}}{9!}}-\cdots $$display$$
Понятно, что бесконечный ряд мы посчитать не можем, кроме случаев, когда есть аналитическая формула бесконечной суммы. Но это не наш случай))) Предположим, что мы хотим посчитать синус в интервале $inline$[0, \frac{\pi}{2}]$inline$. Более подробно работу с интервалами обсудим в части 3. Зная, что $inline$\sin(\frac{\pi}{2})=1$inline$ оценим найдём первый член который можно отбросить исходя из условия, что $inline$\frac{(\pi/2)^n}{n!}<e$inline$, где $inline$e$inline$ это разница между числом 1 и наименьшем числом, которое больше 1. Грубо говоря это последний бит мантиссы (wiki). Решить данное уравнение проще перебором. Для $inline$e \approx 2.22\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;}
Как ускорить программу я думаю, что многие сообразили сразу. Как вы думаете, во сколько раз ваши изменения могут ускорить программу? Оптимизированная версия и ответ на вопрос под спойлером.
Оптимизированная версия программы.
Ускорение больше чем в 10000 раз (GNU C++ v10; -O2)
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;}
Ускорение больше чем в 10000 раз (GNU C++ v10; -O2)
Улучшение точности
Методика
Точность вычисления функции будем определять 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;}
Результаты приведены в табличке.
| Функция | Среднее ошибки | STD | Лучше наша | Лучше libm |
| sin_e1 | -1.28562e-18 | 8.25717e-17 | 0.0588438% | 53.5466% |
| sin_e3 | -3.4074e-21 | 3.39727e-17 | 0.0423% | 10.8049% |
| sin_e4 | 8.79046e-18 | 4.77326e-17 | 0.0686% | 27.6594% |
| sin_e5 | 8.78307e-18 | 3.69995e-17 | 0.0477062% | 13.5105% |
Может показаться, что использование алгоритмов умного суммирования уберёт ошибку практически до 0, но это не так. Конечно эти алгоритмы дадут увеличение точности, но для полного избавления от ошибок требуются ещё и алгоритмы умного умножения. Они существуют, но очень накладны: очень много лишних операций. Применение их здесь не оправдано. Впрочем позднее мы к ним вернёмся в другом контексте.
Осталось совсем немного. Объединить быстрый и точный алгоритмы. Для этого снова вернёмся к ряду Тейлора. Ограничем его для примера 4-мя членами и сделаем следующее преобразование.
$$display$$\sin(x)\approx x(1+x^2(-1/3!+x^2(1/5!+x^2(-1/7!+x^2\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$$\sin(x)\approx x+x \cdot x^2(-1/3!+ \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.
