четверг, 22 октября 2020 г.

Понижение сложности экспоненциальных рядов

Часто используемые элементарные функции вида экспоненты, синусов и косинусов, имеют ряд разложения, в котором члены немного похожи. Если вычислять ряд в лоб, то вычислительная сложность относительно умножения элементов алгебры равна n2. Можно ли её понизить? Попробуем разобраться.

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

Рассмотрим разложение в ряд функции экспоненты: $$ e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\ldots $$ В этом ряду каждый последующий элемент содержит в качестве своей части предыдущие. Чтобы было видно лучше, перепишем чуть иначе: $$ 1+\frac{x}{1}+\frac{x}{1}\frac{x}{2}+\frac{x}{1}\frac{x}{2}\frac{x}{3}+ \frac{x}{1}\frac{x}{2}\frac{x}{3}\frac{x}{4}+ \frac{x}{1}\frac{x}{2}\frac{x}{3}\frac{x}{4}\frac{x}{5}+\ldots $$ В действительности мы не вычисляем бесконечные ряды. Если бы так было, то никакие вычисления не закончились бы. Конечно, мы ограничиваем выисление некоторым членом. Проведем предварительные оценивания, при каком числе членов ряда $n$ достигается разумная точность. В моем эксперименте это оказалось $n=20$. Увеличение $n$ уже не давало существенного увеличения точности, но давало рост времени вычисления.

Ограничим ряд некоторым числом. Правый хвост ряда будет выглядеть в случае ограничения на $n=5$: $$ \frac{x}{1}\frac{x}{2}\frac{x}{3}\frac{x}{4}+ \frac{x}{1}\frac{x}{2}\frac{x}{3}\frac{x}{4}\frac{x}{5} $$ Этот хвост можно сократить до: $$ \frac{x}{1}\frac{x}{2}\frac{x}{3}\frac{x}{4}\left(1+\frac{x}{5}\right) $$ Приложим к этой части хвоста следующую от нее слева часть: $$ \frac{x}{1}\frac{x}{2}\frac{x}{3}+ \frac{x}{1}\frac{x}{2}\frac{x}{3}\frac{x}{4}\left(1+\frac{x}{5}\right) $$ Эту часть хвоста также можно сократить: $$ \frac{x}{1}\frac{x}{2}\frac{x}{3}\left(1+ \frac{x}{4}\left(1+\frac{x}{5}\right)\right) $$ Соответственно, рекурсивно продолжая до самого начала получим итеративное вычисление со сложностью по умножению порядка $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);
};
Таким образом, действительно удалось понизить сложность вычисления экспоненциального ряда с $n^2$ до $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);
};

Комментариев нет:

Отправить комментарий