Дмитpий Hecтepук

Блог о программировании — C#, F#, C++, архитектура, и многое другое

Archive for the ‘Mathematics’ Category

Приемы взятия сложных интегралов

3 комментария

Интегралы, что может быть веселее? Ну, возможно не для всех, но все же, я уже давно ничего не постил такого сугубо математического, так что попробую. Этот пост – про то как брать «сложные» интегралы. Этот пост подразумевает что читатель учился таки в школе и знает тривиальные подходы (например, интегрирование по частям). В посте мы будем обсуждать только интегралы Римана, а не интегралы Лебега-Стилтьеса, Ито, Скорохода и так далее (хотя я бы с удовольствием, чесслово).

Весь этот пост — маленькая выборка рецептов или «паттернов» которые можно взять в копилку и потом применять. Пост рекомендуется читать на high-DPI дисплее дабы предотвратить глазное кровотечение. Я предупредил.

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

Начнем с немного избитого метода — перехода к полярным координатам. Примечательно, что переход к полярным координатам можно применять даже там где, казалось бы, речь о декартовых координатах не идет вообще. Например, неопределенный интеграл Гаусса \textstyle \int e^{-x^2} {\mathrm d}x не имеет аналитического решения, а вот определенный интеграл \textstyle \int_{-\infty}^{\infty} e^{-x^2} {\mathrm d}x = \sqrt{\pi}.

Доказать это можно вот как: сначала, чтобы применить преобразование координат, мы вводим две переменные интегрирования \textstyle x и \textstyle y так что

I = \int_{-\infty}^{\infty} e^{-x^2} {\mathrm d}x = \int_{-\infty}^{\infty} e^{-y^2} {\mathrm d}y

Декартовы координаты можно выразить через полярные \textstyle (r, \theta) вот так:

\begin{align*}
x &= r \cos \theta \\
y &= r \sin \theta \\
r^2 &= x^2 + y^2
\end{align*}

Интегрирование от \textstyle -\infty до \textstyle \infty в декартовой системе координат — это то же, что интегрирование \textstyle r от \textstyle 0 до \textstyle \infty и \textstyle \theta от \textstyle 0 до \textstyle 2\pi.

В результате получим следующее:

\begin{aligned}
I\cdot I &= \int_{-\infty}^{\infty} e^{-x^2} {\mathrm d}x \int_{-\infty}^{\infty} e^{-y^2} {\mathrm d}y \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-x^2-y^2} \;{\mathrm d}x\;{\mathrm d}y \\
&= \int_{0}^{2\pi} {\mathrm d}\theta \int_{0}^{\infty} e^{-r^2} r \;{\mathrm d}r \\
&= 2\pi\int_{0}^{\infty}  e^{-r^2} r \;{\mathrm d}r \\
&= \pi\int_0^{\infty} e^{-r^2} \;{\mathrm d}r^2 = \pi \\
\end{aligned}

\therefore I = \sqrt{\pi}

Этот же подход может применять и в 3-х измерениях с использованим сферических координат \textstyle (x,y,z) \rightarrow (r,\theta,\phi).

Геометрические интерпретации

Вообще, «скатывание в геометрию» порой приносит плоды. Вот например допустим вам надо посчитать

\int_0^\infty \frac{{\mathrm d}x}{1+x^2}

Уверен, многие из вас знают что у этого интеграла есть аналитическое решение \textstyle \tan^{-1}x, поэтому посчитать определенный интеграл не составляет труда. Но на самом деле, этот интеграл можно посчитать даже без этого знания.

Представьте круг с радиусом \textstyle r с центром \textstyle (0,0). Длина дуги этого круга с центральным углом \textstyle \theta равна \textstyle L = r\theta, а если круг единичный – то просто \textstyle \theta. Тогда

L = \theta = \int_0^{\theta} \;{\mathrm d}t

где \textstyle t — это произвольная переменная интегрирования.

При таком раскладе, подынтегральное выражение равно \textstyle 1, но мы можем его усложнить, например

