고유값, 고유벡터 구하기

내가 이걸 왜 만들었지.
아무튼, 복소수로 3차 정사각형 행렬의 eigenvalue와 eigenvector 찾는 프로그램 소스다. 사용법은 각자 공부해 보시길. 정 모르겠으면 댓글로 문의바람. 임의의 n차 행렬이었으면 어딘가에 올렸을지도 모르지만, 그냥 3차 행렬로 고정되어 있음. 포인터같은거 전혀 없이 전부 변수로만 넘겨줌. 언젠가 포인터로 변수 넘겨주는 프로그램으로 바꾸고 싶음.
혹시 버그 발견하게 되면 알려주시기 바람.
산수 공부하는데 도움이 될까 하여, 연구와 관련된 부분을 빼고 올림. 빼는 과정에서 뭔가 이상해질수도 있음.
gcc로 컴파일 할 때는 -lm(link math)옵션이 필요함.
라이센스는 GPL이며, 특별히 문의해도 좋음.
*최적화같은거 전혀 신경쓰지 않았음.

//This source is under GNU Public License.
//Copyright @ Kee-Hwan Nam, all right reserved.
//You can use this source under GPL, or after asking to me.
//email : snowall a_t gmail d_o_t com (a_t = @ / d_o_t = . )

#include<math.h>
#define PI (3.141592653589793)
#define TRUE 1 //True
#define FALSE 0 //False
#define DEBUG

//the type definition of complex number
typedef struct {
    double x;
    double y;
} cnumber;

//the type definition of 3x3 matrix with complex elements
typedef struct {
    cnumber a[3][3];
} matrix;

//the type definition of coefficient for dealing the characteristic equation of the 3x3 matrix
typedef struct {
    cnumber a;
    cnumber b;
    cnumber c;
    cnumber d;
} coeff;

//the type definition of roots of cubic equations
typedef struct {
    cnumber a[3];
//    cnumber b;
//    cnumber c;
} cbeqrt;

//the type definitions of vectors
typedef struct {
    cnumber x;
    cnumber y;
    cnumber z;
} vector;


//function declarifications
cbeqrt eigenvalue_find(matrix mix);
cbeqrt find_cubic(coeff x);
cnumber add_cnum(cnumber z, cnumber w);
cnumber arccos_cnum(cnumber y);
cnumber cbrt_cnum(cnumber z);
cnumber conj_cnum(cnumber z);
double alpha, double beta, int evnumber);
cnumber div_cnum(cnumber z, cnumber w);
cnumber dotproduct(vector a, vector b);
cnumber log_cnum(cnumber y);
cnumber minus_cnum(cnumber z);
cnumber mul_cnum(cnumber z, cnumber w);
cnumber pow_cnum(cnumber z,double k);
cnumber sqrt_cnum(cnumber z);
cnumber sub_cnum(cnumber z, cnumber w);
cnumber det_matrix(matrix m);
coeff character_matrix(matrix m);
double absolute(double a);
double arg_cnum(cnumber z);
double norm_cnum(cnumber z);
double norm_vec(vector a);
int equality(cnumber a,cnumber b);
matrix echelon(matrix);
matrix eigenvectoreq(matrix m,cbeqrt z,int k);
vector add_vec(vector a, vector b);
vector col_vector(matrix m,int i);
vector eigenvector_find(matrix mix,int i,int j);
vector equationsystem(matrix m,int k);
vector normalize(vector);
vector opposite_vec(vector a);
vector row_vector(matrix m,int i);
vector scalar_mul_vec(cnumber a, vector b);
vector sub_vec(vector a, vector b);
void printing_cnum(cnumber c);
void printing_matrix(matrix);
void printing_vec(vector v);


//The absolute value of the given real number a
double absolute(double a){
    if(a>0){
        return a;
    }
    else if(a<0){
        return -a;
    }
    else if(a=0){
        return a;
    }
}

