Продолжаем тему вычислений. Чтобы считать эффективно, мало иметь много процессоров — надо написать хорошо распараллеленный алгоритм, который и сам по себе должен быть эффективен.
Например, пусть нам надо возвести матрицу в какую-то большую степень. По определению, степень матрицы — это матрица, умноженная сама на себя указанное число раз. Умножение матриц довольно дорого, это n³ операций умножения (и много сложений ещё). Если действовать тупо по определению, то сорок вторая степень — это 41 умножение матриц.
А можно разложить число 42 в двоичную систему, получив 101010, и тогда сорок вторая степень есть тридцать вторая плюс восьмая плюс квадрат. Квадрат получим одним матричным умножением, восьмую степень из квадрата — еще двумя, 32-ую из восьмой — еще двумя. В итоге 5 матричных умножений вместо 41: более, чем в восемь раз меньше! В целом же такая техника дает двоичный логарифм степени вместо самой степени, что довольно-таки экономно.
Пусть нам надо посчитать экспоненту через ряд Тейлора. Надо суммировать слагаемые вида xⁿ/n!. Напрашивается вычисление этих слагаемых по формуле, но это неправильно. Надо считать их рекуррентно, умножая предыдущее на x и деля на (n+1). Сами прикиньте, сколько операций делается лишних, повторно, при первом способе. Не говоря уж о том, что сами-то слагаемые невелики, скорее малы, а вот степень и особенно факториал весьма велики, а при делении больших чисел теряется точность.
То же самое происходит при подсчете числа сочетаний. Формула-то есть, но пользоваться ею надо с умом. С(n,m) есть n!/m!/(n-m)!, но там при подсчете много сокращается. Но и не надо это вычислять, а потом сокращать, тем более что при этом и точность теряется. Надо записать это длинно и скучно:
n/m (n-1)/(m-1) (n-2)/(m-2)...(m-1)/1
и перебирать множители по убыванию.
Давайте для примера решим классическую задачку про кодовый замок. Десять кнопок, нажимать надо три одновременно. В итоге разных кодов — число сочетаний из 10 по 3. По формуле, это 10!=3628800 делить на 3!=6 и потом на 7!=5040. Однако 7! и 10! сокращаются, оставляя в числителе 8∙9∙10, а знаменателе остается 3!=1∙2∙3. Можно сократить ещё немного, сведя подсчёт к 4∙3∙5=120.
Потеря точности — тоже важная вещь. Число хранится в виде "целая часть, дробная часть, степень десятки". Или не десятки, а двойки, или ещё чего-нибудь, тогда еще надо следить за конвертацией. Но пусть пока десятки.
Если у вас 16 цифр после запятой, то это точность порядка 10⁻¹⁶. Но если у вас степень десятки сорок вторая, то абсолютная погрешность получается очень большая. Так что вычислить два больших факториала и поделить их друг на друга — идея весьма хреновая.
А целочисленная арифметика такого вообще может не выдержать.
Сам факториал определен рекурсивно: n!=n(n-1)!, 0!=1. Но вычислять его рекурсивно — идея плохая, так как рекурсивный вызов — дорогая вещь. Цикл на порядок быстрее. И проще.
Вот на скриншоте результат теста. На языке Вим, который для вычислений не предназначен (тем чище эксперимент), получилось вот что: факториал 20 циклом вычислен в полтора раза быстрее, чем рекурсией. Но для 30 наступает переполнение и получаем отрицательный (нелепый) результат.
Сделав вычисления вещественными (просто заменой 1 на 1.0) мы получаем возможность вычислить и факториал 50: цикл занимает 0.000415 времени, а рекурсия вдвое больше: 0.000977. При этом результат один и тот же: 3.041409e64.
А вот для 100 рекурсия не справляется, так как глубина рекурсии превышает установленный предел. А циклу норм, он спокойно вычисляет 9.332622e157
На Перле я бы считал факториал так:
perl -E '$n=shift || 5; $S=1;map {$S*=$_} 1..$n; say "F($n)=$S"' 42
Кстати, как я уже рассказывал, факториал вообще необязательно вычислять. Формула Стирлинга дает точность 8% уже для n=1, а потом всё точнее и точнее. После 10 относительная погрешность уже менее одного процента. А по скорости Стирлинг в десять раз быстрее цикла.
Еще о точности. Все мы знаем, что синус представляется рядом Тейлора со слагаемыми вида (-1)ⁿx⁻ⁿ/n! при нечетных n . Знаем и то, что ряд этот сходится, и быстро, на всей комплексной плоскости, так что пригоден для вычисления синуса. Например, взяв шесть слагаемых при x=π/6:
perl -MMath::Trig -E'$pi=acos(-1);$x=1*$pi/6;$f=$x;$S=$x;for(1..5){$f*=-$x*$x/(2*$_)/(2*$_+1);$S+=$f}say"sin(п/".$pi/$x.")=$S"'
получим sin(п/6)=0.499999999999964
Однако взяв 50 слагаемых и поставив множитель 601 при $pi, с удивлением получим нелепый ответ: 1.85579666374312e+92. Это, если кто не знаком с такой нотацией, 1.8 с небольшим, умноженное на 10 в степени 92 (прописью "девяносто два"). Это больше числа электронов в видимой Вселенной. Это не просто плохо, это предельно плохо.
На пяти слагаемых не сильно лучше: степень 19. И это не Перл виноват, это ошибка в методе. Теоретически-то оно сходится, а вот практически - только если достаточно слагаемых И если хватит машинной точности.
То есть да — Тейлор гарантирует сходимость; но есть нюанс. Даже синус от 2п не надо так вычислять, а надо применить школьные формулы. А уж что 100п — это не то же самое, что 2п, наверное, не надо даже начинать.
Вот еще пример. Есть такой прекрасный метод Гаусса для решения систем линейных уравнений. Он без труда программируется, и допустим, что мы используем обычные вещественные числа, одинарной точности. Те, что в Фортране REAL, а в С — float. Матрицу берем случайную, решение тоже, умножаем и получаем правую часть, ну и решаем систему, чтобы восстановить решение. И сравниваем. Для матрицы 10 на 10 получили относительную погрешность в 0.0004% и максимальное отклонение в 3 миллионных.
Берем матрицу 1000 на 1000. Погрешность 4.5%, уклонение 0.07. Берем матрицу порядка 2000 и всё: погрешность в тысячах процентов, уклонение порядка 30 (притом случайные числа были от 0 до 1).
Да, в двойной точности таких эффектов при таких порядках нет, но они появляются тоже, если матрицы больше и если разброс чисел больше. Когда-то в юности я столкнулся с тем, что метод Гаусса на больших матрицах меня подвел, и двойная точность не спасла.
Проблема в том, что преподаватели матанализа часто сами не считают, а счетчики не вникают глубоко в матан.
