C++. Сборник рецептов - Д. Стефенс

template<class Iter_T>

void outputRowOrColumn(Iter_T iter, int n) {

 for (int i=0; i < n; ++i) {

  cout << iter[i] << " ";

 }

 cout << endl;

}

template<class Matrix_T>

void initializeMatrix(Matrix_T& m) {

 int k = 0;

 for (int i=0; i < m.rows(); ++i) {

  for (int j=0; j < m.cols(); ++j) {

   m[i][j] = k++;

  }

 }

}

template<class Matrix_T>

void outputMatrix(Matrix_T& m) {

 for (int i=0; i < m.rows(); ++i) {

  cout << "Row " << i << " = ";

   outputRowOrColumn(m.row(i), m.cols());

 }

 for (int i=0; i < m.cols(); ++i) {

  cout << "Column " << i << " = ";

   outputRowOrColumn(m.col(i), m.rows());

 }

}

int main() {

 kmatrix<int, 2, 4> m;

 initializeMatrix(m); m *= 2;

 outputMatrix(m);

}

Программа примера 11.31 выдает следующий результат.

Row 0 = 0 2 4 6

Row 1 = 8 10 12 14

Column 0 = 0 8

Column 1 = 2 10

Column 2 = 4 12

Column 3 = 6 14

Представленные в примерах 11.30 и 11.31 определение шаблона класса kmatrix и пример его использования очень напоминают шаблон класса matrix из рецепта 11.14. Единственным существенным отличием является то, что при объявлении экземпляра kmatrix приходится передавать размерности матрицы через параметры шаблона, например;

kmatrix<int 5, 6> m; // объявляет матрицу с пятью строками и шестью

                     // столбцами

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

Как и рассмотренный ранее шаблон статического вектора (kvector), шаблон kmatrix особенно эффективен при небольших размерах матрицы.

Рецепты 11.14 и 11.16.

11.16. Умножение матриц

Требуется эффективно выполнить умножение двух матриц.

Пример 11.32 показывает, как можно выполнить умножение матриц, причем эта реализация подходит как для динамических, так и для статических матриц. Формально этот алгоритм реализует соотношение A=A+B*C, которое (возможно, неожиданно) вычисляется более эффективно, чем A=B*C.

Пример 11.32. Умножение матриц

#include "matrix.hpp" // рецепт 11.13

#include "kmatrix.hpp" // рецепт 11.14

#include <iostream>

#include <cassert>

using namespace std;

template<class M1, class M2, class M3>

void matrixMultiply(const M1& m1, const M2& m2, M3& m3) {

 assert(m1.cols() == m2.rows());

 assert(m1.rows() == m3.rows());

 assert(m2.cols() == m3.cols());

 for (int i=m1.rows()-1; i >= 0; --i) {

  for (int j=m2.cols()-1; j >= 0; --j) {

   for (int k = m1.cols()-1; k >= 0; --k) {

    m3[i][j] += m1[i][k] * m2[k][j];

   }

  }

 }

}

int main() {

 matrix<int> m1(2, 1);

 matrix<int> m2(1, 2);

 kmatrix<int, 2, 2> m3;

 m3 = 0;

 m1[0][0] = 1;

 m1[1][0] = 2;

 m2[0][0] = 3;

 m2[0][1] = 4;

 matrixMultlply(m1, m2, m3);

 cout << "(" << m3[0][0] << ", " << m3[0][1] << ")" << endl;

 cout << "(" << m3[1][0] << ", " << m3[1][1 ] << ")" << endl;

}

Программа примера 11.32 выдает следующий результат.

(3, 4)

(6, 8)

При умножении двух матриц число столбцов первой матрицы должно равняться числу строк второй матрицы. Число строк полученной матрицы равно числу строк первой матрицы, а число столбцов равно числу столбцов второй матрицы. Я обеспечиваю эти условия в отладочной версии с помощью макроса assert, определенного в заголовочном файле <cassert>.

