C语言数组数值计算程序返回-nan(ind)错误排查

程序设计逻辑 #

开发目标为不依赖第三方额外库，实现指定长度数组生成及数值运算，逻辑如下：

• 首先生成长度为N的数组M1，之后生成长度为N/2的数组M2
• 对M1数组的每个元素执行除以Pi的操作，随后将结果取三次方
• 对M2数组，将每个元素依次与前一个元素相加，对相加结果执行取正切绝对值的操作
• 对M1和M2相同索引的元素执行幂运算，得到的结果数组使用侏儒排序（dwarf sorting）进行排序
• 最后计算M2数组中满足“元素除以M2数组最小非零元素结果为偶数”的元素的正弦值之和

问题现象 #

程序运行后输出结果为-nan(ind)，初步怀疑为变量数据类型使用错误，但未定位到根因，完整实现代码如下：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
 
const int A = 441;
const double PI = 3.1415926535897931159979635;
 
inline void dwarf_sort(double* array, int size) {
    size_t i = 1;
    while (i < size) {
        if (i == 0) {
            i = 1;
        }
        if (array[i - 1] <= array[i]) {
            ++i;
        }
        else
        {
            long tmp = array[i];
            array[i] = array[i - 1];
            array[i - 1] = tmp;
            --i;
        }
    }
}
 
inline double reduce(double* array, int size) {
 
    size_t i;
    double min = RAND_MAX, sum = 0;
 
    for (i = 0; i < size; ++i) {
        if (array[i] < min && array[i] != 0) {
            min = array[i];
        }
    }
 
    for (i = 0; i < size; ++i) {
        if ((int)(array[i] / min) % 2 == 0) {
            sum += sin(array[i]);
        }
    }
    return sum;
}
 
int main(int argc, char* argv[])
{
    int i, N, j;
    double* M1 = NULL, * M2 = NULL, * M2_copy = NULL;
    double X;
    unsigned int seed = 0;
 
    N = atoi(argv[1]); /* N равен первому параметру командной строки */
    M1 = malloc(N * sizeof(double));
    M2 = malloc(N / 2 * sizeof(double));
    M2_copy = malloc(N / 2 * sizeof(double));
 
    for (i = 0; i < 100; i++) 
    {
        seed = i;
        srand(i); 
        
        /*generate*/
        for (j = 0; j < N; ++j) {
            M1[j] = (rand_r(&seed) % A) + 1;
        }
 
        for (j = 0; j < N / 2; ++j) {
            M2[j] = (rand_r(&seed) % (10 * A)) + 1;
        }
 
        /*map*/
        for (j = 0; j < N; ++j)
        {
            M1[j] = pow(M1[j] / PI, 3);
        }
 
        for (j = 0; j < N / 2; ++j) {
            M2_copy[j] = M2[j];
        }
 
        M2[0] = fabs(tan(M2_copy[0]));
        for (j = 0; j < N / 2; ++j) {
            M2[j] = fabs(tan(M2[j] + M2_copy[j]));
        }
 
        /*merge*/
        for (j = 0; j < N / 2; ++j) {
            M2[j] = pow(M1[j], M2[j]);
        }
 
        /*sort*/
        dwarf_sort(M2, N / 2);
 
        /*sort*/
        X = reduce(M2, N / 2);
    }
    
    printf("\nN=%d.\n", N); 
    printf("X=%f\n", X);
    return 0;
}

根因定位 #

代码存在多处会触发nan的逻辑错误，按影响优先级排序：

1. 侏儒排序临时变量类型不匹配：交换double类型数组元素时，临时变量tmp被定义为long整型，赋值时会直接截断小数部分，造成数据损坏。修复方式：将long tmp修改为double tmp。
2. M2数组累加逻辑不符合需求：需求要求每个元素和前一个处理后的元素相加再取正切绝对值，但现有逻辑错误使用原始数组副本M2_copy参与计算，且循环从j=0开始，既没有实现逐元素累加的要求，还会导致正切计算的输入值无规律放大。同时tan函数在输入值接近π/2 + kπ时返回值会趋近于无穷大，后续将无穷大作为指数传入pow函数时，会直接生成inf或nan值，污染整个数组的运算结果。
3. reduce函数强转触发未定义行为：将浮点数运算结果array[i]/min强转为int时，若结果超过int类型的取值范围（幂运算很容易生成超过int上限的极大值），会直接触发未定义行为；同时直接用!=0判断浮点数、用整型取余判断偶数，对极大浮点数完全失效，一旦出现除零、非法值就会生成nan。修复方式：用fmod函数做浮点数取余，用精度阈值判断结果是否为偶数，浮点数判零也要增加精度阈值判断。
4. 边界逻辑缺失：未校验命令行参数argv[1]是否存在，未传入参数时直接调用atoi会触发段错误；申请的堆内存未在程序退出前释放，存在内存泄漏。

内容的提问来源于stack exchange，提问作者Dmitrii Shurygin