visual studio 调试Expression: (“unexpected input value； log10 failed”, 0)

#include "stdio.h"
#include "math.h"
void my_dft(double* x, double* y, 
    double* a, double* b, int n, int sign)
{
    int i, k;
    double c, d, q, w, s;
    q = 6.28318530718 / (double)n;
    for (k = 0; k < n; k++) {
        w = k * q;
        a[k] = b[k] = 0.0;
        for (i = 0; i < n; i++) {
            d = i * w;
            c = cos(d);
            s = sin(d) * (double)sign;
            a[k] += c * x[i] + s * y[i];
            b[k] += c * y[i] - s * x[i];
            
        }
    }

    if (sign == -1) {
        c = 1.0 / n;
        for (k = 0; k < n; k++) {
            a[k] = c * a[k];
            b[k] = c * b[k];
        }
        
    }
}

static void print_complex(double* a, double* b, int n)
{
    for (int i = 0; i < n / 2; i++) {
        for (int j = 0; j < 2; j++) {
            printf("%10.7f + J %10.7f        ", a[2 * i + j], b[2 * i + j]);
        }
        printf("\n");
    }
}

/*
    设输入序列 x(i)为 x(i) = Q^i, i=0，l，,..., n-1
    其离散傅立叶变换为，X(k)=(1 - Q^n)/(1 - Q*W^k), k=0，l，,..., n-1
    这里 W = e^(-j2pi/n)
    选取Q = 0.9 + j0.3,n=32,
    计算x(i)的离散傅里叶变换X(k)和X(k)的离散傅里叶反变换x(i)，并与理论值进行比较
*/
int my_dft_test(void)
{
    int i,j,n;
    double a1,a2,c,c1,c2,d1,d2,q1,q2,w,w1,w2;
    double x[32] = { 0 }, y[32] = { 0 };
    double a[32] = { 0 }, b[32] = { 0 };

    n = 32;
    a1 = 0.9;
    a2 = 0.3;
    x[0] = 1.0;
    y[0] = 0.0;
    for (i = 1; i < n; i++) {
        x[i] = a1 * x[i - 1] - a2 * y[i - 1];
        y[i] = a2 * x[i - 1] + a1 * y[i - 1];
    }
    
    printf("\n Original Sequence\n");
    print_complex(x, y, n);
    q1 = x[n - 1];
    q2 = y[n - 1];
    my_dft(x, y, a, b, n, 1);

    printf("\n Discrete Fourier Transform\n");
    print_complex(a, b, n);

    for (i = 0; i < n; i++) {
        w = 6.28318530718 / n * i;
        w1 = cos(w);
        w2 = -sin(w);
        c1 = 1.0 - a1 * w1 + a2 * w2;
        c2 = a1 * w2 + a2 * w1;
        c = c1 * c1 + c2 * c2;
        d1 = 1.0 - a1 * q1 + a2 * q2;
        d2 = a1 * q2 + a2 * q1;
        x[i] = (c1 * d1 + c2 * d2) / c;
        y[i] = (c2 * d1 - c1 * d2) / c;
    }
    printf("\n Theoretical Discrete Fourier Transform\n") ;
    print_complex(x, y, n);

    my_dft(a, b, x, y, n, -1);
    printf("\n Inverse Discrete Fourier Transform\n");
    print_complex(x, y, n);
}

int main()
{
    printf("Hello DSP!\n");

    my_dft_test();

    return 0;
}

环境Visual studio 2019，运行上面程序出现错误 Expression: (“unexpected input value; log10 failed”, 0)

进入cfout.cpp里面查看错误来源：

https://github.com/ojdkbuild/tools_toolchain_sdk10_1607/blob/master/Source/10.0.14393.0/ucrt/convert/cfout.cpp

    // k is the decimal exponent, such that the resulting decimal number is of
    // the form 0.mmmmm * 10^k, where mmmm are the mantissa digits we generate.
    // Note that the use of log10 here may not give the correct result:  it may
    // be off by one, e.g. as is the case for powers of ten (log10(100) is two,
    // but the correct value of k for 100 is 3).  We detect off-by-one errors
    // later and correct for them.
    int32_t k = static_cast<int32_t>(ceil(log10(value)));
    if (k == INT32_MAX || k == INT32_MIN)
    {
        _ASSERTE(("unexpected input value; log10 failed", 0));
        k = 0;
    }

原因大概猜测是，printf函数打印浮点数时，进入了整数的浮点优化策略，导致这个assert出现。

 

解决方式：引用头文件math.h，引出cos函数和sin函数，这里会给出double cos(double)的函数，不然的话，visual studio按照默认int cos(int)的方式运算，导致算出来的cos(0) = 0.