//Following definitions are the functions for dealing complex numbers, complex vectors, and complex matrices.
//The addition of two complex numbers z and w
cnumber add_cnum(cnumber z, cnumber w){
    cnumber a={z.x+w.x,z.y+w.y};
#ifdef DEBUG
    printf("\nadd_cnum\n");printing_cnum(z);printing_cnum(w);printing_cnum(a);printf("\n");
#endif
    return a;
}

//The negative number of the given complex number z
cnumber minus_cnum(cnumber z){
    z.x=-z.x;
    z.y=-z.y;
#ifdef DEBUG
    printf("\nminus_cnum\n");printing_cnum(z);printf("\n");
#endif
    return z;
}

//Subtracting w from z
cnumber sub_cnum(cnumber z, cnumber w){
#ifdef DEBUG
    printf("\nsub_cnum\n");printing_cnum(z);printing_cnum(w);printf("\n");
#endif
    return add_cnum(z,minus_cnum(w));
}

//Multiplication of z by w
cnumber mul_cnum(cnumber z, cnumber w){
    cnumber a;
    a.x=z.x*w.x-z.y*w.y;
    a.y=z.x*w.y+z.y*w.x;
#ifdef DEBUG
    printf("\nmul_cnum\n");printing_cnum(z);printing_cnum(w);printing_cnum(a);printf("\n");
#endif
    return a;
}

//Complex conjugate of the given complex number z
cnumber conj_cnum(cnumber z){
    cnumber a;
    a.x=z.x;
    a.y=-z.y;
#ifdef DEBUG
    printf("\nconj_cnum\n");printing_cnum(z);printing_cnum(a);printf("\n");
#endif
    return a;
}

//The norm of the given complex number z
double norm_cnum(cnumber z){
    double a;
    a=sqrt(z.x*z.x+z.y*z.y);
#ifdef DEBUG
    printf("\nnorm_cnum\n");printing_cnum(z);printf("%lf",a);printf("\n");
#endif
    return a;
}

//z is divided by w
cnumber div_cnum(cnumber z, cnumber w){
    cnumber a;
    a=mul_cnum(z,conj_cnum(w));
    a.x=a.x/(norm_cnum(w)*norm_cnum(w));
    a.y=a.y/(norm_cnum(w)*norm_cnum(w));
#ifdef DEBUG
    printf("\ndiv_cnum\n");printing_cnum(z);printing_cnum(w);printing_cnum(a);printf("\n");
#endif
    return a;
}

/*
//Real part of the given complex number z
double re_cnum(cnumber z){
    double a;
    a=z.x;
    return a;
}

//Imaginary part of the given complex number z
double im_cnum(cnumber z){
    double a;
    a=z.y;
    return a;
}


//complex phase transformation of the given complex number z by k
cnumber phase_cnum(cnumber z, double k){
    cnumber x={cos(k),sin(k)};
    cnumber m=mul_cnum(z,x);
    return m;
}
*/

//The argument of the given complex number z
double arg_cnum(cnumber z){
    double x=atan2(z.y,z.x);
    return x;

}

cnumber log_cnum(cnumber y){
    cnumber a;
    a.x=log(norm_cnum(y));
    a.y=arg_cnum(y);
#ifdef DEBUG
    printf("\nlog_cnum\n");printing_cnum(y);printing_cnum(a);printf("%lflf,%lf\n",a.x,a.y);
#endif
    return a;
}

cnumber arccos_cnum(cnumber y){
    cnumber a;
#ifdef DEBUG
    printf("\narccos_cnum_start\n");printing_cnum(y);
#endif
    cnumber minus_image={0.0,-1.0};
    cnumber one={1.0,0.0};
    a=mul_cnum(minus_image,log_cnum(add_cnum(y,sqrt_cnum(sub_cnum(mul_cnum(y,y),one)))));
#ifdef DEBUG
    printf("\narccos_cnum\n");printing_cnum(y);printing_cnum(a);printf("\n");
#endif
    return a;
}

