Каков правильный алгоритм для выполнения двойного деления?

Question

Я следую алгоритмам, предоставленным в этой статье Эндрю Таллом , описывающим алгоритмы для выполнения математических операций с использованием типа данных df64, пары 32-битных чисел с плавающей запятой, используемой для имитации точности 64 -битное число с плавающей запятой. Однако, как представляется, существуют некоторые несоответствия (ошибки?) В том, как они написали свои функции деления и квадрата Root.

Вот как в статье написана функция деления:

float2 df64_div(float2 B, float2 A) {
    float xn = 1.0f / A.x;
    float yn = B.x * xn;
    float diff = (df64_diff(B, df64_mult(A, yn))).x;
    float2 prod = twoProd(xn, diffTerm);

    return df64_add(yn, prodTerm);
}

^{Язык, используемый для написания этого кода, кажется, Cg, для справки, хотя вы должны иметь возможность интерпретировать этот код в C ++, если вы обрабатываете float2, как будто это просто псевдоним для struct float2{float x, y;};, с некоторым дополнительным синтаксисом для поддержки арифметических c операций непосредственно над типом.}

Для справки, это заголовки функций, используемых в этом коде:

float2 df64_add(float2 a, float2 b);
float2 df64_mult(float2 a, float2 b);
float2 df64_diff(/*Not provided...*/);
float2 twoProd(float a, float b);

Таким образом, сразу несколько проблем выделяются:

• diffTerm и prodTerm никогда не определяются. Есть две переменные, diff и prod, которые определены , но нет уверенности, что это термины, которые были предназначены в этом коде.
• Нет объявления df64_diff при условии. Предположительно это предназначено для поддержки вычитания; но опять же, это не ясно.
• df64_mult - это функция, которая не принимает 32-битное значение с плавающей точкой в ​​качестве аргумента; он поддерживает только две пары 32-битных чисел с плавающей точкой в ​​качестве аргументов. Непонятно, каким образом документ ожидает, что этот вызов функции скомпилирует
• То же самое и для df64_add, который также принимает в качестве аргументов только пары 32-битных чисел с плавающей точкой, но здесь он вызывается только с первым аргументом, являющимся один 32-разрядный тип с плавающей точкой.

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

float2 df64_div(float2 B, float2 A) {
    float xn = 1.0f / A.x;
    float yn = B.x * xn;
    float diff = (df64_diff(B, df64_mult(A, float2(yn, 0)))).x;
    float2 prod = twoProd(xn, diff);

    return df64_add(float2(yn, 0), prod);
}

float2 df64_diff(float2 a, float2 b) {
    return df64_add(a, float2(-b.x, -b.y));
}

Поэтому мой вопрос заключается в следующем: точна ли письменная реализация этого алгоритма, как видно из статьи (потому что это зависит о поведении языка Cg, о котором я не знаю?) или нет? И независимо от того, является ли моя интерполяция этого кода правильной реализацией алгоритма деления, описанного в статье?

Примечание: мой целевой язык - C ++, поэтому, хотя различия между языками (для этого вида алгоритма) второстепенный, мой код написан на C ++, и я ищу правильность для языка C ++.

Answer 1

Ответ Xirema обеспечивает точную визуализацию высоко-радикального алгоритма деления длинной руки Талла в C ++. Основываясь на довольно широкомасштабном тестировании по эталону с более высокой точностью, я обнаружил, что его максимальная относительная погрешность составляет порядка 2 ^-45 , при условии, что в промежуточных вычислениях нет недостатков.

Вкл. на платформах, которые обеспечивают операцию смешанного сложения ( FMA ), следующий алгоритм деления на основе Ньютона-Рафсона из-за Nagai et. al. , вероятно, будет более эффективным и достигает идентичной точности в моем тестировании, то есть максимальная относительная погрешность составляет 2 ^-45 .

/*
  T. Nagai, H. Yoshida, H. Kuroda, Y. Kanada, "Fast Quadruple Precision 
  Arithmetic Library on Parallel Computer SR11000/J2." In: Proceedings 
  of the 8th International Conference on Computational Science, ICCS '08, 
  Part I, pp. 446-455.
*/
float2 div_df64 (float2 a, float2 b)
{
    float2 t, c;
    float r, s;
    r = 1.0f / b.x;
    t.x = a.x * r;
    s = fmaf (-b.x, t.x, a.x);
    t.x = fmaf (r, s, t.x);
    t.y = fmaf (-b.x, t.x, a.x);
    t.y = a.y + t.y;
    t.y = fmaf (-b.y, t.x, t.y);
    s = r * t.y;
    t.y = fmaf (-b.x, s, t.y);
    t.y = fmaf (r, t.y, s);
    c.x = t.x + t.y;
    c.y = (t.x - c.x) + t.y;
    return c;
}

