В качестве примера я выбрал тригонометрическую функцию синуса. Это широко распространённая функция, математика которой хорошо известна со школы и университета. В тоже время при её имплементации появятся много ярких примеров «правильной» работы с числами. В качестве числа с плавающей точкой я буду использовать double.
В данном цикле статей планируется много всего начиная от математики, заканчивая машинными кодами и опциями компилятора. Язык написания статьи С++, но без «излишеств». В отличии от языка С, работающие примеры будут более удобочитаемыми даже для людей не знакомых с этим языком и занимать меньше строк.
Статьи будут написаны методом погружения. Будут обсуждаться подзадачи, которые потом соберутся вместе в единое решение проблемы.
Разложение синуса в ряд Тейлора.
Функция синуса раскладывается в бесконечный ряд Тейлора.
Понятно, что бесконечный ряд мы посчитать не можем, кроме случаев, когда есть аналитическая формула бесконечной суммы. Но это не наш случай))) Предположим, что мы хотим посчитать синус в интервале . Более подробно работу с интервалами обсудим в части 3. Зная, что оценим найдём первый член который можно отбросить исходя из условия, что , где это разница между числом 1 и наименьшем числом, которое больше 1. Грубо говоря это последний бит мантиссы (wiki). Решить данное уравнение проще перебором. Для . У меня получилось уже можно отбросить. Правильный выбор количества слагаемых будет обсужден в одной из следующей частей, поэтому на сегодня «перестрахуемся» и возьмём слагаемые до включительно.
Последнее слагаемое приблизительно в 10000 раз меньше, чем .
Простейшее решение
Руки уже чешутся, пишем:
#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, 26> 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;
}
inline
double fsin(double x) {
double result;
asm ("fsin" :"=t" (result) : "0" (x));
return result;
}
#define SIN(a) fsin(a)
//#define SIN(a) sin(a)
//#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());
}
double res = 0;
{
boost::timer::auto_cpu_timer at;
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 << res << endl;
// Количество случаев, когда внутренняя функция дала более правильный результат к общему количеству вычислений.
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;
}
Ускорение больше чем в 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-мя членами и сделаем следующее преобразование.
Можно раскрыть скобки и проверить, что получится исходное выражение. Такое представление очень просто ложится на цикл.
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. Опять же проблема в округлении. Давайте рассмотрим последний шаг цикла и немного преобразуем исходное выражение
И соответствующий код.
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.
firewind1
Круто конечно, но где нужно на практике знать синус с такой потрясающей точностью?
dorex Автор
Хороший вопрос.
Во-первых хочу заметить, что одна из целей статьи показать, что от «перемены мест слагаемых сумма меняется» и показать методы улучшения точности. Если можно вычислить точнее в одно и тоже время, почему так не сделать?
Ну а если говорить по синус, то частый пример — численное дифференцирование функции. Оно требуется в решении дифференциальных уравнений и оптимизации функции. Самый простой пример. cos(x) = (sin(x + dx) — sin(x)/dx. Ясно, что при уменьшении dx, ошибка в последнем знаке будет давать гораздо большую ошибку функции cos. Вычислении численного градиента, например, используется при минимизации функций. Минимизация с неточным градиентом может потребовать больше шагов и/или привести к неточному решению. Вычисление производных более высокого порядка, которое может быть нужно для численного решения дифференциальных уравнений, даст ещё большую ошибку. К тому же для многих уравнений ошибка накапливается и спустя какое-то время становится существенной.
Из самых простых примеров дифференциальных уравнений — это обычный маятник. Только в первом приближении уравнение его движения линейно (Маятник).
Если Вас не интересует точность double, то, конечно, можно использовать числа float. Для них все перечисленные методы работают абсолютно так же. А последний знак там всего лишь 0.000 000 119
dorex Автор
P.S. Большая точность в расчёте sin далеко не всегда. Тему быстрых неточных вычислений я раскрою в одной из следующих статей.
firewind1
Спасибо, буду ждать :)
А насчет точности, просто со времен курсового проекта по термодинамике приучили, что хватает четырех значащих цифр после запятой, вот и проецирую на всё.
maxzhurkin
нужна не всегда?
dorex Автор
Забегу вперёд. В компиляторах GNU есть опция fast-math. При её включении используется код, который считает синус очень быстро, но с точностью 4 ULP. glibc source
Я подробно буду обсуждать это в одной из следующих статей.
jorikdima
Там где может происходить накопление ошибки, ну, например, суммирование таких вот синусов. Там нужно контролировать накопленную ошибку и минимизировать ее.