Можно ли сложить N чисел типа double наиболее точно?

    В предыдущих сериях…


    Прошлая статья рассказала о двух способах сложения двух двоичных чисел с плавающей запятой без потери точности. Чтобы добиться этого, мы представили сумму c=a+b в виде двух чисел (s,t)=a+b, причём таких, что s — наиболее близкое к a+b точно-представимое число, а t=(a+b)-s — это отсекаемая в результате округления часть, составляющая точную погрешность. У читателей был вопрос: а можно ли достаточно точно сложить массив чисел типа double? Оказывается, можно! Но только, вероятно, не всегда и не абсолютно… и не алгоритмом Кэхэна, который тогда вспоминали в комментариях. За подробностями прошу под кат, где мы и найдём приложение тому, о чём я рассказал в прошлый раз.




    Вступление


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

    Традиционно, полное содержание статьи также доступно в форме презентации с закадровым голосом:



    Дело в том, что абсолютно точно сложить N чисел типа double (binary64 в IEEE-754) с сохранением суммы в типе double не получится в силу самой специфики формата с плавающей запятой. Если мы возьмём арифметику с плавающей запятой с бесконечной точностью (или даже воспользуемся рациональными числами), то при достаточном объёме памяти получим абсолютно точную сумму S массива чисел X[N]. Но эту абсолютно точную сумму нужно будет всё равно округлить к ближайшему числу типа double. Назовём эту округлённую сумму S’. Если вспомнить функцию RN(x) из предыдущей статьи, которая выполняет округление произвольного вещественного x «к ближайшему чётному», то можно сказать, что S’=RN(S), поэтому |S-S’| ≤ 0.5ulp (т. е. ошибка составляет не более половины значения младшего бита). Таким образом, об абсолютной точности мы говорить не можем, но мы будем называть число S’ — «наиболее точной суммой» (из возможных). Поэтому переформулируем наш вопрос иначе: можно ли сложить числа X[N] так, чтобы получить наиболее точную сумму S’, не пользуясь длинной арифметикой, а оставаясь лишь в рамках операций сложения с double?

    Да! На случайном наборе чисел этого, вероятно, можно добиться. Я покажу один алгоритм, который мне не удалось «завалить», хотя в теории, наверное, можно придумать какие-то особенные тесты, чтобы он ошибся хотя бы на 1ulp и выдал бы неправильно хотя бы последний бит результата. Всего мы рассмотрим три алгоритма и сравним их между собой через абсолютно точную сумму S, округлённую до наиболее точной суммы S’.

    Нам понадобятся два алгоритма из предыдущей статьи, и сейчас мы должны дать им названия. Я буду пользоваться синтаксисом C++, уверен, это не вызовет трудностей с переводом у любого читателя.

    Первый алгоритм для двух чисел


    Напомню, эти алгоритмы возвращают два числа: наиболее близкую к реальной сумму s=RN(a+b) и точную погрешность t=(a+b)-s. Данный алгоритм работает правильно только при условии |a|≥|b|, если сумма RN(a+b)≠∞.

    s = RN (a+b);
    z = RN (s-a);
    t = RN (b-z);
    return (s, t);

    Я запрограммировал его на C++ вот так:

    double __fastcall fast_two_sum (double &t, double a, double b) {
      double s = a+b;
      t = b-(s-a);
      return s;
    }

    Не забудьте, что при компилировании подобного кода нельзя применять флаги агрессивной оптимизации вычислений с плавающей запятой, потому что переставлять порядок вычислений здесь категорически нельзя. Также мне неизвестно, будут ли работать подобные конструкции на GPU, так как раньше (по крайней мере, на заре их применения) они могли обрабатывать числа с плавающей запятой даже близко не по Стандарту. О текущем состоянии дел мне неизвестно, прошу знающих читателей подсказать в комментариях.

    Второй алгоритм для двух чисел


    Этот алгоритм работает при любом соотношении между $a$ и $b$, но по-прежнему RN(a+b)≠∞. Вместо 3-х операций здесь их 6.

    s = RN (a+b);
    b' = RN (s-a);
    a' = RN (s-b');
    d_b = RN (b-b');
    d_a = RN (a-a');
    t = RN (d_a+d_b);
    return (s, t)

    В моём исполнении на C++ выглядит так:

    double __fastcall two_sum (double &t, double a, double b) {
      double s = a+b;
      double bs = s-a;
      double as = s-bs;
      t = (b-bs) + (a-as);
      return s;
    }
    

    Теперь рассмотрим какие бывают варианты сложения нескольких чисел из массива. Сразу предупреждаю, что их очень много, но я взял два универсальных и наиболее интересных из книги [1], которые действительно хорошо себя показывают. На остальные алгоритмы, если вам интересно, авторы упомянутой книги дают ряд ссылок.

    Обзор алгоритмов сложения


    Задача. Задан массив X[N] типа double, требуется отыскать сумму всех чисел в этом массиве, сумма имеет тот же тип double.

    Ниже по тексту будут представлены мои варианты реализации алгоритмов с описанием (где оно требуется). Все реализации будут опираться на псевдоним для массива чисел типа double:

    using dvect_t = vector<double>;

    Алгоритм 0. Абсолютно точная и наиболее точная сумма


    Этот алгоритм нужен только для проверки остальных, его практический смысл лично для меня почти полностью отсутствует. Чтобы посчитать сумму чисел точно, нужна арифметика с бесконечной точностью. Это может быть длинная арифметика с плавающей запятой, а может быть длинная арифметика рациональных чисел. Я выбрал последнюю. Взял библиотеку MPIR (можно было взять и GMP) и тип данных mpq_t. Написанный в скрытом тексте код выдаёт наиболее точную сумму S’, получаемую из абсолютно точной суммы S путём правильного округления. Пропустите скрытый текст, если вы не знакомы с внутренним устройством формата double, это не помешает вам понять остальной материал.

    Код алгоритма 0 и пояснения
    Если нет цели получить максимальную скорость вычислений, то код может быть написан очень прямолинейно, однако далее возникает проблема. Функция mpq_get_d из библиотеки MPIR округляет число вниз, просто отрезая «лишнее» (а другой функции округления там нет). Нам нужно округлить по правилу «к ближайшему чётному», а для этого приходится совершать «танцы с бубном» над битовым представлением чисел типа double. Поэтому в коде ниже есть UB, будьте внимательны.

    
    using u64 = unsigned long long;
    #define DAU(x) (*(u64*)(&x))  // Это UB, друзья! Перепишите как вам правильнее.
    
        // Наиболее близкая сумма.
        // Суммирование через рациональные дроби с бесконечной точностью.
        // Внимание, здесь танцы с UB! Проверьте, прежде чем запускать.
    double __fastcall sum_exact (const dvect_t &X) {
      mpq_t s, a, b;
      mpq_inits (s, a, b, NULL);
    
      // Сначала нужно всё сложить
      mpq_set_d (s, 0.0); // Общая сумма дробей s = 0
      for (auto x: X) {
        mpq_set_d (a, x);  // Переводим 'x' к дроби 'a'
        mpq_add (s, s, a); // Складываем s += a.
      }
    
      // Теперь получаем точно-округлённый результат.
      double res = mpq_get_d (s); // Здесь округление к нулю! Это ещё не точный результат!
      
      // Достаём "бубен"...
    
      // Нужно на время убрать знак "минус"
      bool is_negative = res < 0.0;
      res = abs(res);
      mpq_abs (s, s);
      u64 res64 = DAU(res); // Получить битовое представление числа типа double.
      u64 res64_next = res64+1;  // Получить следующее за 'res' число в битовой форме
      double res_next;
      DAU(res_next) = res64_next;  // 'res_next' - это число, следующее за 'res'
    
      // Начинаем "танец"...
    
      mpq_set_d (a, res);      // a - точное дробное представление res
      mpq_set_d (b, res_next); // b - точное дробное представление res_next
      // В этот момент наша сумма s где-то внутри отрезка [a, b]=[res, res_next].
      mpq_sub (a, s, a);  // a - расстояние от res до s
      mpq_sub (b, b, s);  // b - расстояние от s до res_next
    
      // Что больше: расстояние от res до s или от s до res_next?
      if (mpq_cmp (a, b) > 0)
        res = res_next; // Если наша s ближе к правой границе отрезка.
    
      if ((res64&1) && mpq_cmp (a, b) == 0)
        res = res_next; // Если s по центру, но младший бит левой границы нечётный.
    
      mpq_clears (s, a, b, NULL);
      return is_negative ? -res : res;
    }
    


    Алгоритм 1. Наивное суммирование


    Очевидный алгоритм и, пожалуй, наиболее быстрый и просторный для оптимизаций компилятора выглядит вот так:

    double __fastcall sum_naive (const dvect_t &X) {
      double s = 0.0;
      for (auto x: X)  s += x;
      return s;
    }
    

    Здесь нечего пояснять.

    Алгоритм 2. Kahan


    double __fastcall sum_kahan (const dvect_t &X) {
      double s=0.0, c=0.0;  
      for (auto x: X) {
        double y = x + c;
        s = fast_two_sum (c, s, y);
      }
      return s;
    }
    

    Схема работы алгоритма основана на первом методе сложения двух чисел без потери точности из предыдущей статьи. «Хвостик» $c$, который получается в результате сложения текущей суммы s и очередного слагаемого, добавляется к следующему слагаемому (у=x+c), как бы компенсируя потери. Иными словами, мы не выбрасываем потерянные биты после операции сложения (s+y) на каждом шаге, а пытаемся «спасти» их, добавляя их на следующем шаге.

    Внимательный читатель вспомнит, что алгоритм s=fast_two_sum(c, s, y) работает корректно только когда |s|≥|y|, что не обязательно будет именно так для беспорядочной последовательности чисел из исходного массива X. Тем не менее, компенсационные способности данного алгоритма перекрывают ошибочное срабатывание для случаев |s|<|y|, по какой причине на практике этот алгоритм чаще всего работает лучше наивного суммирования. «Завалить» его именно на случайных тестах мне не удалось, хотя подобрать специальную последовательность вполне можно (см. последнюю таблицу в статье). Также в единственной книге, на которую я ссылаюсь, приводится пример:

      
      X[0] = 18014398509481984.0; // 2**54
      X[1] = 18014398509481982.0; // 2**54-2
      X[2] = -9007199254740991.0; // -(2**53-1)
      X[3] = -9007199254740991.0; 
      X[4] = -9007199254740991.0;
      X[5] = -9007199254740991.0;
    

    Правильный ответ 2. Алгоритм Кэхэна вернёт 3. Наивное суммирование вернёт 1. Согласитесь, это чудовищная ошибка в 50% (или 251ulp). А вот следующий алгоритм выдаст в этом случае правильный ответ.

    Алгоритм 3. Rump–Ogita–Oishi


    Разумный вопрос: а что если вместо первого алгоритма сложения двух чисел из предыдущей статьи взять второй, который безразличен к порядку этих чисел? Этот подход применяется в алгоритме Rump–Ogita–Oishi и даёт такой код:

    double __fastcall sum_rump (const dvect_t &X) {
      double s=0.0, c=0.0, e;
      for (double x: X) {
        s = two_sum (e, s, x);
        c += e;
      }
      return s+c;
    }
    

    Однако здесь авторы алгоритма решили складывать «хвостики» $e$ в отдельную переменную, что, в общем-то, логично. В конце выполняем сложение итоговой суммы и суммы всех «хвостиков». Данный алгоритм можно представить себе так, как будто мы суммируем числа не в одной переменной типа double, а как бы удваиваем её, получая своего рода «double-double», в котором не 52, а 104 бита дробной части мантиссы. Чтобы «сломать» такой алгоритм, нужно чтобы погрешность «закрыла» 52 бита. Сделать это обычными бытовыми задачами практически нереально, если только вы не складываете триллионы чисел. Впрочем, это лишь домыслы автора текста. Если можете «завалить» алгоритм, прошу показать тест в комментариях (см. этот комментарий, но это не бытовой тест).

    Тестирование


    Замечание о методике и способе вывода погрешности


    Я прошу прощения за многословное описание методики тестирования и преподнесения результата, но если его не сделать, то найдётся очень много скептиков, которые скажут, что я что-то не то протестировал и получил не те данные, какие они получили у себя. А этим объяснением я снимаю с себя всякую ответственность за несовпадения, которые обязательно будут.

    Читатель понимает, что величина N может быть произвольной, так же как произвольными могут быть числа в массиве. Поэтому совершенно ясно, что любой обзор работы алгоритмов на практике будет субъективным и неполным. Какое распределение чисел X[i] выбрать? Какими выбрать верхнюю и нижнюю границы суммируемых чисел? Сколько чисел взять для получения репрезентативного результата? Сколько случайных тестов нужно для вычисления адекватной средней погрешности?

    Видите в какой трудной ситуации находится обзорщик? Стоит что-нибудь сделать неправильно, читатель найдёт к чему придраться «со своей колокольни». Поэтому я заранее предупреждаю, что мой обзор будет обусловлен исключительно моими предпочтениями выбора тестов. Они описаны рядом с таблицами ниже. Если вам нужны другие результаты (а результаты можно получить любые, всё зависит от степени предвзятости экспериментатора), вы можете создать другие тесты :)

    Второе предупреждение: как подавать результат? В каких единицах его будет удобнее всего принять читателям, не нуждающимся в детальном научном исследовании? Ну, допустим, я дам вам абсолютную погрешность, которая вычисляется по формуле

    $|s_{got} – s_{exact}|,$


    где $s_{got}$ — сумма, полученная алгоритмами 1–3, а $s_{exact}$ — наиболее точная сумма, полученная алгоритмом 0.

    Однако эта разность точного и приближённого ответа не будет информативной, потому что сильно зависит от порядка результата (порядок этот, как вы помните, может быть от -324 до +308), вам говорит о чём-нибудь число 1.23456e+224 в качестве погрешности? Вряд ли. Гораздо более информативной будет относительная погрешность:

    $\frac{|s_{got} – s_{exact}|}{|s_{exact}|}\cdot 100\%$


    Но и она будет представлять из себя числа, порядки которых будут очень маленькими и трудными для восприятия через ощущения человека. Смотрите сами: относительная погрешность 1.23456e-10%. Удобно?

    Тогда принимают решение показывать погрешность через ulp. Напомню, что ulp(x) — это, грубо говоря, ценность последнего бита числа x. Например, для чисел x типа double на интервале [1,2) величина ulp(x)=2-52. На интервале [2,4) ulp(x)=2-51 и так далее, вырастает вдвое при каждом увеличении экспоненты числа. Такой способ описания погрешности удобен тем, что показывает относительную погрешность в единицах измерения, равных ценности одного последнего бита результата. То есть зная эту погрешность, вы можете быстро понять, условно выражаясь, сколько битов результат вы «запороли» в ходе расчётов. 1ulp — потеряли 1 бит, 2-3ulp — 2 бита, 4-7ulp — 3 бита и т. д. Здесь логарифмическая зависимость.

    Чтобы получить погрешность в ulp, нужно посчитать выражение

    $\frac{|s_{got} – s_{exact}|}{ulp(s_{exact})}$


    На моём любимом языке UB++ функция ulp(x) для double вычисляется вот так:

        // Получить значение смещённой экспоненты для 'x' типа double
        // Внимание, здесь UB! Перепишите код, если он вам не подходит
    #define GET_EXP(x) (((*(unsigned long long*)(&x))&0x7FFFFFFFFFFFFFFFULL)>>52)
    
        // Получить значение ULP для числа x.
    double get_ulp (double x) { return ldexp (1.0, GET_EXP(x)-1075); }
    


    Таблицы


    Каждый алгоритм прогонялся на 100 тестах. Один тест — это массив X[N] из чисел типа double. При этом получалось два вида погрешностей: средняя по всем 100 тестам и максимальная по ним же. Оба эти числа указаны в каждой ячейке таблицы: сверху — средняя, снизу — максимальная погрешность. Измерение погрешности ведётся в ulp относительно наиболее точной суммы, полученной абсолютно точным алгоритмом.

    Помимо массива X[N], через алгоритм также одновременно прогонялись два других массива, состоящие из этих же чисел. Один из них упорядочен по возрастанию абсолютного значения, а второй — по убыванию.

    Первый набор тестов: равномерно распределённые случайные числа на интервале [1,2) в количестве N=1000. Пояснения к символам «~» «↑» и «↓» даётся ниже.

    Порядок Naive Kahan Rump–Ogita–Oishi
    ~ 4.86
    21
    0
    0
    0
    0
    4.97
    14
    0
    0
    0
    0
    4.50
    19
    0
    0
    0
    0

    Сейчас я поясню как читать такую таблицу. Символ «~» говорит, что числа в массиве расположены хаотично. Символ «↑» — упорядочены по возрастанию. Символ «↓» — по убыванию. Верхнее число в ячейке — средняя погрешность на 100 тестах, нижнее — максимальная на них же. Как понять, например, число 19 в таблице? Оно значит, что если упорядочить 100 разных массивов по убыванию и на каждом запустить алгоритм, то максимальная погрешность на этих тестах составит 19ulp, то есть, грубо говоря, 4-5 битов «потеряли». В десятичных цифрах это будет 1-2 цифры. Если учесть, что число типа double держит почти 16 десятичных цифр, то потерю 2-х цифр в бытовых задачах можно считать несущественной. При этом в среднем на этих 100 тестах погрешность составила 4,5ulp, то есть почти одну десятичную цифру.

    Теперь случайным образом назначим знак «минус» всем нашим случайным числам из интервала [1,2), чтобы положительных чисел и отрицательных стало поровну (в вероятностном смысле). Ещё к нашим символам в первом столбце таблицы мы добавляем «±» с очевидным смыслом.

    Порядок Naive Kahan Rump–Ogita–Oishi
    ±~ 158.96
    7936
    0
    0
    0
    0
    ±↑ 86.35
    2560
    0
    0
    0
    0
    ±↓ 175.70
    11776
    0
    0
    0
    0

    Вот это да! Разрешили числам быть отрицательными и тут же у наивного алгоритма, что называется, почва ушла из под ног. Погрешность в 11776 ulp эквивалентна потере 4-5 десятичных цифр (13-14 битов из 52 битов дробной части мантиссы). Остальные алгоритмы по-прежнему выдают наиболее точную сумму.

    Теперь увеличим объём чисел, установим N=106. Мы ожидаем рост погрешности. Ожидания, в целом, оправдываются.

    Порядок Naive Kahan Rump–Ogita–Oishi
    ~ 143.00
    562
    0
    0
    0
    0
    126.60
    473
    0
    0
    0
    0
    161.91
    482
    0
    0
    0
    0
    ±~ 2015.41
    38889
    0
    0
    0
    0
    ±↑ 1520.84
    33965
    0
    0
    0
    0
    ±↓ 1672.76
    36513
    0
    0
    0
    0

    Перейдём к более показательному интервалу [1e-10, 1e+10) и возьмём сразу N=106 чисел, только генерация чисел выполняется не равномерно по интервалу, а равномерно по множеству всех возможных чисел типа double на этом интервале. На целочисленном интервале [0x3ddb7cdfd9d7bdbb, 0x4202a05f20000000) (это битовая запись чисел 1e-10 и 1e+10 в формате double) выбирается целое число $a$, которое представляет собой битовую последовательность. Эту битовую последовательность мы интерпретируем как число $x$ типа double, и именно это число $x$ считаем случайным в указанном смысле. Поясню, зачем так делается. Если брать равномерное распределение чисел на математическом интервале [1e-10, 1e+10), то в основном будут попадаться числа, близкие к 1e+10, потому что их больше всего. Их в 10 раз больше, чем чисел, близких к 1e+9, а этих, в свою очередь, в 10 раз больше чем чисел, близких к 1e+8 и так далее. Если же брать числа из множества точно-представимых чисел типа double, выборка будет более репрезентативной, потому что в ней будут самые разные числа с самыми разными экспонентами. А мне это и нужно.

    Порядок Naive Kahan Rump–Ogita–Oishi
    ~ 4277.17
    4662
    0
    0
    0
    0
    29.17
    111
    0
    0
    0
    0
    7508.68
    7915
    0
    0
    0
    0
    ±~ 475.68
    8221
    1.09
    21
    0
    0
    ±↑ 65.39
    861
    0.01
    1
    0
    0
    ±↓ 270.34
    1736
    0
    0
    0
    0

    Как видим, ситуация стала хуже, потому что числа очень разные и ясно, что их сумма будет давать гораздо большую погрешность. Здесь даже алгоритм Кэхэна начал ошибаться. А последний алгоритм всё ещё держится.

    Далее у нас экспоненциальное распределение с параметром $\lambda=0.5$ в количестве чисел N=1000. Та же картина: пока числа имеют один знак — всё в целом неплохо. Когда мы назначаем половине всех чисел знак «минус», ситуация меняется.

    Порядок Naive Kahan Rump–Ogita–Oishi
    ~ 4.84
    19
    0.01
    1
    0
    0
    2.81
    11
    0.01
    1
    0
    0
    6.53
    26
    0
    0
    0
    0
    ±~ 33.83
    1650
    1.26
    23
    0
    0
    ±↑ 12.76
    422
    0.31
    6
    0
    0
    ±↓ 18.54
    548
    0
    0
    0
    0

    Напоследок посмотрим на нормальное распределение $\mathcal N(0,1)$ для N=1000 чисел.

    Порядок Naive Kahan Rump–Ogita–Oishi
    ~ 17.77
    483
    1.46
    61
    0
    0
    10.92
    243
    0.34
    19
    0
    0
    19.71
    734
    0
    0
    0
    0

    Мы видим, что если есть отрицательные числа, то наивный алгоритм начинает складывать совсем плохо. Возникает вопрос: что если складывать отдельно отрицательные числа, отдельно положительные, а затем сложить обе суммы? Ничего не будет кроме того, что станет хуже:

    Порядок Всё подряд Отдельно + и –
    ±~ 78.65
    2336
    2121.61
    50048
    ±↑ 76.27
    2496
    2465.55
    66432
    ±↓ 52.74
    480
    2863.61
    99200

    Напоследок я обещал показать в каком случае алгоритм Кэхэна может выдать ответ хуже, чем прямолинейное сложение. Возьмём тест, в котором X[i]=(double)cos(i), i=0..106-1. Убедитесь сами (здесь лишь один тест, поэтому не указана средняя погрешность, только максимальная):

    Порядок Naive Kahan Rump–Ogita–Oishi
    ~ 210 522 0
    410 0 0
    42 0 0

    Такая вот картина, друзья. Я провёл намного больше тестов, пока работал над статьёй, и сделал для себя следующие выводы.

    • Если нужно испортить сумму, получаемую наивным алгоритмом, то можно складывать отдельно положительные и отдельно отрицательные числа. Иногда, но не всегда, сортировка чисел по убыванию модуля также «ломает» наивный алгоритм, лучше сортировать по возрастанию, но и это в ряде случаев даёт более плохой результат, чем в случае хаотичного порядка.
    • Если нужно получить ожидаемо хороший результат, следует взять алгоритм Rump–Ogita–Oishi, я не нашёл случайных тестов, при которых он выдавал бы не наиболее точную сумму. Недостаток алгоритма в том, что работать он будет медленнее наивного сложения. Ускорить такой алгоритм, вероятно, можно, пользуясь битовым представлением числа с плавающей запятой и кое-какими трюками, но у нас такой задачи на этот обзор не было.
    • Алгоритм Кэхэна значительно лучше наивного суммирования и гораздо проще алгоритма Rump–Ogita–Oishi, поэтому наиболее целесообразен на практике, где точность наивного метода вам не подходит. Просто убедитесь заранее, что ваши числа не обладают какими-то экстраординарными свойствами, которые «завалят» данный алгоритм.
    • Наивное суммирование в общем-то не такое страшное, как могло показаться. Ну, подумаешь, потеряли мы 5 десятичных цифр (если N не более миллиона). Много? Если их у вас 16, то вроде бы не страшно, а если 7 (для float), то вроде бы… хотя, друзья, неужели кто-то будет складывать миллион чисел в типе float? Вы серьёзно? Это допустимо только если вам нужна одна правильная цифра, либо если массив обладает какой-то особой спецификой, и то я бы не рискнул.

    Исходный код программы тестирования я не даю по следующим причинам.

    • Это творческий процесс, я создал свою субъективную схему тестирования и не считаю это научно-достоверным знанием, а потому не выношу на суд общественности. Коды алгоритмов я дал выше, если вам очень нужно, то не составит труда написать своё сравнение алгоритмов и получить свои данные.
    • Мой код не автоматизирован, чтобы его запускать для получения разных тестов, нужно ковыряться в программе и менять какие-то параметры, это просто неприлично — давать такой код. А написание полноценной системы тестирования в задачи публикации не входило.

    Прошу отнестись с пониманием.

    Другие приложения


    У алгоритмов точного сложения двух чисел из предыдущей статьи есть и другие приложения, которые я надеюсь рассмотреть в последующих обзорах. Это вычисление скалярного произведения и значения полинома в точке. Хотя, возможно, читатель уже сам видит, как применить туда эти алгоритмы. Проблема возникает только с корректным промежуточным умножением двух чисел. И я бы не сказал, что это прямо-таки простая тема. Моя мотивация писать продолжение зависит от вашего интереса.

    Благодарю за внимание и творческих вам успехов!

    Продолжение: Статья о скалярном произведении.

    Список источников


    [1] Раздел 5.3 книги Jean-Michel Muller, “Handbook of floating-point arithmetic”, 2018.
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 45

      +4

      Помню как переизобрёл алгоритм Кахана :)
      Тот случай когда лучше было бы это знать заранее, не не вышло.

        +7
        Что я запомнил:
        — складывать даблы можно, если точность в каком то там знаке после запятой не важна
        — если важна, есть готовые алгоритмы

        Спасибо, почитать было интересно, особенно то, что алгоритмы в целом простые и логичные, при этом дают довольно качественный результат =)
          +4
          Много воды. Не знай я до прочтения этих двух статей алгоритм Кэхэна – даже первую бы не прочитал, пожав плечами – «чего автору надо?».
          Навскидку, в двух статьях материала на половину от одной, если сильно не ужимать.

          ЗЫ: За «Rump–Ogita–Oishi» спасибо, про него не знал.

          ЗЫ2: В принципе, можно эти алгоритмы рассматривать как обобщение сложения в столбик для случая, когда каждая цифра – double (только из-за специфики floating point перенос идёт в противоположном направлении, ну и вычисляем его с помощью трюка с вычитанием). Тут перешли от «одноразрядного» калькулятора к «двухразрядному» (удвоив количество значащих бит). Я ожидал увидеть переход к произвольной «разрядности» – с гарантированной точностью результата.
          Можно было бы точный результат посчитать, домножив все числа на 2^-e_min (e_min – порядок самого маленького по модулю числа) ив длинные целые, но на FPU интереснее построить цепочку «погрешностей» для алгоритма типа Кэхэна.
            +2
            Я ожидал увидеть переход к произвольной «разрядности» – с гарантированной точностью результата

            Не было такой цели, но если бы она была, то я бы начал с модификации алгоритма Rump–Ogita–Oishi. Там где 2 double, точно также делается и 3 и так далее. Но задача была всё-таки как можно быстрее победить обычную наивную сумму.


            Много воды.

            Каждому своё. Мне не так просто отыскать правильный баланс между сухой научной работой и обучающим материалом. Беру что-то среднее, что примерно в 10 раз короче, чем тот максимальный размах, с которого я вообще начинаю планирование статьи. Я работал в вузе 11 лет и хорошо представляю себе в чём нуждается большинство из тех людей, которых я представляю себе в качестве читателей. Для остальных чтение подобных статей может быть только таким: берём список источников, и не читая мою статью сразу погружаемся в их содержание, не теряя на меня больше ни минуты. Всё, никаких проблем не вижу :), и даже не нужно тратить время на объяснение того, что субъективно не понравилось, потому что ясно же, что всем не угодишь.

            +1
            Если можете «завалить» алгоритм, прошу показать тест в комментариях

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

              Это не то. Нужен конкретный тест суммы чисел, а не более сложных операций. Именно конкретный тест, а не слова.

                +1
                Ну тут либо в тест загнать именно симплекс метод, либо составить последовательность конкретных операций над конкретными значениями.
                Но в жизни ведь использовать нужно не «синтетические» тесты, а прикладные задачи. Если алгоритм хорош в тестах, но бесполезен в прикладных задачах…
                  +3

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

                    +5
                    Ну, не скажу за симплекс-метод, но как добиться – довольно-таки понятно: нужно сделать так, чтобы точная сумма в процессе вычислений была непредставима в виде суммы двух double (минимум трёх). Необходимое (но не достаточное) условие – промежуточная сумма должна иметь более 53*2 значащих бит мантиссы. Причём чтобы выявить эту погрешность после вычисления – старшие значащие биты должны уйти.

                    Т.е. навскидку что-то вроде 1e34 + 1e17 + 1 — 1e34 — 1e17.
                    Сортировка по убыванию модуля, конечно, исправит ситуацию в данном случае, да и вообще пример искусственный.
                      +2

                      Да, верно! Благодарю :) Ошибка алгоритма 3 составила 4503599627370496ulp (52 бита). Однако если числа отсортировать, то алгоритм 3 работает правильно (как по возрастанию, так и по убыванию). Теперь нужно найти именно бытовой тест, чтобы алгоритм перестал работать. Такого я найти не смог.

                        +2
                        Ну, более сотни бит мантиссы, 33 значащих цифры – тут непросто упереться в недостаток знаков :-)

                        Навскидку в реальных задачах проблемы точности вычислений (не при суммировании, при более сложных операциях) возникают, когда работаешь с плохо обусловленными матрицами. Например, расчёт течения почти несжимаемой жидкости.
                        Да банальное решение СЛАУ или обращение матрицы методом Гаусса – уже лажает. Но есть другие методы обращения матриц.
                        Наверняка можно довести задачу до такого состояния, что не только точности double, но и удвоенной точности (как обеспечивают алгоритмы, про которые вы писали) не хватит.
                          +2
                          Да, кстати, насчёт «не хватает длины двух double для представления результата»: если заглянуть в английскую вики про алгоритм Кэхэна – можно увидеть там его обобщение на три double. А где три – там и четыре, и пять…
                          en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements
                            +1

                            Да, благодарю. Вообще ведь очевидно, что в алгоритме Rump–Ogita–Oishi эти погрешности, которые складываются отдельно, можно также складывать этим же алгоритмом, сразу получая 3 double. Вообще нетрудно собрать целый каскад. Получаем аналог длинной арифметики, в которой метод two_sum как будто играет роль переноса, как это делается при сложении длинных целых чисел. Там мы тоже ведь разбиваем сумму на две части: сумма по модулю 2**64-1 и перенос.

                              +3
                              Угу, только перенос в противоположную сторону идёт, забавно.
                              BTW, Rump-Ogita-Oishi, я так полагаю, разработан только для того, чтобы избежать ветвлений в «improved Kahan–Babuška algorithm», так-то они должны давать строго одинаковые результаты.

                              Конечно, эффективней вместо Кэхэна и т.п. было бы просто использовать длинную мантиссу, но FPU такого не умеет.
                                +1
                                Конечно, эффективней вместо Кэхэна и т.п. было бы просто использовать длинную мантиссу, но FPU такого не умеет.

                                Именно так, поэтому я и намекнул в статье, что для ускорения можно попробовать смотреть на исходные числа глубже, то есть как на битовые последовательности. Это позволит взять, например, два целых числа int64 и по полной программе использовать 127-11=116 битов под мантиссу. Это будет даже точнее чем double-double. Но говорить об этом прямо в данной статье я не решился. Потому что у меня был опыт ручной обработки double, когда я решил вместо медленных операций с плавающей запятой перейти к арифметике с фиксированной запятой на целых числах, там как раз было очень удобно, что все экспоненты в моей задаче были равны и по сути можно было избавиться от сдвигов. И вот, возомнив себя самым умным, я получил замедление программы в 8 раз. Точно также однажды решил заменить умножение и сложение на fma, и получил замедление в 4 раза. Тогда я понял, что нельзя просто так предполагать, что хорошо понимаешь как что будет работать с этой плавающей арифметикой, нужны масштабные исследования. И по этой же причине я не стал сравнивать предложенные здесь алгоритмы по скорости.

                +1

                Спасибо! Очень полезная статья. Хотел уточнить есть ли и может ли быть алгоритм, который сработал бы, не хуже описанных, но с опциями компилятора типа --ffast-math?

                  +1
                  Можно поиграться с compiler barriers, чтобы результат каждого шага вычисления всегда выгружался в память.
                    +2

                    …но нет смысла: проще --ffast-math убрать из опций.

                      0
                      Если речь идёт о реальном применении (библиотека для вычислений повышенной точности), а не просто поиграться – то глобально убирать --ffast-math нежелательно. Можно локально – через #pragma или __attribute__ (надо читать доки к своему компилятору)
                    +4
                    Спасибо, автор. Тема злободневная.

                    Недавно делал, казалось бы, простую вещь. Брал 1000 результатов измерений на микроконтроллере (float) и пытался найти их среднее и дисперсию.

                    Со средним более-менее получилось, а вот дисперсия (на основе вычисления суммы квадратов) — с треском провалилась ввиду полной потери точности. Иногда получалась отрицательной! Выяснил, что проблема в алгоритме сложения. Тоже начал искать и нашёл алгоритм Кэхена, а вот за Rump-Ogita-Oishi отдельное спасибо.

                    Микроконтроллер с FPU одинарной точности, так что там материал статьи тем более актуален.
                      0
                      Забавно, получается при переносе данных из аналогового мира в цифровой, не стоит сразу думать, что проблемы с точностью автоматически решены:)
                      +2

                      Ещё большой плюс этих алгоритмов — они подходят для потоковых данных, когда чтение числа необратимо удаляет его из буфера.


                      Для массивов есть ещё вариант "дёшево и сердито", как сделано в языке Julia для суммирования по умолчанию — рекуррентное деление пополам, если длина стала 32 или меньше — наивное суммирование. Точность при этом обычно выше, чем просто при наивном последовательном суммировании, но при этом вычисления хорошо векторизуются, и по сравнению с последовательным суммированием точность практически не страдает.

                        +2
                        Спасибо за статью!

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

                        Как обстоят дела с сортировкой, если все значения одинакового знака? Можно ли считать что в таком случае сортировка точно не ухудшает результат? Если складывать отдельно положительные и отрицательные значения, но предварительно отсортировав их, ситуация лучше?

                        Ну и последний, практический вопрос: я правильно понимаю, что все приведенные алгоритмы не лучше простого перехода к арифметике с более высокой точностью? Скажем, если сравнивать приведенные алгориты на float против использования double в качестве аккамулятора.
                          +1

                          Конечно, я не мог опубликовать все таблицы, которые получил у себя, и может быть что-то не видно, но отчасти поэтому я и предложил читателям самостоятельно провести необходимые исследования. Потому что верить мне — это оказать себе медвежью услугу. Но если вам так нужно моё субъективное мнение, то вот оно...


                          Насколько это ухудшает/улучшает результат?

                          Не знаю, всегда по-разному. Может в два раза, может в 3, а может и не получиться разницы. Проверьте сами и убедитесь, пожалуйста. На экспоненциальных тестах сортировка по возрастанию обычно именно спасает (в два раза), а по убыванию — портит (тоже в два раза). На равномерных — зависит от экспоненты. Поэтому я лично не могу подобрать однозначной рекомендации. Интуиция тут не работает. Например, когда речь идёт о равномерном распределении на интервале [1,2), то вы при использовании интуиции вряд ли учли тот факт, что уже после первого действия сумма S становится наверняка больше любого из слагаемых. После второго — совсем больше и дальше остальные тысячи элементов уже складываются с большой погрешностью. При этом в случае хаотичного порядка получается то же самое! Иными словами, при одинаковом знаке сортировка по возрастанию имеет смысл только когда числа очень разные, с разными экспонентами. Когда экспоненты одинаковые — сортировать может быть вредно. Это мой субъективный вывод.


                          Как обстоят дела с сортировкой, если все значения одинакового знака? Можно ли считать что в таком случае сортировка точно не ухудшает результат?

                          Скорее всего, да, но только если экспоненты очень разные.


                          Если складывать отдельно положительные и отрицательные значения, но предварительно отсортировав их, ситуация лучше?

                          Ответ на этот вопрос есть в одной из таблиц.


                          Ну и последний, практический вопрос: я правильно понимаю, что все приведенные алгоритмы не лучше простого перехода к арифметике с более высокой точностью?

                          С точки зрения точности да, но только при условии, что вы затем не возвращаетесь обратно к арифметике с более низкой точностью, так как в этом случае вы рискуете получить ошибку двойного округления. То есть результат применения алгоритма 3 и прямолинейного применения float128 может дать разбежку в 1 бит из-за того, что вы должны будете выполнить дополнительное действие float128 → double. Но это вряд ли будет происходить на бытовых задачах. Особенно если учесть, что float128 имеет 112 битов против 104 у double-double.


                          Теперь что касается ваших слов не лучше. Что значит не лучше? У вас есть аппаратный float128? Если нет, то может быть есть хороший программный, который будет быстрее double-double? Я сомневаюсь, хотя экспериментов не проводил. Поэтому не знаю, что действительно будет лучше. А если модифицировать алгоритм 3 немного и складывать хвостики этим же алгоритмом, мы получим double-double-double. И вообще можем создать любой каскад по принципу того, как создаётся длинная арифметика с целыми числами. Поэтому опять я не понимаю, что значит не лучше. Я давно убедился, что в профессиональной сфере нельзя сравнивать различные технологии по правилу лучше / не лучше. Каждая технология предназначена для сугубо своей задачи. Для конкретной задачи ДА, она может быть лучше или хуже. Сама по себе она лучше или хуже быть не может в принципе. Такой вот у меня ответ получился :)


                          Но повторюсь, лучше вы эти (или другие) ответы получите сами, проверив всё на своей практике. Это будет куда продуктивнее.

                            +1
                            Когда экспоненты одинаковые — сортировать может быть вредно. Это мой субъективный вывод.

                            Вряд ли именно вредно, скорее просто бесполезно.


                            С точки зрения точности да, но только при условии, что вы затем не возвращаетесь обратно к арифметике с более низкой точностью, так как в этом случае вы рискуете получить ошибку двойного округления.

                            Вы так пишете как будто алгоритм Кэхена два раза не округляет!

                              0

                              Да, можно и так сказать, но я подразумевал ещё потерю времени на сортировку или на искусственное создание заведомо отсортированного порядка на входе в алгоритм.

                                0
                                Вы так пишете как будто алгоритм Кэхена два раза не округляет!

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


                                все приведенные алгоритмы не лучше простого перехода к арифметике с более высокой точностью

                                Я отвечаю, что мы не знаем, что будет лучше, но, вероятно, ответы будут разными. Какой из них точнее — неизвестно.

                                  0

                                  Ну так точного совпадения результатов никто и не ждёт же (хотя у меня есть ощущение, что аккуратно подобрав разрядность и режимы округления в теории можно и его добиться), вопрос лишь в оценке полученной погрешности.


                                  Какой из них точнее — неизвестно.

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

                                    0
                                    А ваш ответ выглядит так, как будто float128 хуже двух double, хотя вы сами же и пишете что у float128 разрядов мантиссы больше.

                                    Может быть и хуже, потому если float128 не аппаратный, то будет работать медленно. Ведь важна не только точность, но приличная скорость. И второе, чем double лучше, так это тем, что из него можно сделать 3-double, а из float128 нельзя. То есть в зависимости от решаемой задачи мы можем сказать, что лучше. А автор вопроса, на который я отвечал, просит дать ответ безотносительно решаемой задачи. А это всегда будет приводить к неопределённым ответам, и, следовательно, к неопределённым замечаниям к этим неопределённым ответам. И так далее. В итоге получится, что правы все, но каждый — по-своему.

                                      0

                                      Извините, но вот этот фрагмент не выглядит как сравнение по скорости:


                                      С точки зрения точности да, но только при условии, что вы затем не возвращаетесь обратно к арифметике с более низкой точностью, так как в этом случае вы рискуете получить ошибку двойного округления. То есть результат применения алгоритма 3 и прямолинейного применения float128 может дать разбежку в 1 бит из-за того, что вы должны будете выполнить дополнительное действие float128 → double.

                                      Он выглядит так, как будто float128 менее точен, чем два double, хотя следующим же предложением вы пишете что это не так.


                                      Кстати, исходный вопрос был не о сравнении float128 против двух double, а о сравнении double против двух float.

                                        0

                                        Здесь всё зависит от субъективного восприятия информации. Если я говорю о возможной погрешности в одном случае и умалчиваю о ней же в другом, то это не значит, что её там нет. Это значит что я о ней не сказал. Разумный читатель сможет подумать и сделать нужные ему выводы, в частности, вот вы же сделали — и хорошо.


                                        Кстати, исходный вопрос был не о сравнении float128 против двух double, а о сравнении double против двух float.

                                        Не совсем. Автор задаёт вопрос вот так:


                                        я правильно понимаю, что все приведенные алгоритмы не лучше простого перехода к арифметике с более высокой точностью?

                                        А уже потом приводит пример (орфография сохранена):


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

                                        На вопрос я ответил так, как посчитал нужным, предупредив об этом заранее. На пример я не посчитал нужным отвечать, автор вопроса в состоянии сделать это самостоятельно. Вся информация для этого у него имеется, а моё время в этом случае дороже.

                            0
                            Этот алгоритм проще было бы показать на числах float, а точный результат посчитать с помощью double. Еще суммирование таких чисел зависит от порядка их следования. Так что если не отсортировать, то точность не будет идеальной. После первой статьи я написал пару тестов на своем процессоре и случайно обнаружил, что «старое» суммирование, но с SIMD подходом, дало лучший результат чем алгоритм Кахана. Просто порядок следования чисел был для SIMD удачнее :)
                              0

                              Ни в коем случае точный результат с помощью double получить нельзя. Потому что разбежка у float по экспоненте может составлять более 250 битов, а double удержит только 52 из них. По поводу порядка я уже написал комментарий другому читателю.

                                0
                                Вы говорите про общий случай. Но для своего теста мне был удобнее именно такой подход. Спасибо за ссылку о сортировке, а то я еще не дочитал, но уже пишу :)
                              0
                              Планирует ли автор писать продолжение про умножение, деление, извлечение квадратного корня и прочих функций повышенной точности?
                                +1

                                К сожалению, моя жизненная позиция запрещает обсуждать планы на будущее и отвечать на подобные вопросы. В конце статьи сказано, что от интереса читателей зависит, напишу ли я здесь про скалярное произведение и вычисление полиномов. Про остальное — не могу сказать.

                                  +1
                                  В таком случае всем заинтересованным рекомендую статью Library for Double-Double and Quad-Double Arithmetic, которая намного более полно и в то же время кратко освещает данный вопрос. Реализацию этих алгоритмов можно найти здесь (qd-2.3.22.tar.gz) — не только сложения, но и умножения/деления с элементарными функциями. Для C# реализации этих алгоритмов есть в библиотеке metanumerics — правда, функции там пока ещё не все реализованы.
                                    +1

                                    Благодарю, вот это действительно конструктивный подход, вместо ожиданий какого-то чуда от заурядного автора на Хабре :)

                                0

                                Хотелось бы продолжениие статьи, особенно про скаляроное произведение и умножение. Это ж вся цифровая фильтрация.

                                  0

                                  Вы возлагаете на меня слишком большие надежды. Я пишу научно-популярный текст, цель которого — положить начало знакомства читателей с обсуждаемой темой. Никоим образом я не хочу чтобы эти материалы именно в таком поверхностном виде ложились бы в основу чьей-то серьёзной работы. Я делаю здесь именно школу, а не науку, показываю верхушку айсберга того, что уже описано в научной литературе, но по какой-то причине плохо описано в литературе популярной. Поэтому если вас интересует судьба всей цифровой фильтрации, то лучше сразу покупать книгу [1], затем те статьи, на которые ссылаются её авторы при обсуждении интересующих вас алгоритмов. Другого пути, увы, нет. Напрасно думать, что моя статья на Хабре может дать хотя бы тысячную долю того знания, которое реально требуется для ответственных проектов. А вся цифровая фильтрация, согласитесь, — это сложная и ответственная область.

                                  0
                                  Ещё такая мысль возникла: а если всегда складывать два самых маленьких числа из имеющихся и заменять их на сумму (кажется, квадратичного времени работы можно избежать с помощью алгоритма heapsort), как это повлияет на точность?
                                    0

                                    Этот алгоритм описан в книге [1], можете сами прочитать. Погрешность в этом случае будет не больше чем значение ulp умножить на сумму всех тех частичных сумм, которые вы будете получать по ходу работы алгоритма. Это довольно плохая теоретическая оценка. А что будет на практике — попробуйте сами.

                                      +1

                                      heapsort — довольно дорогой вариант, проще воспользоваться алгоритмом "двух очередей": первая очередь — это суммируемый массив, сортированный по возрастанию. Во второй очереди будут частичные суммы, тоже по возрастанию, но их сортировать не надо — они автоматически упорядоченными будут. Каждый раз как требуется очередное слагаемое — забираем его из той очереди где "голова" меньше.

                                      0
                                      В студенческую бытность на курсе высшей математики познакомился с рядами, узнал как представить Пи через арктангенс, а последний разложить в виде бесконечной суммы.
                                      Программировал в BorlandC 3.1, использовал float. Каково же было моё удивление, что Пи у меня получается 3.12, причём как ни увеличивай количество членов в ряду, ничего не меняется.
                                      Инета тогда ещё не было.
                                      Откудаж мне было знать, что даже суммировал я неверно. Вот теперь знаю.

                                      Only users with full accounts can post comments. Log in, please.