\begin{align*}
L &= \int_0^{\theta}1 \;{\mathrm d}t \\
&= \int_0^{\theta}\frac{\frac{1}{\cos^2t}}{\frac{1}{\cos^2t}} \;{\mathrm d}t \\
&= \int_0^{\theta}\frac{\frac{1}{\cos^2t}}{\frac{\cos^2t+\sin^2t}{\cos^2t}} \;{\mathrm d}t \\
&= \int_0^{\theta}\frac{\frac{1}{\cos^2t}}{1+\tan^2t} \;{\mathrm d}t \\
\end{align*}

Далее, делаем подстановку

x = \tan t \Rightarrow {\mathrm d}x = \frac{{\mathrm d}t}{\cos^2 t}

Тем самым, получаем

L = \int_0^{\tan \theta}\frac{{\mathrm d}x}{1+x^2}

Допустим что \textstyle \theta = \frac{\pi}{2}. Тогда \textstyle \tan \theta = \tan \frac{\pi}{2} = \infty, а поскольку \textstyle \frac{\pi}{2} отмеряет нам ровно четверть круга (длина всего единичного круга \textstyle 2\pi), мы моментально получаем результат

\frac{\pi}{2}=\int_0^{\infty} \frac{{\mathrm d}x}{1+x^2}

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

\begin{align*}
\frac{\pi}{4} &= \int_0^1 \frac{{\mathrm d}x}{1+x^2} \\
\frac{\pi}{3} &= \int_0^{\sqrt{3}} \frac{{\mathrm d}x}{1+x^2} \\
\end{align*}

и так далее.

Разбиение диапазона интегрирования

Допустим вам надо посчитать

\int_0^{\infty} \frac{\ln x}{1 + x^2} \;{\mathrm d}x

Для взятия этого интеграла, разобъем диапазон интегрирования на два, т.к. \textstyle \int_0^{\infty}=\int_0^1+\int_1^{\infty}.

Займемся сначала первым интегралом, т.е. \textstyle \int_0^1. Сделаем подстановку \textstyle t = 1/x \Rightarrow {\mathrm d}x=-{\mathrm d} t/t^2 . Получим

\begin{align*}
\int_0^1 \frac{\ln x}{1+x^2} \;{\mathrm d}x &= \int_{\infty}^1 \frac{\ln(1/t)}{1+1/(t^2)}\left(-\frac{1}{t^2}\;{\mathrm d}t\right) \\
&= - \int_{\infty}^1 \frac{\ln(1/t)}{t^2+1}\;{\mathrm d}t \\
&= \int_1^{\infty} \frac{\ln(1/t)}{t^2+1}\;{\mathrm d}t \\
&= - \int_1^{\infty} \frac{\ln t}{t^2+1}\;{\mathrm d}t
\end{align*}

То есть внезапно оказалось, что поставленная переменная \textstyle t выполняет такую же функцию что и \textstyle x. Другими словами, \textstyle \int_0^1 = -\int_1^{\infty} а это значит что мы автоматически получаем значение искомого интеграла:

\int_0^{\infty}\frac{\ln x}{1+x^2}\;{\mathrm d}x = 0

Разбиениe на четное и нечетное

Вот нужно вам например посчитать

\int_{-1}^{1} \frac{\cos x}{e^{1/x}+1} \;{\mathrm d}x

Давайте сделаем несколько замен:

\begin{align*}
f(x) &:= e^{1/x} \\
g(x) &:= \frac{\cos x}{f(x)+1}
\end{align*}

Теперь нам нужно посчитать \textstyle \int_{-1}^{1} g(x) \;{\mathrm d}x, и вот тут начинается самое интересное. Мы переписываем \textstyle g(x) как сумму четной и нечетной функции:

g(x) = g_e(x) + g_o(x)

Многие спросят «а так вообще можно?» — на самом деле да, и вот почему. Возьмите и воткните в определение выше \textstyle -x вместо \textstyle x. Вы получите