//complex square-root
cnumber sqrt_cnum(cnumber z){
    cnumber a;
    a.x=sqrt(norm_cnum(z))*cos(arg_cnum(z)/2.0);
    a.y=sqrt(norm_cnum(z))*sin(arg_cnum(z)/2.0);
#ifdef DEBUG
    printf("\nsqrt_cnum\n");printing_cnum(z);printing_cnum(a);printf("\n");
#endif
    return a;
}

cnumber double_sqrt_cnum(cnumber z,int i){
    while(i--){
    z=sqrt_cnum(z);
    }
    return z;
}    
//complex cubic-root
cnumber cbrt_cnum(cnumber z){

#ifdef DEBUG
    printf("\ncbrt_cnum\n");printing_cnum(z);
#endif
    cnumber w;
    z=sqrt_cnum(sqrt_cnum(z));
    z=mul_cnum(z,double_sqrt_cnum(z,2));
    z=mul_cnum(z,double_sqrt_cnum(z,4));
    z=mul_cnum(z,double_sqrt_cnum(z,8));
    z=mul_cnum(z,double_sqrt_cnum(z,16));
    z=mul_cnum(z,double_sqrt_cnum(z,32));
    z=mul_cnum(z,double_sqrt_cnum(z,64));
#ifdef DEBUG
    printf("\ncbrt_cnum\n");printing_cnum(z);
#endif
    return z;
}
/*
//complex power function
cnumber pow_cnum(cnumber z,double k){
    cnumber a;
    a.x=pow(norm_cnum(z),k)*cos(arg_cnum(z)*k);
    a.y=pow(norm_cnum(z),k)*sin(arg_cnum(z)*k);
    return a;
}
*/
/*
//equality
int equality(cnumber a,cnumber b){
    if((a.x==b.x)||(a.y==b.y))
        return TRUE;
    else
        return FALSE;
}
*/
//Complex dot product : (a^*).b
cnumber dotproduct(vector a, vector b){
    cnumber p;
    p=add_cnum(mul_cnum(conj_cnum(a.x),b.x),add_cnum(mul_cnum(conj_cnum(a.y),b.y),mul_cnum(conj_cnum(a.z),b.z)));
#ifdef DEBUG
    printf("\ndotproduct\n");printing_vec(a);printing_vec(b);printing_cnum(p);
#endif
    return p;
}

//The addition of two vectors a and b
vector add_vec(vector a, vector b){
    vector c;
    c.x=add_cnum(a.x,b.x);
    c.y=add_cnum(a.y,b.y);
    c.z=add_cnum(a.z,b.z);
#ifdef DEBUG
    printf("\nadd_vec\n");printing_vec(a);printing_vec(b);printing_vec(c);
#endif    
    return c;

}

//The opposite vector of the given vector a
vector opposite_vec(vector a){
    a.x=minus_cnum(a.x);
    a.y=minus_cnum(a.y);
    a.z=minus_cnum(a.z);
#ifdef DEBUG
    printf("\nopposite_vec\n");printing_vec(a);
#endif
    return a;
}

//The vector subtraction of b from a
vector sub_vec(vector a, vector b){
    vector c;
    c=add_vec(a,opposite_vec(b));
#ifdef DEBUG
    printf("\nsub_vec\n");printing_vec(a);printing_vec(b);printing_vec(c);
#endif
    return c;
}

//The vector scalar multiplication of a by b
vector scalar_mul_vec(cnumber a, vector b){
    vector c;
    c.x=mul_cnum(a,b.x);
    c.y=mul_cnum(a,b.y);
    c.z=mul_cnum(a,b.z);
#ifdef DEBUG
    printf("\nscalar_vec\n");printing_cnum(a);printing_vec(b);printing_vec(c);
#endif
    return c;
}

//The norm of the given vector a
double norm_vec(vector a){
    double b;
    b=sqrt(norm_cnum(dotproduct(a,a)));
#ifdef DEBUG
    printf("\nnorm_vec\n");printing_vec(a);printf("%lf",b);
#endif
    return b;
}