Answer 2

Просмотр алгоритма псевдокода, как написано в книге, похоже, поддерживает реализацию этого алгоритма на C ++, хотя мое незнакомство с Cg означает, что я не могу доказать, что эта реализация верна для Cg.

Таким образом, разбивая эти шаги на простой engli sh:

1. Функция принимает два параметра, каждый из которых [псевдо-] с плавающей запятой двойной точности точечные значения, а второй параметр не равен 0
2. Переменной x _n присваивается арифметическая c величина, обратная компоненту высшего порядка двойного делителя [псевдо-] , рассчитывается с использованием математических операций с плавающей запятой одинарной точности
3. Переменной y _n присваивается произведение компонента высшего порядка [псевдо-] двойного дивиденда и x _n, рассчитывается с использованием математических вычислений с плавающей запятой одинарной точности
4. Произведение [псевдо-] двойного делителя и y _n вычисляется
   • Это Первая сложная часть, потому что в статье не описан алгоритм для [псевдо-] двойного х однократного умножения. Из алгоритма Cg видно, что алгоритм Cg четко отображается на этот шаг 1-к-1, но правила Cg для преобразования скалярного значения в векторное значение неизвестны.
   • Что мы можем сказать, однако является то, что у нас есть функция для умножения двойного на двойное, и одиночный может быть повышен до двойного, заполнив его младший компонент 0, так что мы можем это сделать.
5. Рассчитывается разница между дивидендом и произведением, рассчитанным на шаге 4, и в качестве значения с плавающей запятой одинарной точности сохраняется только компонент более высокого порядка
   • . опишите алгоритм вычитания. Тем не менее, он описывает алгоритм для преобразования двойника [IEEE754-] в [псевдо-] двойное, и мы можем сделать замечание, что отрицательные двойники [IEEE754-] при преобразовании имеют отрицательные значения как для его более высокого порядка, так и для более высокого порядка. компоненты более низкого порядка. Логично, что псевдо-двойник можно отрицать, просто отрицая обе его составляющие. И добавленное отрицательное число математически эквивалентно вычитанию, поэтому мы можем построить алгоритм вычитания, используя эти знания.
6. Произведение x _n и шаг 5 выполнен сохраняя расширенную точность, которая в противном случае была бы потеряна при однократном умножении на единицу.
   • Для этой цели существует функция twoProd.
7. Сумма шага 6 и y _n рассчитывается
   • Опять же, мы можем использовать алгоритм [псевдо-] двойного сложения, если просто повышаем y _n до [псевдо-] двойного, заполняя компонент младшего разряда 0
8. Результатом шага 7 является возвращаемое значение

Таким образом, понимая этот алгоритм, мы можем сопоставить каждый из этих шагов непосредственно с алгоритмом C ++, который я написал:

//(1) Takes two [pseudo-]doubles, returns a [pseudo-]double
float2 df64_div(float2 B, float2 A) {
    //(2) single float divided by single float
    float xn = 1.0f / A.x;
    // (3) single float multiplied by single float
    float yn = B.x * xn;
    //                        (4) double x double multiplication
    //                                       (4a) yn promoted to [pseudo-]double
    //            (5) subtraction                           (5a) only higher order component kept
    float diff = (df64_diff(B, df64_mult(A, float2(yn, 0)))).x;
    // (6) single x single multiplication with extra precision preserved using twoProd
    float2 prod = twoProd(xn, diff);
    // (7) adding higher-order division to lower order division
    //              (7a) yn promoted to [pseudo-]double
    // (8) value is returned
    return df64_add(float2(yn, 0), prod);
}

float2 df64_diff(float2 a, float2 b) {
    //                 (5a) negating both components is a logical negation of the whole number
    return df64_add(a, float2(-b.x, -b.y));
}

Из Таким образом, мы можем сделать вывод, что это правильная реализация алгоритма, описанного в этой статье, подтвержденная некоторыми тестами, которые я провел для проверки того, что выполнение этих операций таким образом дает результаты, которые кажутся правильными.