Решающее значение для эффективной реализации умножения имеет отсутствие избыточных операций по созданию и копированию временных объектов. Так, представленная в примере 11.32 функция умножения матриц передает результат по ссылке. Если бы алгоритм умножения я реализовал впрямую путем перегрузки оператора operator*, это привело бы к лишним операциям распределения, копирования и освобождения памяти, занимаемой временной матрицей. Потенциально такой подход может оказаться очень затратным при работе с большими матрицами.

В примере 11.32 реализуется равенство A=A+B*C, а не A=B*C, для того чтобы избежать лишней инициализации значений матрицы A.

Рецепт 11.17.

11.17. Вычисление быстрого преобразования Фурье

Требуется выполнить эффективный расчет дискретного преобразования Фурье (ДПФ), используя алгоритм быстрого преобразования Фурье (БПФ).

Программный код примера 11.33 обеспечивает базовую реализацию БПФ.

Пример 11.33. Реализация БПФ

#include <iostream>

#include <complex>

#include <cmath>

#include <iterator>

using namespace std;

unsigned int bitReverse(unsigned int x, int log2n) {

 int n = 0;

 int mask = 0x1;

 for (int i=0; i < log2n; i++) {

  n <<= 1;

  n |= (x & 1);

  x >>= 1;

 }

 return n;

}

const double PI = 3.1415926536;

template<class Iter_T>

void fft(Iter_r a, Iter_r b, int log2n) {

 typedef typename iterator_traits<Iter T>::value_type complex;

 const complex J(0, 1);

 int n = 1 << log2n;

 for (unsigned int i=0; i < n; ++i) {

  b[bitReverse(i, log2n)] = a[i];

 }

 for (int s = 1; s <= log2n; ++s) {

  int m = 1 << s;

  int m2 = m >> 1;

  complex w(1, 0);

  complex wm = exp(-J * (PI / m2));

  for (int j=0; j < m2; ++j) {

   for (int k=j; k < n; k += m) {

    complex t = w * b[k + m2];

    complex u = b[k];

    b[k] = u + t;

    b[k + m2] = u - t;

   }

   w *= wm;

  }

 }

}

int main() {

 typedef complex<double> cx;

 cx a[] = { cx(0, 0), cx(1, 1), cx(3, 3), cx(4, 4),

  cx(4, 4), cx(3, 3), cx(1, 1), cx(0, 0) };

 cx b[8];

 fft(a, b, 3);

 for (int i=0; i<8; ++i) cout << b[i] << "n";

}

Программа примера 11.33 выдает следующий результат.

(16,16)

(-4.82843,-11.6569)

(0,0)

(-0.343146,0.828427)

(0.0)

(0.828427,-0.343146)

(0,0)

(-11.6569,-4.82843)

Преобразование Фурье играет важную роль в спектральном анализе и часто используется в технических и научных приложениях. БПФ — это алгоритм вычисления ДПФ, который имеет сложность порядка N log2(N) в отличие от ожидаемой сложности N² для простой реализации ДПФ. Такое впечатляющее ускорение достигается в БПФ благодаря устранению избыточных вычислений.

Очень не просто найти хорошую реализацию БПФ, написанную на «правильном» C++ (т. е. когда программа на C++ не является механическим переложением алгоритмов, написанных на Фортране или С) и которая не была бы защищена сильно ограничивающей лицензией. Представленный в примере 11.33 программный код основан на открытом коде, который можно найти в сетевой конференции Usenet, посвященной цифровой обработке сигналов (comp.dsp). Большим преимуществом реализации БПФ на правильном C++ по сравнению с более распространенным решением в стиле С является то, что стандартная библиотека содержит шаблон complex, который позволяет существенно снизить объем необходимого программного кода. В представленной в примере 11.33 функции fft() основное внимание уделялось простоте, а не эффективности.

11.18. Работа с полярными координатами

Требуется обеспечить представление полярных координат и манипулирование ими.

Шаблон complex из заголовочного файла <complex> содержит функции преобразования в полярные координаты и обратно. Пример 11.34 показывает, как можно использовать класс шаблона complex для представления и манипулирования полярными координатами.

Пример 11.34. Применение шаблонного класса complex для представления полярных координат

#include <complex>

#include <iostream>

using namespace std;

int main() {

 double rho = 3.0; // длина

 double theta = 3.141592 / 2; // угол