Глава 11 Наука и математика

11.0. Введение

Язык программирования C++ хорошо подходит для решения научных и математических задач из-за своей гибкости, выразительности и эффективности. Одно из самых больших преимуществ применения C++ для выполнения численных расчетов связано с тем, что он помогает избегать избыточности.

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

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

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

Многие программисты, использующие С++, все же с недоверием относятся к шаблонам и обобщенному программированию из-за очевидной их сложности. Когда шаблоны впервые были введены в язык, они не были хорошо реализованы, а программисты и разработчики компиляторов не очень хорошо их понимали. В результате многие программисты, включая автора, избегали пользоваться обобщенным программированием на C++ в течение нескольких лет, пока эта технология не достигла зрелости.

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

11.1. Подсчет количества элементов в контейнере

Проблема

Требуется найти количество элементов в контейнере.

Решение

Подсчитать количество элементов в контейнере можно при помощи функции-члена 

size
или функции 
distance
, определенной в заголовочном файле
, как это делается в примере 11.1.

Пример 11.1. Подсчет количества элементов в контейнере

#include 

#include 

#include 


using namespace std;


int main() {

 vector v;

 v.push_back(0);

 v.push_back(1);

 v.push_back(2);

 cout << v.size() << endl;

 cout << distance(v.begin(), v.end()) << endl;

}

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

3

3

Обсуждение

Функция-член 

size
, которая возвращает количество элементов стандартного контейнера, является наилучшим решением в тех случаях, когда доступен объект контейнера. В примере 11.1 я также продемонстрировал применение функции 
distance
, потому что при написании обобщенного программного кода обычно имеешь дело только с парой итераторов. Работая с итераторами, вы часто не знаете тип контейнера и не имеете доступа к его функциям-членам.

Функция 

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

cout << distance::iterator>(v.begin(), v.end()) << endl;

Производительность функции 

distance
зависит от типа используемого итератора. Время ее выполнения будет постоянным, если итератор ввода является итератором с произвольным доступом; в противном случае время ее работы будет линейным. (Концепция итератора рассматривается в рецепте 7.1.)

Смотри также

Рецепт 15.1.

11.2. Поиск наибольшего или наименьшего значения в контейнере

Проблема

Требуется найти максимальное или минимальное значение в контейнере.

Решение

Пример 11.2 показывает, как можно находить максимальные и минимальные элементы контейнера с помощью функций 

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

Пример 11.2. Поиск минимального или максимального элемента контейнера

#include 

#include 

#include 


using namespace std;


int getMaxInt(vector& v) {

 return *max_element(v.begin(), v.end());

}


int getMinInt(vector& v) {

 return *min_element(v.begin(), v.end());

}


int main() {

 vector v;

 for (int i=10; i < 20; ++i) v.push_back(i);

 cout << "min integer = " << getMinInt(v) << endl;

 cout << "max integer = " << getMaxInt(v) << endl;

}

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

min integer = 10

max integer =19

Обсуждение

Вероятно, вы заметили, что выполняется разыменование значения, возвращаемого функциями 

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

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

Пример 11.3. Обобщенная функция, возвращающая минимальное и максимальное значения

#include 

#include 

#include 


using namespace std;


template

void computeMinAndMax(Iter_T first, Iter_T last, Value_T& min, Value_T& max) {

 min = *min_element(first, last);

 max = *max_element(first, last);

}


int main() {

 vector v;

 for (int i=10; i < 20; ++i) v.push_back(i);

 int min = -1;

 int max = -1;

 computeMinAndMax(v.begin(), v.end(), min, max);

 cout << "min integer = " << min << endl;

 cout << "max integer = " << max << endl;

}

В примере 11.3 я написал шаблон функции 

computeMinAndMax
, которая принимает два параметра шаблона: один — это тип итератора, другой — тип минимальных и максимальных значений. Поскольку оба параметра шаблона являются также параметрами функции, компилятор C++ может догадаться, какие два отдельных типа (
Iter_T
и 
Value_T
) используются, как это я продемонстрировал в рецепте 11.1. Это позволяет мне не указывать явно тип параметров шаблона, как это сделано ниже.

compute_min_max::iterator, int>(...)

При выполнении функций 

min_element
и 
max_element
используется оператор 
operator<
для сравнения значений, на которые ссылаются итераторы. Это значит, что, если итератор ссылается на тип, который не поддерживает этот тип сравнения, компилятор выдаст сообщение об ошибке. Однако функции 
min_element
и 
max_element
можно также использовать с функтором сравнения, определенным пользователем, т.е. с указателем на функцию или с объектом-функцией.

Для функций 

min_element
и 
max_element
необходим специальный функтор, принимающий два значения (они имеют тип объектов, на которые ссылается итератор) и возвращающий значение типа 
Boolean
, показывающее, является ли первое значение меньше, чем второе. Функтор, который возвращает значение типа 
Boolean
, называется предикатом. Рассмотрим, например, поиск самого большого элемента в наборе пользовательских типов (пример 11.4).

Пример 11.4. Поиск максимального элемента для пользовательских типов

#include 

#include 

#include 


using namespace std;


struct Chessplayer {

 ChessPlayer(const char* name, int rating)

 : name_(name), rating_(rating) { }

 const char* name_;

 int rating_;

};


struct IsWeakerPlayer {

 bool operator()(const ChessPlayer& x, const ChessPlayer& y) {

 return x.rating_ < y.rating_;

};


int main() {

 ChessPlayer kasparov("Garry Kasparov", 2805);

 ChessPlayer anand("Viswanathan Anand", 2788);

 ChessPlayer topalov("Veselin Topalov", 2788);

 vector v;

 v.push_back(kasparov);

 v.push_back(anand);

 v.push_hack(topalov);

 cout << "the best player is ";

 cout << max_element(v.begin(), v.end(), IsWeakerPlayer())->name_;

 cout << endl;

}

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

the best player is Garry Kasparov (лучший игрок - Гарри Каспаров)

Функторы

Многие STL-алгоритмы в качестве параметров используют определенные пользователем объекты-функции и указатели на функции. И те и другие называются функторами (functors). Иногда в литературе термин «объект-функция» используется как синоним термина «функтор», однако я использую термин «объект-функция» для обозначения только экземпляров класса или структур, которые перегружают 

operator()
. Какой из двух типов функторов лучше использовать? В большинстве случаев объект-функция более эффективен, потому что большинство компиляторов могут легко его реализовать в виде встроенной функции.

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

Наконец, объекты-функции могут определяться внутри другой функции или класса. Указатели на функции приходится объявлять в области видимости пространства имен.

В примере 11.4 я показал, как в функции 

max_element
можно использовать пользовательский предикат. Этот предикат является объектом-функцией 
IsWeakerPlayer
.

Альтернативой пользовательскому предикату, показанному в примере 11.4, является перегрузка оператора 

operator<
для структуры 
ChessPlayer
. Это хорошо работает в определенных случаях, но предполагает, что самой важной является сортировка игроков по рейтингу. Может оказаться, что более распространенной является сортировка по именам. Поскольку в данном случае выбор метода сортировки может быть произвольным, я предпочитаю не определять оператор 
operator<
.

11.3. Вычисление суммы и среднего значения элементов контейнера

Проблема

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

Решение

Для расчета суммы можно использовать функцию 

accumulate
из заголовочного файла
и затем разделить ее на количество элементов, получая среднее значение. Пример 11.5 демонстрирует, как это можно сделать, используя вектор.

Пример 11.5. Вычисление суммы и среднего значения элементов контейнера

#include 

#include 

#include 


using namespace std;


int main() {

 vector v;

 v.push_back(1);

 v.push_back(2);

 v.push_back(3);

 v.push_back(4);

 int sum = accumulate(v.begin(), v.end(), 0);

 double mean = double(sum) / v.size();

 cout << "sum = " << sum << endl;

 cout << "count = " << v.size() << endl;

 cout << "mean = " << mean << endl;

}

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

sum = 10

count = 4

mean = 2.5

Обсуждение

Как правило, функция 

accumulate
обеспечивает самый эффективный и самый простой способ вычисления суммы всех элементов, содержащихся в контейнере.

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

Пример 11.6. Обобщенная функция по расчету среднего значения

template

double computeMean(Iter_T first, Iter_T last) {

 return static_cast(accumulate(first, last, 0.0))

  / distance(first, last);

}

Функция 

computeMean
из примера 11.6 подойдет в большинстве случаев, но она имеет одно ограничение: не работает она с такими итераторами ввода, как 
istream_iterator
.

Итераторы istream_iterator и ostream_iterator

Шаблоны классов 

istream_iterator
и 
ostream_iterator
представляют собой специализированные итераторы, определенные в заголовочном файле
которые позволяют рассматривать потоки как однопроходные контейнеры.

istream_iterator
является итератором ввода, который выступает в роли оболочки такого потока ввода, как 
cin
или 
ifstream
, позволяя использовать его в качестве параметра во многих обобщенных функциях. 
ostream_iterator
является итератором вывода, который позволяет использовать потоки вывода, как будто они являются контейнерами. Использование итераторов 
istream_iterator
и 
ostream_iterator
является хорошей привычкой, так как с их помощью легче создавать повторно используемый программный код

Итератор 

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

Пример 11.7. Более обобщенная функция по расчету среднего значения

#include 

#include 

#include 


using namespace std;


template

Value_T computeMean(Iter_T first, Iter_T last) {

 if (first == last) throw domain_error("mean is undefined");

 Value_T sum;

 int cnt = 0;

 while (first != last) {

  sum += *first++;

  ++cnt;

 }

 return sum / cnt;

)


int main() {

 cout << "please type in several integers separated by newlines" << endl;

 cout << "and terminated by an EOF character (i.e , Ctrl-Z)" << endl;

 double mean = computeMean(

  istream_iterator(cin), istream_iterator());

 cout << "the mean is " << mean << endl;

}

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

istream_iterator
. Кроме того, алгоритмы с единственным проходом часто более эффективны.

Возможно, вас удивляет то, что я решил тип, возвращаемый функцией 

computeMean
из примера 11.7, передать в качестве параметра шаблона, а не выводить его из типа итератора. Это сделано по той причине, что обычно статистические расчеты выполняются с более высокой точностью, чем точность значений, содержащихся в контейнере. Так, в программном коде примера 11.7 возвращаемое среднее значение набора чисел целого типа имеет тип 
double
.

11.4. Фильтрация значений, выпадающих из заданного диапазона

Проблема

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

Решение

Используйте функцию 

remove_copy_if
, определенную в
, как показано в примере 11.8.

Пример 11.8 Удаление из последовательности элементов, значения которых меньше заданного

#include 

#include 

#include 

#include 


using namespace std;


struct OutOfRange {