g(-x)=g_e(-x)+g_o(-x)=g_e(x) - g_o(x)

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

g_e(x)=\frac{g(x)+g(-x)}{2}

и

g_o(x)=\frac{g(x)-g(-x)}{2}

Так-то. Соответственно, наш интеграл можно переписать как

\int_{-1}^{1}g(x) \;{\mathrm d}x = \int_{-1}^{1}g_e(x) \;{\mathrm d}x + \int_{-1}^{1}g_o(x) \;{\mathrm d}x = \int_{-1}^{1}g_e(x) \;{\mathrm d}x

Как видно выше, нечетная функция пропала полностью, осталась только четная сторона, т.к.

\int_{-1}^{1}g_o(x) \;{\mathrm d}x = 0

Ладно, вам уже наверное надоело ждать сути этого примера. Так вот, у нас есть формула \textstyle g_e(x)=\frac{g(x)+g(-x)}{2}, дайвате воткнем в эту формулу \textstyle g(x). Мы получим

g_e(x)=\frac{1}{2}\left(\frac{\cos x}{f(x)+1}+\frac{\cos (-x)}{f(-x)+1}\right)

Но мы-то знаем, что \textstyle \cos x — четная функция, поэтому \textstyle g_e(x) можно переписать как

\begin{align*}
g_e(x) &= \frac{\cos x}{2}\left(\frac{1}{f(x)+1} + \frac{1}{f(-x)+1}\right) \\
&= \frac{\cos x}{2}\left(\frac{f(-x)+1+f(x)+1}{f(x)f(-x)+f(x)+f(-x)+1}\right) \\
&= \frac{\cos x}{2}\left(\frac{2+f(-x)+f(x)}{f(x)f(-x)+f(x)+f(-x)+1}\right) \\
\end{align*}

Это какое-то месиво и непонятно что с ним делать. Но с другой стороны посмотрите, у нас в формуле присутствует \textstyle f(x)f(-x). Давайте вспомним, что \textstyle f(x)=e^{1/x} и мы получим

f(x)f(-x)=e^{1/x}e^{-1/x}=e^0=1

Ну вот и всё — наша страшная дробь выше уже совсем не страшная т.к. числитель и знаменатель равны, а это значит что

g_e(x) = \frac{\cos x}{2}

а сам интеграл теперь легко посчитать:

\int_{-1}^{1} \frac{\cos x}{e^{1/x}+1}\;{\mathrm d}x &= \int_{-1}^{1} \frac{\cos x}{2} \;{\mathrm d}x = \sin(1) = 0.841...

Хотите ещё?

Я на самом деле понял, что по объему для одного поста вполне достаточно. Сорри если что написал не так — я по-русски прочитал ровно нуль математических книг (чего и вам советую), так что терминология может страдать.

Существует еще вагон разных трюков, так что, если интересно, советую глянуть соответствующую литературу. Удачи! ■

Реклама

Written by Dmitri

9 ноября 2016 at 16:57

Опубликовано в Mathematics

Tagged with

Задача про 2008

2 комментария

Этот пост — репост одного моего поста который был на Hexlet до того как Hexlet решил изобрести себя снова и сдвинул мой курс в «песочницу». Он немного технический, и я не хочу его терять. Надеюсь вы поймете. К сожалению, комментарии оригинального поста потеряны, но с этим ничего не поделать.

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

Как с помощью 13 (тринадцати) нулей и простых математических операций получить число 2008?

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

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

Вещественность

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

Самое очевидное решение — это 0^0 = 1 (это кстати спорно и вызвало много дебатов), но это в результате дает нам \left \lfloor{\frac{13}{2}}\right \rfloor = 6 нулей из которых, как вы понимаете, 2008 никаким образом не сделать. Должны ведь быть и другие варианты? Попробуйте угадать, какой оператор делает из нуля единицу?

