Сам интерес понизить сложность вычисления экспоненциальных рядов был вызван возможностью использовать не только действительные числа, но и числа другой природы, других алгебр. Для них можем дать определение операций сложения, умножения, что такое нулевой и единичный элемент. Используя эти операции, можем вычислить и функцию экспоненты. Но при усложнении операции умножения время вычисления функции растет чувствительно. И, если бы нашлась возможность понизить вычислительную сложность экспоненциального ряда, то было бы здорово.
Рассмотрим разложение в ряд функции экспоненты: ex=1+x1!+x22!+x33!+x44!+… В этом ряду каждый последующий элемент содержит в качестве своей части предыдущие. Чтобы было видно лучше, перепишем чуть иначе: 1+x1+x1x2+x1x2x3+x1x2x3x4+x1x2x3x4x5+… В действительности мы не вычисляем бесконечные ряды. Если бы так было, то никакие вычисления не закончились бы. Конечно, мы ограничиваем выисление некоторым членом. Проведем предварительные оценивания, при каком числе членов ряда n достигается разумная точность. В моем эксперименте это оказалось n=20. Увеличение n уже не давало существенного увеличения точности, но давало рост времени вычисления.
Ограничим ряд некоторым числом. Правый хвост ряда будет выглядеть в случае ограничения на n=5: x1x2x3x4+x1x2x3x4x5 Этот хвост можно сократить до: x1x2x3x4(1+x5) Приложим к этой части хвоста следующую от нее слева часть: x1x2x3+x1x2x3x4(1+x5) Эту часть хвоста также можно сократить: x1x2x3(1+x4(1+x5)) Соответственно, рекурсивно продолжая до самого начала получим итеративное вычисление со сложностью по умножению порядка n. Здесь для оценки сложности считаем лишь умножения элементов алгебры x между собой, считая операции деления на действительное число несущественным по сложности.
Итого, если у нас есть определение класса с операциями сложения, умножения, инициализации единичного и нулевого элемента, а также деления на действительное число, то шаблон функции экспоненты может выглядеть примерно так:
template <class T> T exp(const T& x) { T accum(0.0); for(int i = 20; i >= 1; i--) accum = (x / (double)i) * ( T(1.0) + accum); return accum + T(1.0); };Таким образом, действительно удалось понизить сложность вычисления экспоненциального ряда с n2 до n.
Функции вычисления синусов и косинусов гиперболических и тригонометрических также понижаемы, нужно лишь учесть их чередование и знакопеременность.
template <class T> T cosh(const T& x) { T accum(0.0); for(int i = 20; i >= 2; i -= 2) accum = ( x / (double)(i-1) ) * ( x / (double)i ) * ( T(1.0) + accum ); return accum + T(1.0); }; template <class T> T sinh(const T& x) { T accum(0.0); for(int i = 19; i >= 3; i -= 2) accum = ( x / (double)(i-1) ) * ( x / (double)i ) * ( T(1.0) + accum ); return x * ( T(1.0) + accum ); }; template <class T> T sin(const T& x) { T accum(0.0); for(int i = 19; i >= 3; i -= 2) accum = - ( x / (double)(i-1) ) * ( x / (double)i ) * ( T(1.0) + accum ); return x * ( T(1.0) + accum ); }; template <class T> T cos(const T& x) { T accum(0.0); for(int i = 18; i >= 2; i -= 2) accum = - ( x / (double)(i-1) ) * ( x / (double)i ) * ( T(1.0) + accum ); return accum + T(1.0); };
Комментариев нет:
Отправить комментарий