e^x 没有使用 math.h

Question

我试图在不使用math.h的情况下计算e^x。当x大于或小于约±20时，我的代码会给出错误的答案。我尝试将所有double类型更改为long double类型，但在输入上它会产生一些垃圾。

我的代码如下：

#include <stdio.h>

double fabs1(double x) {
    if(x >= 0){
        return x;
    } else {
        return x*(-1);
    }
}

double powerex(double x) {
    double a = 1.0, e = a;
    for (int n = 1; fabs1(a) > 0.001; ++n) {
        a = a * x / n;
        e += a;
    }
    return e;
}

int main(){
    freopen("input.txt", "r", stdin);
    freopen("output.txt", "w", stdout);
    int n;
    scanf("%d", &n);
    for(int i = 0; i<n; i++) {
        double number;
        scanf("%lf", &number);
        double e = powerex(number);
        printf("%0.15g\n", e);
    }
    return 0;
}

输入：

8
0.0
1.0
-1.0
2.0
-2.0
100.0
-100.0
0.189376476361643

我的输出：

1
2.71825396825397
0.367857142857143
7.38899470899471
0.135379188712522
2.68811714181613e+043
-2.91375564689153e+025
1.20849374134639

正确的输出：

1
2.71828182845905
0.367879441171442
7.38905609893065
0.135335283236613
2.68811714181614e+43
3.72007597602084e-44
1.20849583696666

你可以看到我的 e⁻¹⁰⁰ 的答案完全是错误的。为什么我的代码会输出这个？我可以做些什么来改进这个算法？

- user17221096

---

使用调试器运行，或在计算循环中添加中间打印语句，查看变量 a 的值何时开始偏离预期。 - Eugene Sh.

可能会出现int溢出吗？ - Alex Guteniev

还有一个有趣的问题：https://dev59.com/xWw05IYBdhLWcg3w3VeX - mkrieger1

在这里 for (int n = 1; fabs1(a) > 0.001; ++n) {，您强制近似通常不会比0.001更好。 - Damien

OP，请看一下大（正或负）参数的渐近展开。这可能在Abramowitz和Stegun的“数学函数手册”中有所涵盖，我认为你可以在网上找到它。 - Robert Dodier

Answer 1

当x为负数时，每一项的符号会交替变化。这意味着每个连续的求和在值上会发生较大的变化，而不是使用正指数时逐渐增加。这意味着随着连续项的损失，精度损失对结果有很大影响。

为了解决这个问题，在开始时检查x的符号。如果为负数，将x的符号切换以进行计算，然后当您到达循环结尾时反转结果。

此外，您可以通过使用以下反直觉条件来减少迭代次数：

e != e + a

表面上看来，这应该总是正确的。然而，当a的值超出了e值的精度时，条件会变为假，此时将a加到e中不会改变e的值。

double powerex(double x) {
    double a = 1.0, e = a;
    int invert = x<0;
    x = fabs1(x);
    for (int n = 1; e != e + a ; ++n) {
        a = a * x / n;
        e += a;
    }
    return invert ? 1/e : e;
}

我们可以进一步优化，通过将e初始化为0而不是a，并在循环的底部计算下一个项，从而减少一个循环迭代：

double powerex(double x) {
    double a = 1.0, e = 0;
    int invert = x<0;
    x = fabs1(x);
    for (int n = 1; e != e + a ; ++n) {
        e += a;
        a = a * x / n;
    }
    return invert ? 1/e : e;
}

- dbush

Answer 2

对于大于1的x值，建议分别处理整数部分，通过平方计算e的幂次。例如，e^9 = ((e²)²)²，只需要4次乘法。

实际上，泰勒级数展开的一般项 x^n/n! 只有在n>x时才开始减小（每次都乘以x/k），所以求和至少需要 x 个项 。另一方面，e^n 最多可以用2lg(n)次乘法来计算，这更有效率、更准确。

因此，建议：

- 取x的小数部分并使用泰勒级数展开， - 当整数部分为正数时，乘以 e 的该幂次， - 当整数部分为零时，已完成计算， - 当整数部分为负数时，除以 e 的该幂次。

甚至可以通过考虑四分之一来省略更多：在最坏情况下(x=1)，泰勒级数展开需要18个项才能使最后一个项变得不重要。如果考虑从x中减去下界最接近的四分之一的倍数（并通过预先计算的 e 的幂次进行补偿），项数将减少到12。

例如，e^0.8 = e^(3/4+0.05) = 2.1170000166126747 . e^0.05

- user1196549

---

如果您要添加常量，请考虑使用ln2，然后将其分解为x=m*ln2+u，其中abs(u)<=0.5*ln2<0.35，以便结果为2 ^ m * exp（u）。否则，想法是计算y=exp(x/2^n)并通过n次平方操作来恢复exp（x）。 - Lutz Lehmann

@LutzLehmann：由于截断误差，exp(x/2^n)的平方计算不起作用。 - user1196549

但是您建议从“e”的平方计算结果的大部分。 - Lutz Lehmann

@LutzLehmann：是的，如果我没记错的话，这是安全的，不会对像1.000000001这样的数字进行平方。 - user1196549

当然，这个想法是将 x/2^n 保持在一些适度的范围内，例如 0.5 到 1，以便幂级数或多项式逼近可以相对快速地工作，并且平方不会对精度产生太大影响。这种方法确实在 bc 多精度命令行工具的数学库中使用，可以通过提高工作精度来弥补平方“技巧”带来的精度损失。 - Lutz Lehmann