Конечно же это факториал. Факториал — это обычно произведение всех чисел от 1 до n, так что например 3!=3\times2\times1, но у факториала есть одна особенность: 0!=1, а следовательно у нас теперь 13 единиц, а с ними уже можно работать. И кстати, пока не поздно, обращу внимание на то, что в условии задачи я упомянул простые математические операции. Если бы можно было использовать, например, логарифмы, то задача уже решена т.к. любое натуральное число можно выразить через логарифмы, корни, и три двойки (доказательство), а соответственно решение нашей задачи выглядело бы вот так:

\displaystyle 2008 = - \log_2\left(\log_2\left(\underbrace{\sqrt{\sqrt{...\sqrt{2}}}}_{2008}\right)\right)

Для решения выше потребовалось бы всего шесть нулей (т.к. 2 = 0! + 0!) а также 2008 квадратных корней но, как я уже сказал, хочется получить решение с помощью известных операторов, не привлекая в решение всякие сложные функции. Корни, возведение в степень, факториалы — это все ок.

Возведение в степень

Итак, нужно получить 2008. Самое очевидное — это найти близлежащее число 2^n. В данном случае, это 2^{11} и соответственно мы можем подобраться к решению вот так:

\displaystyle \begin{aligned}  2008 &= 2^{11} - 40 \\  &= 2^{2^3+3} - 5\times2^3  \end{aligned}

К сожалению, числа 2, 3 и 5 идут «по себестоимости» (т.е. чтобы сделать 5 нужно 5 нулей), и соответственно в 13 нулей мы ну никак не уложимся: для решения выше нужно аж 20 нулей!

Соответственно нашей первоочередной задачей становится поиск более «дешевых» чисел, и тут на помощь приходят…

Факториалы

Факториалы дают нам не только единицы с которыми можно работать, но также дешевые крупные числа. Например, 3!=6 а это значит что за три нуля мы можем получить не только 3 но также 6=3!, 720=(3!)!, и так далее. Если форма записи X_Y обозначает что число X требует для записи Y нулей, то мы получим следующие значения:

\displaystyle \begin{aligned}  0_1 =& 0 \\  1_1 =& 0! \\  2_2 =& 0! + 0! \\  3_3 =& 0! + 0! + 0! \\  4_4 =& 3 + 0! = 2^2 \\  5_4 =& 3! - 0! \\  6_3 =& 3! \\  7_4 =& 3! + 1 \\  8_5 =& 3! + 0! + 0! = 2^3 \\  9_5 =& 3^2  \end{aligned}

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

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

\displaystyle \begin{aligned}  2008 &= 2^3 \times 251 \\  &= 2^3 \times (216+36-1) \\  &= 2^3 \times (6^3 + 6^2 - 1) \\  &= \underbrace{2^3 \times \left(6^2(6+1)-1)\right)}_{14}  \end{aligned}

Увы, решение выше требует на один ноль больше, чем у нас имеется.

Другой подход к решению задачи — это разложить на множители не 2008 а, скажем 2008\pm10. Вот что нам дает MATLAB:

1998=2   3   3   3  37
1999=1999
2000=2  2  2  2  5  5  5
2001=3  23  29
2002=2   7  11  13
2003=2003
2004=2    2    3  167
2005=5  401
2006=2  17  59
2007=3    3  223
2008=2    2    2  251
2009=7   7  41
2010=2   3   5  67
2011=2011
2012=2    2  503
2013=3  11  61
2014=2  19  53
2015=5  13  31
2016=2  2  2  2  2  3  3  7
2017=2017
2018=2  1009

Из интересных чисел тут, пожалуй, 2000 и 2016. Также соблазнительно число 1998 (у него есть множитель 37=3!^2+1). Возьмем сначала 2000:

\displaystyle \begin{aligned}  2000 =& \underbrace{2^4\times5^3}_{12} \\  =& \underbrace{2\times\left(2^{3}+2\right)^3}_{11}  \end{aligned}

Это решение нам явно не подходит т.к. 2008 = 2000+8_5 т.е. для решения нужно 16 единиц. Попробуем 2016:

