How to become an author
Search
Write a publication
Pull to refresh

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

Reading time6 min
Views28K
Algorithms*Mathematics*Programming*
From sandbox
Внимательно прочитал очень хорошие статьи от ArtemKaravaev по сложению чисел с плавающей точкой. Тема очень интересная и хочется её продолжить и показать на примерах, как работать с числами с плавающей точкой на практике. В качестве эталона возьмём библиотеку GNU glibc (libm). А чтобы статья не была уж скучной, добавим соревновательную составляющую: попробуем не только повторить, но и улучшить код библиотеки, сделав его более быстрым/точным.

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

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

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

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


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

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


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

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


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

Полный текст программы для тестирования 
#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-мя членами и сделаем следующее преобразование.

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



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

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. Опять же проблема в округлении. Давайте рассмотрим последний шаг цикла и немного преобразуем исходное выражение

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



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

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.
Tags:
Hubs:
Total votes 40: ↑39 and ↓1+59
Comments58
17
Karma
0
Rating
@dorex

User

Send message

Articles

Your account

Sections

Information

Services

Support
© 2006–2024, Habr