 OutOfRange(int min, int max) :

 min_(min), max_(max) {}


 bool operator()(int x) {

 return (x < min_) || (x > max_);

 }

 int min_;

 int max_;

};


int main() {

 vector v;

 v.push_back(6);

 v.push_back(12);

 v.push_back(10);

 v.push_back(24);

 v.push_back(30);

 remove_copy_if(v.begin(), v.end(),

 ostream_iterator(cout, "\n"), OutOfRange(10, 25));

}

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

12

18

24

Обсуждение

Функция 

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

Для функции 

remove_copy_if
требуется унарный предикат (функтор, который принимает один аргумент и возвращает значение типа 
boolean
), который возвращает значение «истина», когда элемент не должен копироваться. В примере 11.8 предикатом является объект-функция 
OutOfRange
. Конструктор 
OutOfRange
принимает нижнюю и верхнюю границу и перегружает оператор 
operator()
. Функция 
operator()
принимает параметр целого типа и возвращает значение «истина», если переданный аргумент меньше, чем нижняя граница, или больше, чем верхняя граница.

11.5. Вычисление дисперсии, стандартного отклонения и других статистических функций

Проблема

Требуется рассчитать значение одной или нескольких обычных статистических функций, например дисперсии (variance), стандартного отклонения (standard deviation), коэффициента асимметрии (skew) и эксцесса (kurtosis) для последовательности чисел.

Решение

Функцию 

accumulate
из заголовочного файла
можно использовать для расчета многих статистических параметров, а не только для суммирования пользовательских объектов-функций. Пример 11.9 показывает, как можно вычислить значения некоторых важных статистические функций при помощи 
accumulate
.

Пример 11.9. Статистические функции

#include 

#include 

#include 

#include 

#include 

#include 


using namespace std;


template

T nthPnwer(T x) {

 T ret = x;

 for (int i=1; i < N; ++i) {

  ret *= x;

 }

 return ret;

}


template

struct SumDiffNthPower {

 SumDiffNthPower(T x) : mean_(x) {};

 T operator()(T sum, T current) {

  return sum + nthPower(current - mean_);

 }

 T mean_;

};


template

T nthMoment(Iter_T first, Iter_T last, T mean) {

 size_t cnt = distance(first, last);

 return accumulate(first, last, T(), SumDiffNthPower(mean)) / cnt;

}


template

T computeVariance(Iter_T first, Iter_T last, T mean) {

 return nthMoment(first, last, mean);

}


template

T computeStdDev(Iter_T first, Iter_T last, T mean) {

 return sqrt(computeVariance(first, last, mean));

}


template

T computeSkew(Iter_T begin, Iter_T end, T mean) {

 T m3 = nthMoment(begin, end, mean);

 T m2 = nthMoment(begin, end, mean);

 return m3 / (m2 * sqrt(m2));

}


template

T computeKurtosisExcess(Iter_T begin, Iter_T end, T mean) {

 T m4 = nthMoment(begin, end, mean);

 T m2 = nthMoment(begin, end, mean);

 return m4 / (m2 * m2) - 3;

}


template

void computeStats(Iter_T first, Iter_T last, T& sum, T& mean,

 T& var, T& std_dev, T& skew, T& kurt) {

 size_t cnt = distance(first, last);

 sum = accumulate(first, last, T());

 mean = sum / cnt;

 var = computeVariance(first, last, mean);

 std_dev = sort(var);

 skew = computeSkew(first, last, mean);

 kurt = computeKurtosisExcess(first, last, mean);

}


int main() {

 vector v;

 v.push_back(2);

 v.push_back(4);

 v.push_back(8);

 v.push_back(10);

 v.push_back(99);

 v.push_back(1);

 double sum, mean, var, dev, skew, kurt;

 computeStats(v.begin(), v.end(), sum, mean, var, dev, skew, kurt);

 cout << "count = " << v.size() << "\n";

 cout << "sum = " << sum << "\n";

 cout << "mean = " << mean << "\n";

 cout << "variance = " << var << "\n";

 cout << "standard deviation = " << dev << "\n";

 cout << "skew = " << skew << "\n";

 cout << "kurtosis excess = " << kurt << "\n";

 cout << endl;

}

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

count = 6

sum = 124

mean = 20.6667

variance = 1237.22

standard deviation = 35.1742

skew = 1.75664

kurtosis excess = 1.14171

Обсуждение

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


Табл. 11.1. Определения статистических функций

Статистическая функция Формула
n-й центральный момент (μn) ∑(xi-mean)n
Дисперсия μ2
Стандартное отклонение √μ2
Коэффициент асимметрии μ2/μ33/2
Эксцесс (μ4/μ2²)-3

Момент характеризует последовательность чисел. Другими словами, он определяет некий способ математического описания последовательности чисел. Моменты являются основой для расчета нескольких важных статистических функций, например дисперсии, стандартного отклонения, коэффициента асимметрии и эксцесса. Центральный момент — это момент, рассчитанный относительно среднего значения, а не нуля. Выборочный момент — это момент, рассчитанный для дискретного набора числовых значений, а не для всех значений функции. Нормализованный момент — это момент, поделенный на некоторую степень стандартного отклонения (стандартное отклонение рассчитывается как квадратный корень второго момента).

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

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

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

Пример 11.10 иллюстрирует работу функции 

accumulate
, показывая, как предоставленный пользователем функтор вызывается для каждого элемента последовательности.

Пример 11.10. Пример реализации функции accumulate

template

Iter_T accumulate(Iter_T begin, Iter_T end, Value_T value, BinOp_T op) {

 while (begin != end) {

  value = op(value, *begin++)

 }

 return value;

}

11.6. Генерация случайных чисел

Проблема

Требуется сгенерировать несколько случайных чисел в формате с плавающей точкой в интервале значений 

[0.0, 1.0)
при равномерном их распределении.

Решение

Стандарт C++ предусматривает наличие C-функции библиотеки этапа исполнения 

rand
, определенной в заголовочном файле
, которая возвращает случайное число в диапазоне от 0 до 
RAND_MAX
включительно. Макропеременная 
RAND_MAX
представляет собой максимальное значение, которое может быть возвращено функцией 
rand
. Пример 11.11 демонстрирует применение функции 
rand
для генерации случайных чисел с плавающей точкой.

Пример 11.11. Генерация случайных чисел функцией rand

#include 

#include 

#include 


using namespace std;


double doubleRand() {

 return double(rand()) / (double(RAND_MAX) + 1.0);

}


int main() {

 srand(static_cast(clock()));

 cout << "expect 5 numbers within the interval [0.0, 1.0)" << endl;

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

  cout << doubleRand() << "\n";

 }