\displaystyle \begin{aligned}  2016 =& \underbrace{2^3\times6^2\times7}_{14} \\  =& \underbrace{6^4+(3!)!}_{10}  \end{aligned}

Но даже так, у нас остается три нуля, а 2008=2016-8_5, то есть для этого решения нам всяко не хватает двух нулей. Что же, подход с обычными (намек!) факториалами почти себя исчерпал, давайте еще попробуем что-нибудь сделать с числом 1998:

\displaystyle \begin{aligned}  1998 &= \underbrace{3!^2\times\left(3!^2+1\right)}_{11}  \end{aligned}

Это достаточно бессмысленное занятие т.к. двумя нулями лишние 10 не набрать, но, по крайней мере, мы выловили тот факт, что возможно стоит попробовать 3!^4:

\displaystyle \begin{aligned}  2008 &= 3!^4 + 712 \\  &= \underbrace{3!^4 + (3!)! - 2^3}_{15}  \end{aligned}

Увы и ах. Как я уже говорил в самом начале, число 2008 выбрано неслучайно — это число вставляет максимальное число палок в колеса. Нужен какой-то трюк чтобы снизить «стоимость» простых чисел вроде 8. К счастью, такой трюк имеется, и это…

Двойной факториал

Двойной факториал n!! — это тоже произведение всех чисел от 1 до n, но с шагом 2, т.к. например 5!!=5\times3\times1=15. Само по себе это кардинально меняет картину, и я хочу обратить ваше внимание на два равенства:

  • 8_4 = 4!!, то есть восьмерка стала «дешевле» на один ноль.

  • 48_3=(3!)!!

Итак, делим 2008 на 48 и получаем 41.8(3), а поскольку 42=6\times7, мы наконец-то можем попробовать получить ответ:

\displaystyle \begin{aligned}  2008 &= 48\times42-8 \\  &= \underbrace{(3!)!!\times\left(3!\times(3!+1)\right)-4!!}_{14}  \end{aligned}

Но что-то тут не так: 42 нам далось адским трудом, мы заплатили за него аж 7 нулей. На один ноль поменьше и все получится. На самом же деле, 42=48-6=(3!)!!-3! и вот, ура, у нас готово первое решение:

\displaystyle 2008 = \underbrace{(3!)!!\times\left((3!)!!-3!\right)-4!!}_{13}

Субфакториал

Разных факториалов бывает много и субфакториалы — это такой особый тип факториала, который определяет количество беспорядков порядка n. Высчитывается он так

\displaystyle !n=n!\sum_{k=0}^n \frac{\left(-1\right)^k}{k!}

Прелесть этого факториала в том, что он в очередной раз дает нам другие равенства, например 265_3 = !(3!). А это в свою очередь заставляет нас в очередной раз посмотреть на разложение 2008 на множители:

\displaystyle \begin{aligned}  2008 &= 4!! \times 251 \\  &= 4!! \times (!(3!) - 14) \\  &= \underbrace{4!! \times \left(!(3!)-5!!+1\right)}_{13}  \end{aligned}

Оптимизация задачи

Кому-то может показаться, что 13 нулей — это предел мечтания, но нет — вот например пара решений, где используются только 12 нулей:

\displaystyle \begin{aligned}  2008 &= \underbrace{(!5)^2 + \sqrt{\left((3!)!+1\right)} + 1}_{12} \\  &= \underbrace{\left(!(4!!)-1\right)\times2+5!}_{12}  \end{aligned}

Предлагаю игру в стиле code golf — оставляйте в комментариях решения, которые используют 11 нулей и меньше, посмотрим чье кунг‑фу лучше. No cheating! ■

Written by Dmitri

15 февраля 2016 at 0:39

Опубликовано в Mathematics

Tagged with

Генерация псевдослучайных чисел

2 комментария

