Эффективный способ вычисления математической константы e

Question

Стандартное представление константы e в виде суммы бесконечного ряда очень неэффективно для вычислений из-за множества операций деления.Так есть ли альтернативные способы эффективного вычисления константы?

Спасибо!

РЕДАКТИРОВАТЬ

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

Answer 1

Поскольку невозможно вычислить каждую цифру 'e', ​​вам нужно будет выбрать точку остановки.

двойная точность: 16 десятичных цифр

Для практических применений "64-битное значение с плавающей запятой двойной точности, максимально близкое к истинному значению 'e' - приблизительно 16 десятичных цифр" более чем адекватно.

As KennyTMсказал, что значение уже было предварительно рассчитано для вас в математической библиотеке.Если вы хотите рассчитать это сами, как отметил Ханс Пассант, факториал уже очень быстро растет.Первые 22 слагаемых в серии уже излишни для вычисления с такой точностью - добавление дополнительных слагаемых из ряда не изменит результат, если он будет сохранен в 64-битной переменной с плавающей запятой двойной точности.Я думаю, что мигание займет у вас больше времени, чем вашему компьютеру, чтобы сделать 22 деления.Поэтому я не вижу причин для дальнейшей оптимизации.

тысячи, миллионы или миллиарды десятичных цифр

Как отметил Матье М., это значение уже вычислено, и выможете загрузить его с веб-сайта Йи.

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

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

Одна страница очень кратко описывает алгоритмы, которые Йи использует для вычисления математических констант .

Статья Wikipedia "двоичное разбиение" содержит гораздо более подробные сведения.Я думаю, что часть, которую вы ищете, представляет собой представление чисел: вместо внутреннего хранения всех чисел в виде длинного ряда цифр до и после десятичной точки (или двоичной точки), Yee сохраняет каждый член и каждую частичную сумму как рациональное число.- как два целых числа, каждое из которых представляет собой длинный ряд цифр.Например, скажем, одному из рабочих процессоров была назначена частичная сумма:

... 1/4! + 1/5! + 1/6! + ... .

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

// extended to a million digits
1/24 + 1/120 + 1/720 => 0.0416666 + 0.0083333 + 0.00138888

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

// faster
1/24 + 1/120 + 1/720 => 1/24 + 840/86400 => 106560/2073600

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

Не забудьте избегать PrematureOptimization и всегда ProfileBeforeOptimizing .

Answer 2

Если вы используете double или float, в math.h уже есть константа M_E.

#define M_E         2.71828182845904523536028747135266250   /* e */

Есть и другие представления e в http://en.wikipedia.org/wiki/Representations_of_e#As_an_infinite_series;, все они будут связаны с делением.

Answer 3

Я не знаю каких-либо «более быстрых» вычислений, чем разложение в ряд Тейлора, т. Е.

e = 1/0!+ 1/1!+ 1/2!+ ...

или

1 / e = 1/0!- 1/1!+ 1/2!- 1/3!+ ...

Учитывая, что они были использованы А. Йи, который вычислил первые 500 миллиардов цифр e , я полагаю, что не так много оптимизаций длясделать (или лучше, это можно было бы оптимизировать, но никто не нашел пути, AFAIK)

РЕДАКТИРОВАТЬ

Очень грубая реализация

#include <iostream>
#include <iomanip>

using namespace std;

double gete(int nsteps)
{
  // Let's skip the first two terms
  double res = 2.0;
  double fact = 1;

  for (int i=2; i<nsteps; i++)
  {
    fact *= i;
    res += 1/fact;
  }

  return res;
}

int main()
{
  cout << setprecision(50) << gete(10) << endl;
  cout << setprecision(50) << gete(50) << endl;
}

Выходы

2.71828152557319224769116772222332656383514404296875
2.71828182845904553488480814849026501178741455078125

Answer 4

Эта страница содержит краткое изложение различных методов расчета.

Это крошечная C-программа от Xavier Gourdon, которая вычисляет 9000 десятичных цифр e на вашем компьютере. Программа такого же типа существует для π и некоторых других констант, определенных средним гипергеометрического ряда.

