Числовая точность для разности квадратов

Наиболее точный способ с числами с плавающей запятой, вероятно, заключается в вычислении как sin, так и cos с помощью одной инструкции x87, fsincos.

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

В статье Что должен знать каждый компьютерный ученый об арифметике с плавающей запятой примечания:

Выражение x² - y² — это еще одна формула, демонстрирующая катастрофическую отмену. Более точно вычислить его как (x - y)(x + y).

Изменить: это сложнее, чем это. Хотя приведенное выше в целом верно, (x - y)(x + y) немного менее точно, когда x и y имеют очень разные величины, как поясняется в сноске к утверждению:

В этом случае (x - y)(x + y) имеет три ошибки округления, а x² - y² имеет только две, так как ошибка округления допущена при вычислении меньшее из x² и y² не влияет на окончательное вычитание.

Другими словами, взятие x - y, x + y и произведения (x - y)(x + y) приводит к ошибкам округления (3 шага ошибки округления). x², y² и вычитание x² - y² также приводят к ошибкам округления, но ошибка округления, полученная путем возведения в квадрат относительно небольшого числа (меньшего из x и y), настолько незначительна, что фактически существует только два шага ошибки округления, что делает разницу квадратов более точной.

Так что вариант 1 на самом деле будет более точным. Это подтверждается тестом Java от dev.brutus.

Я написал небольшой тест. Он вычисляет ожидаемое значение с двойной точностью. Затем он вычисляет ошибку с вашими параметрами. Первый вариант лучше:

Algorithm: FloatTest$1
option 1 error = 3.802792362162126
option 2 error = 4.333273185303996
Algorithm: FloatTest$2
option 1 error = 3.802792362167937
option 2 error = 4.333273185305868

Java-код:

import org.junit.Test;

public class FloatTest {

    @Test
    public void test() {
        testImpl(new ExpectedAlgorithm() {
            public double te(double cos_theta) {
                return Math.sqrt(1.0f - cos_theta * cos_theta);
            }
        });
        testImpl(new ExpectedAlgorithm() {
            public double te(double cos_theta) {
                return Math.sqrt((1.0f + cos_theta) * (1.0f - cos_theta));
            }
        });
    }

    public void testImpl(ExpectedAlgorithm ea) {
        double delta1 = 0;
        double delta2 = 0;
        for (double cos_theta = -1; cos_theta <= 1; cos_theta += 1e-8) {
            double[] delta = delta(cos_theta, ea);
            delta1 += delta[0];
            delta2 += delta[1];
        }

        System.out.println("Algorithm: " + ea.getClass().getName());
        System.out.println("option 1 error = " + delta1);
        System.out.println("option 2 error = " + delta2);
    }

    private double[] delta(double cos_theta, ExpectedAlgorithm ea) {
        double expected = ea.te(cos_theta);
        double delta1 = Math.abs(expected - t1((float) cos_theta));
        double delta2 = Math.abs(expected - t2((float) cos_theta));

        return new double[]{delta1, delta2};
    }

    private double t1(float cos_theta) {
        return Math.sqrt(1.0f - cos_theta * cos_theta);
    }

    private double t2(float cos_theta) {
        return Math.sqrt((1.0f + cos_theta) * (1.0f - cos_theta));
    }

    interface ExpectedAlgorithm {
        double te(double cos_theta);
    }

}

Правильный способ рассуждать о числовой точности некоторого выражения состоит в следующем:

1. Измерьте несоответствие результата относительно правильного значения в ULP (единица на последнем месте), представленная в 1960 г. У. Х. Кахан. Вы можете найти реализации C, Python и Mathematica здесь и узнать больше о здесь.
2. Различайте два или более выражений на основе наихудшего случая, который они производят, а не по средней абсолютной ошибке, как это делается в других ответах, или по какой-либо другой произвольной метрике. Вот как строятся полиномы численной аппроксимации (алгоритм Ремеза), как анализируются реализации стандартных библиотечных методов (например, Intel atan2) , так далее...

Имея это в виду, version_1: sqrt(1 - x * x) и version_2: sqrt((1 - x) * (1 + x)) дают существенно разные результаты. Как показано на графике ниже, версия_1 демонстрирует катастрофическую производительность для x, близкого к 1, с ошибкой > 1_000_000 ulps, в то время как, с другой стороны, ошибка версии_2 ведет себя хорошо.

Вот почему я всегда рекомендую использовать версию_2, то есть использовать формулу квадратной разницы.

Код Python 3.6, создающий файл square_diff_error.csv:

from fractions import Fraction
from math import exp, fabs, sqrt
from random import random
from struct import pack, unpack


def ulp(x):
    """
    Computing ULP of input double precision number x exploiting
    lexicographic ordering property of positive IEEE-754 numbers.

    The implementation correctly handles the special cases:
      - ulp(NaN) = NaN
      - ulp(-Inf) = Inf
      - ulp(Inf) = Inf

    Author: Hrvoje Abraham
    Date: 11.12.2015
    Revisions: 15.08.2017
               26.11.2017
    MIT License https://opensource.org/licenses/MIT

    :param x: (float) float ULP will be calculated for
    :returns: (float) the input float number ULP value
    """

    # setting sign bit to 0, e.g. -0.0 becomes 0.0
    t = abs(x)

    # converting IEEE-754 64-bit format bit content to unsigned integer
    ll = unpack('Q', pack('d', t))[0]

    # computing first smaller integer, bigger in a case of ll=0 (t=0.0)
    near_ll = abs(ll - 1)

    # converting back to float, its value will be float nearest to t
    near_t = unpack('d', pack('Q', near_ll))[0]

    # abs takes care of case t=0.0
    return abs(t - near_t)


