Разминка для мозгов: вычисление экспоненты через ряд Маклорена

Вычисление экспоненциальной функции с помощью разложения в ряд Маклорена

Для тех, кто не знает, что такое ряд Маклорена: понимания математического аппарата для решения задачи не требуется.

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

Известно, что экспоненциальная функция раскладывается в следующий ряд (см. по ссылке первую формулу).

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

Шаблон программы прилагается. Изменять шаблон программы нежелательно: ваш код должен быть только там, где указано комментарием (в начало программы можно добавить заголовочные файлы, если это совершенно необходимо). Для удобства отладки под Visual Studio в шаблон программы вставлены операторы условной компиляции: добавляется pause перед закрытием окна программы. Постарайтесь написать наиболее эффективный код.

Публикуйте, пожалуйста, только свой вариант функции ex.

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

#include <iostream>
#include <math.h>
#include <iomanip>
#ifdef _MSC_VER
#include <cstdlib>
#endif

using namespace std;

const double epsilon = 1e-10;

double ex(double x, double eps) {

    // здесь должен быть ваш код

}

int main() {

    double arg;
    cout << "x = ";
    cin >> arg;

    cout << "exp(" << arg << ")\n"
        << "Calculated using Maclaurin series: " 
        << setiosflags(ios::fixed) << setprecision(9) 
        << ex(arg, epsilon) << endl
        << "Calculated using library function: " 
        << exp(arg) << endl;

#ifdef _MSC_VER
    system("pause");
#endif
    return 0;
}

Дополнительная информация (если кто ещё не знает или уже забыл):

• Факториал
• Возведение в степень
• Модуль

Ну и чего все притихли? Сами же просили подкинуть какую-нибудь интересную задачку.

Если решать «в лоб», то задача совсем ПРОСТАЯ. Чуть сложнее, чем посчитать сумму чисел от 1 до 100 ))

А если немного подумать, то можно решить ещё проще )))

PS. Cranium — это по латыни череп.

Череп, Just to me I've guessed that cranium means Череп :) but now I haven't so much time to solve this problem, becouse protecting my diploma design. But I have one solution which I will try put at the blog tomorrow morning.

I will try put at the blog tomorrow morning.

Sorry for the mistake. I was be hurry :)

I will try to put at the blog tomorrow morning. It will be right.
By the way it's MY new authorized nick :)

Череп,

#ifdef _MSC_VER
#include <cstdlib>
#endif
// ...
#ifdef _MSC_VER
system("pause");
#endif

Я бы сделал проверку через #ifdef _WIN32. В Dev C++ тоже system("pause") может пригодиться.

А <cstdlib> можно и глобально включить, т.к. там есть другие полезные функции, которые могут быть использованы в решении.

Ну и чего все притихли? Сами же просили подкинуть какую-нибудь интересную задачку.

Череп, с моим школьным уровнем знаний алгебры( только 8 класс )? я задание-то не совсем ещё понял :)

selevit, у меня при прогоне программы под Dev C++ (Windows) IDE сама тормозит закрытие окна. A MS VS окно программы закрывает. Поэтому поставил именно такое условие. Вообще-то это не принципиально. Тот, кто знает зачем нужно system("pause"), сам внесёт в код такое изменение, с которым ему будет комфортно отладить программу ))

Кстати, а почему на сайте нельзя регистрироваться под русскими никами?

Юрий, точнее говоря, Alf :) а с какого бодуна ты вдруг на английский перешёл? Над дипломом перетрудился? ;-)

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

Попробую переформулировать задание:

Экспоненциальную функцию от x можно посчитать по формуле:

f(x) = 1 + x/1! + x^2/2! + x^3/3! + ... +x^n/n!

здесь ^ — операция возведения в степень,
n! — факториал от n (n! = 1 * 2 * 3 * ... * n).

Собственно говоря, нужно запрограммировать именно это выражение.

Поскольку ряд (см. выражение выше) бесконечный (n возрастает до бесконечности), вводится ограничение по точности вычисления: последнее вычисленное слагаемое (член ряда) должно быть по модулю меньше заданного предела точности (epsilon), который передаётся в подпрограмму вторым аргументом.

Так понятнее?

