Nan ind c как исправить

Давайте вспомним математику и выполним небольшое преобразование:

Теперь нужно пройтись по массиву только один раз… Запишем код для double:

double Var(double* arr, size_t N)
{
    double x2 = 0, x = 0;
    for(size_t i = 0; i < N; ++i)
    {
        x2 += arr[i]*arr[i];
        x += arr[i];
    }
    return sqrt(N*x2-x*x)/N;
}

Теперь просто расширим его в шаблон для float и double с использованием концепта floating_point:

template<floating_point Double>
Double Var(Double* arr, size_t N)
{
    Double x2 = 0, x = 0;
    for(size_t i = 0; i < N; ++i)
    {
        x2 += arr[i]*arr[i];
        x += arr[i];
    }
    return sqrt(N*x2-x*x)/N;
}

Собственно, это всё…

Впрочем, можно еще один вариант, если передавать именно массив, то аргументом можно сделать ссылку на массив:

template<floating_point Double, size_t N>
Double Var(Double(&arr)[N])
{
    Double x2 = 0, x = 0;
    for(size_t i = 0; i < N; ++i)
    {
        x2 += arr[i]*arr[i];
        x += arr[i];
    }
    return sqrt(N*x2-x*x)/N;
}

и тогда при передаче массива не нужно передавать его размер…

P.S. Ах да, ваш код…

template <typename T>  // Ой. Сюда можно передать что угодно, например, int...
T Sigma_Array(const T* array, int count)
{
    T sum = 0;
    for (int i = 0; i < count; count++) // ЗАЧЕМ ВЫ УВЕЛИЧИВАЕТЕ COUNT?!!!
    {
        sum += array[i];
    }
    T average = sum / count;
    // Что вы делаете дальше?

    T square_s = 1. / count;  // Откуда это?

    for (int i = 0; i < count; count++) // ЗАЧЕМ ВЫ УВЕЛИЧИВАЕТЕ COUNT?!!!
    {
        T difference = array[i] - average;
        square_s *= pow(difference, 2); // Какое произведение?!
    }
    return sqrt(square_s);
}

Хотите в два прохода – так делайте это правильно:

template <typename T>
T Sigma_Array(const T* array, int count)
{
    T sum = 0;
    for (int i = 0; i < count; ++i) sum += array[i];
    T average = sum / count;
    T square_s = 0;
    for (int i = 0; i < count; ++i)
    {
        T difference = array[i] - average;
        square_s += difference*difference;
    }
    return sqrt(square_s/count);
}

By using complex math it is possible … The problem is that the pow is using ln operation which is on complex domain multi-valued which means it has infinite number of valid result (as imaginary part of result is added with k*2*PI where k is any integer). So when the used result is not added with correct k then the result will be complex and not the one you want.

However for rooting (1.0/exponent -> integer) the k can be obtained directly as this:

int k=int(floor((1.0/exponent)+0.5))/2;

Which I discovered just now empirically. When I rewrite my complex math cpow to use this I got this C++ code:

//---------------------------------------------------------------------------
vec2 cadd(vec2 a,vec2 b)    // a+b
    {
    return a+b;
    }
vec2 csub(vec2 a,vec2 b)    // a-b
    {
    return a-b;
    }
vec2 cmul(vec2 a,vec2 b)    // a*b
    {
    return vec2((a.x*b.x)-(a.y*b.y),(a.x*b.y)+(a.y*b.x));
    }
vec2 cdiv(vec2 a,vec2 b)    // a/b
    {
    float an=atan2(-a.y,-a.x)-atan2(-b.y,-b.x);
    float  r=length(a)/length(b);
    return r*vec2(cos(an),sin(an));
    }
vec2 csqr(vec2 a)           // a^2
    {
    return cmul(a,a);
    }
vec2 cexp(vec2 a)           // e^a
    {
    //  e^(x+y*i)= e^x * e^(y*i) = e^x * ( cos(y) + i*sin(y) )
    return exp(a.x)*vec2(cos(a.y),sin(a.y));
    }