//The normalzation : the unit vector with the same direction to the given vector a
vector normalize(vector a){
#ifdef DEBUG
    printf("\nnormalize_vec1\n");printing_vec(a);
#endif
    cnumber norm;
    norm.x=1.0/norm_vec(a);
    norm.y=0.0;
    a=scalar_mul_vec(norm,a);
#ifdef DEBUG
    printf("\nnormalize_vec\n");printing_vec(a);printing_cnum(norm);
#endif
    return a;
}

//For a given matrix m having 3 rows and 3 columns, finding the coefficients of the characteristic equation of m
//k.a, k.b, k.c, k.d are the coefficients of cubic term, square term, linear term, constant term respectively.
coeff character_matrix(matrix m){
    cnumber minus_one={-1.0,0.0};
    
    coeff k;
    k.d=det_matrix(m);
    k.b=add_cnum(add_cnum(m.a[0][0],m.a[1][1]),m.a[2][2]); //trace of matrix m
    k.c=minus_cnum(mul_cnum(conj_cnum(k.b),k.d));
    k.a=minus_one;
#ifdef DEBUG
  printf("\ncharacter_matrix\n");printing_matrix(m);printing_cnum(k.a);printing_cnum(k.b);printing_cnum(k.c);printing_cnum(k.d);
#endif
    return k;
}

cnumber det_matrix(matrix m){
    cnumber a[6],b;
    a[0]=mul_cnum(mul_cnum(m.a[0][0],m.a[1][1]),m.a[2][2]);
    a[1]=mul_cnum(mul_cnum(m.a[1][0],m.a[2][1]),m.a[0][2]);
    a[2]=mul_cnum(mul_cnum(m.a[0][1],m.a[1][2]),m.a[2][0]);
    a[3]=minus_cnum(mul_cnum(mul_cnum(m.a[2][0],m.a[1][1]),m.a[0][2]));
    a[4]=minus_cnum(mul_cnum(mul_cnum(m.a[1][0],m.a[0][1]),m.a[2][2]));
    a[5]=minus_cnum(mul_cnum(mul_cnum(m.a[2][1],m.a[1][2]),m.a[0][0]));
    int i;
    b=add_cnum(add_cnum(add_cnum(a[0],a[1]),add_cnum(a[2],a[3])),add_cnum(a[4],a[5]));
#ifdef DEBUG
    printf("\ndet_matrix\n");printing_matrix(m);printing_cnum(b);
#endif
    return b;
}

