C++实现R语言向量化运算（向量类：c 矩阵类：matrix）2015.9.11
2016-07-16 21:49
716 查看
类源代码：
#include
using std::cout;
using std::endl;
using std::cin;
using std::istream;
using std::ios_base;
#include
using std::vector;
#include
using std::initializer_list;
#include
using std::for_each;
using std::sort;
using std::copy;
using std::find;
#include
using std::runtime_error;
#include
using std::srand;
using std::rand;
using std::system;
#include
using std::time;
#include
using std::sqrt;
#include
using std::pair;
using std::make_pair;
#include
using std::setw;
#include
using std::multiset;
using std::set;
#include
using std::unordered_map;
using std::unordered_multimap;
#include
using std::multimap;
using std::map;
class matrix;
class c
{
private:
vector c_val;
size_t dim;
public:
c():c_val(vector()),dim(0){}
c(initializer_listlst);
c(const vector&c_v):c_val(c_v),dim(c_v.size()){}
c(const matrix&m);
void out(size_t n=1)const{cout<<"The vector is:\n";for(size_t i=0;i<<setw(n)<<c_val[i]<<" ";cout<<endl<<endl;}
double operator[](size_t i)const{return c_val[i];}
double& operator[](size_t i){return c_val[i];}
size_t re_dim()const{return dim;}
friend c operator+(const
4000
c&a,const c&b);
friend c operator-(const c&a);
friend c operator-(const c&a,const c&b);
friend c rep(const double&d,size_t n);
friend c operator+(const c&a,const double&d);
friend c operator-(const c&a,const double&d);
friend c operator+(const double&d,const c&a);
friend c operator-(const double&d,const c&a);
friend c operator*(const double&d,const c&a);
friend c operator*(const c&a,const double&d);
friend c operator/(const double&d,const c&a);
friend c operator/(const c&a,const double&d);
friend c operator*(const c&a,const c&b);
friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}
friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}
vector const_re_c_val()const{return c_val;}
};
class matrix
{
private:
vector m_val;
size_t rdim;
size_t cdim;
public:
matrix():m_val(vector()),rdim(0),cdim(0){}
matrix(const vector&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}
matrix(size_t m,size_t n,const double&d):m_val(vector(m,vector(n,d))),rdim(m),cdim(n){}
matrix(const c&v,size_t r,size_t c);
matrix(const c&v);
c operator[](size_t i)const{return m_val[i];}
c& operator[](size_t i){return m_val[i];}
void out(size_t n=3)const;
c col(size_t i)const;
void col_change(size_t i,const c&a);
friend matrix operator+(const matrix&a,const matrix&b);
friend matrix operator-(const matrix&a);
friend matrix operator-(const matrix&a,const matrix&b);
friend matrix operator+(const matrix&a,const double&d);
friend matrix operator+(const double&d,const matrix&a);
friend matrix operator-(const matrix&a,const double&d);
friend matrix operator-(const double&b,const matrix&a);
friend matrix operator*(const matrix&a,const double&d);
friend matrix operator*(const double&d,const matrix&a);
friend matrix operator/(const matrix&a,const double&d);
friend matrix operator/(const double&d,const matrix&a);
void Swap(size_t i,size_t j)
{c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}
void col_Swap(size_t i,size_t j)
{
auto temp=this->col(i);
this->col_change(i,this->col(j));
this->col_change(j,temp);
}
friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}
friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}
pair re_dim()const{return make_pair(rdim,cdim);}
friend c operator%(const matrix&m,const c&v);
friend matrix operator%(const matrix&a,const matrix&b);
friend matrix operator%(const c&v,const matrix&m);
};
//类c的有关定义
c::c(initializer_listlst)
{
vectorout;
for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});
c_val=out;
dim=out.size();
}
c operator+(const c&a,const c&b)
{
size_t MAX=b.re_dim();
if(a.re_dim()>b.re_dim())
MAX=a.re_dim();
vectortemp0(MAX,0);
vectortemp1(MAX,0);
vectortemp2(MAX,0);
copy(a.c_val.begin(),a.c_val.end(),temp0.begin());
copy(b.c_val.begin(),b.c_val.end(),temp1.begin());
for(size_t i=0;i
temp2[i]=temp0[i]+temp1[i];
return c(temp2);
}
c operator-(const c&a)
{
c out(a);
for(size_t i=0;i
out[i]=-out[i];
return out;
}
c operator-(const c&a,const c&b)
{
return (a+(-b));
}
c rep(const double&d,size_t n)
{
return c(vector(n,d));
}
c operator+(const c&a,const double&d)
{
return a+rep(d,a.re_dim());
}
c operator-(const c&a,const double&d)
{
return (a+(-d));
}
c operator+(const double&d,const c&a)
{
return a+d;
}
c operator-(const double&d,const c&a)
{
return (-a)+d;
}
c operator*(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i
out[i]*=d;
return out;
}
c operator*(const c&a,const double&d)
{
return d*a;
}
c operator/(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i
out[i]=d/a[i];
return out;
}
c operator/(const c&a,const double&d)
{
c out(a);
for(size_t i=0;i
out[i]=a[i]/d;
return out;
}
c::c(const matrix&m)
{
vectortemp;
for(size_t j=0;j
for(size_t i=0;i
temp.push_back(m[i][j]);
c_val=temp;
dim=temp.size();
}
c operator*(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
c out(a);
for(size_t i=0;i
out[i]*=b[i];
return out;
}
//类matrix的有关定义
void matrix::out(size_t n)const
{
cout<<"The matrix is:\n";
for(size_t i=0;i
{
for(size_t j=0;j
cout<<setw(n)<<(*this)[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
c matrix::col(size_t i)const
{
c out(rep(0,rdim));
for(size_t j=0;j
out[j]=(*this)[j][i];
return out;
}
void matrix::col_change(size_t i,const c&a)
{
for(size_t j=0;j
(*this)[j][i]=a[j];
}
matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)
{
matrix out(r,c,0);
for(size_t i=0;i
out[i%r][i/r]=v[i];
(this->m_val)=out.m_val;
}
matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)
{
matrix out(v.re_dim(),1,0);
for(size_t i=0;i
out[i][0]=v[i];
this->m_val=out.m_val;
}
matrix operator+(const matrix&a,const matrix&b)
{
size_t r_MAX=a.rdim,c_MAX=a.cdim;
if(a.rdim
r_MAX=b.rdim;
if(a.cdim
c_MAX=b.cdim;
matrix tempa(r_MAX,c_MAX,0);
for(size_t i=0;i
for(size_t j=0;j
tempa[i][j]=a[i][j];
matrix tempb(r_MAX,c_MAX,0);
for(size_t i=0;i
for(size_t j=0;j
tempb[i][j]=b[i][j];
matrix tempout(r_MAX,c_MAX,0);
for(size_t i=0;i
for(size_t j=0;j
tempout[i][j]=tempa[i][j]+tempb[i][j];
return tempout;
}
matrix operator-(const matrix&a)
{
matrix out(a);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=-out[i][j];
return out;
}
matrix operator-(const matrix&a,const matrix&b)
{
return ((-b)+a);
}
matrix operator+(const matrix&a,const double&d)
{
return a+matrix(a.rdim,a.cdim,d);
}
matrix operator+(const double&d,const matrix&a)
{
return a+d;
}
matrix operator-(const matrix&a,const double&d)
{
return a+(-d);
}
matrix operator-(const double&b,const matrix&a)
{
return (-a)+b;
}
matrix operator*(const matrix&a,const double&d)
{
matrix out(a);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]*=d;
return out;
}
matrix operator*(const double&d,const matrix&a)
{
return a*d;
}
matrix operator/(const matrix&a,const double&d)
{
return a*(1/d);
}
matrix operator/(const double&d,const matrix&a)
{
matrix out(a.rdim,a.cdim,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=d/a[i][j];
return out;
}
//基本运算函数
matrix diag(const c&v)
{
matrix out(v.re_dim(),v.re_dim(),0);
for(size_t i=0;i
out[i][i]=v[i];
return out;
}
matrix diag(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i
out[i][i]=m[i][i];
return out;
}
double inter_product(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
double out=0;
for(size_t i=0;i
out+=(a[i]*b[i]);
return out;
}
c operator%(const matrix&m,const c&v)
{
if(m.re_dim().second!=v.re_dim())
throw runtime_error("dim error!\n");
c out(rep(0,m.re_dim().first));
for(size_t i=0;i
out[i]=inter_product(m[i],v);
return out;
}
//同下
double prod(const c&v)
{
double out=1;
for(size_t i=0;i
out*=v[i];
return out;
}
//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)
double sum(const c&v)
{
double out=v[0];
for(size_t i=1;i
out+=v[i];
return out;
}
matrix operator%(const c&v,const matrix&m)
{
if(v.re_dim()!=m.re_dim().first)
throw runtime_error("dim error!\n");
matrix out(1,m.re_dim().second,0);
for(size_t i=0;i
out[0][i]=inter_product(v,m.col(i));
return out;
}
matrix outer_product(const c&a,const c&b)
{
matrix out(a.re_dim(),b.re_dim(),0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=a[i]*b[j];
return out;
}
size_t length(const c&v)
{
return v.re_dim();
}
c sort(const c&v)
{
vectorout;
for(size_t i=0;i
out.push_back(v[i]);
sort(out.begin(),out.end());
return c(out);
}
double mean(const c&v)
{
return sum(v)/length(v);
}
//其转换版本返回矩阵转为向量的方差
double var(const c&v)
{
return mean((v-mean(v))*(v-mean(v)));
}
c connect(initializer_list lst)
{
vectortemp;
for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});
return c(temp);
}
//标准差
template
double STD(const any&a)
{
return sqrt(var(a)*(length(a))/(length(a)-1));
}
matrix rbind(initializer_listlst)
{
size_t r=0;
auto a=lst.begin();
while(a!=lst.end())
{
a++;
r++;
}
size_t co=0;
for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});
matrix out(r,co,0);
size_t i=0;
a=lst.begin();
while(a!=lst.end())
{
vectortemp(co,0);
for(size_t i=0;iconst_re_c_val().size();i++)
temp[i]=a->operator[](i);
out[i]=c(temp);
i++;
a++;
}
return out;
}
size_t dim(const c&v)
{
return v.re_dim();
}
pair dim(const matrix&m)
{
return m.re_dim();
}
double min(const c&v)
{
double temp=v[0];
for(size_t i=1;i
if(v[i]
temp=v[i];
return temp;
}
double max(const c&v)
{
double temp=v[0];
for(size_t i=1;i
if(v[i]>temp)
temp=v[i];
return temp;
}
c pmin(initializer_listlst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i
{
temp[i]=min(*a);
a++;
}
return temp;
}
c pmax(initializer_listlst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i
{
temp[i]=max(*a);
a++;
}
return temp;
}
template
pair range(const any&v)
{
cout<<"MIN: "<<min(v)<<" MAX: "<<max(v)<<endl;
return make_pair(min(v),max(v));
}
double median(const c&v)
{
return sort(v)[length(v)/2];
}
double var(const c&a,const c&b)
{
if(length(a)!=length(b))
throw runtime_error("dim error!\n");
c t1(a-mean(a)),t2(b-mean(b));
return inter_product(t1,t2)/(length(a)-1);
}
double cov(const c&a,const c&b)
{
return var(a,b);
}
double cor(const c&a,const c&b)
{
return var(a,b)/STD(a)/STD(b);
}
matrix var(const matrix&a,const matrix&b)
{
if(a.re_dim().first!=b.re_dim().first)
throw runtime_error("dim error!\n");
matrix temp(a.re_dim().second,b.re_dim().second,0);
for(size_t i=0;i
for(size_t j=0;j
temp[i][j]=var(a.col(i),b.col(j));
return temp;
}
matrix cov(const matrix&a,const matrix&b)
{
return var(a,b);
}
matrix cor(const matrix&a,const matrix&b)
{
matrix temp=var(a,b);
for(size_t i=0;i
for(size_t j=0;j
temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));
return temp;
}
matrix cov(const matrix&a)
{
return var(a,a);
}
matrix cor(const matrix&a)
{
return cor(a,a);
}
//由于关联容器的删除方法利用误差排序的方法
//利用顺序容器迭代器的删除方法有困难
c rank(const c&v)
{
vectortemp00(v.const_re_c_val());
for(size_t i=0;i
temp00[i]+=(0.00000001*i);
vectortemp0(temp00);
multisettemp2(temp0.begin(),temp0.end());
vectortemp1(v.re_dim());
size_t i=0,j=0;
double pre_val=min(c(vector(temp2.begin(),temp2.end())));
vector::iterator s0;
while(i
{
double val=min(c(vector(temp2.begin(),temp2.end())));
s0=find(temp0.begin(),temp0.end(),val);
if(sqrt((pre_val-val)*(pre_val-val))>0.000001)
j++;
temp1[s0-temp0.begin()]=j;
temp2.erase(val);
i++;
pre_val=val;
}
return c(temp1);
}
c rev(const c&v)
{
vector temp;
auto a=v.const_re_c_val().rbegin();
while(a!=v.const_re_c_val().rend())
{
temp.push_back(*a);
a++;
}
return c(temp);
}
//向量先排序再输出排序前的下标值
c order(const c&v)
{
unordered_multimapm0;
for(size_t i=0;i
m0.insert(make_pair(v[i],i));
multimapm1(m0.begin(),m0.end());
c temp(rep(0,length(v)));
auto a=m1.begin();
size_t i=0;
while(a!=m1.end())
{
temp[i]=a->second;
i++;
a++;
}
return temp;
}
c cumsum(const c&v)
{
c temp(v);
for(size_t i=1;i
temp[i]+=temp[i-1];
return temp;
}
c cumprod(const c&v)
{
c temp(v);
for(size_t i=1;i
temp[i]*=temp[i-1];
return temp;
}
matrix cumsum(const matrix&m)
{
c temp(m);
for(size_t i=1;i
temp[i]+=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
matrix cumprod(const matrix&m)
{
c temp(m);
for(size_t i=1;i
temp[i]*=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
//矩阵运算函数
//初等行变换化上三角求解行列式
double det(const matrix&m)
{
auto f=[](const matrix&mm)->pair
{
if(mm.re_dim().first!=mm.re_dim().second)
throw runtime_error("dim error!\n");
size_t MIN=mm.re_dim().first;
double mult=1;
size_t col=0;
size_t bar=0;
size_t b_c=0;
matrix t(mm);
while(b_c
{
for(size_t i=bar+1;i
{
if(t[i][col]==0)
i++;
else
{
t.Swap(i,bar);
mult*=(-1);
break;
}
}
if(bar!=MIN+1)
{
if(t[bar][col]!=0)
{
for(size_t i=bar+1;i
{
t[bar][i]=t[bar][i]/t[bar][col];
if(i==bar+1)
mult*=t[bar][col];
}
if(col!=MIN-1)
t[bar][col]=1;
}
for(size_t i=bar+1;i
{
t[i]=t[i]-t[i][col]*t[ba
ee37
r];
}
}
col++;
bar++;
b_c=bar;
if(col>bar);
b_c=col;
}
return make_pair(t,mult);
};
pair tt=f(m);
return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));
}
//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c solve(const matrix&m,const c&v)
{
c temp(rep(0,v.re_dim()));
double low=det(m);
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
if(low==0)
throw runtime_error("singular!\n");
for(size_t i=0;i
{
matrix temp0(m);
temp0.col_change(i,v);
temp[i]=det(temp0)/low;
}
return temp;
}
double tr(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));
}
matrix T(const matrix&m)
{
matrix out(m.re_dim().second,m.re_dim().first,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=m[j][i];
return out;
}
matrix cbind(initializer_listlst)
{
return T(rbind(lst));
}
matrix operator%(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first,a.re_dim().second,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=inter_product(a[i],b.col(j));
return out;
}
double F_norm(const matrix&m)
{
return sqrt(tr(T(m)%m));
}
//a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善
//只能求实值
//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中
//可用于求解但效率低
c eigen(const matrix&m,const double&a,const double& b,const double &d=0.0001,const double &e=1)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
vectortemp;
double pre_beg=100000000;
for(double beg=a;beg
{
matrix mm0(diag(rep(beg,m.re_dim().first))-m);
double s0=det(mm0);
if(sqrt(s0*s0)
{
if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)
{
temp.push_back(beg);
pre_beg=beg;
cout<<"Done!\n";
cout<<beg<<endl;//eigenvalue
}
}
}
return c(temp);
}
//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）
c eigen(const matrix&m,const double&d=0.1)
{
double dd=d;
double ra=F_norm(m);
double trr=tr(m);
c temp;
double r=0;
do
{
dd*=0.1;
temp=eigen(m,-ra,ra,dd);
if(temp.re_dim()==m.re_dim().first)
break;
r=sum(temp)-trr;
}
while(sqrt((r*r)>1));
//与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略
//此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握
//只适用于所有特征值为实值的情况
//如有复特征值，只输出实数的并循环不能结束
return temp;
}
matrix M(const matrix&m,size_t k,size_t l)
{
matrix temp(m);
set s1;
for(size_t i=0;i
s1.insert(i);
set s2;
for(size_t i=0;i
s2.insert(i);
s1.erase(k);
s2.erase(l);
vectortemp00;
for(auto j:s2)
for(auto i:s1)
temp00.push_back(temp[i][j]);
return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);
}
matrix A_accompany(const matrix&m)
{
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);
//代数余子式
return T(out);
}
//求方阵的逆矩阵（利用伴随矩阵方法）
matrix solve(const matrix&m)
{
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return A_accompany(m)/de;
}
//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆
c solve_eq(const matrix&m,const c&v)
{
c out;
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return solve(m)%v;
}
测试程序：
int main()
{
//测试程序：
c c0({35.1,12});
vectorv0;
for(size_t i=0;i<10;i++)
v0.push_back(i);
c c1(v0);
c0.out();
c1.out();
cout<<c0[1]<<endl;
c0[1]=10;
cout<<c0[1]<<endl;
(c0+c1).out();
(-c0).out();
(c0-c1).out();
rep(1,10).out();
(c0+rep(1,10)).out();
(c0-rep(1,10)).out();
(rep(1,10)+1).out();
(rep(1,10)-1).out();
(10/rep(1,10)).out();
(rep(1,10)/10).out();
cout<<((10/rep(1,10))==(rep(1,10)/10))<<endl;
matrix m0=matrix(3,4,10);
m0[1]=rep(4,4);
m0.out();
m0.col(2).out();
m0.col_change(0,c({1,2,3}));
m0.out();
matrix m1=(matrix(c({1,5,7,8,9,0,3,4,5}),3,3)+matrix(c({3,5,7,9}),2,2));
m1.out();
(-m1).out();
(m1-matrix(rep(1,10),5,2)).out();
(m1+1).out();
(m1*0.6).out();
(m1-1).out();
(1-m1).out();
(0.6*m1).out();
(1/m1).out();
(m1/10).out();
m1.out();
m1.Swap(1,2);
m1.out();
m1.col_Swap(1,2);
m1.out();
matrix m2=m1;
m1.Swap(0,1);
cout<<(m1==m2)<<endl;
cout<<m1.re_dim().first<<" "<<m1.re_dim().second<<endl;
matrix A(c({1,2,1,4,7,1,5,2,1}),3,3);
cout<<det((A+diag(rep(10,A.re_dim().first))))<<endl;
solve(matrix(c({30,7.3,4.2,0,0,6,8,7,0}),3,3),c({11,12,13})).out();
m1.out();
cout<<F_norm(m1)<<endl;
eigen(matrix(c({50,1,1,9,7,1,4,1,13}),3,3)).out();
solve(matrix(c({50,1,1,9,7,1,4,1,13}),3,3),c({1,1,1})).out();
solve(matrix(c({50,1,1,9,7,1,4,1,13}),3,3)).out();
solve_eq(matrix(c({50,1,1,9,7,1,4,1,13}),3,3),c({1,1,1})).out();
m2.out();
M(m2,0,0).out();
cout<<sum(m2)<<endl;
cout<<prod(m2)<<endl;
(c({1,1})%matrix(c({5,6,7,8,9,10}),2,3)).out();
outer_product(c({1,1}),c({5,6,7})).out();
diag(m2).out();
diag(c(diag(m2))).out();
cout<<length(m2)<<endl;
cout<<length(m2[0])<<endl;
sort(c({3,2,1})).out();
cout<<mean(m2)<<endl;
cout<<mean(m2.col(0))<<endl;
cout<<var(matrix(c({1,2,3,4,5,6,7,8,9}),3,3))<<endl;
connect({c({1,1,1}),c({-1,1,-1})}).out();
cout<<STD(matrix(c({1,2,3,4,5,6,7,8,9}),3,3))<<endl;
matrix m3=rbind({c({1,2}),c({3,4,5})});
cbind({c({1,2}),c({3,4,5})}).out();
cout<<dim(m3).first<<" "<<dim(m3).second<<endl;
m3.out();
cout<<dim(m3[0])<<endl;
cout<<min(m3)<<endl;
cout<<min(m3[0])<<endl;
cout<<max(m3)<<endl;
cout<<max(m3[0])<<endl;
c(m3).out();
pmin({m3[0],m3[1]}).out();
pmax({m3[0],m3[1]}).out();
range(m3);
cout<<median(m3)<<endl;
var(m3,m3).out();
cov(m3).out();
cor(m3).out();
rank(c({0,3,4,4,5,6,7,6,6,3})).out();
rev(c({1,5,7,6,9})).out();
order(c({0,3,4,4,5,6,7,6,3})).out();
m3.out();
cumsum(m3[0]).out();
cumsum(m3).out();
cumprod(m3[0]).out();
cumprod(m3).out();
cout<<det(matrix(c({50,7,8,9,10,4,3,2,9,4,8,7,1,1,1,10}),4,4))<<endl;
return 0;
}
运行结果：
The vector is:
35.1 12
The vector is:
0 1 2 3 4 5 6 7 8 9
12
10
The vector is:
35.1 11 2 3 4 5 6 7 8 9
The vector is:
-35.1 -10
The vector is:
35.1 9 -2 -3 -4 -5 -6 -7 -8 -9
The vector is:
1 1 1 1 1 1 1 1 1 1
The vector is:
36.1 11 1 1 1 1 1 1 1 1
The vector is:
34.1 9 -1 -1 -1 -1 -1 -1 -1 -1
The vector is:
2 2 2 2 2 2 2 2 2 2
The vector is:
0 0 0 0 0 0 0 0 0 0
The vector is:
10 10 10 10 10 10 10 10 10 10
The vector is:
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
0
The matrix is:
10 10 10 10
4 4 4 4
10 10 10 10
The vector is:
10 4 10
The matrix is:
1 10 10 10
2 4 4 4
3 10 10 10
The matrix is:
4 15 3
10 18 4
7 0 5
The matrix is:
-4 -15 -3
-10 -18 -4
-7 -0 -5
The matrix is:
3 14 3
9 17 4
6 -1 5
-1 -1 0
-1 -1 0
The matrix is:
5 16 4
11 19 5
8 1 6
The matrix is:
2.4 9 1.8
6 10.8 2.4
4.2 0 3
The matrix is:
3 14 2
9 17 3
6 -1 4
The matrix is:
-3 -14 -2
-9 -17 -3
-6 1 -4
The matrix is:
2.4 9 1.8
6 10.8 2.4
4.2 0 3
The matrix is:
0.25 0.0666667 0.333333
0.1 0.0555556 0.25
0.142857 inf 0.2
The matrix is:
0.4 1.5 0.3
1 1.8 0.4
0.7 0 0.5
The matrix is:
4 15 3
10 18 4
7 0 5
The matrix is:
4 15 3
7 0 5
10 18 4
The matrix is:
4 3 15
7 5 0
10 4 18
0
3 3
1880
The vector is:
-0.12533 2.2544 1.84499
The matrix is:
7 5 0
4 3 15
10 4 18
27.6405
Done!
6.67574
Done!
6.67374
Done!
12.9967
Done!
50.3227
The vector is:
6.67374 12.9967 50.3227
The vector is:
-0.00961538 0.134615 0.0673077
The matrix is:
0.0206044 -0.02587 -0.00434982
-0.00274725 0.147894 -0.0105311
-0.00137363 -0.00938645 0.0780678
The vector is:
-0.00961538 0.134615 0.0673077
The matrix is:
4 3 15
7 5 0
10 4 18
The matrix is:
5 0
4 18
66
0
The matrix is:
11 15 19
The matrix is:
5 6 7
5 6 7
The matrix is:
4 0 0
0 5 0
0 0 18
The matrix is:
4 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 5 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 18
9
3
The vector is:
1 2 3
7.33333
7
6.66667
The vector is:
1 1 1 -1 1 -1
2.73861
The matrix is:
1 3
2 4
0 5
2 3
The matrix is:
1 2 0
3 4 5
3
0
0
5
2
The vector is:
1 3 2 4 0 5
The vector is:
0 3
The vector is:
2 5
MIN: 0 MAX: 5
3
The matrix is:
2 2 5
2 2 5
5 5 12.5
The matrix is:
2 2 5
2 2 5
5 5 12.5
The matrix is:
1 1 1
1 1 1
1 1 1
The vector is:
0 1 2 2 3 4 5 4 4 1
The vector is:
9 6 7 5 1
The vector is:
0 8 1 3 2 4 7 5 6
The matrix is:
1 2 0
3 4 5
The vector is:
1 3 3
The matrix is:
1 6 10
4 10 15
The vector is:
1 2 0
The matrix is:
1 6 0
3 24 0
6194
Process returned 0 (0x0) execution time : 4.891 s
Press any key to continue.
#include
using std::cout;
using std::endl;
using std::cin;
using std::istream;
using std::ios_base;
#include
using std::vector;
#include
using std::initializer_list;
#include
using std::for_each;
using std::sort;
using std::copy;
using std::find;
#include
using std::runtime_error;
#include
using std::srand;
using std::rand;
using std::system;
#include
using std::time;
#include
using std::sqrt;
#include
using std::pair;
using std::make_pair;
#include
using std::setw;
#include
using std::multiset;
using std::set;
#include
using std::unordered_map;
using std::unordered_multimap;
#include
using std::multimap;
using std::map;
class matrix;
class c
{
private:
vector c_val;
size_t dim;
public:
c():c_val(vector()),dim(0){}
c(initializer_listlst);
c(const vector&c_v):c_val(c_v),dim(c_v.size()){}
c(const matrix&m);
void out(size_t n=1)const{cout<<"The vector is:\n";for(size_t i=0;i<<setw(n)<<c_val[i]<<" ";cout<<endl<<endl;}
double operator[](size_t i)const{return c_val[i];}
double& operator[](size_t i){return c_val[i];}
size_t re_dim()const{return dim;}
friend c operator+(const
4000
c&a,const c&b);
friend c operator-(const c&a);
friend c operator-(const c&a,const c&b);
friend c rep(const double&d,size_t n);
friend c operator+(const c&a,const double&d);
friend c operator-(const c&a,const double&d);
friend c operator+(const double&d,const c&a);
friend c operator-(const double&d,const c&a);
friend c operator*(const double&d,const c&a);
friend c operator*(const c&a,const double&d);
friend c operator/(const double&d,const c&a);
friend c operator/(const c&a,const double&d);
friend c operator*(const c&a,const c&b);
friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}
friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}
vector const_re_c_val()const{return c_val;}
};
class matrix
{
private:
vector m_val;
size_t rdim;
size_t cdim;
public:
matrix():m_val(vector()),rdim(0),cdim(0){}
matrix(const vector&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}
matrix(size_t m,size_t n,const double&d):m_val(vector(m,vector(n,d))),rdim(m),cdim(n){}
matrix(const c&v,size_t r,size_t c);
matrix(const c&v);
c operator[](size_t i)const{return m_val[i];}
c& operator[](size_t i){return m_val[i];}
void out(size_t n=3)const;
c col(size_t i)const;
void col_change(size_t i,const c&a);
friend matrix operator+(const matrix&a,const matrix&b);
friend matrix operator-(const matrix&a);
friend matrix operator-(const matrix&a,const matrix&b);
friend matrix operator+(const matrix&a,const double&d);
friend matrix operator+(const double&d,const matrix&a);
friend matrix operator-(const matrix&a,const double&d);
friend matrix operator-(const double&b,const matrix&a);
friend matrix operator*(const matrix&a,const double&d);
friend matrix operator*(const double&d,const matrix&a);
friend matrix operator/(const matrix&a,const double&d);
friend matrix operator/(const double&d,const matrix&a);
void Swap(size_t i,size_t j)
{c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}
void col_Swap(size_t i,size_t j)
{
auto temp=this->col(i);
this->col_change(i,this->col(j));
this->col_change(j,temp);
}
friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}
friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}
pair re_dim()const{return make_pair(rdim,cdim);}
friend c operator%(const matrix&m,const c&v);
friend matrix operator%(const matrix&a,const matrix&b);
friend matrix operator%(const c&v,const matrix&m);
};
//类c的有关定义
c::c(initializer_listlst)
{
vectorout;
for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});
c_val=out;
dim=out.size();
}
c operator+(const c&a,const c&b)
{
size_t MAX=b.re_dim();
if(a.re_dim()>b.re_dim())
MAX=a.re_dim();
vectortemp0(MAX,0);
vectortemp1(MAX,0);
vectortemp2(MAX,0);
copy(a.c_val.begin(),a.c_val.end(),temp0.begin());
copy(b.c_val.begin(),b.c_val.end(),temp1.begin());
for(size_t i=0;i
temp2[i]=temp0[i]+temp1[i];
return c(temp2);
}
c operator-(const c&a)
{
c out(a);
for(size_t i=0;i
out[i]=-out[i];
return out;
}
c operator-(const c&a,const c&b)
{
return (a+(-b));
}
c rep(const double&d,size_t n)
{
return c(vector(n,d));
}
c operator+(const c&a,const double&d)
{
return a+rep(d,a.re_dim());
}
c operator-(const c&a,const double&d)
{
return (a+(-d));
}
c operator+(const double&d,const c&a)
{
return a+d;
}
c operator-(const double&d,const c&a)
{
return (-a)+d;
}
c operator*(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i
out[i]*=d;
return out;
}
c operator*(const c&a,const double&d)
{
return d*a;
}
c operator/(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i
out[i]=d/a[i];
return out;
}
c operator/(const c&a,const double&d)
{
c out(a);
for(size_t i=0;i
out[i]=a[i]/d;
return out;
}
c::c(const matrix&m)
{
vectortemp;
for(size_t j=0;j
for(size_t i=0;i
temp.push_back(m[i][j]);
c_val=temp;
dim=temp.size();
}
c operator*(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
c out(a);
for(size_t i=0;i
out[i]*=b[i];
return out;
}
//类matrix的有关定义
void matrix::out(size_t n)const
{
cout<<"The matrix is:\n";
for(size_t i=0;i
{
for(size_t j=0;j
cout<<setw(n)<<(*this)[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
c matrix::col(size_t i)const
{
c out(rep(0,rdim));
for(size_t j=0;j
out[j]=(*this)[j][i];
return out;
}
void matrix::col_change(size_t i,const c&a)
{
for(size_t j=0;j
(*this)[j][i]=a[j];
}
matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)
{
matrix out(r,c,0);
for(size_t i=0;i
out[i%r][i/r]=v[i];
(this->m_val)=out.m_val;
}
matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)
{
matrix out(v.re_dim(),1,0);
for(size_t i=0;i
out[i][0]=v[i];
this->m_val=out.m_val;
}
matrix operator+(const matrix&a,const matrix&b)
{
size_t r_MAX=a.rdim,c_MAX=a.cdim;
if(a.rdim
r_MAX=b.rdim;
if(a.cdim
c_MAX=b.cdim;
matrix tempa(r_MAX,c_MAX,0);
for(size_t i=0;i
for(size_t j=0;j
tempa[i][j]=a[i][j];
matrix tempb(r_MAX,c_MAX,0);
for(size_t i=0;i
for(size_t j=0;j
tempb[i][j]=b[i][j];
matrix tempout(r_MAX,c_MAX,0);
for(size_t i=0;i
for(size_t j=0;j
tempout[i][j]=tempa[i][j]+tempb[i][j];
return tempout;
}
matrix operator-(const matrix&a)
{
matrix out(a);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=-out[i][j];
return out;
}
matrix operator-(const matrix&a,const matrix&b)
{
return ((-b)+a);
}
matrix operator+(const matrix&a,const double&d)
{
return a+matrix(a.rdim,a.cdim,d);
}
matrix operator+(const double&d,const matrix&a)
{
return a+d;
}
matrix operator-(const matrix&a,const double&d)
{
return a+(-d);
}
matrix operator-(const double&b,const matrix&a)
{
return (-a)+b;
}
matrix operator*(const matrix&a,const double&d)
{
matrix out(a);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]*=d;
return out;
}
matrix operator*(const double&d,const matrix&a)
{
return a*d;
}
matrix operator/(const matrix&a,const double&d)
{
return a*(1/d);
}
matrix operator/(const double&d,const matrix&a)
{
matrix out(a.rdim,a.cdim,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=d/a[i][j];
return out;
}
//基本运算函数
matrix diag(const c&v)
{
matrix out(v.re_dim(),v.re_dim(),0);
for(size_t i=0;i
out[i][i]=v[i];
return out;
}
matrix diag(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i
out[i][i]=m[i][i];
return out;
}
double inter_product(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
double out=0;
for(size_t i=0;i
out+=(a[i]*b[i]);
return out;
}
c operator%(const matrix&m,const c&v)
{
if(m.re_dim().second!=v.re_dim())
throw runtime_error("dim error!\n");
c out(rep(0,m.re_dim().first));
for(size_t i=0;i
out[i]=inter_product(m[i],v);
return out;
}
//同下
double prod(const c&v)
{
double out=1;
for(size_t i=0;i
out*=v[i];
return out;
}
//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)
double sum(const c&v)
{
double out=v[0];
for(size_t i=1;i
out+=v[i];
return out;
}
matrix operator%(const c&v,const matrix&m)
{
if(v.re_dim()!=m.re_dim().first)
throw runtime_error("dim error!\n");
matrix out(1,m.re_dim().second,0);
for(size_t i=0;i
out[0][i]=inter_product(v,m.col(i));
return out;
}
matrix outer_product(const c&a,const c&b)
{
matrix out(a.re_dim(),b.re_dim(),0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=a[i]*b[j];
return out;
}
size_t length(const c&v)
{
return v.re_dim();
}
c sort(const c&v)
{
vectorout;
for(size_t i=0;i
out.push_back(v[i]);
sort(out.begin(),out.end());
return c(out);
}
double mean(const c&v)
{
return sum(v)/length(v);
}
//其转换版本返回矩阵转为向量的方差
double var(const c&v)
{
return mean((v-mean(v))*(v-mean(v)));
}
c connect(initializer_list lst)
{
vectortemp;
for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});
return c(temp);
}
//标准差
template
double STD(const any&a)
{
return sqrt(var(a)*(length(a))/(length(a)-1));
}
matrix rbind(initializer_listlst)
{
size_t r=0;
auto a=lst.begin();
while(a!=lst.end())
{
a++;
r++;
}
size_t co=0;
for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});
matrix out(r,co,0);
size_t i=0;
a=lst.begin();
while(a!=lst.end())
{
vectortemp(co,0);
for(size_t i=0;iconst_re_c_val().size();i++)
temp[i]=a->operator[](i);
out[i]=c(temp);
i++;
a++;
}
return out;
}
size_t dim(const c&v)
{
return v.re_dim();
}
pair dim(const matrix&m)
{
return m.re_dim();
}
double min(const c&v)
{
double temp=v[0];
for(size_t i=1;i
if(v[i]
temp=v[i];
return temp;
}
double max(const c&v)
{
double temp=v[0];
for(size_t i=1;i
if(v[i]>temp)
temp=v[i];
return temp;
}
c pmin(initializer_listlst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i
{
temp[i]=min(*a);
a++;
}
return temp;
}
c pmax(initializer_listlst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i
{
temp[i]=max(*a);
a++;
}
return temp;
}
template
pair range(const any&v)
{
cout<<"MIN: "<<min(v)<<" MAX: "<<max(v)<<endl;
return make_pair(min(v),max(v));
}
double median(const c&v)
{
return sort(v)[length(v)/2];
}
double var(const c&a,const c&b)
{
if(length(a)!=length(b))
throw runtime_error("dim error!\n");
c t1(a-mean(a)),t2(b-mean(b));
return inter_product(t1,t2)/(length(a)-1);
}
double cov(const c&a,const c&b)
{
return var(a,b);
}
double cor(const c&a,const c&b)
{
return var(a,b)/STD(a)/STD(b);
}
matrix var(const matrix&a,const matrix&b)
{
if(a.re_dim().first!=b.re_dim().first)
throw runtime_error("dim error!\n");
matrix temp(a.re_dim().second,b.re_dim().second,0);
for(size_t i=0;i
for(size_t j=0;j
temp[i][j]=var(a.col(i),b.col(j));
return temp;
}
matrix cov(const matrix&a,const matrix&b)
{
return var(a,b);
}
matrix cor(const matrix&a,const matrix&b)
{
matrix temp=var(a,b);
for(size_t i=0;i
for(size_t j=0;j
temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));
return temp;
}
matrix cov(const matrix&a)
{
return var(a,a);
}
matrix cor(const matrix&a)
{
return cor(a,a);
}
//由于关联容器的删除方法利用误差排序的方法
//利用顺序容器迭代器的删除方法有困难
c rank(const c&v)
{
vectortemp00(v.const_re_c_val());
for(size_t i=0;i
temp00[i]+=(0.00000001*i);
vectortemp0(temp00);
multisettemp2(temp0.begin(),temp0.end());
vectortemp1(v.re_dim());
size_t i=0,j=0;
double pre_val=min(c(vector(temp2.begin(),temp2.end())));
vector::iterator s0;
while(i
{
double val=min(c(vector(temp2.begin(),temp2.end())));
s0=find(temp0.begin(),temp0.end(),val);
if(sqrt((pre_val-val)*(pre_val-val))>0.000001)
j++;
temp1[s0-temp0.begin()]=j;
temp2.erase(val);
i++;
pre_val=val;
}
return c(temp1);
}
c rev(const c&v)
{
vector temp;
auto a=v.const_re_c_val().rbegin();
while(a!=v.const_re_c_val().rend())
{
temp.push_back(*a);
a++;
}
return c(temp);
}
//向量先排序再输出排序前的下标值
c order(const c&v)
{
unordered_multimapm0;
for(size_t i=0;i
m0.insert(make_pair(v[i],i));
multimapm1(m0.begin(),m0.end());
c temp(rep(0,length(v)));
auto a=m1.begin();
size_t i=0;
while(a!=m1.end())
{
temp[i]=a->second;
i++;
a++;
}
return temp;
}
c cumsum(const c&v)
{
c temp(v);
for(size_t i=1;i
temp[i]+=temp[i-1];
return temp;
}
c cumprod(const c&v)
{
c temp(v);
for(size_t i=1;i
temp[i]*=temp[i-1];
return temp;
}
matrix cumsum(const matrix&m)
{
c temp(m);
for(size_t i=1;i
temp[i]+=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
matrix cumprod(const matrix&m)
{
c temp(m);
for(size_t i=1;i
temp[i]*=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
//矩阵运算函数
//初等行变换化上三角求解行列式
double det(const matrix&m)
{
auto f=[](const matrix&mm)->pair
{
if(mm.re_dim().first!=mm.re_dim().second)
throw runtime_error("dim error!\n");
size_t MIN=mm.re_dim().first;
double mult=1;
size_t col=0;
size_t bar=0;
size_t b_c=0;
matrix t(mm);
while(b_c
{
for(size_t i=bar+1;i
{
if(t[i][col]==0)
i++;
else
{
t.Swap(i,bar);
mult*=(-1);
break;
}
}
if(bar!=MIN+1)
{
if(t[bar][col]!=0)
{
for(size_t i=bar+1;i
{
t[bar][i]=t[bar][i]/t[bar][col];
if(i==bar+1)
mult*=t[bar][col];
}
if(col!=MIN-1)
t[bar][col]=1;
}
for(size_t i=bar+1;i
{
t[i]=t[i]-t[i][col]*t[ba
ee37
r];
}
}
col++;
bar++;
b_c=bar;
if(col>bar);
b_c=col;
}
return make_pair(t,mult);
};
pair tt=f(m);
return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));
}
//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c solve(const matrix&m,const c&v)
{
c temp(rep(0,v.re_dim()));
double low=det(m);
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
if(low==0)
throw runtime_error("singular!\n");
for(size_t i=0;i
{
matrix temp0(m);
temp0.col_change(i,v);
temp[i]=det(temp0)/low;
}
return temp;
}
double tr(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));
}
matrix T(const matrix&m)
{
matrix out(m.re_dim().second,m.re_dim().first,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=m[j][i];
return out;
}
matrix cbind(initializer_listlst)
{
return T(rbind(lst));
}
matrix operator%(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first,a.re_dim().second,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=inter_product(a[i],b.col(j));
return out;
}
double F_norm(const matrix&m)
{
return sqrt(tr(T(m)%m));
}
//a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善
//只能求实值
//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中
//可用于求解但效率低
c eigen(const matrix&m,const double&a,const double& b,const double &d=0.0001,const double &e=1)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
vectortemp;
double pre_beg=100000000;
for(double beg=a;beg
{
matrix mm0(diag(rep(beg,m.re_dim().first))-m);
double s0=det(mm0);
if(sqrt(s0*s0)
{
if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)
{
temp.push_back(beg);
pre_beg=beg;
cout<<"Done!\n";
cout<<beg<<endl;//eigenvalue
}
}
}
return c(temp);
}
//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）
c eigen(const matrix&m,const double&d=0.1)
{
double dd=d;
double ra=F_norm(m);
double trr=tr(m);
c temp;
double r=0;
do
{
dd*=0.1;
temp=eigen(m,-ra,ra,dd);
if(temp.re_dim()==m.re_dim().first)
break;
r=sum(temp)-trr;
}
while(sqrt((r*r)>1));
//与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略
//此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握
//只适用于所有特征值为实值的情况
//如有复特征值，只输出实数的并循环不能结束
return temp;
}
matrix M(const matrix&m,size_t k,size_t l)
{
matrix temp(m);
set s1;
for(size_t i=0;i
s1.insert(i);
set s2;
for(size_t i=0;i
s2.insert(i);
s1.erase(k);
s2.erase(l);
vectortemp00;
for(auto j:s2)
for(auto i:s1)
temp00.push_back(temp[i][j]);
return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);
}
matrix A_accompany(const matrix&m)
{
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i
for(size_t j=0;j
out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);
//代数余子式
return T(out);
}
//求方阵的逆矩阵（利用伴随矩阵方法）
matrix solve(const matrix&m)
{
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return A_accompany(m)/de;
}
//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆
c solve_eq(const matrix&m,const c&v)
{
c out;
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return solve(m)%v;
}
测试程序：
int main()
{
//测试程序：
c c0({35.1,12});
vectorv0;
for(size_t i=0;i<10;i++)
v0.push_back(i);
c c1(v0);
c0.out();
c1.out();
cout<<c0[1]<<endl;
c0[1]=10;
cout<<c0[1]<<endl;
(c0+c1).out();
(-c0).out();
(c0-c1).out();
rep(1,10).out();
(c0+rep(1,10)).out();
(c0-rep(1,10)).out();
(rep(1,10)+1).out();
(rep(1,10)-1).out();
(10/rep(1,10)).out();
(rep(1,10)/10).out();
cout<<((10/rep(1,10))==(rep(1,10)/10))<<endl;
matrix m0=matrix(3,4,10);
m0[1]=rep(4,4);
m0.out();
m0.col(2).out();
m0.col_change(0,c({1,2,3}));
m0.out();
matrix m1=(matrix(c({1,5,7,8,9,0,3,4,5}),3,3)+matrix(c({3,5,7,9}),2,2));
m1.out();
(-m1).out();
(m1-matrix(rep(1,10),5,2)).out();
(m1+1).out();
(m1*0.6).out();
(m1-1).out();
(1-m1).out();
(0.6*m1).out();
(1/m1).out();
(m1/10).out();
m1.out();
m1.Swap(1,2);
m1.out();
m1.col_Swap(1,2);
m1.out();
matrix m2=m1;
m1.Swap(0,1);
cout<<(m1==m2)<<endl;
cout<<m1.re_dim().first<<" "<<m1.re_dim().second<<endl;
matrix A(c({1,2,1,4,7,1,5,2,1}),3,3);
cout<<det((A+diag(rep(10,A.re_dim().first))))<<endl;
solve(matrix(c({30,7.3,4.2,0,0,6,8,7,0}),3,3),c({11,12,13})).out();
m1.out();
cout<<F_norm(m1)<<endl;
eigen(matrix(c({50,1,1,9,7,1,4,1,13}),3,3)).out();
solve(matrix(c({50,1,1,9,7,1,4,1,13}),3,3),c({1,1,1})).out();
solve(matrix(c({50,1,1,9,7,1,4,1,13}),3,3)).out();
solve_eq(matrix(c({50,1,1,9,7,1,4,1,13}),3,3),c({1,1,1})).out();
m2.out();
M(m2,0,0).out();
cout<<sum(m2)<<endl;
cout<<prod(m2)<<endl;
(c({1,1})%matrix(c({5,6,7,8,9,10}),2,3)).out();
outer_product(c({1,1}),c({5,6,7})).out();
diag(m2).out();
diag(c(diag(m2))).out();
cout<<length(m2)<<endl;
cout<<length(m2[0])<<endl;
sort(c({3,2,1})).out();
cout<<mean(m2)<<endl;
cout<<mean(m2.col(0))<<endl;
cout<<var(matrix(c({1,2,3,4,5,6,7,8,9}),3,3))<<endl;
connect({c({1,1,1}),c({-1,1,-1})}).out();
cout<<STD(matrix(c({1,2,3,4,5,6,7,8,9}),3,3))<<endl;
matrix m3=rbind({c({1,2}),c({3,4,5})});
cbind({c({1,2}),c({3,4,5})}).out();
cout<<dim(m3).first<<" "<<dim(m3).second<<endl;
m3.out();
cout<<dim(m3[0])<<endl;
cout<<min(m3)<<endl;
cout<<min(m3[0])<<endl;
cout<<max(m3)<<endl;
cout<<max(m3[0])<<endl;
c(m3).out();
pmin({m3[0],m3[1]}).out();
pmax({m3[0],m3[1]}).out();
range(m3);
cout<<median(m3)<<endl;
var(m3,m3).out();
cov(m3).out();
cor(m3).out();
rank(c({0,3,4,4,5,6,7,6,6,3})).out();
rev(c({1,5,7,6,9})).out();
order(c({0,3,4,4,5,6,7,6,3})).out();
m3.out();
cumsum(m3[0]).out();
cumsum(m3).out();
cumprod(m3[0]).out();
cumprod(m3).out();
cout<<det(matrix(c({50,7,8,9,10,4,3,2,9,4,8,7,1,1,1,10}),4,4))<<endl;
return 0;
}
运行结果：
The vector is:
35.1 12
The vector is:
0 1 2 3 4 5 6 7 8 9
12
10
The vector is:
35.1 11 2 3 4 5 6 7 8 9
The vector is:
-35.1 -10
The vector is:
35.1 9 -2 -3 -4 -5 -6 -7 -8 -9
The vector is:
1 1 1 1 1 1 1 1 1 1
The vector is:
36.1 11 1 1 1 1 1 1 1 1
The vector is:
34.1 9 -1 -1 -1 -1 -1 -1 -1 -1
The vector is:
2 2 2 2 2 2 2 2 2 2
The vector is:
0 0 0 0 0 0 0 0 0 0
The vector is:
10 10 10 10 10 10 10 10 10 10
The vector is:
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
0
The matrix is:
10 10 10 10
4 4 4 4
10 10 10 10
The vector is:
10 4 10
The matrix is:
1 10 10 10
2 4 4 4
3 10 10 10
The matrix is:
4 15 3
10 18 4
7 0 5
The matrix is:
-4 -15 -3
-10 -18 -4
-7 -0 -5
The matrix is:
3 14 3
9 17 4
6 -1 5
-1 -1 0
-1 -1 0
The matrix is:
5 16 4
11 19 5
8 1 6
The matrix is:
2.4 9 1.8
6 10.8 2.4
4.2 0 3
The matrix is:
3 14 2
9 17 3
6 -1 4
The matrix is:
-3 -14 -2
-9 -17 -3
-6 1 -4
The matrix is:
2.4 9 1.8
6 10.8 2.4
4.2 0 3
The matrix is:
0.25 0.0666667 0.333333
0.1 0.0555556 0.25
0.142857 inf 0.2
The matrix is:
0.4 1.5 0.3
1 1.8 0.4
0.7 0 0.5
The matrix is:
4 15 3
10 18 4
7 0 5
The matrix is:
4 15 3
7 0 5
10 18 4
The matrix is:
4 3 15
7 5 0
10 4 18
0
3 3
1880
The vector is:
-0.12533 2.2544 1.84499
The matrix is:
7 5 0
4 3 15
10 4 18
27.6405
Done!
6.67574
Done!
6.67374
Done!
12.9967
Done!
50.3227
The vector is:
6.67374 12.9967 50.3227
The vector is:
-0.00961538 0.134615 0.0673077
The matrix is:
0.0206044 -0.02587 -0.00434982
-0.00274725 0.147894 -0.0105311
-0.00137363 -0.00938645 0.0780678
The vector is:
-0.00961538 0.134615 0.0673077
The matrix is:
4 3 15
7 5 0
10 4 18
The matrix is:
5 0
4 18
66
0
The matrix is:
11 15 19
The matrix is:
5 6 7
5 6 7
The matrix is:
4 0 0
0 5 0
0 0 18
The matrix is:
4 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 5 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 18
9
3
The vector is:
1 2 3
7.33333
7
6.66667
The vector is:
1 1 1 -1 1 -1
2.73861
The matrix is:
1 3
2 4
0 5
2 3
The matrix is:
1 2 0
3 4 5
3
0
0
5
2
The vector is:
1 3 2 4 0 5
The vector is:
0 3
The vector is:
2 5
MIN: 0 MAX: 5
3
The matrix is:
2 2 5
2 2 5
5 5 12.5
The matrix is:
2 2 5
2 2 5
5 5 12.5
The matrix is:
1 1 1
1 1 1
1 1 1
The vector is:
0 1 2 2 3 4 5 4 4 1
The vector is:
9 6 7 5 1
The vector is:
0 8 1 3 2 4 7 5 6
The matrix is:
1 2 0
3 4 5
The vector is:
1 3 3
The matrix is:
1 6 10
4 10 15
The vector is:
1 2 0
The matrix is:
1 6 0
3 24 0
6194
Process returned 0 (0x0) execution time : 4.891 s
Press any key to continue.
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- [c/c++]指针（4）
- Ugly Number II
- SCPPO（九）:性能优化之停工装置按照开工装置更新价格
- c++ string的实现。
- #include<head.h>和#include "head.h"有什么区别
- C++与C语言的不同
- PAT 1054. 求平均值
- C++判断字符串是否所有字符全都不同
- 【C++优先队列】
- C++ PP Chapter IX 内存模型和名称空间
- vector使用需要注意的一些问题
- [C-C++]获取系统时间
- 一起学习c++11——c++11中的新增的容器
- c和指针 前三章知识点以及常见的问题和程序分析
- interbase C++Builder 简单例子
- [C-C++]控制台用定时器
- 华为OJ——字符串匹配
- 华为OJ——整形数组合并
- 【C++专题】static_cast, dynamic_cast, const_cast探讨
- 华为OJ——计算字符串的相似度