计算级数求和S时出现-nan(ind)错误，请求排查问题原因

问题背景 #

在区间[0, π/2]内以π/20为步长，依据等式：
$$\cos^2x = 1 - \sum_{k=1}^{\infty} \frac{(-1)^{k+1} \cdot 2^{2k-1} \cdot x^{2k}}{(2k)!}$$
计算S时，x取π/20及后续值均返回-nan(ind)，相关代码与输出如下：

相关C代码 #

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

#define pi       3.1415926535f
#define start    0
#define end      pi / 2
#define delta    pi / 20
#define accuracy powf(10, -6)

float factorial(float number)
{
    if (number == 0)
    {
        return 1;
    }
    int c = 0;
    float factorial = 1;
    for (c = 1; c <= number; c++)
    {
        factorial *= c;
    }
    return(factorial);
}

int main()
{
    setlocale(LC_ALL, "Russian");
    int k = 0;
    float x = 0, y = 0, s = 0;
    printf("1 - Аргумент (x)\n2 - Значение функции y(x)\n3 - Сумма (s)\n4 - Количество итераций (k)\n\n+-------------+-------------+-------------+-------------+\n|      1      |      2      |      3      |      4      |\n+-------------+-------------+-------------+-------------+\n");
    for (x = start; x <= end; x += delta)
    {
        do
        {
            k++;
            s += powf(-1, k + 1) * ((powf(2, 2 * k - 1) * powf(x, 2 * k)) / (factorial(2 * k)));
        } while (fabs(s) >= accuracy);
        s = 1 - s;
        y = powf(cos(x), 2);
        printf("| %11.6f | %11.6f | %11.6f | %11i |\n+-------------+-------------+-------------+-------------+\n", x, y, s, k);
        s = 0;
        k = 1;
    }
    return 0;
}

控制台输出 #

1 - Аргумент (x)
2 - Значение функции y(x)
3 - Сумма (s)
4 - Количество итераций (k)

+-------------+-------------+-------------+-------------+
|      1      |      2      |      3      |      4      |
+-------------+-------------+-------------+-------------+
|    0,000000 |    1,000000 |    1,000000 |           1 |
+-------------+-------------+-------------+-------------+
|    0,157080 |    0,975528 |   -nan(ind) |          65 |
+-------------+-------------+-------------+-------------+
|    0,314159 |    0,904509 |   -nan(ind) |          65 |
+-------------+-------------+-------------+-------------+
|    0,471239 |    0,793893 |   -nan(ind) |          65 |
+-------------+-------------+-------------+-------------+
|    0,628319 |    0,654509 |   -nan(ind) |          65 |
+-------------+-------------+-------------+-------------+
|    0,785398 |    0,500000 |   -nan(ind) |          65 |
+-------------+-------------+-------------+-------------+
|    0,942478 |    0,345491 |   -nan(ind) |          65 |
+-------------+-------------+-------------+-------------+
|    1,099557 |    0,206107 |   -nan(ind) |          57 |
+-------------+-------------+-------------+-------------+
|    1,256637 |    0,095491 |   -nan(ind) |          49 |
+-------------+-------------+-------------+-------------+
|    1,413717 |    0,024472 |   -nan(ind) |          44 |
+-------------+-------------+-------------+

问题原因分析 #

• 阶乘计算溢出
  factorial函数使用float类型存储结果，而float的最大值约为3.4×10³⁸。当2k达到34时，34!≈2.95×10³⁸，已接近float上限；2k超过34后，阶乘结果溢出为inf（无穷大）。后续计算中出现inf/inf这类非法运算，最终生成NaN。

• 循环终止条件错误
  当前循环用累加和s的绝对值与精度比较（fabs(s) >= accuracy），这会导致循环无法及时终止：随着累加s逐渐收敛，该条件会让循环持续迭代到阶乘溢出，触发NaN。正确的终止条件应该是判断当前迭代项的绝对值小于精度要求，而非累加和。

• 数值精度与迭代逻辑缺陷
  使用float进行高精度级数计算本身就容易出现精度丢失，加上每次迭代重复计算高次幂和阶乘，进一步加剧了数值溢出的风险。

修复建议 #

1. 改用double类型：替换所有float为double，提升数值范围与精度，延迟阶乘溢出的发生；
2. 修正终止条件：计算当前项的绝对值，当该项小于accuracy时终止循环；
3. 递推计算项值：利用前一项推导当前项，避免重复计算幂和阶乘，比如：

   // 初始项（k=1时）
   double term = (2 * x * x) / 2;
   s += term;
   k = 2;
   do {
       // 递推计算下一项
       term = term * (-1) * 4 * x * x / ((2*k-1)*(2*k));
       s += term;
       k++;
   } while (fabs(term) >= accuracy);
   

   这种方式无需单独计算阶乘和高次幂，效率更高且不易溢出。

内容的提问来源于stack exchange，提问作者Ivanov Vladislav