//finding the roots of the given cubic equation
cbeqrt find_cubic(coeff x){
    cnumber one={1.0,0.0};
    cnumber two={2.0,0.0};
    cnumber three={3.0,0.0};
    cnumber two_over_27={2.0/27.0,0.0};
    cnumber twentyseven={27.0,0.0};
    cnumber imaginary={0.0,1.0};
    cnumber twelve={12.0,0.0};
    cnumber eightyone={81.0,0.0};
    cnumber fiftyfour={54.0,0.0};
    cnumber thirtysix={36.0,0.0};
    cnumber eight={8.0,0.0};
    cnumber one08={108.0,0.0};
    cnumber six={6.0,0.0};
    cbeqrt z;
    cnumber p,q,D;
    cnumber u,v;
    cnumber y1,y2,y3;
    
//    printing_cnum(x.a);printing_cnum(x.b);printing_cnum(x.c);printing_cnum(x.d);    
  D=add_cnum(add_cnum(add_cnum(mul_cnum(mul_cnum(mul_cnum(x.c,x.c),x.c),minus_cnum(twelve)),mul_cnum(mul_cnum(mul_cnum(x.b,x.b),mul_cnum(x.c,x.c)),minus_cnum(three))),add_cnum(mul_cnum(mul_cnum(x.c,mul_cnum(x.b,x.d)),fiftyfour),mul_cnum(mul_cnum(mul_cnum(x.b,x.b),x.b),mul_cnum(twelve,x.d)))),mul_cnum(eightyone,mul_cnum(x.d,x.d)));
  u=div_cnum(cbrt_cnum(add_cnum(add_cnum(mul_cnum(mul_cnum(x.c,six),mul_cnum(x.b,six)),mul_cnum(one08,x.d)),add_cnum(mul_cnum(eight,mul_cnum(mul_cnum(x.b,x.b),x.b)),mul_cnum(twelve,sqrt_cnum(D))))),six);
    v=div_cnum(add_cnum(div_cnum(x.c,three),mul_cnum(div_cnum(x.b,three),div_cnum(x.b,three))),u);
    
    y1=add_cnum(u,v);
    y2=add_cnum(minus_cnum(div_cnum(y1,two)),mul_cnum(mul_cnum(imaginary,div_cnum(sub_cnum(u,v),two)),sqrt_cnum(three)));
    y3=sub_cnum(minus_cnum(div_cnum(y1,two)),mul_cnum(mul_cnum(imaginary,div_cnum(sub_cnum(u,v),two)),sqrt_cnum(three)));
    
    z.a[0]=add_cnum(y1,div_cnum(x.b,three));
    z.a[1]=add_cnum(y2,div_cnum(x.b,three));
    z.a[2]=add_cnum(y3,div_cnum(x.b,three));
#ifdef DEBUG
  printf("\nfind_cubic\n");printing_cnum(x.a);printing_cnum(x.b);printing_cnum(x.c);printing_cnum(x.d);printing_cnum(D);printing_cnum(u);printing_cnum(v);printing_cnum(z.a[0]);printing_cnum(z.a[1]);printing_cnum(z.a[2]);
#endif    
    return z;
//    printf("-------------------------\n");
}

//finding the eivencector equation : m -> m-ik
matrix eigenvectoreq(matrix m,cbeqrt z,int k){
//    cnumber one={1.0,0.0};
    m.a[0][0]=sub_cnum(m.a[0][0],z.a[k]);
    m.a[1][1]=sub_cnum(m.a[1][1],z.a[k]);
    m.a[2][2]=sub_cnum(m.a[2][2],z.a[k]);
#ifdef DEBUG
    printf("\neigenvectoreq\n");printing_matrix(m);
#endif
    return m;
}

//making the given matrix echelon form
matrix echelon(matrix m){
    int i,j;
    cnumber k;
    for(i=1;i<3;i++){
        k=m.a[i][0];
        for(j=0;j<3;j++){
            m.a[i][j]=sub_cnum(m.a[i][j],div_cnum(mul_cnum(m.a[0][j],k),m.a[0][0]));
        }
    }
    for(i=2;i<3;i++){
        k=m.a[i][1];
        for(j=1;j<3;j++){
            m.a[i][j]=sub_cnum(m.a[i][j],div_cnum(mul_cnum(m.a[1][j],k),m.a[1][1]));
        }
    }
#ifdef DEBUG
    printf("\nechelon\n");printing_matrix(m);
#endif

    return m;
}

//solving the matrix equation mx=0
vector equationsystem(matrix m,int k){
//printing_matrix(m);
    cnumber zero={0.0,0.0};
    cnumber one={1.0,0.0};
    cnumber image={0.0,1.0};
    cnumber a,b,c,d,e,f,g,h,i,j;
    vector v;
    
//    m=echelon(m);
            a=m.a[0][0];
            b=m.a[0][1];
            c=m.a[0][2];
            
            d=m.a[1][0];
            e=m.a[1][1];
            f=m.a[1][2];
            
            h=m.a[2][0];
            j=m.a[2][1];
            i=m.a[2][2];

/*    if(norm_cnum(e)==norm_cnum(zero)) //degenerated case of 1 and 2
        {
            v.z=zero;
            v.y=one;
            v.x=minus_cnum(div_cnum(b,a));
        }
    else//degenerated case of 1 and 3
        {
            v.z=one;
            v.y=minus_cnum(div_cnum(f,e));
            v.x=minus_cnum(div_cnum(add_cnum(mul_cnum(v.y,b),c),a));
        }
*/
v.x=sub_cnum(mul_cnum(m.a[1][2],m.a[0][1]),mul_cnum(m.a[1][1],m.a[0][2]));
v.y=sub_cnum(mul_cnum(m.a[1][0],m.a[0][2]),mul_cnum(m.a[0][0],m.a[1][2]));
v.z=sub_cnum(mul_cnum(m.a[0][0],m.a[1][1]),mul_cnum(m.a[0][1],m.a[1][0]));

//printing_matrix(m);
#ifdef DEBUG
    printf("\nequationsystem\n");printing_vec(v);
#endif
    return normalize(v);
}

