线性代数实现（部分）

1.数学能力决定程序员发展上限。

2. C语言编程也是可以使用部分面向对象思想的，关键在人，而不在工具。

3. 如下编程风格大量出现在linux kernel& 其他开源代码中，很明显使用的面向对象编程思想。

think about it !

--------------------------------------------------------------------------------------------------

 

linear_algebra.h

#ifndef _LINEAR_ALGEBRA_H_
#define _LINEAR_ALGEBRA_H_

static void matrix_show(unsigned char *fmt,...);//for test use;
/*
static int det(int *A,int n);
static void matrix_transposed(int *m,int r,int c);
static void matrix_add(int *A,int *B,int r,int c);
static void matrix_sub(int *A,int *B,int r,int c);
static void matrix_multi_lambda(int *A,int r,int c ,int lambda);
static void matrix_div_lambda (double *A,int r,int c,double lambda);
static int  *matrix_multi_AB(int *A,int *B,int Ar,int Ac,int Br,int Bc);
static void matrix_power(int *A,int n,int m);
static int matrix_adjugate(int *A,int n);
static double  *matrix_inverse(int *A,int n);
static int matrix_rank(int *A,int n);
*/

struct namespace_matrix_func{
    void    (*show)(unsigned char *fmt,...);
    int     (*det)(int *A,int n);
    void    (*transposed)(int *m,int r,int c);
    void    (*add)(int *A,int *B,int r,int c);
    void    (*sub)(int *A,int *B,int r,int c);
    void    (*multi_lambda)(int *A,int r,int c ,int lambda);
    void    (*div_lambda) (double *A,int r,int c,double lambda);
    int    *(*multi)(int *A,int *B,int Ar,int Ac,int Br,int Bc);
    void    (*power)(int *A,int n,int m);
    int     (*adjugate)(int *A,int n);
    double *(*inverse)(int *A,int n);
    int     (*rank)(int *A,int n);
} ;

void namespace_matrix_func_init(struct namespace_matrix_func *mf); 

#endif

linear_algebra.c 主体部分

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <math.h>
#include <assert.h>


#include "linear_algebra.h"
static void matrix_show_double(void *M,int r,int c)
{
   int i,j;
   double *A=(double *)M;             
   for(i=0;i<r;i++){                    
         for(j=0;j<c;j++)                   
           printf("%7.2f ",A[i*c+j]);       
         printf("\n");                      
   } 
} 
static void matrix_show_int(void *M,int r,int c)                  
{                                                           
    int i,j;                                                                              
    int *A=(int *)M;               
    for(i=0;i<r;i++){                     
         for(j=0;j<c;j++)                   
           printf("%6d ",A[i*c+j]);         
         printf("\n");                      
    }                                                                        
}
static int fmt_transfer(unsigned char *fmt)
{
  unsigned char c=1,*p;
  
  int ret=0;
  if (fmt==NULL)return 0;
  
  while(c!='\0') 
  {
    c=*fmt++;
    if(c=='%'){
      p=fmt-1;

      c=*fmt++;
      while(c!='\0'&&c!='f'&&c!='d')c=*fmt++;
      
      if(c=='d')ret=sizeof(int);
      if(c=='f')ret=sizeof(double);
      
      if(c!='\0')
         strcpy(p,fmt);
      break;
    }
  
  }

   return ret;

}
static void matrix_show(unsigned char *fmt,...)
{
  void *a;
  int r,c,ret;
  va_list ap;
  
  ret=fmt_transfer(fmt);
  printf("%s\n",fmt);
  
  va_start(ap,fmt); 
   
  a=va_arg(ap,void*); 
  r=va_arg(ap,int);
  c=va_arg(ap,int);
  
  if(ret==sizeof(double))
           matrix_show_double(a,r,c);
  else if(ret==sizeof(int))
           matrix_show_int(a,r,c);
  else     printf("unkown format\n");
  
  
  va_end(ap);
}