 cout << endl;

}

Программа примера 11.11 должна выдать результат, подобный следующему.

expect 5 numbers within the interval [0.0, 1.0)

0.010437

0.740997

0.34906

0.369293

0.544373

Обсуждение

Необходимо уточнить, что функции, генерирующие случайные числа (в том числе 

rand
), возвращают псевдослучайные числа, а не реальные случайные числа, поэтому там, где я говорю «случайное число», я на самом деле имею в виду псевдослучайное число.

Перед применением функции 

rand
вы должны «посеять» (т.е. инициализировать) генератор случайных чисел с помощью вызова функции 
srand
. Это обеспечивает генерацию последующими вызовами 
rand
разных последовательностей чисел при каждом новом исполнении программы. Проще всего инициализировать генератор случайных чисел путем передачи ему результата вызова функции 
clock
из заголовочного файла
, имеющего тип 
unsigned int
. Повторная инициализация генератора случайных чисел приводит к тому, что генерируемые числа становятся менее случайными.

Функция 

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

Значительно более изощренную альтернативу 

rand
представляет написанная Джензом Маурером (Jens Maurer) библиотека Boost Random; она была инспирирована предложениями по генерации случайных чисел, представленными в TR1.

TR1 означает «Technical Report One» и представляет собой официальный проект по расширению стандартной библиотеки C++98.

Библиотека Boost Random содержит несколько высококачественных функций по генерации случайных чисел как для целых типов, так и для типов с плавающей точкой, причем с поддержкой многочисленных распределений. Пример 11.12 показывает, как можно сгенерировать случайные числа с плавающей точкой в интервале значений 

[0,1)
.

Пример 11.12. Использование библиотеки Boost Random

#include 

#include 

#include 


using namespace std;

using namespace boost;


typedef boost::mt19937 BaseGenerator;

typedef boost::uniform_real Distribution;

typedef boost::variate_generator Generator;


double boostDoubleRand() {

 static BaseGenerator base;

 static Distribution dist;

 static Generator rng(base, dist);

 return rng();

}


int main() {

 cout << "expect 5 numbers within the interval [0.1)" << endl;

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

  cout << boostDoubleRand() << "\n";

 }

 cout << endl;

}

Основное преимущество библиотеки Boost Random в том, что алгоритм генерации псевдослучайных чисел обеспечивает гарантированные и воспроизводимые свойства случайных последовательностей, зависящих от выбранного алгоритма. В примере 11.12 я использую генератор Mersenne Twister (

mt19937
), потому что он дает хорошее сочетание производительности и качества последовательности случайных чисел.

11.7. Инициализация контейнера случайными числами

Проблема

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

Решение

Можно использовать функции 

generate
и 
generate_n
из заголовочного файла
совместно с функтором, возвращающим случайные числа. Пример 11.13 показывает, как это можно сделать.

Пример 11.13. Инициализация контейнеров случайными числами

#include 

#include 

#include 

#include 

#include 


using namespace std;


struct RndIntGen {

 RndIntGen(int l, int h) : low(l), high(h) {}

 int operator()() const {

  return low + (rand() % ((high - low) + 1));

 }

private:

 int low;

 int high;

};


int main() {

 srand(static_cast(clock()));

 vector v(5);

 generate(v.begin(), v.end(), RndIntGen(1, 6));

 copy(v.begin(), v.end(), ostream_iterator(cout, "\n"));

}

Программа примера 11.13 должна выдать результат, подобный следующему.

3

1

2

6

4

Обсуждение

Стандартная библиотека C++ содержит функции 

generate
и 
generate_n
, специально предназначенные для заполнения контейнеров значениями, полученными функцией генератора случайных чисел. Эти функции принимают нуль-арный функтор (указатель на функцию или объект-функцию без аргументов), результат которого присваивается соседним элементам контейнера. Пример реализации функции 
generate
и 
generate_n
показан в примере 11.14.

Пример 11.14. Пример реализации функций generate и generate_n

template

void generate(Iter_T first, Iter_T last, Fxn_T f) {

 while (first != last) *first++ = f();

}


template

void generate_n(Iter_T first, int n, Fxn_T f) {

 for (int i=0; i < n; ++i) *first++ = f();

}

11.8. Представление динамического числового вектора

Проблема

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

Решение

Вы можете использовать шаблон 

valarray
из заголовочного файла
. Пример 11.15 показывает, как можно использовать шаблон 
valarray
.

Пример 11.15. Применение шаблона valarray

#include 

#include 


using namespace std;


int main() {

 valarray v(3);

 v[0] = 1;

 v[1] = 2;

 v[2] = 3;

 cout << v[0] << ", " << v[1] << ", " << v[2] << endl;

 v = v + v;

 cout << v[0] << ", " << v[1] << ", " << v[2] << endl;

 v /= 2;

 cout << v[0] << ", " << v[1] << ", " << v[2] << endl;

}

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

1, 2, 3

2, 4, 6

1, 2, 3

Обсуждение

Вопреки своему названию тип 

vector
не предназначен для использования в качестве числового вектора, для этой цели используется шаблон 
valarray
. Этот шаблон написан так, чтобы в конкретных реализациях С++, особенно на высокопроизводительных машинах, можно было применить к нему специальную векторную оптимизацию. Другое большое преимущество 
valarray
состоит в наличии многочисленных перегруженных операторов, предназначенных для работы с числовыми векторами. Эти операторы обеспечивают выполнение таких операций, как сложение и скалярное умножение векторов.

Шаблон 

valarray
может также использоваться в стандартных алгоритмах, работающих с массивами, представленными в C-стиле. Пример 11.16 показывает, как можно создавать итераторы, ссылающиеся на начальный элемент 
valarray
и на элемент, следующий за последним.

Пример 11.16. Получение итераторов для valarray

template

T* valarray_begin(valarray& x) {

 return &x[0];

}


template T* valarray_end(valarray& x) {

 return valarray_begin(x) + x.size();

}

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

valarray
, используя выражение 
&x[х.size()]
. Если это сработает, то только случайно, поскольку индексация 
valarray
, выходящая за допустимый индексный диапазон, приводит к непредсказуемому результату.

Отсутствие в 

valarray
функций-членов 
begin
и 
end
, несомненно, противоречит стилю STL. Отсутствие этих функций подчеркивает то, что в 
valarray
реализуется модель, отличная от концепции контейнера STL. Несмотря на это, вы можете использовать 
valarray
в любом обобщенном алгоритме, где требуется итератор с произвольным доступом.

11.9. Представление числового вектора фиксированного размера

Проблема

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

Решение

В программном обеспечении обычного типа часто более эффектный результат по сравнению с 

valarray
дает применение специальной реализации вектора, когда его размер заранее известен на этапе компиляции. Пример 11.17 показывает, как можно реализовать шаблон вектора фиксированного размера, названный здесь 
kvector
.

Пример 11.17. kvector.hpp

