Программирование — страница 44 из 57

Числа

“Любая сложная проблема имеет ясное, простое

и при этом неправильное решение”.

Г.Л. Менкен (H.L. Mencken)

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

24.1. Введение

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

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

• Массивы, как в стиле языка С, так и класс из библиотека Matrix, который лучше подходит для числовых расчетов.

• Введение в случайные числа.

• Стандартные математические функции из библиотеки.

• Комплексные числа.

Основное внимание уделено многомерным массивам в стиле языка С и библиотеке N-мерных матриц Matrix, которая позволяет упростить работу с матрицами (многомерными массивами).

24.2. Размер, точность и переполнение

 Когда вы используете встроенные типы и обычные методы вычислений, числа хранятся в областях памяти фиксированного размера; иначе говоря, целочисленные типы (int, long и др.) представляют собой лишь приближение целых чисел, а числа с плавающей точкой (float, double и др.) являются лишь приближением действительных чисел. Отсюда следует, что с математической точки зрения некоторые вычисления являются неточными или неправильными. Рассмотрим пример.

float x = 1.0/333;
float sum = 0;
for (int i=0; i<333; ++i) sum+=x;
cout << setprecision(15) << sum << "\n";

Выполнив эту программы, мы получим не единицу, а

0.999999463558197

Мы ожидали чего-то подобного. Число с плавающей точкой состоит только из фиксированного количества битов, поэтому мы всегда можем “испортить” его, выполнив вычисление, результат которого состоит из большего количества битов, чем допускает аппаратное обеспечение. Например, рациональное число 1/3 невозможно представить точно как десятичное число (однако можно использовать много цифр его десятичного разложения). Точно так же невозможно точно представить число 1/333, поэтому, когда мы складываем 333 копии числа x (наилучшее машинное приближение числа 1/333 с помощью типа float), то получим число, немного отличающееся от единицы. При интенсивном использовании чисел с плавающей точкой возникает ошибка округления; остается лишь оценить, насколько сильно она влияет на результат.

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

ПОПРОБУЙТЕ

Замените в примере число 333 числом 10 и снова выполните программу. Какой результат следовало ожидать? Какой результат вы получили? А ведь мы предупреждали!

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

Рассмотрим целочисленную задачу.

short int y = 40000;
int i = 1000000;
cout << y << " " << i*i << "\n";

Выполнив эту программу, получим следующий результат:

–25536 –727379968

Этого следовало ожидать. Здесь мы видим эффект переполнения. Целочисленные типы позволяют представить лишь относительно небольшие целые числа. Нам просто не хватит битов, чтобы точно представить каждое целое число, поэтому нам необходим способ, позволяющий выполнять эффективные вычисления. В данном случае двухбайтовое число типа short не может представить число 40 000, а четырехбайтовое число типа int не может представить число 1 000 000 000 000. Точные размеры встроенных типов в языке C++ (см. раздел A.8) зависят от аппаратного обеспечения и компилятора; размер переменной x или типа x в байтах можно определить с помощью оператора sizeof(x). По определению sizeof(char)==1. Это можно проиллюстрировать следующим образом.

 Эти размеры характерны для операционной системы Windows и компилятора компании Microsoft. В языке С++ есть много способов представить целые числа и числа с плавающей точкой, используя разные размеры, но при отсутствии важных причин лучше придерживаться типов char, int и double. В большинстве программ (но, разумеется, не во всех) остальные типы целых чисел и чисел с плавающей точкой вызывают больше проблем, чем хотелось бы.

Целое число можно присвоить переменной, имеющей тип числа с плавающей точкой. Если целое число окажется больше, чем может представить тип числа с плавающей точкой, произойдет потеря точности. Рассмотрим пример.

cout << "размеры: " << sizeof(int) << ' ' << sizeof(float) << '\n';
int x = 2100000009; // большое целое число
float f = x;
cout << x << ' ' << f << endl;
cout << setprecision(15) << x << ' ' << f << '\n';

На нашем компьютере мы получили следующий результат:

Sizes: 4 4
2100000009 2.1e+009
2100000009 2100000000

Типы float и int занимают одинаковое количество памяти (4 байта). Тип float состоит из мантиссы (как правило, числа от нуля до единицы) и показателя степени (т.е. мантисса*10 показатель степени), поэтому он не может точно выразить самое большое число int. (Если бы мы попытались сделать это, то не смогли бы выделить достаточно памяти для мантиссы после размещения в памяти показателя степени.) Как и следовало ожидать, переменная f представляет число 2100000009 настолько точно, насколько это возможно. Однако последняя цифра 9 вносит слишком большую ошибку, — именно поэтому мы выбрали это число для иллюстрации.

 С другой стороны, когда мы присваиваем число с плавающей точкой перемен- ной целочисленного типа, происходит усечение; иначе говоря, дробная часть — цифры после десятичной точки — просто отбрасываются. Рассмотрим пример.

float f = 2.8;
int x = f;
cout << x << ' ' << f << '\n';

Значение переменной x будет равно 2. Оно не будет равным 3, как вы могли подумать, если применили “правило округления 4/5”. В языке C++ преобразование типа float в тип int сопровождается усечением, а не округлением.

 При вычислениях следует опасаться возможного переполнения и усечения.

Язык C++ не решит эту проблему за вас. Рассмотрим пример.

void f(int i, double fpd)
{
 char c = i;      // да: тип char действительно представляет
                  // очень маленькие целые числа
 short s = i;     // опасно: переменная типа int может
                  // не поместиться
                  // в памяти, выделенной для переменной
                  // типа short
 i = i+1;         // что, если число i станет максимальным?
 long lg = i*i;   // опасно: переменная типа long не может
                  // вместить результат
 float fps = fpd; // опасно: большее число типа large может
                  // не поместиться в типе float
 i = fpd;         // усечение: например, 5.7 –> 5
 fps = i;         // можно потерять точность (при очень
                  // больших целых)
}

void g()
{
  char ch = 0;
  for (int i = 0; i<500; ++i)
    cout << int(ch++) << '\t';
}

Если сомневаетесь, поэкспериментируйте! Не следует отчаиваться и в то же время нельзя просто читать документацию. Без экспериментирования вы можете не понять содержание весьма сложной документации, связанной с числовыми типами.

ПОПРОБУЙТЕ

Выполните функцию g(). Модифицируйте функцию f() так, чтобы она выводила на печать переменные c, s, i и т.д. Протестируйте программу на разных значениях.

 Представление целых чисел и их преобразование еще будет рассматриваться в разделе 25.5.3. По возможности ограничивайтесь немногими типами данных, чтобы минимизировать вероятность ошибок. Например, используя только тип double и избегая типа float, мы минимизируем вероятность возникновения проблем, связанных с преобразованием double — float. Например, мы предпочитаем использовать только типы int, double и complex (см. раздел 24.9) для вычислений, char — для символов и bool — для логических сущностей. Остальные арифметические типы мы используем только при крайней необходимости.