Должен признаться, что меня никогда толком не волновли механизмы генерации случайных чисел. Ну есть себе rand() или Random, и пускай живет. Но в последнее время, особенно когда внезапно захотелось генерировать числа с нормальным распределением, пришлось немного вникнуть в тему. Поэтому этот пост – про случайные числа.

Равномерное Распределение

Когда вы вызываете какой-нибудь rand() – вы получаете псевдослучайное число с равномерным распределением. Давайте рассмотрим, почему псевдо, и почему с равномерным.

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

Кому-то может показаться, что псевдослучайно – значит плохо. Это вовсе не так. Например, если у вас есть генератор случайных картинок, вы можете сгенерировать много картинок с небольшим разрешением, а потом воссоздать понравившуюся вам картинку с большим разрешением. Главное запомнить seed value каждой картинки что, собственно, не сложно.

Теперь насчет равномерности распределения. Вообще, идея достаточно простая – когда вы генерируете число от 0 до 1, идея в том что вероятность любого числа – одинакова. Иначе говоря, если \theta_1 < \theta_2, то случайная величина Y имеет непрерывное равномерное вероятностное распределение на интервале (\theta_1, \theta_2) если фунция плотности Y выглядит вот так:

\displaystyle f(x)=\frac{1}{\theta_{1}-\theta_{2}}, \theta_{1} \leq y \leq \theta_{2}

Вообщем-то, самая примитивная реализация этой функции это нечто вроде x_n = ax_{n-1}+b \ (\mathrm{mod} \ m), но на практике используются чуть более сложные варианты которые дают более «качественные» случайные числа. И кому-то из вас наверное хочется закричать что «это все булшит» и «какая разница?», но когда вы делаете, например, Монте-Карло симуляции в которых случайных чисел нужно реально много, важно иметь хорошие случайные числа чтобы ваша модель ближе отражала полный хаос который творится во вселенной. (Правильное слово энтропия, кстати.)

Есть целый ряд статистических тестов для проверки качества случайных чисел. Мой подход к этому вопросу – брать лучшее, что можно получить в тех или иных случаях, например рутину ran4 из Numeical Recipes in C. Хотя, по правде сказать, даже дефолтная реализация случайных чисел в .NET – сделана «по Кнуту», и является вполне юзабельной.

Что еще сказать про равномерно распределенные случайные числа? Ну, то что среднее значение – это среднее всего диапазона, т.е.

\displaystyle E[Y]=\frac{\theta_1+\theta_2}{2}

А дисперсия равна

\displaystyle V[Y]=\frac{(\theta_2-\theta_1)^2}{12}

Последний результат особенно интересен. Попробуйте угадать, что будет если сложить 12 случайных чисел и отнять 6!

Нормальное Распределение

Нормальное распределение само по себе выглядит чуть иначе и имеет функцию плотности

\displaystyle f(y)=\frac{1}{\sigma\sqrt{2\pi}}\exp(\frac{-(y-\mu)^2}{2\sigma^2})

Соответственно, это распределение параметризовано средним и дисперсией, т.е. E[Y]=\mu и V[Y]=\sigma^2.

Грубым, но действительно работающим приближением к получению случайного нормально распределенного значения является подход, в котором мы берем, скажем, 12 случайных равномерно распределенных значений и вычитаем 6. Почему? Это возможно благодаря центральной предельной теореме, котороая говорит нам, что сумма большого количества независимых случайных величин (мы это называем i.i.d. – independent and identically distributed) имеет распределение, близкое к нормальному. Соответственно, количество случайных чисел (12) выбрано не случайно, а так чтобы «погасить» дисперсию.

Случай выше подходит разве что для простых вычислений в Excel, где нет встроенной функции RAND() с нормальным распределением. Если же говорить о программировании, то тут скорее полезно преобразование Бокса-Мюллера. Идея в том, что мы берем два случайных числа v_1 и v_2 и проверяем что они находятся в единичном круге (т.е. v_1^2+v_2^2<1). Если они подходят, то берем и возвращаем одно из этих двух чисел, наприемер v_1, но не само по себе а как

