Ошибка nan ind

Потому что вы делите на 0.

Запишите

assert(a[k][i] != 0);

перед

float s = a[k + 1][i] / a[k][i];

и убедитесь в этом сами…

Что и неудивительно — ведь матрица

10 23 17  44
15 35 26  69
25 57 42 108
30 69 51 133

вырожденная, ее определитель равен 0…

P.S. Что-то странный у вас вложенный цикл

for (int i = 0; i < n-1; i++) {
    for (int i = 0; i < n - 1; i++) {

нет? i во вложенном цикле скрывает i из внешнего, да и дальнейшее вызывает сомнения… int z = a[i][y]; тоже как-то не комильфо — зачем double до int обрезать?
Но главное — какой смысл в проверках до преобразований? Ведь именно в их процессе очередной опорный элемент может стать нулевым! Возьмите простую матрицу, типа

1  1
1  1

и вручную убедитесь, что первая ваша проверка работает, но уже после первой итерации преобразований вы получите нулевую вторую строку — и никаких проверок такой неприятности у вас уже нет.

Написал функцию поворота прямой.

// Line - y = _k * x + _b
void Line::rotation(double angle) {
    double temp = _k;
    _k = (_k + tan(angle))/(1 - _k*tan(angle));
    _b *= _k/temp;
}

double Line::rotationAngle(const Line &l) { // counterclock-wise from l to this
    return atan((_k - l._k)/(1. + _k*l._k));
}

При вводе следующих данных

Line Line1(0, 0), Line2(2, 1);

    cout << Line1;
    cout << Line2;

    Line1.rotation(Line2.rotationAngle(Line1));

    cout << Line1;

Почему-то выдает следующее

y = 0x + 0
y = 2x + 1
y = 2x + -nan(ind)

Коэффициент повернутой прямой программа посчитала правильно (_k = 2), но почему-то потом не может поделить 2 на 2:
_b *= _k/temp;
Хотя если temp заменить 2, то все ок. Выдаст: y = 2x + 0
Вроде бы никаких корней и делений на ноль. В чем причина?

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