#include 

#include 


template

class kvector {

public:

 // открытые поля

 Value_T m[N];


 // открытые имена, вводимые typedef

 typedef Value_T value_type;

 typedef Value_T* iterator;

 typedef const Value_T* const_iterator;

 typedef Value_T& reference;

 typedef const Value_T& const_reference;

 typedef size_t size_type;


 // определение более короткого синонима для kvector

 typedef kvector self;


 // функции-члены

 template

 void copy(Iter_T first, Iter_T last) {

  copy(first, last, begin());

 }

 iterator begin() { return m; }

 iterator end() { return m + N; }

 const_iterator begin() const { return m; }

 const_iterator end() const { return m + N; }

 reference operator[](size_type n) { return m[n]; }

 const_reference operator[](size_type n) const { return m[n]; }

 static size_type size() { return N; }


 // векторные операции

 self& operator+=(const self& x) {

  for (int i=0; i

  return *this;

 }

 self& operator-=(const self& x) {

  for (int i=0; i

  return *this;

 }


 // скалярные операции

 self& operator=(value_type x) {

  std::fill(begin(), end(), x);

  return *this;

 }

 self& operator+=(value_type x) {

  for (int i=0; i

  return *this;

 }

 self& operator-=(value_type x) {

  for (int i=0; i

  return *this;

 }

 self& operator*=(value_type x) {

  for (int i=0; i

  return *this;

 }

 self& operator/=(value_type x) {

  for (int i=0; i

  return *this;

 }

 self& operator%=(value_type x) {

  for (int i=0; i

  return *this;

 }

 self operator-() {

 self x;

  for (int i=n; i

  return x;

 }


 // дружественные операторы

 friend self operator+(self x, const self& y) { return x += у; }

 friend self operator-(self x, const self& y) { return x -= y; }

 friend self operator+(self x, value_type y) { return x += y; }

 friend self operator-(self x, value_type y) { return x -= y; }

 friend self operator*(self x, value_type y) { return x *= y; }

 friend self operator/(self x, value_type y) { return x /= y; }

 friend self operator%(self x, value type y) { return x %= y; }

};

Пример 11.18 показывает, как можно применять шаблон класса 

kvector
.

Пример 11.18. Применение вектора kvector

#include "kvector.hpp"

#include 

#include 

#include 


using namespace std;


int main() {

 kvector v = { 1, 2, 3, 4 };

 cout << "sum = " << accumulate(v.begin(), v.end(), 0) << endl;

 v *= 3;

 cout << "sum = " << accumulated.begin(), v.end(), 0) << endl;

 v += 1;

 cout << "sum = " << accumulate(v.begin(), v.end(), 0) << endl;

}

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

sum = 10

sum = 30

sum = 34

Обсуждение

Представленный в примере 11.17 шаблон 

kvector
является гибридом 
valarray
и шаблона массива, предложенного в TR1. Как и 
valarray
, вектор 
kvector
представляет собой последовательность значений заданного числового типа, однако подобно массиву 
TR1::array
его размер известен на этапе компиляции.

Характерной особенностью шаблона 

kvector
является то, что для его инициализации может использоваться синтаксис, применяемый для массивов, и то, что он имеет функции-члены 
begin
и 
end
. Фактически 
kvector
можно рассматривать как псевдоконтейнер, т.е. он удовлетворяет некоторым, но не всем требованиям концепции стандартного контейнера. Следствие этого — более легкое применение 
kvector
в стандартных алгоритмах по сравнению с 
valarray
.

Другое преимущество шаблонного класса 

kvector
состоит в том, что он поддерживает синтаксис, используемый при инициализации массивов.

int x;

kvector k = { x = 1, x+2, 5}

Этот синтаксис возможен только потому, что 

kvector
является агрегатом. Агрегат (aggregate) — это массив или класс, который не имеет объявленных пользователем конструкторов, закрытых или защищенных данных-членов, базового класса и виртуальных функций. Следует отметить, что все же можно при объявлении 
kvector
его заполнить значениями по умолчанию.

kvector k = {};

В результате этот вектор будет заполнен нулями.

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

array
, удовлетворяющего требованиям TR1.

Возможно, самое большое преимущество 

kvector
над реализациями динамического вектора проявляется в его высокой производительности. По двум причинам шаблон kvector значительно эффективнее, чем большинство реализаций динамических 
векторов
: компиляторы очень хорошо справляются с оптимизацией циклов фиксированною размера, и здесь нет динамического распределения памяти. Различия в производительности особенно проявляются при работе с небольшими матрицами (например, 2×2 или 3×3), которые часто встречаются во многих приложениях.

Что означает имя «self», введенное оператором typedef?

Введенное с помощью typedef имя 

self
я использую в примере 11.17 и в последующих примерах; оно представляет собой удобное краткое имя, которое я использую для ссылки на тип текущего класса. Программу значительно легче писать и воспринимать при использовании self вместо имени класса.

11.10. Вычисление скалярного произведения

Проблема

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

Решение

Пример 11.19 показывает, как можно вычислить скалярное произведение, используя функцию 

inner_product
из заголовочного файла
.

Пример 11.19. Расчет скалярного произведения

#include 

#include 

#include 


using namespace std;


int main() {

 int v1[] = { 1, 2, 3 };

 int v2[] = { 4, 6, 8 };

 cout << "the dot product of (1,2,3) and (4,6,8) is ";

 cout << inner_product(v1, v1 + 3, v2, 0) << endl;

}

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

the dot product of (1,2,3) and (4,6,8) is 40

Обсуждение

Скалярное произведение (dot product) является одной из форм обобщенного скалярного произведения (inner product), называемой евклидовым скалярным произведением (Euclidean Inner Product). Функция 

inner_product
объявляется следующим образом.

template

T inner_product(In first, In last, In2 first2, T init);

template

T inner_product(In first, In last, In2 first2, T init, BinOp op, BinOp2 op2);

Первый вариант функции 

inner_product
суммирует произведения соответствующих элементов двух контейнеров. Второй вариант функции 
inner_product
позволяет вам самому предоставить операцию над парой чисел и функцию суммирования. В примере 11.20 продемонстрирована простая реализация функции 
inner_product
.

Пример 11.20. Пример реализации функции inner_product()

template

T inner_product(In first, In last, In2 first2, T init, BinOp op, Binop2 op2) {

 while (first != last) {

  BinOp(init, BinOp2(*first++, *first2++));

 }

 return init;

}

Благодаря гибкости реализации функции 

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

Смотри также

Рецепты 11.11 и 11.12.

11.11. Вычисление нормы вектора

Проблема

Требуется найти норму (т. е. длину) числового вектора.

Решение

Можно использовать функцию 

inner_product
из заголовочного файла
для умножения вектора на самого себя, как показано в примере 11.21.

Пример 11.21. Вычисление нормы вектора

#include 

#include 

#include 

#include 


using namespace std;


template

long double vectorNorm(Iter_T first, Iter_T last) {

 return sqrt(inner_product(first, last, first, 0.0L));

}


int main() {

 int v[] = { 3, 4 };

 cout << "The length of the vector (3.4) is ";

 cout << vectorNorm(v, v + 2) << endl;

}

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

The length of the vector (3,4) is 5

Обсуждение

В примере 11.21 функция 

inner_product
из заголовочного файла
используется для вычисления скалярного произведения числового вектора на самого себя. Квадратный корень полученного значения, как известно, является нормой вектора, или длиной вектора.

Вместо того чтобы в функции 

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

11.12. Вычисление расстояния между векторами

Проблема

Требуется найти евклидово расстояние между векторами.

Решение

Евклидово расстояние между векторами определяется как квадратный корень суммы квадратов разностей соответствующих элементов. Рассчитать его можно так, как показано в примере 11.22.

Пример 11.22. Расчет расстояния между двумя векторами

#include 

#include 


using namespace std;


template

double vectorDistance(Iter_T first, Iter_T last, Iter2_T first2) {

 double ret = 0.0;

 while (first != last) {

  double dist = (*first++) - (*first2++);

  ret += dist * dist;

 }

 return ret > 0.0 ? sqrt(ret) : 0.0;

}


int main() {

 int v1[] = { 1, 5 };

 int v2[] = { 4, 9 };

 cout << "distance between vectors (1,5) and (4,9) is ";

 cout << vectorDistance(v1, v1 + 2, v2) << endl;

}

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

distance between vectors (1,5) and (4,9) is 5

Обсуждение

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

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

Пример 11.23. Расчет расстояния между векторами с использованием функции inner_product

#include 

#include 

#include 

#include 


using namespace std;


template

struct DiffSquared {

 Value_T operator()(Value_T x, Value_T y) const {

  return (x - y) * (x - y);

 }

};


template

double vectorDistance(Iter_T first, Iter_T last, Iter2_T first2) {

 double ret = inner_product(first, last, first2, 0.0L,

 plus(), DiffSquared());

 return ret > 0.0 ? sqrt(ret) : 0.0;

}


int main() {

 int v1[] = { 1, 5 };

 int v2[] = { 4, 9 };

 cout << "distance between vectors (1,5) and (4,9) is ";

 cout << vectorDistance(v1, v1 + 2, v2) << endl;

}

Поскольку реализация функции 

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

11.13. Реализация итератора с шагом

Проблема

Имеются смежные числовые ряды и требуется обеспечить доступ к каждому n-му элементу.

Решение

В примере 11.24 представлен заголовочный файл, реализующий класс итератора с шагом.

Пример 11.24. stride_iter.hpp

#ifndef STRIDE_ITER_HPP

#define STRIDE_ITER_HPP


#include 

#include 


template

class stride_iter {

public:


 // открытые имена, вводимые typedef

 typedef typename std::iterator_traits::value_type value_type;

 typedef typename std::iterator_traits::reference reference;

 typedef typename std::iterator_traits::difference_type

  difference_type;

 typedef typename std::iterator_traits::pointer pointer;

 typedef std::random_access_iterator_tag iterator_category;

 typedef stride_iter self;


 // конструкторы

 stride_iter() : m(NULL), step(0) {};

 stride_iter(const self& x) : m(x.m), step(x.step) {}

 stride_iter(Iter_T x, difference_type n) : m(x), step(n) {}


 // операторы

 self& operator++() { m += step; return *this; }

 self operator++(int) { self tmp = *this; m += step; return tmp; }

 self& operator+=(difference_type x) { m += x * step; return *this; }

 self& operator--() { m -= step; return *this; }

 self operator--(int) { self tmp = *this; m -= step; return trap; }

 self& operator--(difference type x) { m -= x + step; return *this; }

 reference operator[](difference_type n) { return m[n * step]; }

 reference operator*() { return *m; }


 // дружественные операторы

 friend bool operator==(const self& x, const self& y) {

  assert(x.step == y.step);

  return x.m == y.m;

 }

 friend bool operator!=(const self& x, const self& y) {

  assert(x.step == y.step);

  return x.m != y.m;

 }

 friend bool operator<(const self& x, const self& y) {

  assert(x.step == y.step);

  return x.m < y.m;

 }

 friend difference type operator-(const self& x, const self& y) {

  assert(x.step == y.step);

  return (x.m - y.m) / x.step;

 }

 friend self operator+(const self& x, difference_type y) {

  assert(x.step == y.step);

  return x += y * x.step;

 }

 friend self operator+(difference_type x, const self& y) {

  assert(x.step == y.step);

  return y += x * x.step;

 }

private:

 Iter_T m;

 difference_type step;

};


#endif

Пример 11.25 показывает, как можно использовать итератор 

stride_iter
из примера 11.24 для получения доступа к каждому второму элементу последовательности.

Пример 11.25. Применение итератора stride_iter

#include "stride_iter.hpp"

#include 

#include 

#include 


using namespace std;


int main() {

 int a[] = { 0, 1, 2, 3, 4, 5, 6, 7 };

 stride_iter first(a, 2);

 stride_iter last(a + 8, 2);

 copy(first, last, ostream_iterator(cout, "\n"));

}

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

0

2

4

6

Обсуждение

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

Я хотел сделать итератор с шагом совместимым с STL, поэтому пришлось выбрать подходящий тип стандартного итератора и удовлетворить его требования. Представленный в примере 11.24 итератор с шагом сделан по образцу итератора с произвольным доступом.

В примере 11.26 я отдельно привел реализацию итератора с шагом (названную 

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

Пример 11.26. kstride_iter.hpp

#ifndef KSTRIDE_ITER_HPP

#define KSTRIDE_ITER_HPP


#include 


template

class kstride_iter {

public:

 // открытые имена, вводимые typedef

 typedef typename std::iterator_traits::value_type value_type;

 typedef typename std::iterator_traits::reference reference;

 typedef typename std::iterator_traits::difference_type

  difference_type;

 typedef typename std::iterator_traits::pointer pointer;

 typedef std::random_access_iterator_tag iterator_category;

 typedef kstride_iter self;


 // конструкторы

 kstride_iter() : m(NULL) {} kstride_iter(const self& x) : m(x.m) {}

 explicit kstride_iter(Iter_T x) : m(x) {}


 // операторы

 self& operator++() { m += Step_N; return *this; }

 self operator++(int) { self tmp = *this; m += Step_N; return tmp; }

 self& operator+=(difference_type x) { m += x * Step_N; return *this; }

 self& operator--() { m -= Step_N; return *this; }

 self operator--(int) { self tmp = *this; m -= Step_N; return tmp; }

 self& operator--(difference_type x) { m -= x * Step_N; return *this; }

 reference operator[](difference_type n) { return m[n * Step_N]; }

 reference operator*() { return *m; }


 // дружественные операторы

 friend bool operator==(self x, self y) { return x.m == y.m; }

 friend bool operator!=(self x, self y) { return x.m != y.m; }

 friend bool operator<(self x, self y) { return x.m < y.m; }

 friend difference_type operator-(self x, self y) {

 return (x.m - y.m) / Step_N;

 }


 friend self operator+(self x, difference_type y) { return x += y * Step_N; }

 friend self operator+(difference_type x, self y) { return y += x * Step_N; }

private:

 Iter_T m;

};


#endif

Пример 11.27 показывает, как можно использовать итератор 

kstride_iter
.

Пример 11.27. Применение итератора kstride_iter

#include "kstride_iter.hpp"

#include 

#include 

#include 


using namespace std;


int main() {

 int a[] = { 0, 1, 2, 3, 4, 5, 6, 7 };

 kstride_iter first(a);

 kstride_iter last(a + 8);

 copy(first, last, ostream_iterator(cout, "\n"));

}

11.14. Реализация динамической матрицы

Проблема

Требуется реализовать числовые матрицы, размерности которых (количество строк и столбцов) неизвестны на этапе компиляции.

Решение

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

valarray
.

Пример 11.28. matrix.hpp

#ifndef MATRIX_HPP

#define MATRIX_HPP


#include "stride_iter.hpp" // см. рецепт 11.12

#include 

#include 

#include 


template

class matrix {

public:

 // открытые имена, вводимые typedef

 typedef Value_T value_type;

 typedef matrix self;

 typedef value_type* iterator;

 typedef const value_type* const_iterator;

 typedef Value_T* row_type;

 typedef stride_iter col_type;

 typedef const value_type* const_row_type;

 typedef stride_iter const_col_type;


 // конструкторы

 matrix() : nrows(0), ncols(0), m() {}

 matrix(int r, int c) : nrows(r), ncols(c), m(r * с) {}

 matrix(const self& x) : m(x.m), nrows(x.nrows), ncols(x.ncols) {}

 template

 explicit matrix(const valarray& x)

  : m(x.size() + 1), nrows(x.size()), ncols(1) {

 for (int i=0; i

 }

 // позволить конструирование из матриц других типов

 template explicit matrix(const matrix& x)

  : m(x.size() + 1), nrows(x.nrows), ncols(x.ncols) {

  copy(x.begin(), x.end(), m.begin());

 }


// открытые функции

 int rows() const { return nrows; }

 int cols() const { return ncols; }

 int size() const { return nrows * ncols; }


 // доступ к элементам

 row_type row begin(int n) { return &m[n * cols()]; }

 row_type row_end(int n) { return row_begin() + cols(); }

 col_type col_begin(int n) { return col_type(&m[n], cols()); }

 col_type col_end(int n) { return col_begin(n) + cols(); }

 const_row_type row_begin(int n) const { return &m[n * cols()]; }

 const_row_type row_end(int n) const { return row_begin() + cols(); }

 const_col_type col_begin(int n) const { return col_type(&m[n], cols()); }

 const_col_type col_end(int n) const { return col_begin() + cols(); }

 iterator begin() { return &m[0]; }

 iterator end() { return begin() + size(); }

 const_iterator begin() const { return &m[0]; }

 const_iterator end() const { return begin() + size(); }


 // операторы

 self& operator=(const self& x) {

  m = x.m;

  nrows = x.nrows;

  ncols = x.ncols;

  return *this;

 }

 self& operator=(value_type x) { m = x; return *this; }

 row_type operator[](int n) { return row_begin(n); }

 const_row_type operator[](int n) const { return row_begin(n); }

 self& operator+=(const self& x) { m += x.m; return *this; }

 self& operator-=(const self& x) { m -= x.m; return *this; }

 self& operator+=(value_type x) { m += x; return *this; }

 self& operator-=(value_type x) { m -= x; return *this; }

 self& operator*=(value_type x) { m *= x; return *this; }

 self& operator/=(value_type x) { m /= x; return *this; }

 self& operator%=(value_type x) { m %= x; return *this; }

 self operator-() { return -m; }

 self operator+() { return +m; }

 self operator!() { return !m; }

 self operator~() { return ~m; }


 // дружественные операторы

 friend self operator+(const self& x, const self& y) { return self(x) += y; }

 friend self operator-(const self& x, const self& y) { return self(x) -= y; }

 friend self operator+(const self& x, value_type y) { return self(x) += y; }

 friend self operator-(const self& x, value_type y) { return self(x) -= y; }

 friend self operator*(const self& x, value type y) { return self(x) *= y; }

 friend self operator/(const self& x, value_type y) { return self(x) /= y; }

 friend self operator%(const self& x, value_type y) { return self(x) %= y; }

private:

 mutable valarray m;

 int nrows;

 int ncols;

};

#endif

Пример 11.29 показывает, как можно использовать шаблонный класс 

matrix
.

Пример 11.29. Применение шаблона matrix

#include "matrix.hpp"

#include 


using namespace std;


int main() {

 matrix m(2,2);

 m = 0;

 m[0][0] = 1;

 m[1][1] = 1;

 m *= 2;

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

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

}

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

(2,0)

(0,2)

Обсуждение

Проект шаблона матрицы, представленный в примере 11.28, в значительной степени инспирирован шаблоном матрицы Бьерна Страуструпа (Bjarne Stroustrup) из его книги «The C++ Programming Language», 3-е издание (издательство «Addison Wesley»). Реализация Страуструпа отличается тем, что его итератор использует класс 

slice
и указатель на 
valarray
для индексации. Реализованная в примере 11.27 матрица использует вместо них итератор с шагом из рецепта 11.12, что делает итераторы более компактными и при некоторых реализациях более эффективными.

Шаблонный класс 

matrix
позволяет индексировать элемент i-й строки и j-го столбца, используя операцию двойной индексации. Например:

matrix m(100,100);

cout << "the element at row 24 and column 42 is " << m[24][42] << endl;

Шаблонный класс 

matrix
также имеет функции-члены 
begin
и 
end
, т.е. его легко можно использовать в различных алгоритмах STL.

Пример 11.28 содержит строку, которая, возможно, вызывает у вас некоторое удивление. Имеется в виду следующее объявление.

mutable valarray m;

Объявление поля-члена 

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

Смотри также

Рецепты 11.15 и 11.16.

11.15. Реализация статической матрицы

Проблема

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

Решение

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

Пример 11.30. kmatrix.hpp

#ifndef KMATRIX_HPP

#define KMATRIX_HPP


#include "kvector.hpp"

#include "kstride_iter.hpp"


template

class kmatrix {

public:

 // открытые имена, вводимые typedef

 typedef Value_T value_type;

 typedef kmatrix self;

 typedef Value_T* iterator;

 typedef const Value_T* const_iterator;

 typedef kstride_iter row_type;

 typedef kstride_iter col_type;

 typedef kstride_iter const_row_type;

 typedef kstride_iter const_col_type;


 // открытые константы

 static const int nRows = Rows_N;

 static const int nCols = Cols_N;


 // конструкторы

 kmatrix() { m = Value_T(); }

 kmatrix(const self& x) { m = x.m; }

 explicit kmatrix(Value_T& x) { m = x.m; }


 // открытые функции

 static int rows() { return Rows_N; }

 static int cols() { return Cols_N; }

 row_type row(int n) { return row_type(begin() * (n * Cols_N)); }

 col_type col(int n) { return col_type(begin() + n); }

 const_row_type row(int n) const {

  return const_row_type(begin() + (n * Cols_N));

 }

 const_col_type col(int n) const {

  return const_col_type(begin() + n);

 }

 iterator begin() { return m.begin(); }

 iterator end() { return m.begin() + size(); }

 const_iterator begin() const { return m; }

 const_iterator end() const { return m + size(); }

 static int size() { return Rows_N * Cols_N; }


 // операторы

 row_type operator[](int n) { return row(n); }

 const_row_type operator[](int n) const { return row(n); }


 // операции присваивания

 self& operator=(const self& x) { m = x.m; return *this; }

 self& operator=(value_type x) { m = x; return *this; }

 self& operator+=(const self& x) { m += x.m; return *this; }

 self& operator-=(const self& x) { m -= x.m; return *this; }

 self& operator+={value_type x) { m += x; return *this; }

 self& operator-=(value_type x) { m -= x; return *this; }

 self& operator*=(value_type x) { m *= x; return *this; }

 self& operator/=(value_type x) { m /= x; return *this; }

 self operator-() { return self(-m); }


 // друзья

 friend self operator+(self x, const self& у) { return x += y; }

 friend self operator-(self x, const self& y) { return x -= y; }

 friend self operator+(self x, value_type y) { return x += y; }

 friend self operator-(self x, value type y) { return x -= y; }

 friend self operator*(self x, value_type y) { return x *= y; }

 friend self operator/(self x, value_type y) { return x /= y; }

 friend bool operator==(const self& x, const self& y) { return x == y; }

 friend bool operator!=(const self& x, const self& y) { return x.m != y.m; }

private:

 kvector m;

};


#endif

В примере 11.31 приведена программа, демонстрирующая применение шаблонного класса 

kmatrix
.

Пример 11.31. Применение kmatrix

#include "kmatrix.hpp"

#include 


using namespace std;


template

void outputRowOrColumn(Iter_T iter, int n) {

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

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

 }

 cout << endl;

}


template

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

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 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 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 

#include 


using namespace std;


template

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 m1(2, 1);

 matrix m2(1, 2);

 kmatrix 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
, определенного в заголовочном файле
.

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

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

В примере 11.32 реализуется равенство 

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

Смотри также

Рецепт 11.17.

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

Проблема

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

Решение

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

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

#include 

#include 

#include 

#include 


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

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

 typedef typename iterator_traits::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 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
из заголовочного файла
содержит функции преобразования в полярные координаты и обратно. Пример 11.34 показывает, как можно использовать класс шаблона complex для представления и манипулирования полярными координатами.

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

#include 

#include 


using namespace std;


int main() {

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

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

 complex coord = polar(rho, theta);

 cout << "rho = " << abs(coord) << ", theta = " << arg(coord) << endl;

 coord += polar(4.0, 0.0);

 cout << "rho = " << abs(coord) << ", theta = " << arg(coord) << endl;

}

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

rho = 3, theta = 1.5708

rho = 5, theta = 0.643501

Обсуждение

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

complex
для представления полярных координат не является элегантным решением, я предусмотрел приведенный в примере 11.25 класс полярных координат, допускающий более естественное применение.

Пример 11.35. Класс полярных координат

#include 

#include 


using namespace std;


template

struct BasicPolar {

 public typedef BasicPolar self;


 // конструкторы

 BasicPolar() : m() {} BasicPolar(const self& x) : m(x.m) {}

 BasicPolar(const T& rho, const T& theta) : m(polar(rho, theta)) {}


 // операторы присваивания

 self operator-() { return Polar(-m); }

 self& operator+=(const self& x) { m += x.m; return *this; }

 self& operator-=(const self& x) { m -= x.m; return *this; }

 self& operator*=(const self& x) { m *= x.m; return *this; }

 self& operator/=(const self& x) { m /= x.m; return *this; }

 operator complex() const { return m; }


 // открытые функции-члены

 T rho() const { return abs(m); }

 T theta() const { return arg(m); }


 // бинарные операции

 friend self operator+(self x, const self& y) { return x += y; }

 friend self operator-(self x, const self& y) { return x -= y; }

 friend self operator*(self x, const self& y) { return x *= y; }

 friend self operator/(self x, const self& y) { return x /= y; }


 // операторы сравнения

 friend bool operator==(const self& x, const self& y) { return x.m == y.m; }

 friend bool operator!=(const self& x, const self& y) { return x.m ! = y.m; }

private:

 complex m;

};


typedef BasicPolar Polar;


int main() {

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

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

 Polar coord(rho, theta);

 cout << "rho = " << coord.rho() << ", theta = " << coord.theta() << endl;

 coord += Polar(4.0, 0.0);

 cout << "rho = " << coord.rho() << ", theta = " << coord.theta() << endl;

 system("pause");

}

В примере 11.35 с помощью 

typedef
я определил тип 
Polar
как специализацию шаблона 
BasicPolar
. Так удобно определять используемый по умолчанию тип, однако вы можете при необходимости специализировать шаблон 
BasicPolar
другим числовым типом. Такой подход используется в стандартной библиотеке в отношении классе 
string
, который является специализацией шаблона 
basic_string
.

11.19. Выполнение операций с битовыми наборами

Проблема

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

Решение

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

bitset
из заголовочного файла
, рассматривая его как целый тип без знака.

Пример 11.36. bitset_arithmetic.hpp

#include 

#include 


bool fullAdder(bool b1, bool b2, bool& carry) {

 bool sum = (b1 ^ b2) ^ carry;

 carry = (b1 && b2) || (b1 && carry) || (b2 && carry);

 return sum;

}


bool fullSubtractor(bool b1, bool b2, bool& borrow) {

 bool diff;

 if (borrow) {

  diff = !(b1 ^ b2);

  borrow = !b1 || (b1 && b2);

 } else {

  diff = b1 ^ b2;

  borrow = !b1 && b2;

 }

 return diff;

}


template

bool bitsetLtEq(const std::bitset& x, const std::bitset& y) {

 for (int i=N-1; i >= 0; i--) {

  if (x[i] && !y[i]) return false;

  if (!x[i] && y[i]) return true;

 }

 return true;

}


template

bool bitsetLt(const std::bitset& x, const std::bitset& y) {

 for (int i=N-1; i >= 0, i--) {

  if (x[i] && !y[i]) return false;

  if (!x[i] && y[i]) return true;

 }

 return false;

}


template

bool bitsetGtEq(const std::bitset& x, const std::bitset& y) {

 for (int i=N-1; i >= 0; i--) {

  if (x[i] && !y[i]) return true;

  if (!x[i] && y[i]) return false;

 }

 return true;

}


template

bool bitsetGt(const std::bitset& x, const std::bitset& y) {

 for (int i=N-1; i >= 0; i--) {

  if (x[i] && !y[i]) return true;

  if (!x[i] && y[i]) return false;

 }

 return false;

}


template

void bitsetAdd(std::bitset& x, const std::bitset& y) {

 bool carry = false;

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

  x[i] = fullAdder(x[i], y[x], carry);

 }

}


template

void bitsetSubtract(std::bitset& x, const std::bitset& y) {

 bool borrow = false;

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

  if (borrow) {

  if (x[i]) {

   x[i] = y[i];

   borrow = y[i];

  } else {

   x[i] = !y[i];

   borrow = true;

  }

  } else {

  if (x[i]) {

   x[i] = !y[i];

   borrow = false;

  } else {

   x[i] = y[i];

   borrow = y[i];

  }

  }

 }

}


template

void bitsetMultiply(std::bitset& x, const std::bitset& y) {

 std::bitset tmp = x;

 x.reset();

 // мы хотим минимизировать количество операций сдвига и сложения

 if (tmp.count() < y.count()) {

  for (int i=0; i < N; i++) if (tmp[i]) bitsetAdd(x, у << i);

 } else {

  for (int i=0; i < N; i++) if (y[i]) bitsetAdd(x, tmp << i);

 }

}


template

void bitsetDivide(std::bitset x, std::bitset y,

 std::bitset& q, std::bitset& r) {

 if (y.none()) {

  throw std::domain_error("division by zero undefined");

 }

 q.reset();

 r.reset();

 if (x.none()) {

  return;

 }

 if (x == y) {

  q[0] = 1;

  return;

 }

 r = x;

 if (bitsetLt(x, y)) {

  return;

 }

 // подсчитать количество значащих цифр в делителе и делимом

 unsigned int sig_x;

 for (int i=N-1; i>=0; i--) {

  sig_x = i;

  if (x[i]) break;

 }

 unsigned int sig_y;

 for (int i=N-1; i>=0; i--) {

  sig_y = i;

  if (y[i]) break;

 }

 // выровнять делитель по отношению к делимому

 unsigned int n = (sig_x — sig_y);

 y <<= n;

 // обеспечить правильное число шагов цикла

 n += 1;

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

 while (n--) {

  // сдвинуть частное влево

  if (bitsetLtEq(y, r)) {

  // добавить новую цифру к частному

  q[n] = true;

  bitset.Subtract(r, y);

  }

  // сдвинуть делитель вправо

  y >>= 1;

 }

}

Пример 11.37 показывает, как можно использовать заголовочный файл bitset_arithmetic.hpp.

Пример 11.37. Применение функций bitset_arithmetic.hpp

#include "bitset_arithmetic.hpp"

#include 

#include 

#include 


using namespace std;


int main() {

 bitset<10> bits1(string("100010001"));

 bitset<10> bits2(string("000000011"));

 bitsetAdd(bits1, bits2);

 cout << bits1.to_string, allocator >() << endl;

}

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

0100010100

Обсуждение

Шаблон класса 

bitset
содержит основные операции по манипулированию битовыми наборами, но не обеспечивает арифметические операции и операции сравнения. Это объясняется тем, что в библиотеке нельзя заранее точно предвидеть, какой числовой тип будет использоваться для представления произвольного битового набора согласно ожиданиям программиста.

В функциях примера 11.36 считается, что 

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

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

Смотри также

Рецепт 11.20.

11.20. Представление больших чисел фиксированного размера

Проблема

Требуется выполнить операции с числами, размер которых превышает размер типа 

long int
.

Решение

Шаблон 

BigInt
в примере 11.38 использует 
bitset
из заголовочного файла
для того, чтобы можно было представить целые числа без знака в виде набора бит фиксированного размера, причем количество бит определяется параметром шаблона.

Пример 11.38. big_int.hpp

#ifndef BIG_INT_HPP

#define BIG_INT_HPP


#include 

#include "bitset_arithmetic.hpp" // Рецепт 11.20


template

class BigInt {

 typedef BigInt self;

public:

 BigInt() : bits() {}

 BigInt(const self& x) : bits(x.bits) {}

 BigInt(unsigned long x) {

 int n = 0;

  while (x) {

  bits[n++] = x & 0x1;

  x >>= 1;

  }

 }

 explicit BigInt(const std::bitset& x) bits(x) {}


 // открытые функции

 bool operator[](int n) const { return bits[n]; }

 unsigned long toUlong() const { return bits.to_ulong(); }


 // операторы

 self& operator<<=(unsigned int n) {

 bits <<= n;

  return *this;

 }

 self& operator>>=(unsigned int n) {

  bits >>= n;

  return *this;

 }

 self operator++(int) {

  self i = *this;

  operator++();

  return i;

 }

 self operator--(int) {

  self i = *this;

  operator--();

  return i;

 }

 self& operator++() {

  bool carry = false;

  bits[0] = fullAdder(bits[0], 1, carry);

  for (int i = 1; i < N; i++) {

  bits[i] = fullAdder(bits[i], 0, carry);

  }

  return *this;

 }

 self& operator--() {

  bool borrow = false;

  bits[0] = fullSubtractor(bits[0], 1, borrow);

  for (int i = 1; i < N; i++) {

  bits[i] = fullSubtractor(bits[i], 0, borrow);

  }

  return *this;

 }

 self& operator+=(const self& x) {

  bitsetAdd(bits, x.bits);

  return *this;

 }

 self& operator-=(const self& x) {

  bitsetSubtract(bits, x.bits);

  return *this;

 }

 self& operator*=(const self& x) {

  bitsetMultiply(bits, x.bits);

  return *this;

 }

 self& operator/=(const self& x) {

  std::bitset tmp;

  bitsetDivide(bits, x.bits, bits, tmp);

  return *this;

 }

 self& operator%=(const self& x) {

  std::bitset tmp;

  bitsetDivide(bits, x.bits, tmp, bits);

  return *this;

 }

 self operator~() const { return ~bits; }

 self& operator&=(self x) { bits x.bits; return *this; }

 self& operator|=(self x) { bits x.bits; return *this; }

 self& operator~=(self x) { bits ~= x.bits; return *this; }


 // дружественные функции

 friend self operator<<(self x, unsigned int n) { return x <<= n; }

 friend self operator>>(self x, unsigned int n) { return x >>= n; }

 friend self operator+(self x, const self& y) { return x += y; }

 friend self operator-(self x, const self& y) { return x -= y; }

 friend self operator*(self x, const self& y) { return x *= y; }

 friend self operator/(self x, const self& y) { return x /= y; }

 friend self operator%(self x, const self& y) { return x %= y; }

 friend self operator^(self x, const self& y) { return x ^= y; }

 friend self operator&(self x, const self& y) { return x &= y; }

 friend self operator|(self x, const self& y) { return x |= y; }


 // операторы сравнения

 friend bool operator==(const self& x, const self& y) {

  return x.bits == y.bits;

 }

 friend bool operator!=(const self& x, const self& y) {

  return x.bits ! = y.bits;

 }

 friend bool operator>(const self& x, const self& y) {

  return bitsetGt(x.bits, y.bits);

 }

 friend bool operator<(const self& x, const self& y) {

  return bitsetLt(x.bits, y.bits);

 }

 friend bool operator>=(const self& x, const self& y) {

  return bitsetGtEq(x.bits, y.bits);

 }

 friend bool operator<=(const self& x, const self& y) {

  return bitsetLtEq(x bits, y.bits);

 }

private:

 std::bitset bits;

};


#endif

Шаблонный класс 

BigInt
можно использовать для вычисления факториалов, как показано в примере 11.39.

Пример 11.39. Применение класса big_int

#include "big_int.hpp"

#include 

#include 

#include 

#include 


using namespace std;


void outputBigInt(BigInt<1024> x) {

 vector v;

 if (x == 0) {

  cout << 0;

  return;

 }

 while (x > 0) {

  v.push_back((x % 10).to_ulong());

  x /= 10;

 }

 copy(v.rbegin(), v.rend(), ostream_iterator(cout, ""));

 cout << endl;

}


int main() {

 BigInt<1024> n(1); // вычислить факториал числа 32

 for (int 1=1; i <= 32; ++i) {

  n *= i;

 }

 outputBigInt(n);

}

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

263130836933693530167218012160000000

Обсуждение

Большие целые числа часто встречаются во многих приложениях. Например, в криптографии нередки числа, которые представляются 1000 и более битами. Однако современный стандарт C++ позволяет работать как максимум с типом 

long int
.

Число бит типа 

long int
зависит от реализации, но оно не может быть меньше 32. И едва ли это число будет больше 1000. Следует помнить, что один из этих битов используется в качестве знака.

Ожидается, что новая версия стандарта (C++0x) последует за стандартом C99 и предусмотрит тип 

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

Представленная здесь реализация основана на двоичном представлении чисел при помощи класса 

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

Смотри также

Рецепт 11.19.

11.21. Реализация чисел с фиксированной точкой

Проблема

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

Решение

В примере 11.40 представлена реализация вещественного числа с фиксированной точкой, когда количество двоичных позиций справа от точки задается параметром шаблона. Например, тип 

basic_fixed_real<10>
имеет 10 двоичных цифр справа от двоичной точки, что позволяет представлять числа с точностью до 1/1024.

Пример 11.40. Представление вещественных чисел, используя формат с фиксированной точкой

#include 


using namespace std;


template

struct BasicFixedReal {

 typedef BasicFixedReal self;

 static const int factor = 1 << (E - 1);

 BasicFixedReal() : m(0) {}

 BasicFixedReal(double d) : m(static_cast(d * factor)) {}

 self& operator+=(const self& x) { m += x.m; return *this; }

 self& operator-=(const self& x) { m -= x.m; return *this; }

 self& operator*=(const self& x) { m *= x.m; m >>=E; return *this; }

 self& operator/=(const self& x) { m /= x.m; m *= factor; return *this; }

 self& operator*=(int x) { m *= x; return *this; }

 self& operator/=(int x) { m /= x; return *this; }

 self operator-() { return self(-m); }

 double toDouble() const { return double(m) / factor; }


 // дружественные функции

 friend self operator+(self x, const self& v) { return x += y; }

 friend self operator-(self x, const self& y) { return x -= y; }

 friend self operator-(self x, const self& y) { return x *= y; }

 friend self operator/(self x, const self& y) { return x /= y; }


 // операторы сравнения

 friend bool operator==(const self& x, const self& y) { return x.m == y.m; }

 friend bool operator!=(const self& x, const self& y) { return x.m != y.m; }

 friend bool operator>(const self& x, const self& y) { return x.m > y.m; }

 friend bool operator<(const self& x, const self& y) { return x.m < y.m; }

 friend bool operator>=(const self& x, const self& y) { return x.m >= y.m; }

 friend bool operator<=(const self& x, const self& y) { return x.m <= y.m; }

private:

 int m;

};


typedef BasicFixedReal<10> FixedReal;


int main() {

 FixedReal x(0);

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

  x += FixedReal(0.0625);

 }

 cout << x.toDouble() << endl;

}

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

6.25

Обсуждение

Число с фиксированной точкой, как и число с плавающей точкой, является приблизительным представлением вещественного числа. Число с плавающей точкой имеет мантиссу (m) и экспоненту (е), обеспечивая значение, вычисляемое по формуле m*bе, где b — некоторая константа.

Число с фиксированной точкой имеет почти такой же формат, но здесь экспонента также фиксирована. Эта константа в примере 11.40 передается шаблону 

basic_fixed_real
в качестве его параметра.

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

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

fixed_real
позволяет представлять значения только в диапазоне от -2 097 151 до +2 097 151 с точностью 1/1024.

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

ТЕЛЕГРАМ

Канал с обзорами, анонсами новинок и книжными подборками

Книжный Вестник

Бот для удобного поиска книг (если не нашлось на сайте)

Поиск книг

Свежие любовные романы в удобных форматах

Любовные романы

О психологии, саморазвитии и личностном росте

Саморазвитие

Детективы и триллеры, все новинки

Детективы

Фантастика и фэнтези, все новинки

Фантастика

Отборные классические книги

Классика
ВКОНТАКТЕ

Цитаты, афоризмы, стихи, книжные подборки, обсуждения и многое другое

Книжный Вестник
БИБЛИОТЕКИ

Библиотека с любовными романами, которая наверняка придётся по вкусу женской части аудитории

Любовные романы

Библиотека с фантастикой и фэнтези, а также смежных жанров

Фантастика

Самые популярные книги в формате фб2

Топ фб2 книги