\displaystyle   r = v_1^2+v_2^2 \\  f = \sqrt{-2 \frac{\log{r}}{r}} \\  \mathrm{return} \ v_1 \times f

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

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

Written by Dmitri

16 апреля 2013 at 22:47

Опубликовано в Mathematics

Tagged with

Оптимизация математических выражений

leave a comment »

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

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

Умножение

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

Например, если изначально у нас есть формула

\displaystyle y=ax^2+bx+c

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

\displaystyle y=x(ax+b)+c

В первой формуле – 3 умножения и 2 сложения. Во втором – 2 умножения и 2 сложения. Так или иначе мы сэкономили стоимость одного умножения, т.е. срезали как минимум 20% стоимости этого вычисления.

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

  • a^2-b^2=(a+b)(a-b) имеет смысл при условии что умножение дороже сложения.

  • x^2+10x+16=(x+8)(x+2) содержит одно умножение вместо 3-х (!) но подразумевает решение уравнения. Такую оптимизацию очень сложно обобщить, и в общем случае она подразумевает решение практически бесконечного количества уравнений.

Возведение в степень

Функции возведения в степень обычно заточены под нецелые степени, т.е. степень это значение типа float или double и соответственно для вычисления используются ряды Тейлора. Это значит что цена вычисления x^2 в десятки (!) раз больше чем цена x \cdot x.

В долгосрочной перспективе, любое вычисление x^n где n \ge 2, n \in \mathbb{Z} подразумевает \left\lceil \log_{2}n\right\rceil умножений. В принципе, некоторые языки (например С++) имеют перегрузки метода возведения в степень, где параметр степени – целое число (т.е. int). В этом случае никаких дополнительных усилий не требуется. Если же такой функции не существует (например, в .NET), то есть два варианта:

  • Для n=2,3 имеет смысл просто заменить возведение в степень на умножение, т.к. это эффективней и понятнее чем использование сторонней функции.

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

Конечно, существует и просто неэффективное использование функции возведения в степень: например, не имеет смысла считать x^{0.5} когда вместо этого можно намного быстрее посчитать \sqrt{x}. Ну и не будем забывать про совсем уж специализированные оптимизации расчета функций, например всем известный механизм расчета 1/\sqrt{x} Джона Кармака.

Деление

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

\displaystyle \frac{a}{x}+\frac{b}{x}=\frac{a+b}{x}

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

\displaystyle \frac{a}{b}+\frac{c}{d}=\frac{ad+bc}{bd}

Два деления вместо одного имеют смысл, особенно когда деление в 6 раз медленнее умножения.

Введение временных переменных

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

\displaystyle ax^6+bx^3=x^3(ax^3+b)

мы продолжаем считать x^3 дважды, что влечет за собой 4 умножения вместо 2-х. Гораздо проще внедрить временную переменную y=x^3 и потом просто посчитать y(ay+b).

Заключение

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

Ну и не забывайте, что компилятор тоже не дремлет, и некоторые вещи все же оптимизирует, например:

\displaystyle x/4=x\gg 2

Так что оптимизируйте с удовольствием, или используйте программы, которые оптимизируют за вас. ■

Written by Dmitri

4 декабря 2012 at 21:42

Опубликовано в Mathematics

Снова про оптимизацию математических выражений

11 комментариев

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

Инлайнинг возведения в степень

Если ваша программа реализует f(x)=ax^2+bx+c как

double f(double a, double b, double c, double x)
{
  return a * std::pow(x, 2.0) + b*x + c;
}

то вы оказываете пользователю медвежью услугу. Для начала, следует помнить, что возведение в целочисленную степень является достаточно бессмысленной операцией, т.к. простое умножение “инлайн” или умножение в цикле всегда будет быстрее. Причина этого проста – функции вроде std::pow() расчитаны на подсчет степени для любых показателей, в т.ч. нецелочисленных. Как только в экспоненту ставится число вроде 2.0, строятся ряды и производительность идет лесом.

