Безопасный и быстрый FFT

Вдохновленный Herb Sutter убедительной лекцией Не ваш отец С++, я решил еще раз взглянуть на последнюю версию С++ с помощью Microsoft Visual Studio 2010. Мне было особенно интересно по утверждению Herb, что С++ "безопасен и быстр", потому что я пишу много критически важных для работы кода.

В качестве эталона я решил попытаться написать один и тот же простой алгоритм БПФ на разных языках.

Я придумал следующий код С++ 11, который использует встроенный тип complex и vector:

#include <complex>
#include <vector>

using namespace std;

// Must provide type or MSVC++ barfs with "ambiguous call to overloaded function"
double pi = 4 * atan(1.0);

void fft(int sign, vector<complex<double>> &zs) {
    unsigned int j=0;
    // Warning about signed vs unsigned comparison
    for(unsigned int i=0; i<zs.size()-1; ++i) {
        if (i < j) {
            auto t = zs.at(i);
            zs.at(i) = zs.at(j);
            zs.at(j) = t;
        }
        int m=zs.size()/2;
        j^=m;
        while ((j & m) == 0) { m/=2; j^=m; }
    }
    for(unsigned int j=1; j<zs.size(); j*=2)
        for(unsigned int m=0; m<j; ++m) {
            auto t = pi * sign * m / j;
            auto w = complex<double>(cos(t), sin(t));
            for(unsigned int i = m; i<zs.size(); i+=2*j) {
                complex<double> zi = zs.at(i), t = w * zs.at(i + j);
                zs.at(i) = zi + t;
                zs.at(i + j) = zi - t;
            }
        }
}

Обратите внимание, что эта функция работает только для векторов n -элементов, где n является интегральной степенью двух. Любой, кто ищет быстрый код FFT, который работает для любого n, должен смотреть на FFTW.

Как я понимаю, традиционный синтаксис xs[i] из C для индексации a vector не выполняет проверку границ и, следовательно, не является безопасным для памяти и может быть источником ошибок памяти, таких как детерминированное повреждение и нарушения доступа к памяти. Поэтому я использовал xs.at(i).

Теперь я хочу, чтобы этот код был "безопасным и быстрым", но я не эксперт на С++ 11, поэтому я хотел бы попросить улучшения этого кода, которые сделают его более идиоматичным или эффективным?

Ответы

Ответ 1

Я думаю, что вы чрезмерно "безопасны" в своем использовании в(). В большинстве ваших случаев используемый индекс тривиально проверяется как ограниченный размером контейнера в цикле for.

например.

  for(unsigned int i=0; i<zs.size()-1; ++i) { 
    ...
    auto t = zs.at(i);

Единственными, которые я оставил бы как(), являются (i + j) s. Не сразу видно, будут ли они всегда ограничены (хотя, если бы я был действительно неуверен, я бы, вероятно, проверил вручную, но я не знаком с БПФ, чтобы иметь мнение в этом случае).

Есть также некоторые фиксированные вычисления, повторяющиеся для каждой итерации цикла:

int m=zs.size()/2;
pi * sign
2*j

И zs.at(i + j) вычисляется дважды.

Возможно, что оптимизатор может их поймать, но если вы рассматриваете это как критическую для производительности, и у вас есть ваши таймеры, проверяющие его, я бы вытащил их из циклов (или, в случае zs.at( я + j), просто возьмите ссылку на первое использование) и посмотрите, влияет ли это на таймер.

Говоря о втором угадывании оптимизатора: я уверен, что вызовы .size() будут привязаны, по крайней мере, к прямому вызову внутренней переменной-члену, но учитывая, сколько раз вы это называете, d также экспериментировать с введением локальных переменных для zs.size() и zs.size() - 1 upfront. Они, скорее всего, тоже будут помещены в регистры.

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

Что касается идиоматического, то единственным моим комментарием является, что size() возвращает std:: size_t (обычно это typedef для unsigned int, но более идиоматично использовать этот тип). Если вы хотите использовать auto, но избегаете предупреждения, вы можете попробовать добавить суффикс ul к 0 - не уверен, что я сказал бы, что это идиоматично. Я полагаю, вы уже идиотичны, не используя итераторы здесь, но я понимаю, почему вы не можете этого сделать (легко).

Обновление

Я дал все свои предложения попробовать, и все они имели измеримое улучшение производительности - за исключением я + j и 2 * j precalcs - они фактически вызвали небольшое замедление! Я предполагаю, что они либо предотвратили оптимизацию компилятора, либо помешали ему использовать регистры для некоторых вещей.

В целом я получил > 10% perf. улучшение с этими предложениями. Я подозреваю, что больше можно было бы, если бы второй блок циклов был реорганизован немного, чтобы избежать прыжков - и, сделав это, включение набора инструкций SSE2 может дать значительный импульс (я попробовал это как есть и увидел небольшое замедление).

Я думаю, что рефакторинг, наряду с использованием чего-то вроде MKL для вызовов cos и sin, должен давать большие и менее хрупкие улучшения. И ни одна из этих вещей не зависит от языка (я знаю, что это первоначально сравнивалось с реализацией F #).

Обновление 2

Я забыл упомянуть, что предварительное вычисление zs.size() действительно имело значение.

Обновление 3

Также забыл сказать (до тех пор, пока не будет напомнено @xeo в комментарии к OP), что блок, следующий за я < j check можно свернуть до std:: swap. Это более идиоматично и, по крайней мере, как исполнитель - в худшем случае должен встроить тот же код, что и написанный. Действительно, когда я это делал, я не видел изменений в производительности. В других случаях это может привести к усилению производительности, если доступны конструкторы перемещения.