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

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

В качестве примера я выбрал тригонометрическую функцию синуса. Это широко распространённая функция, математика которой хорошо известна со школы и университета. В тоже время при её имплементации появятся много ярких примеров «правильной» работы с числами. В качестве числа с плавающей точкой я буду использовать 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.22times10^{-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$$sin(x)approx x(1+x^2(-1/3!+x^2(1/5!+x^2(-1/7!+x^2cdot1/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.

Автор: dorex

Источник