vec2 cln(vec2 a)            // ln(a) + i*pi2*k where k={ ...-1,0,+1,... }
    {
    return vec2(log(length(a)),atan2(a.y,a.x));
    }
vec2 cpow(vec2 a,vec2 b)    // a^b
    {
    return cexp(cmul(cln(a),b));
    }
vec2 ctet(vec2 a,int b)     // a^^b
    {
    vec2 c=vec2(1.0,0.0);
    for (;b>0;b--) c=cpow(a,c);
    return c;
    }
//-------------------------------------------------------------------------
vec2 cln(vec2 a,int k)          // ln(a) + i*pi2*k where k={ ...-1,0,+1,... }
    {
    return vec2(log(length(a)),atan2(a.y,a.x)+float(k+k)*M_PI);
    }
float mypow(float a,float b)        // a^b
    {
    if (b<0.0) return 1.0/mypow(a,-b);
    int k=0;
    if ((a<0.0)&&(b<1.0))       // rooting with negative base
        {
        k=floor((1.0/b)+0.5);
        k/=2;
        }
    return cexp(cmul(cln(vec2(a,0.0),k),vec2(b,0.0))).x;
    }
//-------------------------------------------------------------------------

I am using GLSL like math vec2 you can easily rewrite to x,y components instead for example to this:

//-------------------------------------------------------------------------
float mypow(float a,float b)        // a^b
    {
    if (b<0.0) return 1.0/mypow(a,-b);
    int k=0;
    if ((a<0.0)&&(b<1.0))           // rooting with negative base
        {
        k=floor((1.0/b)+0.5);
        k/=2;
        }
    float x,y;
    x=log(fabs(a));
    y=atan2(0.0,a)+(float(k+k)*M_PI);
    x*=b; y*=b;
//  if (fabs(exp(x)*sin(y))>1e-6) throw domain error;  // abs imaginary part is not zero
    return exp(x)*cos(y);                              // real part of result
    }
//-------------------------------------------------------------------------

Which returns real part of complex domain pow and also select the correct multi-valued cln sub-result. I also added domain test as a comment in the code in case you need/want to implement it.

Here result for your case:

mypow(-12.411202,0.200000) = -1.654866

Have tested more number and 1/odd number exponents and looks like it is working.

[Edit1] handling mixed exponents

Now if we have exponents in form a^(b0+1/b1) where b0,b1 are integers and a<0 we need to dissect the pow to 2 parts:

a^(b0+1/b1) = a^b0 * a^(1/b1)

so when added to above function:

//-------------------------------------------------------------------------
float mypow(float a,float b)            // a^b
    {
    if (b<0.0) return 1.0/mypow(a,-b);  // handle negative exponents
    int k;
    float x0,x,y,e,b0;
    k=0;                                // for normal cases k=0
    b0=floor(b);                        // integer part of exponent
    b-=b0;                              // decimal part of exponent
    // integer exponent (real domain power |a|^b0 )
    x0=pow(fabs(a),b0);
    if ((a<0.0)&&((int(b0)&1))==1) x0=-x0; // just add sign if odd exponent and negative base
    // decimal exponent (complex domain rooting a^b )
    if ((a<0.0)&&(b>0.0))               // rooting with negative base
        {
        k=floor((1.0/b)+0.5);
        k&=0xFFFFFFFE;
        }
    x=b*(log(fabs(a))); e=exp(x);
    y=b*(atan2(0.0,a)+(float(k)*M_PI));
    x=e*cos(y);
//  y=e*sin(y); if (fabs(y)>1e-6) throw domain error;
    return x*x0; // full complex result is x0*(x+i*y)
    }
//-------------------------------------------------------------------------

and sample output:

mypow(-12.41120243,3.200000) = 3163.76635742
mypow(3163.76635742,0.312500) = 12.41120243