with open('e:/square_diff_error.csv', 'w') as f:
    for _ in range(100_000):
        # nonlinear distribution of x in [0, 1] to produce more cases close to 1
        k = 10
        x = (exp(k) - exp(k * random())) / (exp(k) - 1)

        fx = Fraction(x)
        correct = sqrt(float(Fraction(1) - fx * fx))

        version1 = sqrt(1.0 - x * x)
        version2 = sqrt((1.0 - x) * (1.0 + x))

        err1 = fabs(version1 - correct) / ulp(correct)
        err2 = fabs(version2 - correct) / ulp(correct)

        f.write(f'{x},{err1},{err2}\n')

Код Mathematica, который создает окончательный график:

data = Import["e:/square_diff_error.csv"];

err1 = {1 - #[[1]], #[[2]]} & /@ data;
err2 = {1 - #[[1]], #[[3]]} & /@ data;

ListLogLogPlot[{err1, err2}, PlotRange -> All, Axes -> False, Frame -> True,
    FrameLabel -> {"1-x", "error [ULPs]"}, LabelStyle -> {FontSize -> 20}]

Кроме того, у вас всегда будут проблемы, когда тета мала, потому что косинус плоский вокруг тета = 0. Если тета находится в диапазоне от -0,0001 до 0,0001, тогда cos(theta) в float ровно единице, поэтому ваш sin_theta будет точно равен нуль.

Чтобы ответить на ваш вопрос, когда cos_theta близок к единице (соответствует малому тета), ваше второе вычисление явно более точное. Это показано в следующей программе, в которой перечислены абсолютные и относительные ошибки для обоих вычислений для различных значений cos_theta. Ошибки вычисляются путем сравнения со значением, которое вычисляется с точностью до 200 бит с использованием библиотеки GNU MP, а затем преобразуется в число с плавающей запятой.

#include <math.h>
#include <stdio.h>
#include <gmp.h>

int main() 
{
  int i;
  printf("cos_theta       abs (1)    rel (1)       abs (2)    rel (2)\n\n");
  for (i = -14; i < 0; ++i) {
    float x = 1 - pow(10, i/2.0);
    float approx1 = sqrt(1 - x * x);
    float approx2 = sqrt((1 - x) * (1 + x));

    /* Use GNU MultiPrecision Library to get 'exact' answer */
    mpf_t tmp1, tmp2;
    mpf_init2(tmp1, 200);  /* use 200 bits precision */
    mpf_init2(tmp2, 200);
    mpf_set_d(tmp1, x);
    mpf_mul(tmp2, tmp1, tmp1);  /* tmp2 = x * x */
    mpf_neg(tmp1, tmp2);        /* tmp1 = -x * x */
    mpf_add_ui(tmp2, tmp1, 1);  /* tmp2 = 1 - x * x */
    mpf_sqrt(tmp1, tmp2);       /* tmp1 = sqrt(1 - x * x) */
    float exact = mpf_get_d(tmp1);

    printf("%.8f     %.3e  %.3e     %.3e  %.3e\n", x,
           fabs(approx1 - exact), fabs((approx1 - exact) / exact),
           fabs(approx2 - exact), fabs((approx2 - exact) / exact));
    /* printf("%.10f  %.8f  %.8f  %.8f\n", x, exact, approx1, approx2); */
  }
  return 0;
}

Выход:

cos_theta       abs (1)    rel (1)       abs (2)    rel (2)

0.99999988     2.910e-11  5.960e-08     0.000e+00  0.000e+00
0.99999970     5.821e-11  7.539e-08     0.000e+00  0.000e+00
0.99999899     3.492e-10  2.453e-07     1.164e-10  8.178e-08
0.99999684     2.095e-09  8.337e-07     0.000e+00  0.000e+00
0.99998999     1.118e-08  2.497e-06     0.000e+00  0.000e+00
0.99996835     6.240e-08  7.843e-06     9.313e-10  1.171e-07
0.99989998     3.530e-07  2.496e-05     0.000e+00  0.000e+00
0.99968380     3.818e-07  1.519e-05     0.000e+00  0.000e+00
0.99900001     1.490e-07  3.333e-06     0.000e+00  0.000e+00
0.99683774     8.941e-08  1.125e-06     7.451e-09  9.376e-08
0.99000001     5.960e-08  4.225e-07     0.000e+00  0.000e+00
0.96837723     1.490e-08  5.973e-08     0.000e+00  0.000e+00
0.89999998     2.980e-08  6.837e-08     0.000e+00  0.000e+00
0.68377221     5.960e-08  8.168e-08     5.960e-08  8.168e-08

Когда cos_theta не близок к единице, точность обоих методов очень близка друг к другу и к ошибке округления.

[Отредактировано для серьезных размышлений] Мне кажется, что вариант 2 будет лучше, потому что для числа, такого как 0.000001, например, вариант 1 вернет синус как 1, а вариант вернет число чуть меньше 1.

Нет разницы в моем варианте, поскольку (1-x) сохраняет точность, не влияя на переносимый бит. Тогда для (1+x) то же верно. Тогда единственное, что влияет на точность бита переноса, — это умножение. Таким образом, в обоих случаях есть одно единственное умножение, поэтому они оба с одинаковой вероятностью дадут одну и ту же ошибку бита переноса.