24.2.1. Пределы числовых диапазонов

 Каждая реализация языка C++ определяет свойства встроенных типов в заголовках , и , чтобы программисты могли проверить пределы диапазонов, установить сигнальные метки и т.д. Эти значения перечислены в разделе Б.9.1. Они играют очень важную роль для создания низкоуровневых инструментов. Если они вам нужны, значит, вы работаете непосредственно с аппаратным обеспечением, хотя существуют и другие приложения. Например, довольно часто возникают вопросы о тонкостях реализации языка, например: “Насколько большим является тип int?” или “Имеет ли знак тип char?” Найти определенные и правильные ответы в системной документации бывает трудно, а в стандарте указаны только минимальные требования. Однако можно легко написать программу, находящую ответы на эти вопросы.

cout << "количество байтов в типе int: " << sizeof(int) << '\n';
cout << "наибольшее число типа int: " << INT_MAX << endl;
cout << "наименьшее число типа int: " << numeric_limits::min()
<< '\n';

if (numeric_limits::is_signed)
  cout << "тип char имеет знак n";
else
  cout << "тип char не имеет знака\n";

cout << "char с минимальным значением: "
<< numeric_limits::min() <<'\n';
cout << "минимальное значение типа char: "
<< int(numeric_limits::min()) << '\n';

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

Эти пределы также могут быть полезными для выявления переполнения.

24.3. Массивы

Массив (array) — это последовательность, в которой доступ к каждому элементу осуществляется с помощью его индекса (позиции). Синонимом этого понятия является вектор (vector). В этом разделе мы уделим внимание многомерным массивам, элементами которых являются тоже массивы. Обычно многомерный массив называют матрицей (matrix). Разнообразие синонимов свидетельствует о популярности и полезности этого общего понятия. Стандартные классы vector (см. раздел Б.4), array (см. раздел 20.9), а также встроенный массив (см. раздел A.8.2) являются одномерными. А что если нам нужен двумерный массив (например, матрица)? А если нам нужны семь измерений? Проиллюстрировать одно- и двухмерные массивы можно так.

Массивы имеют фундаментальное значение в большинстве вычислений, связанных с так называемым “перемалыванием чисел” (“number crunching”). Наиболее интересные научные, технические, статистические и финансовые вычисления тесно связаны с массивами.

 Часто говорят, что массив состоит из строки столбцов.

Столбец — это последовательность элементов, имеющих одинаковые первые координаты (х-координаты). Строка — это множество элементов, имеющих одинаковые вторые координаты (y-координаты).

24.4. Многомерные массивы в стиле языка С

В качестве многомерного массива можно использовать встроенный массив в языке С++ . В этом случае многомерный массив интерпретируется как массив массивов, т.е. массив, элементами которого являются массивы. Рассмотрим пример.

int ai[4];        // 1-мерный массив
double ad[3][4];  // 2-мерный массив
char ac[3][4][5]; // 3-мерный массив
ai[1] = 7;
ad[2][3] = 7.2;
ac[2][3][4] = 'c';

 Этот подход наследует все преимущества и недостатки одномерного массива.

• Преимущества

• Непосредственное отображение с помощью аппаратного обеспечения.

• Эффективные низкоуровневые операции.

• Непосредственная языковая поддержка.

• Проблемы

• Многомерные массивы в стиле языка являются массивами массивов(см. ниже).

• Фиксированные размеры (например, фиксированные на этапе компиляции). Если хотите определять размер массива на этапе выполнения программы, то должны использовать свободную память.

• Массивы невозможно передать аккуратно. Массив превращается в указатель на свой первый элемент при малейшей возможности.

• Нет проверки диапазона. Как обычно, массив не знает своего размера.

• Нет операций над массивами, даже присваивания (копирования).

Встроенные массивы широко используются в числовых расчетах. Они также являются основным источником ошибок и сложностей. Создание и отладка таких программ у большинства людей вызывают головную боль. Если вы вынуждены использовать встроенные массивы, почитайте учебники (например, The C++ Programming Language, Appendix C, p. 836–840). К сожалению, язык C++ унаследовал многомерные массивы от языка C, поэтому они до сих пор используются во многих программах.

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

void f1(int a[3][5]);     // имеет смысл только в матрице [3][5]

void f2(int [ ][5], int dim1);  // первая размерность может быть
                                // переменной

void f3(int [5 ][ ], int dim2); // ошибка: вторая размерность
                                // не может быть переменной

void f4(int[ ][ ], int dim1, int dim2); // ошибка (совсем
                                        // не работает)