Первое, перехожу на русский :)))
Второе, как и обещал — свой вариант кода.

double ex(double x, double eps) {

    int i = 1, j = 2;
    long double fact[1000] = {1, 1};
    long double e = 1.0, temp = 0;

    do{
        e += temp;
        temp =(pow(x, i) / fact[i]);
        fact[j] = j * fact[i];
        ++i; ++j;
    }while (fabs(temp) > fabs(eps));


    return e;
}

А вот работа программы при разных значениях:

    x = -66
    exp(-66)
    Calculated using Maclaurin series: 1.#INF00000
    Calculated using library function: 0.000000000

    x = -60
    exp(-60)
    Calculated using Maclaurin series: -850380742.421332480
    Calculated using library function: 0.000000000

    x = -17
    exp(-17)
    Calculated using Maclaurin series: 0.000000042
    Calculated using library function: 0.000000041

    x = 1
    exp(1)
    Calculated using Maclaurin series: 2.718281828
    Calculated using library function: 2.718281828

    x = 15
    exp(15)
    Calculated using Maclaurin series: 3269017.372472111
    Calculated using library function: 3269017.372472111

    x = 17
    exp(17)
    Calculated using Maclaurin series: 24154952.753575306
    Calculated using library function: 24154952.753575299

    x = 66
    exp(66)
    Calculated using Maclaurin series: 1.#INF00000
    Calculated using library function: 46071866343312918000000000000.000000000

В общем диапазон вычислений, который по точности совпадает с библиотечной функцией, от -16 до 16. При дальнейшем увеличении или уменьшении значений появляется погрешность вычислений, при чем для положительных значений погрешность не велика (последний обрабатываемый аргумент 65), а при вычислении отрицательных — вычисление аргумента после -17 дает увеличение результата, а при значении -27 — результат со знаком минус.
Пока грешу на типы данных, которые я использовал в функции.
Я посмотрел на работу библиотечной функции и на данном этапе не могу даже предположить с каким типом данных она работает, так как она стабильно выдает результат, при заданной точности, в диапазоне аргумента от -22 до 709 (а это более 10^300, которым располагает long double).
Конечно может быть ошибка в коде? Но тогда как понимать точные вычисления от -16 до 16?
Череп, объясняй в чем соль... :)

Протестировал программу все работает. Погрешность дает факториал с типом данных long double, у которого сначала не хватает точности при вычислении в диапазонах от -20 до -65 и от 20 до 65, а потом и вовсе не хватает диапазона, так как значения выходят за 10^306.
С какими же типами данных работает функция exp()?

Так понятнее?

Это я и сразу понял) Я только сразу не понял что такое точность, хотя и щас толком не понял зачем она нужна) ( всё это время я пытался постичь суть всего этого, пока не понял, что это бессмысленное занятие) )

Вот мой код:

double ex(double x, double eps)
{
    double ret = 1, t = 1, f = 1;
    int k;
    for ( int n = 1; t > eps; n++ )
    {
        t = pow( x, n ) / f;
        ret += t;
        f = k = n+1;
        while ( k > 1 )
              f *= --k;           
    }
    return ret;
}

Не-е-ет, брат porshe, не пойдет. Если ты протестируешь свой код с отрицательными значениями, то увидишь, что код не рабочий:)
В проверочном условии проверку нужно делать по модулю.

Вот еще один вариант. Есть пару моментов — я думаю, что тип double здесь маловат, поэтому я заменил его на long double и добавил в функцию main() цикл для удобства сравнения значений.

#include <iostream>
#include <math.h>
#include <iomanip>

using namespace std;

const double epsilon = 1e-10;

long double ex(long double x, double eps)
{
  int i = 1;
  long double e = 1.0, temp = 0, fact = 1.0;

  while (fabs(temp = pow(x, i) / (fact *= i)) > fabs(eps)){
      e += temp;
      ++i;
  }
  return e;
}

int main() {

    long double arg;

    ios_base::fmtflags old = cout.setf(ios_base::fixed, ios_base::floatfield);

    cout << "x = ";

    while (cin >> arg){
        cout.setf(old, ios_base::floatfield);
        cout << "exp(" << arg << ")\n"
            << "Calculated using Maclaurin series: " 
            << setiosflags(ios::fixed) << setprecision(9) 
            << ex(arg, epsilon) << endl
            << "Calculated using library function: " 
            << exp(arg) << endl << endl;
        cout << "x = ";
    }
}