[версия разложена из https://codereview.stackexchange.com/a/33019]

#include <stdio.h>
int main() {
      int N = 9009, a[9009], x;
      for (int n = N - 1; n > 0; --n) {
          a[n] = 1;
      }
      a[1] = 2;
      while (N > 9) {
          int n = N--;
          while (--n) {
              a[n] = x % n;
              x = 10 * a[n-1] + x/n;
          }
          printf("%d", x);
      }
      return 0;
  }

Эта программа [при игре в гольф] имеет 117 символов. Его можно изменить, чтобы вычислить больше цифр (изменить значение 9009 на большее) и сделать его более быстрым (измените постоянную 10 на другую степень 10 и команду printf). Не столь очевидный вопрос - найти используемый алгоритм.

Answer 5

Я дал этот ответ в CodeReviews на вопрос , касающийся вычисления e по его определению через ряд Тейлора (поэтому другие методы не были возможны).Кросс-пост здесь был предложен в комментариях.Я удалил свои замечания, относящиеся к этой другой теме;Те, кто интересуется дальнейшими объяснениями migth, хотят проверить исходное сообщение.

---

Решение на C (должно быть достаточно простым для адаптации к C ++):

#include <stdio.h>
#include <math.h>

int main ()
{
    long double n = 0, f = 1;
    int i;
    for (i = 28; i >= 1; i--) {
        f *= i;  // f = 28*27*...*i = 28! / (i-1)!
        n += f;  // n = 28 + 28*27 + ... + 28! / (i-1)!
    }  // n = 28! * (1/0! + 1/1! + ... + 1/28!), f = 28!
    n /= f;
    printf("%.64llf\n", n);
    printf("%.64llf\n", expl(1));
    printf("%llg\n", n - expl(1));
    printf("%d\n", n == expl(1));
}

Вывод:

2.7182818284590452354281681079939403389289509505033493041992187500
2.7182818284590452354281681079939403389289509505033493041992187500
0
1

Есть два важных момента:

1. Этот код не вычисляет 1, 1 * 2, 1 * 2 * 3, ... которыйO (n ^ 2), но вычисляет 1 * 2 * 3 * ... за один проход (а это O (n)).

2. Он начинается с меньших чисел.Если бы мы попытались вычислить

   1/1 + 1/2 + 1/6 + ... + 1/20!

   и попытались добавить его 1/21 !, мы быдобавляем

   1/21!= 1/51090942171709440000 = 2E-20,

   до 2.что-то, что не влияет на результат (double содержит около 16 значащих цифр).Этот эффект называется underflow .

   Однако, когда мы начнем с этих чисел, т. Е. Если мы вычислим 1/32! +1/31! + ... все они имеют некоторое влияние.

Это решение выглядит в соответствии с тем, что C вычисляет с его функцией expl на моей 64-битной машине, скомпилированной с gcc 4.7.2 20120921.

Answer 6

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

e = 1 + 1/1! + 1/2! + 1/3! ...  

Расширяя уравнение:

e = 1 + 1/(1 * 1) + 1/(1 * 1 * 2) + 1/(1 * 2 * 3) ...

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

Answer 7

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

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

(e - 1) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)...))

Мы можем определить ряд функций f1 (x) ... fn(x) таким образом:

f1(x) = 1 + (1/2)x
f2(x) = 1 + (1/3)x
f3(x) = 1 + (1/4)x
...

Значение e находится из состава всех этих функций:

(e-1) = f1(f2(f3(...fn(x))))

Мы можем наблюдать, что значение x в каждой функцииопределяется следующей функцией, и что каждое из этих значений ограничено диапазоном [1,2], то есть для любой из этих функций значение x будет 1 <= x <= 2

, так как этоВ этом случае мы можем установить нижнюю и верхнюю границы для e, используя значения 1 и 2 для x соответственно:

lower(e-1) = f1(1) = 1 + (1/2)*1 = 3/2 = 1.5
upper(e-1) = f1(2) = 1 + (1/2)*2 = 2