void f5(int* m, int dim1, int dim2) // странно, но работает
{
  for (int i=0; i
  for (int j = 0; j
}

Здесь мы передаем массив m как указатель int*, даже если он является двумерным. Поскольку вторая переменная должна быть переменной (параметром), у нас нет никакой возможности сообщить компилятору, что массив m является массивом (dim1, dim2), поэтому мы просто передаем указатель на первую его ячейку. Выражение m[i*dim2+j] на самом деле означает m[i,j], но, поскольку компилятор не знает, что переменная m — это двумерный массив, мы должны сначала вычислить позицию элемента m[i,j] в памяти.

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

24.5. Библиотека Matrix

 Каково основное предназначение массива (матрицы) в численных расчетах?

• “Мой код должен выглядеть очень похожим на описание массивов, изложенное в большинстве учебников по математике”.

• Это относится также к векторам, матрицам и тензорам.

• Проверка на этапах компиляции и выполнения программы.

• Массивы любой размерности.

• Массивы с произвольным количеством элементов в любой размерности.

• Массивы являются полноценными переменными/объектами.

• Их можно передавать куда угодно.

• Обычные операции над массивами.

• Индексирование: ().

• Срезка: [].

• Присваивание: =.

• Операции пересчета (+=, –=, *=, %= и т.д.).

• Встроенные векторные операции (например, res[i] = a[i]*c+b[2]).

• Скалярное произведение (res = сумма элементов a[i]*b[i]; известна также как inner_product).

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

• Массивы при необходимости можно увеличивать (при их реализации не используются “магические” числа).

Библиотека Matrix делает это и только это. Если вы хотите большего, то должны самостоятельно написать сложные функции обработки массивов, разреженных массивов, управления распределением памяти и так далее или использовать другую библиотеку, которая лучше соответствует вашим потребностям. Однако многие эти потребности можно удовлетворить с помощью алгоритмов и структур данных, надстроенных над библиотекой Matrix. Библиотека Matrix не является частью стандарта ISO C++. Вы можете найти ее описание на сайте в заголовке Matrix.h. Свои возможности она определяет в пространстве имен Numeric_lib. Мы выбрали слово Matrix, потому что слова “вектор” и “массив” перегружены в библиотеках языка C++. Реализация библиотеки Matrix основана на сложных методах, которые здесь не описываются. 

24.5.1. Размерности и доступ

Рассмотрим простой пример.

#include "Matrix.h"
using namespace Numeric_lib;

void f(int n1, int n2, int n3)
{
 Matrix ad1(n1); // элементы типа double;
                           // одна размерность
 Matrix ai1(n1);    // элементы типа int;
                           // одна размерность
 ad1(7) = 0; // индексирование ( ) в стиле языка Fortran
 ad1[7] = 8; // индексирование [ ] в стиле языка C

 Matrix ad2(n1,n2);    // двумерный
 Matrix ad3(n1,n2,n3); // трехмерный
 ad2(3,4) = 7.5;                 // истинное многомерное
                                 // индексирование
 ad3(3,4,5) = 9.2;
}

 Итак, определяя переменную типа Matrix (объект класса Matrix), вы должны указать тип элемента и количество размерностей. Очевидно, что класс Matrix является шаблонным, а тип элементов и количество размерностей представляют собой шаблонные параметры. В результате, передав пару шаблонных параметров классу Matrix (например, Matrix), получаем тип (класс), с помощью которого можно определить объекты, указав аргументы (например, Matrixad2(n1,n2)); эти аргументы задают размерности. Итак, переменная ad2 является двумерным массивом с размерностями n1 и n2, которую также называют матрицей n1 на n2. Для того чтобы получить элемент объявленного типа из одномерного объекта класса Matrix, следует указать один индекс. Для того чтобы получить элемент объявленного типа из двумерного объекта класса Matrix, следует указать два индекса.

Как и во встроенных массивах и объектах класса vector, элементы в объекте класса Matrix индексируются с нуля (а не с единицы, как в языке Fortran); иначе говоря, элементы объекта класса Matrix нумеруются в диапазоне [0,max], где max — количество элементов.

 Это просто и взято прямо из учебника. Если у вас возникнут проблемы, нужно лишь обратиться к нужному учебнику по математике, а не к руководству по программированию. Единственная тонкость здесь заключается в том, что мы не указали количество размерностей в объекте класса Matrix: по умолчанию он является одномерным. Обратите внимание также на то, что мы можем использовать как индексирование с помощью оператора [] (в стиле языков C и C++), так и с помощью оператора () (в стиле языка Fortran).

Это позволяет нам лучше справляться с большим количеством размерностей. Индекс [x] всегда означает отдельный индекс, выделяя отдельную строку в объекте класса Matrix; если переменная a является n мерным объектом класса Matrix, то a[x] — это (n–1)-размерный объект класса Matrix. Обозначение (x,y,z) подразумевает использование нескольких индексов, выделяя соответствующий элемент объекта класса Matrix; количество индексов должно равняться количеству размерностей.

Посмотрим, что произойдет, если мы сделаем ошибку.

void f(int n1,int n2,int n3)
{
 Matrix ai0;   // ошибка: 0-размерных матриц не бывает
 Matrix ad1(5);
 Matrix ai(5);

 Matrix ad11(7);
 ad1(7) = 0;               // исключение Matrix_error
                           // (7 — за пределами диапазона)
 ad1 = ai;                 // ошибка: разные типы элементов
 ad1 = ad11;               // исключение Matrix_error
                           // (разные размерности)
Matrix ad2(n1);  // ошибка: пропущена длина 2-й
                           // размерности
 ad2(3) = 7.5;             // ошибка: неправильное количество
                           // индексов
 ad2(1,2,3) = 7.5;         // ошибка: неправильное количество
                           // индексов

 Matrix ad3(n1,n2,n3);
 Matrix ad33(n1,n2,n3);
 ad3 = ad33;               // OK: одинаковые типы элементов,
                           // одинаковые размерности
}

Несоответствия между объявленным количеством размерностей и их использованием обнаруживается на этапе компиляции. Выход за пределы диапазона перехватывается на этапе выполнения программы; при этом генерируется исключение Matrix_error.

 Первая размерность матрицы — это строка, а вторая — столбец, поэтому индекс — это двумерная матрица (двумерный массив), имеющая вид (строка,столбец). Можно также использовать обозначение [строка][столбец], так как индексирование двумерной матрицы с помощью одномерного индекса порождает одномерную матрицу — строку. Эту ситуацию можно проиллюстрировать следующим образом.

Этот объект класса Matrix размещается в памяти построчно.

Класс Matrix знает свою размерность, поэтому его элементы можно очень просто передавать как аргумент,

void init(Matrix& a) // инициализация каждого элемента
                            // характеристическим значением
{
  for (int i=0; i
    for (int j = 0; j
      a(i,j) = 10*i+j;
}

void print(const Matrix& a) // вывод элементов построчно
{
  for (int i=0; i
    for (int j = 0; j
      cout << a(i,j) <<'\t';
    cout << '\n';
  }
}

 Итак, dim1() — это количество элементов в первой размерности, dim2() — количество элементов во второй размерности и т.д. Тип элементов и количество размерностей являются частью класса Matrix, поэтому невозможно написать функцию, получающую объект класса Matrix как аргумент (но можно написать шаблон).

void init(Matrix& a); // ошибка: пропущены тип элементов
                      // и количество размерностей

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

24.5.2. Одномерный объект класса Matrix

Что можно сделать с простейшим объектом класса Matrix — одномерной матрицей?

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

Matrix a1(8); // a1 — это одномерная матрица целых чисел
Matrix a(8);    // т.е. Matrix a(8);

Таким образом, объекты a и a1 имеют одинаковый тип (Matrix). У каждого объекта класса Matrix можно запросить общее количество элементов и количество элементов в определенном измерении. У одномерного объекта класса Matrix эти параметры совпадают.

a.size(); // количество элементов в объекте класса Matrix
a.dim1(); // количество элементов в первом измерении

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

int* p = a.data(); // извлекаем данные с помощью указателя на массив

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

a(i);   // i-й элемент (в стиле языка Fortran) с проверкой
        // диапазона
a[i];   // i-й элемент (в стиле языка C) с проверкой диапазона
a(1,2); // ошибка: a — одномерный объект класса Matrix

Многие алгоритмы обращаются к части объекта класса Matrix. Эта часть называется срезкой и создается функцией slice() (часть объекта класса Matrix или диапазон элементов). В классе Matrix есть два варианта этой функции.

a.slice(i); // элементы, начиная с a[i] и заканчивая последним
a.slice(i,n); // n элементов, начиная с a[i] и заканчивая a[i+n–1]

Индексы и срезки можно использовать как в левой части оператора присваивания, так и в правой. Они ссылаются на элементы объекта класса Matrix, не создавая их копии. Рассмотрим пример.

a.slice(4,4) = a.slice(0,4); // присваиваем первую половину матрицы
                             // второй

Например, если объект a вначале выглядел так:

{ 1 2 3 4 5 6 7 8 }

то получим

{ 1 2 3 4 1 2 3 4 }

Обратите внимание на то, что чаще всего срезки задаются начальными и последними элементами объекта класса Matrix; т.е. a.slice(0,j) — это диапазон [0:j], а a.slice(j) — диапазон [j:a.size()]. В частности, приведенный выше пример можно легко переписать:

a.slice(4) = a.slice(0,4); // присваиваем первую половину матрицы
                           // второй

Иначе говоря, обозначения — дело вкуса. Вы можете указать такие индексы i и n, так что a.slice(i,n) выйдет за пределы диапазона матрицы a. Однако полученная срезка будет содержать только те элементы, которые действительно принадлежат объекту a. Например, срезка a.slice(i,a.size()) означает диапазон [i:a.size()], а a.slice(a.size()) и a.slice(a.size(),2) — это пустые объекты класса Matrix. Это оказывается полезным во многих алгоритмах. Мы подсмотрели это обозначение в математических текстах. Очевидно, что срезка a.slice(i,0) является пустым объектом класса Matrix. Нам не следовало бы писать это намеренно, но существуют алгоритмы, которые становятся проще, если срезка a.slice(i,n) при параметре n, равном 0, является пустой матрицей (это позволяет избежать ошибки).

 Копирование всех элементов выполняется как обычно.

Matrix a2 = a;  // копирующая инициализация
a = a2;              // копирующее присваивание

 К каждому элементу объекта класса Matrix можно применять встроенные операции.

a *= 7;   // пересчет: a[i]*=7 для каждого i (кроме того, +=, –=, /=
          // и т.д.)
a = 7;    // a[i]=7 для каждого i

Это относится к каждому оператору присваивания и каждому составному оператору присваивания (=, +=, –=, /=, *=, %=, ^=, &=, |=, >>=, <<=) при условии, что тип элемента поддерживает соответствующий оператор. Кроме того, к каждому элементу объекта класса Matrix можно применять функции.

a.apply(f);    // a[i]=f(a[i]) для каждого элемента a[i]
a.apply(f,7);  // a[i]=f(a[i],7) для каждого элемента a[i]

Составные операторы присваивания и функция apply() модифицируют свои аргументы типа Matrix. Если же мы захотим создать новый объект класса Matrix, то можем выполнить следующую инструкцию:

b = apply(abs,a); // создаем новый объект класса Matrix
                  // с условием b(i)==abs(a(i))

Функция abs — это стандартная функция вычисления абсолютной величины (раздел 24.8). По существу, функция apply(f,x) связана с функцией x.apply(f) точно так же, как оператор + связан с оператором +=. Рассмотрим пример.

b = a*7;        // b[i] = a[i]*7 для каждого i
a *= 7;         // a[i] = a[i]*7 для каждого i
y = apply(f,x); // y[i] = f(x[i]) для каждого i
x.apply(f);     // x[i] = f(x[i]) для каждого i

В результате a==b и x==y.

 В языке Fortran второй вариант функции apply называется функцией пересылки (“broadcast” function). В этом языке чаще пишут вызов f(x), а не apply(f,x). Для того чтобы эта возможность стала доступной для каждой функции f (а не только для отдельных функций, как в языке Fortran), мы должны присвоить операции пересылки конкретное имя, поэтому (повторно) использовали имя apply.

Кроме того, для того чтобы обеспечить соответствие с вариантом функции-члена apply, имеющим вид a.apply(f,x), мы пишем

b = apply(f,a,x); // b[i]=f(a[i],x) для каждого i

Рассмотрим пример.

double scale(double d, double s) { return d*s; }
b = apply(scale,a,7); // b[i] = a[i]*7 для каждого i

Обратите внимание на то, что “автономная” функция apply() принимает в качестве аргумента функцию, вычисляющую результат по ее аргументам, а затем использует этот результат для инициализации итогового объекта класса Matrix. Как правило, это не приводит к изменению объекта класса Matrix, к которому эта функция применяется. В то же время функция-член apply() отличается тем, что принимает в качестве аргумента функцию, модифицирующую ее аргументы; иначе говоря, она модифицирует элементы объекта класса Matrix, к которому применяется. Рассмотрим пример.

void scale_in_place(double& d, double s) { d *= s; }
b.apply(scale_in_place,7); // b[i] *= 7 для каждого i

В классе Matrix предусмотрено также много полезных функций из традиционных математических библиотек.

Matrix a3 = scale_and_add(a,8,a2); // объединенное умножение
                                        // и сложение
int r = dot_product(a3,a);              // скалярное произведение

 Операцию scale_and_add() часто называют объединенным умножением и сложением (fused multiply-add), или просто fma; ее определение выглядит так: result(i)=arg1(i)*arg2+arg3(i) для каждого i в объекте класса Matrix. Скалярное произведение также известно под именем inner_product и описано в разделе 21.5.3; ее определение выглядит так: result+=arg1(i)*arg2(i) для каждого i в объекте класса Matrix, где накопление объекта result начинается с нуля.

Одномерные массивы очень широко распространены; их можно представить как в виде встроенного массива, так и с помощью классов vector и Matrix. Класс Matrix следует применять тогда, когда необходимо выполнять матричные операции, такие как *=, или когда объект класса Matrix должен взаимодействовать с другими объектами этого класса, имеющими более высокую размерность.

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

Класс Matrix имеет два конструктора для копирования данных из встроенных массивов в объект класса Matrix. Рассмотрим пример.

void some_function(double* p, int n)
{
  double val[] = { 1.2, 2.3, 3.4, 4.5 };
  Matrix data(p,n);
  Matrix constants(val);
  // ...
}

Это часто бывает полезным, когда мы получаем данные в виде обычных массивов или векторов, созданных в других частях программы, не использующих объекты класса Matrix.

Обратите внимание на то, что компилятор может самостоятельно определить количество элементов в инициализированном массиве, поэтому это число при определении объекта constants указывать не обязательно — оно равно — 4. С другой стороны, если элементы заданы всего лишь указателем, то компилятор не знает их количества, поэтому при определении объекта data мы должны задать как указатель p, так и количество элементов n.

24.5.3. Двумерный объект класса Matrix

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

Matrix a(3,4);
int s = a.size();  // количество элементов
int d1 = a.dim1(); // количество элементов в строке
int d2 = a.dim2(); // количество элементов в столбце
int* p = a.data(); // извлекаем данные с помощью указателя в стиле
                   // языка С

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

Мы можем использовать индексы.

a(i,j);  // (i,j)-й элемент (в стиле языка Fortran) с проверкой
         // диапазона
a[i];    // i-я строка (в стиле языка C) с проверкой диапазона
a[i][j]; // (i,j)-й элемент (в стиле языка C)

 В двумерном объекте класса Matrix индексирование с помощью конструкции [i] создает одномерный объект класса Matrix, представляющий собой i-ю строку. Это значит, что мы можем извлекать строки и передавать их операторам и функциям, получающим в качестве аргументов одномерные объекты класса Matrix и даже встроенные массивы (a[i].data()). Обратите внимание на то, что индексирование вида a(i,j) может оказаться быстрее, чем индексирование вида a[i][j], хотя это сильно зависит от компилятора и оптимизатора.

Мы можем получить срезки.

a.slice(i);   // строки от a[i] до последней
a.slice(i,n); // строки от a[i] до a[i+n–1]

Срезка двумерного объекта класса Matrix сама является двумерным объектом этого класса (возможно, с меньшим количеством строк). Распределенные операции над двумерными матрицами такие же, как и над одномерными. Этим операциям неважно, как именно хранятся элементы; они просто применяются ко всем элементам в порядке их следования в памяти.

Matrix a2 = a; // копирующая инициализация
a = a2;          // копирующее присваивание
a *= 7;          // пересчет (и +=, –=, /= и т.д.)
a.apply(f);      // a(i,j)=f(a(i,j)) для каждого элемента a(i,j)
a.apply(f,7);    // a(i,j)=f(a(i,j),7) для каждого элемента a(i,j)
b=apply(f,a);    // создаем новую матрицу с b(i,j)==f(a(i,j))
b=apply(f,a,7);  // создаем новую матрицу с b(i,j)==f(a(i,j),7)

Оказывается, что перестановка строк также полезна, поэтому мы предусмотрим и ее.

a.swap_rows(1,2); // перестановка строк a[1] <–> a[2]

 Перестановки столбцов swap_columns() не существует. Если она вам потребуется, то вы сможете написать ее самостоятельно (см. упр. 11). Из-за построчной схемы хранения матриц в памяти строки и столбцы не совсем равноправны. Эта асимметрия проявляется также в том, что оператор [i] возвращает только строку (а для столбцов аналогичный оператор не предусмотрен). Итак, в тройке (i,j) первый индекс i выбирает строку. Эта асимметрия имеет глубокие математические корни.

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

enum Piece { none, pawn, knight, queen, king, bishop, rook };
Matrix board(8,8); // шахматная доска

const int white_start_row = 0;
const int black_start_row = 7;

Piece init_pos[] = {rook,knight,bishop, queen,king,bishop,knight,rook};
Matrix start_row(init_pos); // инициализация элементов из
                                   // init_pos
Matrix clear_row(8);        // 8 элементов со значениями
                                   // по умолчанию

Инициализация объекта clear_row использует возможность задать условие none==0 и то, что эти элементы по умолчанию инициализируются нулем. Мы могли бы предпочесть другой код.

Matrix start_row
 = {rook,knight,bishop,queen,king,bishop,knight,rook};

Однако он не работает (по крайней мере, пока не появится новая версия языка C++ (C++0x)), поэтому пока приходится прибегать к трюкам с инициализацией массива (в данном случае init_pos) и использовать его для инициализации объектов класса Matrix. Мы можем использовать объекты start_row и clear_row следующим образом:

board[white_start_row] = start_row;             // расставить белые фигуры
for (int i = 1; i<7; ++i) board[i] = clear_row; // очистить середину
                                                // доски
board[black_start_row] = start_row;             // расставить черные фигуры

Обратите внимание на то, что когда мы извлекли строку, используя выражение [i], мы получили значение lvalue (см. раздел 4.3); иначе говоря, мы можем присвоить результат элементу board[i].

24.5.4. Ввод-вывод объектов класса Matrix

Библиотека Matrix предоставляет очень простые средства для ввода и вывода одно- и двухмерных объектов класса Matrix:

Matrix a(4);
cin >> a;
cout << a;

Этот фрагмент кода прочитает четыре разделенные пробелами числа типа double, заключенных в фигурные скобки; например:

{ 1.2 3.4 5.6 7.8 }

Вывод очень прост, поэтому мы просто можем увидеть то, что ввели. Механизм ввода-вывода двумерных объектов класса Matrix просто считывает и записывает последовательности одномерных объектов класса Matrix, заключенные в квадратные скобки. Рассмотрим пример.

Matrix m(2,2);
cin >> m;
cout << m;

Он прочитает запись

{
  { 1 2 }
  { 3 4 }
}

Вывод очень похож.

Операторы <<и >> из класса Matrix позволяют писать простые программы. В более сложных ситуациях нам потребуется заменить их своими операторами. По этой причине определение операторов << и >> из класса Matrix помещены в заголовок MatrixIO.h (а не Matrix.h), поэтому, для того чтобы использовать матрицы в своей программе, вам не обязательно включать заголовок MatrixIO.h.

24.5.5. Трехмерный объект класса Matrix

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

Matrix a(10,20,30);
a.size();             // количество элементов
a.dim1();             // количество элементов в размерности 1
a.dim2();             // количество элементов в размерности 2
a.dim3();             // количество элементов в размерности 3
int* p = a.data();    // извлекает данные по указателю (в стиле языка С)
a(i,j,k);             // (i,j,k)-й элемент (в стиле языка Fortran)
                      // с проверкой диапазона
a[i];                 // i-я строка (в стиле языка C)
                      // с проверкой диапазона
a[i][j][k];           // (i,j,k)-й элемент (в стиле языка С)
a.slice(i);           // строки от i-й до последней
a.slice(i,j);         // строки от i-й до j-й
Matrix a2 = a; // копирующая инициализация
a = a2;               // копирующее присваивание
a *= 7;               // пересчет (и +=, –=, /= и т.д.)
a.apply(f);           // a(i,j,k)=f(a(i,j,k)) для каждого элемента a(i,j,k)
a.apply(f,7);         // a(i,j,k)=f(a(i,j,k),7) для каждого элемента a(i,j,k)
b=apply(f,a);         // создает новую матрицу с условием b(i,j,k)==f(a(i,j,k))
b=apply(f,a,7);       // создает новую матрицу с условием b(i,j,k)==f(a(i,j,k),7)
a.swap_rows(7,9);     // переставляет строки a[7] <–> a[9]

Если вы умеете работать с двумерными объектами класса Matrix, то сможете работать и с трехмерными. Например, здесь a — трехмерная матрица, поэтому a[i] — двумерная (при условии, что индекс i не выходит за пределы допустимого диапазона); a[i][j] — одномерная (при условии, что индекс j не выходит за пределы допустимого диапазона); a[i][j][k] — элемент типа int (при условии, что индекс k не выходит за пределы допустимого диапазона).

Поскольку мы видим трехмерный мир, при моделировании чаще используются трехмерные объекты класса Matrix (например, в физическом моделировании в декартовой системе координат).

int grid_nx; // разрешение сетки; задается вначале
int grid_ny;
int grid_nz;
Matrix cube(grid_nx, grid_ny, grid_nz);

Если добавить время в качестве четвертого измерения, то получим четырехмерное пространство, в котором необходимы четырехмерные объекты класса Matrix. И так далее. 

24.6. Пример: решение систем линейных уравнений

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

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

a_1,1x₁ + ... + a_1,nx_n = b₁

...

a_n,1x

1 + ... + a_n,nx_n = b_n

где буквы x обозначают n неизвестных, а буквы a и b — константы. Для простоты предполагаем, что неизвестные и константы являются числами с плавающей точкой.

Наша цель — найти неизвестные, которые одновременно удовлетворяют указанные n уравнений. Эти уравнения можно компактно выразить с помощью матрицы и двух векторов.

Ax = b

где A — квадратная матрица n×n коэффициентов:

Векторы x и b векторы неизвестных и константа соответственно.

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

Алгоритм несложен. Для того чтобы элемент в позиции (i,j) стал равным нулю, необходимо умножить строку i на константу, чтобы элемент в позиции (i,j) стал равным другому элементу в столбце j, например a(k, j). После этого просто вычтем одно уравнение из другого и получим a(i,j)==0. При этом все остальные значения в строке i изменятся соответственно.

Если все диагональные элементы окажутся ненулевыми, то система имеет единственное решение, которое можно найти в ходе обратной подстановки. Сначала решим последнее уравнение (это просто).

a_n,nx_n = b_n

Очевидно, что x[n] равен b[n]/a(n,n). Теперь исключим строку n из системы, найдем значение x[n–1] и будем продолжать процесс, пока не вычислим значение x[1].

При каждом значении n выполняем деление на a(n,n), поэтому диагональные значения должны быть ненулевыми. Если это условие не выполняется, то обратная подстановка завершится неудачей. Это значит, что система либо не имеет решения, либо имеет бесконечно много решений. 

24.6.1. Классическое исключение Гаусса

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

typedef Numeric_lib::Matrix Matrix;
typedef Numeric_lib::Matrix Vector;

Затем выразим сам алгоритм.

Vector classical_gaussian_elimination(Matrix A,Vector b)
{
  classical_elimination(A, b);
  return back_substitution(A, b);
}

Иначе говоря, мы создаем копии входных матрицы A и вектора b (используя механизм передачи аргументов по значению), вызываем функцию для решения системы, а затем вычисляем результат с помощью обратной подстановки. Такое разделение задачи на части и система обозначений приняты во всех учебниках. Для того чтобы завершить программу, мы должны реализовать функции classical_elimination() и back_substitution(). Решение также можно найти в учебнике.

void classical_elimination(Matrix& A,Vector& b)
{
  const Index n = A.dim1();
  // проходим от первого столбца до последнего,
  // обнуляя элементы, стоящие ниже диагонали:
  for (Index j = 0; j
    const double pivot = A(j, j);
    if (pivot == 0) throw Elim_failure(j);

    // обнуляем элементы, стоящие ниже диагонали в строке i
    for (Index i = j+1; i
      const double mult = A(i, j) / pivot;
      A[i].slice(j) = scale_and_add(A[j].slice(j),
      –mult, A[i].slice(j));
      b(i) –= mult * b(j); // изменяем вектор b
    }
  }
}

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

Vector back_substitution(const Matrix& A, const Vector& b)
{
  const Index n = A.dim1();
  Vector x(n);
  for (Index i = n – 1; i >= 0; ––i) {
    double s = b(i)–dot_product(A[i].slice(i+1),x.slice(i+1));
    if (double m = A(i, i))
      x(i) = s / m;
    else
      throw Back_subst_failure(i);
  }
  return x;
} 

24.6.2. Выбор ведущего элемента

Для того чтобы избежать проблем с нулевыми диагональными элементами и повысить устойчивость алгоритма, можно переставить строки так, чтобы нули и малые величины на диагонали не стояли. Говоря “повысить устойчивость”, мы имеем в виду понижение чувствительности к ошибкам округления. Однако по мере выполнения алгоритма элементы матрицы будут изменяться, поэтому перестановку строк приходится делать постоянно (иначе говоря, мы не можем лишь один раз переупорядочить матрицу, а затем применить классический алгоритм).

void elim_with_partial_pivot(Matrix& A, Vector& b)
{
  const Index n = A.dim1();
  for (Index j = 0; j < n; ++j) {
    Index pivot_row = j;

    // ищем подходящий опорный элемент:
    for (Index k = j + 1; k < n; ++k)
      if (abs(A(k, j)) > abs(A(pivot_row, j))) pivot_row = k;

    // переставляем строки, если найдется лучший опорный
    // элемент
    if (pivot_row != j) {
      A.swap_rows(j, pivot_row);
      std::swap(b(j), b(pivot_row));
    }

    // исключение:
    for (Index i = j + 1; i < n; ++i) {
      const double pivot = A(j, j);
      if (pivot==0) error("Решения нет: pivot==0");
        onst double mult = A(i, j)/pivot;
      A[i].slice(j) = scale_and_add(A[j].slice(j),
      –mult, A[i].slice(j));
      b(i) –= mult * b(j);
    }
  }
}

Для того чтобы не писать циклы явно и привести код в более традиционный вид, мы используем функции swap_rows() и scale_and_multiply().

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

Очевидно, что мы должны протестировать нашу программу. К счастью, это сделать несложно.

void solve_random_system(Index n)
{
  Matrix A = random_matrix(n); // см. раздел 24.7
  Vector b = random_vector(n);

  cout << "A = " << A << endl;
  cout << "b = " << b << endl;

  try {
    Vector x = classical_gaussian_elimination(A, b);
    cout << "Решение методом Гаусса x = " << x << endl;
    Vector v = A * x;
    cout << " A * x = " << v << endl;
  }
  catch(const exception& e) {
    cerr << e.what() << std::endl;
  }
}

Существуют три причины, из-за которых можно попасть в раздел catch.

• Ошибка в программе (однако, будучи оптимистами, будем считать, что этого никогда не произойдет).

• Входные данные, приводящие к краху алгоритма classical_elimination (целесообразно использовать функцию elim_with_partial_pivot).

• Ошибки округления.

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

Для того чтобы проверить наше решение, выводим на экране произведение A*x, которое должно быть равно вектору b (или достаточно близким к нему с учетом ошибок округления). Из-за вероятных ошибок округления мы не можем просто ограничиться инструкцией

if (A*x!=b) error("Неправильное решение");

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

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

Vector operator*(const Matrix& m,const Vector& u)
{
  const Index n = m.dim1();
  Vector v(n);
  for (Index i = 0; i < n; ++i) v(i) = dot_product(m[i], u);
  return v;
}

И вновь простая операция над объектом класса Matrix делает за нас большую часть работы. Как указывалось в разделе 24.5.3, операции вывода объектов класса Matrix описаны в заголовке MatrixIO.h. Функции random_matrix() и random_vector() просто используют случайные числа (раздел 24.7). Читатели могут написать эти функции в качестве упражнения. Имя Index является синонимом типа индекса, используемого в библиотеке Matrix, и определено с помощью оператора typedef (раздел A.15). Этот тип включается в программу с помощью объявления using.

using Numeric_lib::Index;

24.7. Случайные числа

 Если вы попросите любого человека назвать случайное число, то они назовут 7 или 17, потому что эти числа считаются самыми случайными. Люди практически никогда не называют число нуль, так как оно кажется таким идеально круглым числом, что не воспринимается как случайное, и поэтому его считают наименее случайным числом. С математической точки зрения это полная бессмыслица: ни одно отдельно взятое число нельзя назвать случайным. То, что мы часто называем случайными числами — это последовательность чисел, которые подчиняются определенному закону распределения и которые невозможно предсказать, зная предыдущие числа. Такие числа очень полезны при тестировании программ (они позволяют генерировать множество тестов), в играх (это один из способов гарантировать, что следующий шаг в игре не совпадет с предыдущим) и в моделировании (мы можем моделировать сущность, которая ведет себя случайно в пределах изменения своих параметров).

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

int rand(); // возвращает числа из диапазона
            // [0:RAND_MAX]
RAND_MAX    // наибольшее число, которое генерирует
            // датчик rand()
void srand(unsigned int); // начальное значение датчика
                          // случайных чисел

Повторяющиеся вызовы функции rand() генерируют последовательность чисел типа int, равномерно распределенных в диапазоне [0:RAND_MAX]. Эта последовательность чисел называется псевдослучайной, потому что она генерируется с помощью математической формулы и с определенного места начинает повторяться (т.е. становится предсказуемой и не случайной). В частности, если мы много раз вызовем функцию rand() в программе, то при каждом запуске программы получим одинаковые последовательности. Это чрезвычайно полезно для отладки. Если же мы хотим получать разные последовательности, то должны вызывать функцию srand() с разными значениями. При каждом новом аргументе функции srand() функция rand() будет порождать разные последовательности.

Например, рассмотрим функцию random_vector(), упомянутую в разделе 24.6.3. Вызов функции random_vector(n) порождает объект класса Matrix, содержащий n элементов, представляющих собой случайные числа в диапазоне от [0:n]:

Vector random_vector(Index n)
{
  Vector v(n);
  for (Index i = 0; i < n; ++i)
    v(i) = 1.0 * n * rand() / RAND_MAX;
  return v;
}

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

Сложнее получить целое число из заданного диапазона, например [0:max]. Большинство людей сразу предлагают следующее решение:

int val = rand()%max;

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

int randint(int max) { return rand()%max; }
int randint(int min, int max) { return randint(max–min)+min; }

 Таким образом, мы можем скрыть определение функции randint(), если окажется, что реализация функции rand() является неудовлетворительной. В промышленных программных системах, а также в приложениях, где требуются неравномерные распределения, обычно используются качественные и широко доступные библиотеки случайных чисел, например Boost::random. Для того чтобы получить представление о качестве вашего датчика случайных чисел, выполните упр. 10.

24.8. Стандартные математические функции

В стандартной библиотеке есть стандартные математические функции (cos, sin, log и т.д.). Их объявления можно найти в заголовке .

Стандартные математические функции могут иметь аргументы типов float, double, long double и complex (раздел 24.9). Эти функции очень полезны при вычислениях с плавающей точкой. Более подробная информация содержится в широко доступной документации, а для начала можно обратиться к документации, размещенной в веб.

 Если стандартная математическая функция не может дать корректного результата, она устанавливает флажок errno. Рассмотрим пример.

errno = 0;
double s2 = sqrt(–1);
if (errno) cerr << "Что-то где-то пошло не так, как надо";
if (errno == EDOM) // ошибка из-за выхода аргумента
                   // за пределы области определения
  cerr << " фунция sqrt() для отрицательных аргументов 
не определена.";
pow(very_large,2); // плохая идея
if (errno==ERANGE) // ошибка из-за выхода за пределы допустимого
                   // диапазона
  cerr << "pow(" << very_large
<< ",2) слишком большое число для double";

Если вы выполняете серьезные математические вычисления, то всегда должны проверять значение errno, чтобы убедиться, что после возвращения результата оно по-прежнему равно 0. Если нет, то что-то пошло не так, как надо. Для того чтобы узнать, какие математические функции могут устанавливать флажок errno и чему он может быть равен, обратитесь к документации.

 Как показано в примере, ненулевое значение флажка errno просто означает, что что-то пошло не так. Функции, не входящие в стандартную библиотеку, довольно часто также устанавливают флажок errno при выявлении ошибок, поэтому следует точнее анализировать разные значения переменной errno, чтобы понять, что именно случилось. В данном примере до вызова стандартной библиотечной функции переменная errno была равна нулю, а проверка значения errno сразу после выполнения функции может обнаружить, например, константы EDOM и ERANGE. Константа EDOM означает ошибку, возникшую из-за выхода аргумента за пределы области определения функции (domain error). Константа ERANGE означает выход за пределы допустимого диапазона значений (range error).

Обработка ошибок с помощью переменной errno носит довольно примитивный характер. Она уходит корнями в первую версию (выпуска 1975 года) математических функций языка C. 

24.9. Комплексные числа

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

template class complex {
  // комплексное число — это пара скалярных величин,
  // по существу, пара координат
  Scalar re, im;
public:
  complex(const Scalar & r, const Scalar & i) :re(r), im(i) { }
  complex(const Scalar & r) :re(r),im(Scalar ()) { }
  complex() :re(Scalar ()), im(Scalar ()) { }
  Scalar real() { return re; } // действительная часть
  Scalar imag() { return im; } // мнимая часть
  // операторы : = += –= *= /=
};

Стандартная библиотека complex поддерживает типы скалярных величин float, double и long double. Кроме членов класса complex и стандартных математических функций (раздел 24.8), заголовок содержит множество полезных операций.

Примечание: в классе complex нет операций < и %.

Класс complex используется так же, как любой другой встроенный тип, например double. Рассмотрим пример.

typedef complex dcmplx; // иногда выражение complex
                                // является слишком громоздким
void f(dcmplx z, vector& vc)
{
  dcmplx z2 = pow(z,2);
  dcmplx z3 = z2*9.3+vc[3];
  dcmplx sum = accumulate(vc.begin(), vc.end(), dcmplx());
  // ...
}

Помните, что не все операции над числами типов int и double определены для класса complex. Рассмотрим пример.

if (z2

Обратите внимание на то, что представление (схема) комплексных чисел в стандартной библиотеке языка С++ сопоставима с соответствующими типами в языках C и Fortran.

24.10. Ссылки

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

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

Архив MacTutor History of Mathematics, размещенный на веб-странице http://www-gap.dcs.st-and.ac.uk/~history.

• Отличная веб-страница для всех, кто любит математику или просто хочет ее применять.

• Отличная веб-страница для всех, кто хочет увидеть гуманитарный аспект математики; например, кто из крупных математиков выиграл Олимпийские игры?

• Знаменитые математики: биографии, достижения.

• Курьезы.

• Знаменитые кривые.

• Известные задачи.

• Математические темы.

• Алгебра.

• Анализ.

• Теория чисел.

• Геометрия и топология.

• Математическая физика.

• Математическая астрономия.

• История математики.

• Многое другое

Freeman, T. L., and Chris Phillips. Parallel Numerical Algorithms. Prentice Hall, 1992.

Gullberg, Jan. Mathematics — From the Birth of Numbers. W. W. Norton, 1996. ISBN 039304002X. Одна из наиболее интересных книг об основах и пользе математики, которую можно читать одновременно и с пользой (например, о матрицах), и с удовольствием.

Knuth, Donald E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms, Third Edition. Addison-Wesley, 1998. ISBN: 0202496842.

Stewart, G. W. Matrix Algorithms, Volume I: Basic Decompositions. SIAM, 1998. ISBN 0898714141.

Wood, Alistair. Introduction to Numerical Analysis. Addison-Wesley, 1999. ISBN 020194291X.

Задание

1. Выведите на экран размеры типов char, short, int, long, float, double, int* и double* (используйте оператор sizeof, а не заголовок ).

2. Используя оператор sizeof, выведите на экран размеры объектов Matrix a(10), Matrix b(10), Matrix c(10), Matrix d(10,10), Matrix e(10, 10,10).

3. Выведите на печать количество элементов в каждом из объектов, перечисленных в упр. 2.

4. Напишите программу, вводящую числа типа int из потока cin и результат применения функции sqrt() к каждому из этих чисел int. Если функцию sqrt(x) нельзя применять к некоторым значениям x, выведите на экран сообщение “корень квадратный не существует” (т.е. проверяйте значения, возвращаемые функцией sqrt()).

5. Считайте десять чисел с плавающей точкой из потока ввода и запишите их в объект типа Matrix. Класс Matrix не имеет функции push_back(), поэтому будьте осторожны и предусмотрите реакцию на попытку ввести неверное количество чисел типа double. Выведите этот объект класса Matrix на экран.

6. Вычислите таблицу умножения [0,n]*[0,m] и представьте ее в виде двумерного объекта класса Matrix. Введите числа n и m из потока cin и аккуратно выведите на экран полученную таблицу (предполагается, что число m достаточно мало, чтобы результаты поместились в одной строке).

7. Введите из потока cin десять объектов класса complex (да, класс cin поддерживает оператор >> для типа complex) и поместите его в объект класса Matrix. Вычислите и выведите на экран сумму десяти комплексных матриц.

8. Запишите шесть чисел типа int в объект класса Matrix m(2,3) и выведите их на экран.

Контрольные вопросы

1. Кто выполняет числовые расчеты?

2. Что такое точность?

3. Что такое переполнение?

4. Каковы обычные размеры типов double и int?

5. Как обнаружить переполнение?

6. Как определить пределы изменения чисел, например наибольшее число типа int?

7. Что такое массив, строка и столбец?

8. Что такое многомерный массив в стиле языка C?

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

10. Что такое размерность матрицы?

11. Сколько размерностей может иметь матрица?

12. Что такое срезка?

13. Что такое пересылка? Приведите пример.

14. В чем заключается разница между индексированием в стиле языков Fortran и C?

15. Как применить операцию к каждому элементу матрицы? Приведите примеры.

16. Что такое объединенное умножение и сложение (fused operation)?

17. Дайте определение скалярного произведения.

18. Что такое линейная алгебра?

19. Опишите метод исключения Гаусса.

20. Что такое опорный элемент (в линейной алгебре и реальной жизни)?

21. Что делает число случайным?

22. Что такое равномерное распределение?

23. Где найти стандартные математические функции? Для каких типов аргументов они определены?

24. Что такое мнимая часть комплексного числа?

25. Чему равен корень квадратный из –1?

Термины

Упражнения

1. Аргументы функции f в выражениях a.apply(f) и apply(f,a) являются разными. Напишите функцию triple() для каждого варианта и примените их для удвоения элементов массива { 1 2 3 4 5 }. Определите отдельную функцию triple(), которую можно было бы использовать как в выражении a.apply(triple), так и в выражении apply(triple,a). Объясните, почему нецелесообразно писать все функции именно так для использования в качестве аргумента функции apply().

2. Повторите упр. 1, используя не функции, а объекты-функции. Подсказка: примеры можно найти в заголовке Matrix.h.

3. Только для экспертов (средствами, описанными в книге эту задачу решить невозможно). Напишите функцию apply(f,a), принимающую в качестве аргумента функции void (T&), T (const T&), а также эквивалентные им объекты-функции. Подсказка: Boost::bind.

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

5. Примените программу исключения методом Гаусса к системе A=={{0 1}{1 0}} и b=={5 6} и убедитесь, что программа завершится крахом. Затем попробуйте вызвать функцию elim_with_partial_pivot().

6. Замените циклами векторные операции dot_product() и scale_and_add() в программе исключения методом Гаусса. Протестируйте и прокомментируйте эту программу.

7. Перепишите программу исключения методом Гаусса без помощи библиотеки Matrix. Иначе говоря, используйте встроенные массивы или класс vector, а не класс Matrix.

8. Проиллюстрируйте метод исключения методом Гаусса.

9. Перепишите функцию apply(), не являющуюся членом класса Matrix, так, чтобы она возвращала объект класса Matrix, содержащий объекты, имеющие тип примененной функции. Иначе говоря, функция apply(f,a) должна возвращать объект класса Matrix, где R — тип значения, возвращаемого функцией f. Предупреждение: это решение требует информации о шаблонах, которая не излагалась в этой книге.

10. Насколько случайной является функция rand()? Напишите программу, принимающую два целых числа n и d из потока ввода, d раз вызывающую функцию randint(n) и записывающую результат. Выведите на экран количество выпавших чисел из каждого диапазона [0:n] и оцените, насколько постоянным является их количество. Выполните программу с небольшими значениями n и небольшими значениями d, чтобы убедиться в том, что очевидные смещения возникают только при небольшом количестве испытаний.

11. Напишите функцию swap_columns(), аналогичную функции swap_rows() из раздела 24.5.3. Очевидно, что для этого необходимо изучить код библиотеки Matrix. Не беспокойтесь об эффективности своей программы: быстродействие функции swap_columns() в принципе не может превышать быстродействие функции swap_rows().

12. Реализуйте операторы

Matrix operator*(Matrix&, Matrix&);

и

Matrix operator+(Matrix&, Matrix&).

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

Послесловие

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

Глава 25