//finding the eigenvalues of the given matrix m
//now, it is set that output is always (1,1,1)
cbeqrt eigenvalue_find(matrix m){
    coeff abcd=character_matrix(m);
    cbeqrt z=find_cubic(abcd);
#ifdef DEBUG
    printf("\neigenvalue_find\n");printing_cnum(z.a[0]);printing_cnum(z.a[1]);printing_cnum(z.a[2]);
#endif
    return z;

}

//finding the i-th eigenvectors of the given matrix m
vector eigenvector_find(matrix m,int i,int j){
//    printing_matrix(m);
    cbeqrt z=eigenvalue_find(m);
    vector v;
    v=equationsystem(eigenvectoreq(m,z,i),j);
//    printing_vec(v);
#ifdef DEBUG
    printf("\neigenvector_find\n");
    printing_cnum(z.a[i]);printing_vec(v);
#endif
    return normalize(v);
}

//the i-th row vector of the given matrix m
vector row_vector(matrix m,int i){
    vector v;
    v.x=m.a[i][0];
    v.y=m.a[i][1];
    v.z=m.a[i][2];
#ifdef DEBUG
    printf("\nrow_vector\n");printing_matrix(m);printing_vec(v);
#endif
    return v;
}

//the i-th column vector of the given matrix m
vector col_vector(matrix m,int i){
    vector v;
    v.x=m.a[0][i];
    v.y=m.a[1][i];
    v.z=m.a[2][i];
#ifdef DEBUG
    printf("\ncol_vector\n");printing_matrix(m);printing_vec(v);
#endif
    return v;
}

//Printing functions (for debugging)
void printing_vec(vector v){
    printf("\nprinting_vec\n");
    printf("%lf+%lf*I,%lf+%lf*I,%lf+%lf*I\n",v.x.x,v.x.y,v.y.x,v.y.y,v.z.x,v.z.y);
}

void printing_cnum(cnumber c){
    printf("%lf+%lf*I\n ",c.x,c.y);
}

void printing_matrix(matrix mix)
{
printf("[[");
            printf("%lf+",mix.a[0][0].x);
            printf("%lf*I",mix.a[0][0].y);
            printf(",%lf+",mix.a[0][1].x);
            printf("%lf*I",mix.a[0][1].y);
            printf(",%lf+",mix.a[0][2].x);
            printf("%lf*I",mix.a[0][2].y);
printf("],\n[");            
            printf("%lf+",mix.a[1][0].x);
            printf("%lf*I",mix.a[1][0].y);
            printf(",%lf+",mix.a[1][1].x);
            printf("%lf*I",mix.a[1][1].y);
            printf(",%lf+",mix.a[1][2].x);
            printf("%lf*I",mix.a[1][2].y);
printf("],\n[");
            printf("%lf+",mix.a[2][0].x);
            printf("%lf*I",mix.a[2][0].y);
            printf(",%lf+",mix.a[2][1].x);
            printf("%lf*I",mix.a[2][1].y);
            printf(",%lf+",mix.a[2][2].x);
            printf("%lf*I",mix.a[2][2].y);
printf("]]");
            printf("\n");
}