Или так

long double ex(long double x, double eps)
{
  long double e = 1.0, temp = 0, fact = 1.0;

  for (int i = 1; fabs(temp = pow(x, i) / (fact *= i)) > fabs(eps); ++i)
      e += temp;
  return e;
}

Для начала мой вариант функции:

double ex(double x, double eps) {

    long double sum = 1.0;
    long double cur = 1.0;
    unsigned int n = 1;

    while (fabs(cur) > eps) {
        sum += (cur *= ((long double)x) / n++);
    }

    return sum;
}

Весь цимес задачи был в том, что бы заметить, что следующий член ряда получается путём умножения предыдущего члена на x/n. Это даёт следующие преимущества:

1. Скорость. Это очевидно. Сразу избавляемся от вызова функции pow() и, не дай Бог, от вычисления факториала на каждой итерации.
2. На значительно большем диапазоне входных значений отсутствует переполнение. Функции f(n) =x^n и f(n) = n! очень быстрорастущие (здесь x — константа!). При их вычислении возможности разрядной сетки исчерпываются очень быстро исчерпываются. Функция f(n) = x/n — убывающая, асимптотически приближающаяся к нулю. В данном случае это очень хорошо. Этот вариант функции считает до x = 709 включительно (так же, как и библиотечная функция).
3. Точность вычисления. При вычислении больших чисел, например функцией pow(), теряются значащие цифры: в MS VC++ из 308 возможных десятичных цифр (диапазон double 1.7E +/- 308) работают только 15. Выражение x/n считается с точностью 13 знаков минимум. Отсюда абсолютное совпадение вычисленных значений с библиотечной функцией в диапазоне -5..709.

Работал с Dev-C++ 5.5.3/TDM-GCC 4.7.1 64-bit. Это важно. Потому что с MS VC++ 2012 результаты другие.

В MS VC++ 2012 тип long double совпадает с double: длина 8 байт, диапазон 1.7E +/- 308, точность 15 десятичных знаков.

В TDM-GCC 4.7.1 64-bit тип double совпадает с майкрософтовским. А long double имеет длину 16 байт, из которых используются 10 (80 бит, x86 Extended Precision Format), диапазон 1.18973E +/- 4932, точность 19 десятичных знаков.

Понятно, что по точности результаты MS VC++ сильно уступают результатам TDM-GCC. (Вообще-то я был неприятно удивлён таким фокусом от MS.)

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

Череп, так у тебя функция тоже не работает со значениями ниже -25 ???

А оно и не должно считать. Точность задана 1e-10, а значение функции exp(-25) = 1.3887943864964E-11, т.е. меньше заданной точности.

Ради интереса я немножко поэкспериментировал: использовал библиотеку для вычислений с произвольной точностью. У меня получилось перекрыть весь диапазон библиотечной функции: от -708 до 709... при epsilon = 1e-1200 и использовании 3000 значащих бит (против 64 для double и 80 для long double). Это, так сказать, решение «в лоб»: путём тупого увеличения количества бит и уменьшения порога точности.

Видимо библиотечная функция использует другой алгоритм ))

Но это всё выходит далеко за рамки изначальной задачи.

Череп, у меня есть два вопроса:
1. В консоли если ввести число типа double с запятой вместо точки, в качестве разделителя, то следующий ввод будет с ошибкой, так как запятая и цифры идущие за ней остаются в потоке — в чем проблема и как с ней «бороться»?
2. Можно ли на одну машину установить два компилятора? Не «подерутся» ли они?

1. В деталях я не помню, но десятичный разделитель в консольном вводе-выводе настраивается. Почитай описание setlocale() и, если не поможет, консольных функций.
2. Можно, если осторожно. Опасность есть только в отладчиках: перехват на отладку упавших программ, удалённая отладка и пр. Надо аккуратно настраивать. У меня на одной машине стоят MS VS 2012, Dev-C++ и Digital Mars — всё отлично уживается.