Мы можем повысить точность, составив функции, определенные выше, и когда цифра maв нижней и верхней границах мы знаем, что наше вычисленное значение e точно соответствует этой цифре:

lower(e-1) = f1(f2(f3(1))) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)*1)) = 41/24 = 1.708333
upper(e-1) = f1(f2(f3(2))) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)*2)) = 7/4 = 1.75

Поскольку цифры 1 и 10 совпадают, можно сказать, что приближение (e-1) с точностью до десятых составляет 1,7.Когда первая цифра совпадает между верхней и нижней границами, мы вычитаем ее, а затем умножаем на 10 - таким образом, рассматриваемая цифра всегда находится на месте 1, где точность с плавающей запятой высока.

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

Другое объяснение метода: http://www.hulver.com/scoop/story/2004/7/22/153549/352

Бумага, котораяописывает алгоритм: http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf

Краткое введение в выполнение линейных преобразований с помощью матричной арифметики: https://people.math.gatech.edu/~cain/notes/cal6.pdf

Примечание: этот алгоритм использует преобразования Мёбиуса, которые представляют собой тип линейного преобразования, кратко описанногов бумаге Гиббонов.

Answer 8

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

3-sqrt(5/63)
2.7182819

Если вам нужно точное значение:

e = (-1)^(1/(j*pi))

, где j - мнимая единица, а piобщеизвестная математическая константа (тождество Эйлера)

Answer 9

С моей точки зрения, наиболее эффективный способ вычислить e с заданной точностью - это использовать следующее представление:

e := lim (n -> inf): (1 + (1/n))^n

Особенно, если вы выберете n = 2^xВы можете вычислить потенцию только с умножением x, так как:

a^n = (a^2)^(n/2), if n % 2 = 0

Answer 10

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

#include <iostream>
#include <iomanip>

template<int NUMER = 0, int DENOM = 1>
struct Rational
{
    enum {NUMERATOR = NUMER};
    enum {DENOMINATOR = DENOM};

    static double value;
};

template<int NUMER, int DENOM>
double Rational<NUMER, DENOM>::value = static_cast<double> (NUMER) / DENOM;

template<int ITERS, class APPROX = Rational<2, 1>, int I = 2>
struct CalcE
{
    typedef Rational<APPROX::NUMERATOR * I + 1, APPROX::DENOMINATOR * I> NewApprox;
    typedef typename CalcE<ITERS, NewApprox, I + 1>::Result Result;
};

template<int ITERS, class APPROX>
struct CalcE<ITERS, APPROX, ITERS>
{
    typedef APPROX Result;
};

int test (int argc, char* argv[])
{
    std::cout << std::setprecision (9);

    // ExpType is the type containing our approximation to e.
    typedef CalcE<13>::Result ExpType;

    // Call result() to produce the double value.
    std::cout << "e ~ " << ExpType::value << std::endl;

    return 0;
}

Другое (не метапрограммное) изменение во время компиляции вычислит двойное приближение e.У этого нет ограничения на число итераций.

#include <iostream>
#include <iomanip>

template<int ITERS, long long NUMERATOR = 2, long long DENOMINATOR = 1, int I = 2>
struct CalcE
{
    static double result ()
    {
        return CalcE<ITERS, NUMERATOR * I + 1, DENOMINATOR * I, I + 1>::result ();
    }
};

template<int ITERS, long long NUMERATOR, long long DENOMINATOR>
struct CalcE<ITERS, NUMERATOR, DENOMINATOR, ITERS>
{
    static double result ()
    {
        return (double)NUMERATOR / DENOMINATOR;
    }
};

int main (int argc, char* argv[])
{
    std::cout << std::setprecision (16);

    std::cout << "e ~ " <<  CalcE<16>::result () << std::endl;

    return 0;
}

В оптимизированной сборке выражение CalcE<16>::result () будет заменено фактическим двойным значением.

Обе, возможно, обадовольно эффективно, так как они вычисляют e во время компиляции: -)