static int det(int *A,int n)
{//|A| 矩阵A所有元素构成的行列式   
    int f=0;
    int *B=(int*)malloc(sizeof(int)*(n-1)*(n-1));
    int i,j,k;
    if(n>2){
       for(j=0;j<n;j++){ 
          for(i=0;i<n;i++){
             if(i<j){
               for(k=0;k<n-1;k++){
                  B[k*(n-1)+i]=A[(k+1)*n+i];
               }
             }else if(i>j){
               for(k=0;k<n-1;k++){
                  B[k*(n-1)+i-1]=A[(k+1)*n+i];
               }
             }        
          }    
          if(j%2==1)f+=A[j]*det(B,n-1)*(-1);
          else      f+=A[j]*det(B,n-1); 
               
       }
          
    }
    else  f=A[0]*A[3]-A[1]*A[2];
  
    free(B);
    return f;
}
static void matrix_transposed(int *m,int r,int c)
{//A^T  转置矩阵
  int i,j,k;
  int *t=(int *)malloc(sizeof(int)*r*c);

  memcpy(t,m,sizeof(int)*r*c);

  for(i=0;i<c;i++)
    for(j=0;j<r;j++)
      m[i*r+j]=t[j*c+i] ;

  free(t);
}
static void matrix_add(int *A,int *B,int r,int c)
{  //A=A+B;矩阵加法
   r=r*c;
   for(c=0;c<r;c++){
     A[c]=A[c]+B[c] ; 
   }
}
static void matrix_sub(int *A,int *B,int r,int c)
{  //A=A-B;矩阵减法
   r=r*c;
   for(c=0;c<r;c++){
     A[c]=A[c]-B[c] ; 
   }
}
static void matrix_multi_lambda(int *A,int r,int c ,int lambda)
{//  λA 常数*矩阵
   int i;  
   for(i=0;i<r*c;i++) {  
      A[i]= lambda*A[i];  
   }
}
static void matrix_div_lambda (double *A,int r,int c,double lambda)
{  //  A/λ  矩阵/常数
   int i;  
   for(i=0;i<r*c;i++) {  
      A[i]= A[i]/lambda;  
   }
}
//                                      A行   A列    B行    B列
static int  *matrix_multi_AB(int *A,int *B,int Ar,int Ac,int Br,int Bc)
{ //AB   矩阵相乘
  int i,j,k,t;
  
  if(Ac!=Br)return NULL;
  int *M=(int *)malloc(Ar*Bc*sizeof(int));

  for(i=0;i<Ar;i++) {
      for(j=0;j<Bc;j++){
         t=0;
         for(k=0;k<Ac;k++)
                t+=A[i*Ac+k]*B[j+k*Bc];                   
         M[i*Bc+j]= t;      
      } 
  }
  return M;
}

static void matrix_power(int *A,int n,int m)
{ //A^n 方阵的n次幂
  int *M=NULL,i;
  for(i=0;i<n;i++){
      M=matrix_multi_AB(A,A,m,m,m,m) ;
      memcpy(A,M,m*m*sizeof(int));  
      free(M);
  }
}

static int matrix_adjugate(int *A,int n)
{ //A*  //伴随方阵  
  if(n<=2)return -1;
  int i,j,count;
  int *B=(int *)malloc(sizeof(int)*(n-1)*(n-1));
  int *T=(int *)malloc(sizeof(int)*(n)*(n));
  
  for(j=0;j<n*n;j++){
     count=0;
     for(i=0;i<n*n;i++){
     
        if(i/n==j/n)continue;
        if(i%n==j%n)continue;
     
        B[count]=A[i];
        count++;
     
     }
     
     //求矩阵余子式 
     if(j%2==1) T[j]=-det(B,n-1) ;
     else       T[j]=det(B,n-1) ;  
  }
  
  matrix_transposed(T,n,n);
  
  memcpy(A,T,sizeof(int)*n*n) ;
          
  free(B);
  free(T);
  return 0;      
}

static double  *matrix_inverse(int *A,int n)
{//A^-1  逆矩阵
   int lambda,i;
   lambda=det(A,n);
   if(lambda==0)return NULL;
   
   
   double *m=(double *)malloc(sizeof(double)*n*n);
   
   matrix_adjugate(A,n);
   
   for(i=0;i<n*n;i++){  
        m[i]=(double )A[i]; 
   }

   matrix_div_lambda (m,n,n,(double )lambda);

   return m;
}
static int matrix_rank(int *A,int n)
{//rank(A);矩阵的秩


  return 0;
}

void namespace_matrix_func_init(struct namespace_matrix_func *mf)
{
  mf->show=matrix_show;
  mf->det=  det;
  mf->transposed= matrix_transposed;
  mf->add= matrix_add;
  mf->sub= matrix_sub;
  mf->multi_lambda=matrix_multi_lambda;
  mf->div_lambda=matrix_div_lambda;
  mf->multi= matrix_multi_AB;
  mf->power=matrix_power;
  mf->adjugate= matrix_adjugate;
  mf->inverse=  matrix_inverse;
  mf->rank=matrix_rank;  
}  

matrix_test.c  测试代码

#include <stdio.h>
#include <stdlib.h>
#include "linear_algebra.h"
int main(int argc,char **argv)
{
  int i,*m;
  int a[12]={1,2,0,-3,  3,-1,1,0, -2,0,4,-3};

  int det_A[]={5,1,1,0, 1,1,-1,1, 2,2,-1,1, 3,1,2,3} ;
  struct namespace_matrix_func *mf=(struct namespace_matrix_func*)malloc(sizeof(struct namespace_matrix_func));
  namespace_matrix_func_init(mf);
  
  //私有函数测试(static 代替 private)，编译器会报错，表示测试成功 。
  //matrix_show("matrix_show test:%d",a,3,4);
  
  mf->transposed(a,3,4);
    
  mf->show("matrix_transposed result:%d",a,4,3);
  
  i=mf->det(det_A,4);
  printf("detA=%d\n",i);
  
  
  
  int A[]={-2,4, 1,-2};
  int B[]={2,4,-3, -6} ;
  
  
  m=mf->multi(A,B,2,2,2,2);
  
  mf->show("matrix_multi_AB result:%d",A,2,2);
  
  
  int A1[]={1,2,3, 2,2,1, 3,4,3};
  double *mt; 
  //matrix_adjugate(A1,3);
  
  mt=mf->inverse(A1,3);
  
  mf->show("matrix_inverse result:%f",mt,3,3);

  system("pause");
  free (mf) ;
  return EXIT_SUCCESS;
}