Решение этой проблемы весьма простое. Во-первых, для всех простых показателей вроде x^2, x^3 и так далее лучше все инлайнить, т.е. писать x*x, x*x*x и т.д. Что касается более крупных степеней то тут тоже все просто – надо чтобы степень была целым числом, т.е. имела тип int. Тогда будет вызвана пегрузка std::pow(), которая в свою очередь вызовет функцию _Pow_int() (по крайней мере в <math.h> от Microsoft), и все будет хорошо.

В .NET все не так ажурно т.к. Math.Pow() берет два параметра типа double. Очень досадное упущение, но это можно устранить самостоятельно:

public static long IntPower(int x, short power)
{
  if (power == 0) return 1;
  if (power == 1) return x;
  int n = 15;
  while ((power <<= 1) >= 0) n--;
  long tmp = x;
  while (--n > 0)
    tmp = tmp * tmp *
         (((power <<= 1) < 0) ? x : 1);
  return tmp;
}

Пример выше – для типа short, т.к. его обычно вполне достаточно. Естественно, следует также понимать что есть некий порог, после которого дешевле использовать канонический pow(double,double) вместо итеративного умножения.

Сокращение количества умножений

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

double f(double a, double b, double c, double x)
{
  return a*x*x + b*x + c;
}

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

double f(double a, double b, double c, double c)
{
  return x*(a*x + b) + c;
}

На самом деле, факторизация выражений, т.е. вынесение за скобки общих членов — задача не из легких, т.к. требует чтобы система знала уйму всяких разных правил. Например, как записать f(x,y)=x^2-y^2? Вы действительно думаете, что без двух умножений не обойтись? Ничего подобного!

double f(double x, double y)
{
  return (x + y) * (x - y);
}

И если равенство a^2-b^2=(a+b)(a-b) мы еще можем запомнить, то какие-то более сложные равенства уже не получится, и придется жить без них. Хотя у пользователя всегда есть вариант использовать какой-нибудь полноценный математический пакет и собственноручно оптимизировать вычисление.

Предварительные вычисления

Понятное дело что выражение \frac{1}{3}x превратится во что-то вроде 0.333333333333 * x будучи оптимизированным компилятором. Но если взять выражение вроде \frac{2x}{3y}, то тут уже все неоднозначно. До какой степени компилятор сможет оптимизировать и “инлайнить” результаты вычислений? Я думаю что в репертуар типичного компилятора не входит глубокий анализ и оптимизация вычислений.

Выделение переменных

Одна из проблем с генерацией кода для вычислений (к пр. на основе Excel) заключается в том, что есть большой шанс выполнить одну и ту же дорогую операцию несколько раз. Например, если напрямую конвертировать выражение 2cos^3x-cos^2x, то в результате можно получить весьма предсказуемый и неэффективный код:

double f(double x)
{
  return 2.0*pow(cos(x),3.0) - pow(cos(x), 2.0);
}

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

double f(double x)
{
  double u = cos(x);
  return (2*u - 1)*u*u;
}

То, что приведено выше – это так называемая схема Горнера. Примечательно, что эта задача — анализ выражений на поиск общих кусков кода — будучи созданной ускорять код, сама по себе имеет огромную стоимость, т.к. подразумевает поиск любого поддерева с определенной стоимостью в остальных частях дерева и выделения их для последующей замены. В общем случае эта задача не решается. К счастью, при трансляции Excel→С++ ее можно решить хотя бы частично: результаты промежуточных вычислений часто хранятся в отдельных клетках, а детектировать их повторное использование легче т.к. это всего лишь сравнение cell reference’ов.

Заключение

Для создания вменяемых математических моделей недостаточно просто конвертировать код “один в один”. В конце концов, главная цель подобных конверсий из Excel в С++ — получить многократный прирост производительности. А сколько эта конверсия займет времени никого толком не волнует.

Written by Dmitri

12 февраля 2012 at 22:26

Опубликовано в Mathematics