C++实现R语言向量化运算(向量类:c 矩阵类:matrix)2015.9.11

来源:互联网 发布:python难找工作 编辑:程序博客网 时间:2024/05/22 12:45

类源代码:

[cpp] view plain copy
print?在CODE上查看代码片派生到我的代码片
  1. #include  
  2. using std::cout;  
  3. using std::endl;  
  4. using std::cin;  
  5. using std::istream;  
  6. using std::ios_base;  
  7. #include  
  8. using std::vector;  
  9. #include  
  10. using std::initializer_list;  
  11. #include  
  12. using std::for_each;  
  13. using std::sort;  
  14. using std::copy;  
  15. using std::find;  
  16. #include  
  17. using std::runtime_error;  
  18. #include  
  19. using std::srand;  
  20. using std::rand;  
  21. using std::system;  
  22. #include  
  23. using std::time;  
  24. #include  
  25. using std::sqrt;  
  26. #include  
  27. using std::pair;  
  28. using std::make_pair;  
  29. #include  
  30. using std::setw;  
  31. #include  
  32. using std::multiset;  
  33. using std::set;  
  34. #include  
  35. using std::unordered_map;  
  36. using std::unordered_multimap;  
  37. #include  
  38. using std::multimap;  
  39. using std::map;   
  40. class matrix;  
  41. class c  
  42. {  
  43. private:  
  44.     vector c_val;  
  45.     size_t dim;  
  46. public:  
  47.     c():c_val(vector()),dim(0){}  
  48.     c(initializer_listlst);  
  49.     c(const vector&c_v):c_val(c_v),dim(c_v.size()){}  
  50.     c(const matrix&m);  
  51.     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;}  
  52.     double operator[](size_t i)const{return c_val[i];}  
  53.     double& operator[](size_t i){return c_val[i];}  
  54.     size_t re_dim()const{return dim;}  
  55.     friend c operator+(const c&a,const c&b);  
  56.     friend c operator-(const c&a);  
  57.     friend c operator-(const c&a,const c&b);  
  58.     friend c rep(const double&d,size_t n);  
  59.     friend c operator+(const c&a,const double&d);  
  60.     friend c operator-(const c&a,const double&d);  
  61.     friend c operator+(const double&d,const c&a);  
  62.     friend c operator-(const double&d,const c&a);  
  63.     friend c operator*(const double&d,const c&a);  
  64.     friend c operator*(const c&a,const double&d);  
  65.     friend c operator/(const double&d,const c&a);  
  66.     friend c operator/(const c&a,const double&d);  
  67.     friend c operator*(const c&a,const c&b);  
  68.     friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}  
  69.     friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}  
  70.     vector const_re_c_val()const{return c_val;}  
  71. };  
  72. class matrix  
  73. {  
  74. private:  
  75.     vector m_val;  
  76.     size_t rdim;  
  77.     size_t cdim;  
  78. public:  
  79.     matrix():m_val(vector()),rdim(0),cdim(0){}  
  80.     matrix(const vector&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}  
  81.     matrix(size_t m,size_t n,const double&d):m_val(vector(m,vector(n,d))),rdim(m),cdim(n){}  
  82.     matrix(const c&v,size_t r,size_t c);  
  83.     matrix(const c&v);  
  84.     c operator[](size_t i)const{return m_val[i];}  
  85.     c& operator[](size_t i){return m_val[i];}  
  86.     void out(size_t n=3)const;  
  87.     c col(size_t i)const;  
  88.     void col_change(size_t i,const c&a);  
  89.     friend matrix operator+(const matrix&a,const matrix&b);  
  90.     friend matrix operator-(const matrix&a);  
  91.     friend matrix operator-(const matrix&a,const matrix&b);  
  92.     friend matrix operator+(const matrix&a,const double&d);  
  93.     friend matrix operator+(const double&d,const matrix&a);  
  94.     friend matrix operator-(const matrix&a,const double&d);  
  95.     friend matrix operator-(const double&b,const matrix&a);  
  96.     friend matrix operator*(const matrix&a,const double&d);  
  97.     friend matrix operator*(const double&d,const matrix&a);  
  98.     friend matrix operator/(const matrix&a,const double&d);  
  99.     friend matrix operator/(const double&d,const matrix&a);  
  100.     void Swap(size_t i,size_t j)  
  101.     {c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}  
  102.     void col_Swap(size_t i,size_t j)  
  103.     {  
  104.         auto temp=this->col(i);  
  105.         this->col_change(i,this->col(j));  
  106.         this->col_change(j,temp);  
  107.     }  
  108.     friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}  
  109.     friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}  
  110.     pair re_dim()const{return make_pair(rdim,cdim);}  
  111.     friend c operator%(const matrix&m,const c&v);  
  112.     friend matrix operator%(const matrix&a,const matrix&b);  
  113.     friend matrix operator%(const c&v,const matrix&m);  
  114. };  
  115.    
  116. //类c的有关定义  
  117. c::c(initializer_listlst)  
  118. {  
  119.     vectorout;  
  120.     for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});  
  121.     c_val=out;  
  122.     dim=out.size();  
  123. }  
  124. c operator+(const c&a,const c&b)  
  125. {  
  126.     size_t MAX=b.re_dim();  
  127.     if(a.re_dim()>b.re_dim())  
  128.         MAX=a.re_dim();  
  129.     vectortemp0(MAX,0);  
  130.     vectortemp1(MAX,0);  
  131.     vectortemp2(MAX,0);  
  132.     copy(a.c_val.begin(),a.c_val.end(),temp0.begin());  
  133.     copy(b.c_val.begin(),b.c_val.end(),temp1.begin());  
  134.     for(size_t i=0;i  
  135.         temp2[i]=temp0[i]+temp1[i];  
  136.     return c(temp2);  
  137. }  
  138. c operator-(const c&a)  
  139. {  
  140.     c out(a);  
  141.     for(size_t i=0;i  
  142.         out[i]=-out[i];  
  143.     return out;  
  144. }  
  145. c operator-(const c&a,const c&b)  
  146. {  
  147.     return (a+(-b));  
  148. }  
  149. c rep(const double&d,size_t n)  
  150. {  
  151.     return c(vector(n,d));  
  152. }  
  153. c operator+(const c&a,const double&d)  
  154. {  
  155.     return a+rep(d,a.re_dim());  
  156. }  
  157. c operator-(const c&a,const double&d)  
  158. {  
  159.     return (a+(-d));  
  160. }  
  161. c operator+(const double&d,const c&a)  
  162. {  
  163.     return a+d;  
  164. }  
  165. c operator-(const double&d,const c&a)  
  166. {  
  167.     return (-a)+d;  
  168. }  
  169. c operator*(const double&d,const c&a)  
  170. {  
  171.     c out(a);  
  172.     for(size_t i=0;i  
  173.         out[i]*=d;  
  174.     return out;  
  175. }  
  176. c operator*(const c&a,const double&d)  
  177. {  
  178.     return d*a;  
  179. }  
  180. c operator/(const double&d,const c&a)  
  181. {  
  182.     c out(a);  
  183.     for(size_t i=0;i  
  184.         out[i]=d/a[i];  
  185.     return out;  
  186. }  
  187. c operator/(const c&a,const double&d)  
  188. {  
  189.     c out(a);  
  190.     for(size_t i=0;i  
  191.         out[i]=a[i]/d;  
  192.     return out;  
  193. }  
  194. c::c(const matrix&m)  
  195. {  
  196.     vectortemp;  
  197.   
  198.         for(size_t j=0;j  
  199.                 for(size_t i=0;i  
  200.         temp.push_back(m[i][j]);  
  201.     c_val=temp;  
  202.     dim=temp.size();  
  203. }  
  204. c operator*(const c&a,const c&b)  
  205. {  
  206.     if(a.re_dim()!=b.re_dim())  
  207.         throw runtime_error(“dim error!\n”);  
  208.     c out(a);  
  209.     for(size_t i=0;i  
  210.         out[i]*=b[i];  
  211.     return out;  
  212. }  
  213.    
  214. //类matrix的有关定义  
  215. void matrix::out(size_t n)const  
  216. {  
  217.     cout<<”The matrix is:\n”;  
  218.     for(size_t i=0;i  
  219.     {  
  220.         for(size_t j=0;j  
  221.             cout<<setw(n)<<(*this)[i][j]<<“ ”;  
  222.         cout<<endl;  
  223.     }  
  224.     cout<<endl;  
  225. }  
  226. c matrix::col(size_t i)const  
  227. {  
  228.     c out(rep(0,rdim));  
  229.     for(size_t j=0;j  
  230.         out[j]=(*this)[j][i];  
  231.     return out;  
  232. }  
  233. void matrix::col_change(size_t i,const c&a)  
  234. {  
  235.     for(size_t j=0;j  
  236.         (*this)[j][i]=a[j];  
  237. }  
  238. matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)  
  239. {  
  240.     matrix out(r,c,0);  
  241.     for(size_t i=0;i  
  242.         out[i%r][i/r]=v[i];  
  243.     (this->m_val)=out.m_val;  
  244. }  
  245. matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)  
  246. {  
  247.     matrix out(v.re_dim(),1,0);  
  248.     for(size_t i=0;i  
  249.         out[i][0]=v[i];  
  250.     this->m_val=out.m_val;  
  251. }  
  252. matrix operator+(const matrix&a,const matrix&b)  
  253. {  
  254.     size_t r_MAX=a.rdim,c_MAX=a.cdim;  
  255.     if(a.rdim  
  256.         r_MAX=b.rdim;  
  257.     if(a.cdim  
  258.         c_MAX=b.cdim;  
  259.     matrix tempa(r_MAX,c_MAX,0);  
  260.     for(size_t i=0;i  
  261.         for(size_t j=0;j  
  262.         tempa[i][j]=a[i][j];  
  263.     matrix tempb(r_MAX,c_MAX,0);  
  264.     for(size_t i=0;i  
  265.         for(size_t j=0;j  
  266.         tempb[i][j]=b[i][j];  
  267.     matrix tempout(r_MAX,c_MAX,0);  
  268.     for(size_t i=0;i  
  269.         for(size_t j=0;j  
  270.         tempout[i][j]=tempa[i][j]+tempb[i][j];  
  271.     return tempout;  
  272. }  
  273. matrix operator-(const matrix&a)  
  274. {  
  275.     matrix out(a);  
  276.     for(size_t i=0;i  
  277.         for(size_t j=0;j  
  278.         out[i][j]=-out[i][j];  
  279.     return out;  
  280. }  
  281. matrix operator-(const matrix&a,const matrix&b)  
  282. {  
  283.     return ((-b)+a);  
  284. }  
  285. matrix operator+(const matrix&a,const double&d)  
  286. {  
  287.     return a+matrix(a.rdim,a.cdim,d);  
  288. }  
  289. matrix operator+(const double&d,const matrix&a)  
  290. {  
  291.     return a+d;  
  292. }  
  293. matrix operator-(const matrix&a,const double&d)  
  294. {  
  295.     return a+(-d);  
  296. }  
  297. matrix operator-(const double&b,const matrix&a)  
  298. {  
  299.     return (-a)+b;  
  300. }  
  301. matrix operator*(const matrix&a,const double&d)  
  302. {  
  303.     matrix out(a);  
  304.     for(size_t i=0;i  
  305.         for(size_t j=0;j  
  306.         out[i][j]*=d;  
  307.     return out;  
  308. }  
  309. matrix operator*(const double&d,const matrix&a)  
  310. {  
  311.     return a*d;  
  312. }  
  313. matrix operator/(const matrix&a,const double&d)  
  314. {  
  315.     return a*(1/d);  
  316. }  
  317. matrix operator/(const double&d,const matrix&a)  
  318. {  
  319.     matrix out(a.rdim,a.cdim,0);  
  320.     for(size_t i=0;i  
  321.         for(size_t j=0;j  
  322.         out[i][j]=d/a[i][j];  
  323.     return out;  
  324. }  
  325.    
  326. //基本运算函数  
  327. matrix diag(const c&v)  
  328. {  
  329.     matrix out(v.re_dim(),v.re_dim(),0);  
  330.     for(size_t i=0;i  
  331.         out[i][i]=v[i];  
  332.     return out;  
  333. }  
  334. matrix diag(const matrix&m)  
  335. {  
  336.     if(m.re_dim().first!=m.re_dim().second)  
  337.         throw runtime_error(“dim error!\n”);  
  338.     matrix out(m.re_dim().first,m.re_dim().second,0);  
  339.     for(size_t i=0;i  
  340.         out[i][i]=m[i][i];  
  341.     return out;  
  342. }  
  343. double inter_product(const c&a,const c&b)  
  344. {  
  345.     if(a.re_dim()!=b.re_dim())  
  346.         throw runtime_error(“dim error!\n”);  
  347.     double out=0;  
  348.     for(size_t i=0;i  
  349.     out+=(a[i]*b[i]);  
  350.     return out;  
  351. }  
  352. c operator%(const matrix&m,const c&v)  
  353. {  
  354.     if(m.re_dim().second!=v.re_dim())  
  355.         throw runtime_error(“dim error!\n”);  
  356.     c out(rep(0,m.re_dim().first));  
  357.     for(size_t i=0;i  
  358.         out[i]=inter_product(m[i],v);  
  359.     return out;  
  360. }  
  361. //同下  
  362. double prod(const c&v)  
  363. {  
  364.     double out=1;  
  365.     for(size_t i=0;i  
  366.         out*=v[i];  
  367.     return out;  
  368. }  
  369. //对于matrix的求和由c复制构造函数的重载函数定义:c(const matrix&m)  
  370. double sum(const c&v)  
  371. {  
  372.     double out=v[0];  
  373.     for(size_t i=1;i  
  374.         out+=v[i];  
  375.     return out;  
  376. }  
  377. matrix operator%(const c&v,const matrix&m)  
  378. {  
  379.     if(v.re_dim()!=m.re_dim().first)  
  380.         throw runtime_error(“dim error!\n”);  
  381.     matrix out(1,m.re_dim().second,0);  
  382.     for(size_t i=0;i  
  383.         out[0][i]=inter_product(v,m.col(i));  
  384.     return out;  
  385. }  
  386. matrix outer_product(const c&a,const c&b)  
  387. {  
  388.     matrix out(a.re_dim(),b.re_dim(),0);  
  389.     for(size_t i=0;i  
  390.         for(size_t j=0;j  
  391.         out[i][j]=a[i]*b[j];  
  392.     return out;  
  393. }  
  394. size_t length(const c&v)  
  395. {  
  396.     return v.re_dim();  
  397. }  
  398. c sort(const c&v)  
  399. {  
  400.     vectorout;  
  401.     for(size_t i=0;i  
  402.         out.push_back(v[i]);  
  403.     sort(out.begin(),out.end());  
  404.     return c(out);  
  405. }  
  406. double mean(const c&v)  
  407. {  
  408.     return sum(v)/length(v);  
  409. }  
  410. //其转换版本返回矩阵转为向量的方差  
  411. double var(const c&v)  
  412. {  
  413.     return mean((v-mean(v))*(v-mean(v)));  
  414. }  
  415. c connect(initializer_list lst)  
  416. {  
  417.     vectortemp;  
  418.     for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});  
  419.     return c(temp);  
  420. }  
  421. //标准差  
  422. template  
  423. double STD(const any&a)  
  424. {  
  425.     return sqrt(var(a)*(length(a))/(length(a)-1));  
  426. }  
  427. matrix rbind(initializer_listlst)  
  428. {  
  429.     size_t r=0;  
  430.     auto a=lst.begin();  
  431.     while(a!=lst.end())  
  432.     {  
  433.         a++;  
  434.         r++;  
  435.     }  
  436.     size_t co=0;  
  437.     for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});  
  438.     matrix out(r,co,0);  
  439.     size_t i=0;  
  440.     a=lst.begin();  
  441.     while(a!=lst.end())  
  442.     {  
  443.         vectortemp(co,0);  
  444.         for(size_t i=0;iconst_re_c_val().size();i++)  
  445.             temp[i]=a->operator[](i);  
  446.         out[i]=c(temp);  
  447.         i++;  
  448.         a++;  
  449.     }  
  450.   
  451.     return out;  
  452. }  
  453. size_t dim(const c&v)  
  454. {  
  455.     return v.re_dim();  
  456. }  
  457. pair dim(const matrix&m)  
  458. {  
  459.     return m.re_dim();  
  460. }  
  461.   
  462. double min(const c&v)  
  463. {  
  464.     double temp=v[0];  
  465.     for(size_t i=1;i  
  466.         if(v[i]  
  467.         temp=v[i];  
  468.     return temp;  
  469. }  
  470. double max(const c&v)  
  471. {  
  472.     double temp=v[0];  
  473.     for(size_t i=1;i  
  474.         if(v[i]>temp)  
  475.         temp=v[i];  
  476.     return temp;  
  477. }  
  478. c pmin(initializer_listlst)  
  479. {  
  480.     c temp=rep(0,lst.size());  
  481.     auto a=lst.begin();  
  482.     for(size_t i=0;i  
  483.     {  
  484.         temp[i]=min(*a);  
  485.         a++;  
  486.     }  
  487.     return temp;  
  488. }  
  489. c pmax(initializer_listlst)  
  490. {  
  491.     c temp=rep(0,lst.size());  
  492.     auto a=lst.begin();  
  493.     for(size_t i=0;i  
  494.     {  
  495.         temp[i]=max(*a);  
  496.         a++;  
  497.     }  
  498.     return temp;  
  499. }  
  500. template  
  501. pair range(const any&v)  
  502. {  
  503.     cout<<”MIN: ”<<min(v)<<“ MAX: ”<<max(v)<<endl;  
  504.     return make_pair(min(v),max(v));  
  505. }  
  506. double median(const c&v)  
  507. {  
  508.     return sort(v)[length(v)/2];  
  509. }  
  510. double var(const c&a,const c&b)  
  511. {  
  512.     if(length(a)!=length(b))  
  513.         throw runtime_error(“dim error!\n”);  
  514.     c t1(a-mean(a)),t2(b-mean(b));  
  515.     return inter_product(t1,t2)/(length(a)-1);  
  516. }  
  517. double cov(const c&a,const c&b)  
  518. {  
  519.     return var(a,b);  
  520. }  
  521. double cor(const c&a,const c&b)  
  522. {  
  523.     return var(a,b)/STD(a)/STD(b);  
  524. }  
  525. matrix var(const matrix&a,const matrix&b)  
  526. {  
  527.     if(a.re_dim().first!=b.re_dim().first)  
  528.         throw runtime_error(“dim error!\n”);  
  529.     matrix temp(a.re_dim().second,b.re_dim().second,0);  
  530.     for(size_t i=0;i  
  531.         for(size_t j=0;j  
  532.         temp[i][j]=var(a.col(i),b.col(j));  
  533.     return temp;  
  534. }  
  535. matrix cov(const matrix&a,const matrix&b)  
  536. {  
  537.     return var(a,b);  
  538. }  
  539. matrix cor(const matrix&a,const matrix&b)  
  540. {  
  541.     matrix temp=var(a,b);  
  542.     for(size_t i=0;i  
  543.         for(size_t j=0;j  
  544.         temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));  
  545.     return temp;  
  546. }  
  547. matrix cov(const matrix&a)  
  548. {  
  549.     return var(a,a);  
  550. }  
  551. matrix cor(const matrix&a)  
  552. {  
  553.     return cor(a,a);  
  554. }  
  555. //由于关联容器的删除方法利用误差排序的方法  
  556. //利用顺序容器迭代器的删除方法有困难  
  557. c rank(const c&v)  
  558. {  
  559.     vectortemp00(v.const_re_c_val());  
  560.     for(size_t i=0;i  
  561.         temp00[i]+=(0.00000001*i);  
  562.     vectortemp0(temp00);  
  563.     multisettemp2(temp0.begin(),temp0.end());  
  564.     vectortemp1(v.re_dim());  
  565.     size_t i=0,j=0;  
  566.     double pre_val=min(c(vector(temp2.begin(),temp2.end())));  
  567.     vector::iterator s0;  
  568.     while(i  
  569.     {  
  570.     double val=min(c(vector(temp2.begin(),temp2.end())));  
  571.     s0=find(temp0.begin(),temp0.end(),val);  
  572.     if(sqrt((pre_val-val)*(pre_val-val))>0.000001)  
  573.         j++;  
  574.     temp1[s0-temp0.begin()]=j;  
  575.     temp2.erase(val);  
  576.     i++;  
  577.     pre_val=val;  
  578.     }  
  579.     return  c(temp1);  
  580. }  
  581. c rev(const c&v)  
  582. {  
  583.     vector temp;  
  584.     auto a=v.const_re_c_val().rbegin();  
  585.     while(a!=v.const_re_c_val().rend())  
  586.     {  
  587.         temp.push_back(*a);  
  588.         a++;  
  589.     }  
  590.     return c(temp);  
  591. }  
  592. //向量先排序再输出排序前的下标值  
  593. c order(const c&v)  
  594. {  
  595.     unordered_multimapm0;  
  596.     for(size_t i=0;i  
  597.         m0.insert(make_pair(v[i],i));  
  598.     multimapm1(m0.begin(),m0.end());  
  599.     c temp(rep(0,length(v)));  
  600.     auto a=m1.begin();  
  601.     size_t i=0;  
  602.     while(a!=m1.end())  
  603.     {  
  604.         temp[i]=a->second;  
  605.         i++;  
  606.         a++;  
  607.     }  
  608.     return temp;  
  609. }  
  610. c cumsum(const c&v)  
  611. {  
  612.     c temp(v);  
  613.     for(size_t i=1;i  
  614.         temp[i]+=temp[i-1];  
  615.     return temp;  
  616. }  
  617. c cumprod(const c&v)  
  618. {  
  619.     c temp(v);  
  620.     for(size_t i=1;i  
  621.         temp[i]*=temp[i-1];  
  622.     return temp;  
  623. }  
  624. matrix cumsum(const matrix&m)  
  625. {  
  626.     c temp(m);  
  627.     for(size_t i=1;i  
  628.         temp[i]+=temp[i-1];  
  629.     return matrix(temp,m.re_dim().first,m.re_dim().second);  
  630. }  
  631. matrix cumprod(const matrix&m)  
  632. {  
  633.     c temp(m);  
  634.     for(size_t i=1;i  
  635.         temp[i]*=temp[i-1];  
  636.     return matrix(temp,m.re_dim().first,m.re_dim().second);  
  637. }  
  638.    
  639. //矩阵运算函数  
  640. //初等行变换化上三角求解行列式  
  641. double det(const matrix&m)  
  642. {  
  643. auto f=[](const matrix&mm)->pair  
  644. {  
  645.     if(mm.re_dim().first!=mm.re_dim().second)  
  646.        throw runtime_error(“dim error!\n”);  
  647.     size_t MIN=mm.re_dim().first;  
  648.   
  649.     double mult=1;  
  650.     size_t col=0;  
  651.     size_t bar=0;  
  652.     size_t b_c=0;  
  653.     matrix t(mm);  
  654.   
  655.     while(b_c  
  656.     {  
  657.     for(size_t i=bar+1;i  
  658.     {  
  659.         if(t[i][col]==0)  
  660.             i++;  
  661.         else  
  662.         {  
  663.            t.Swap(i,bar);  
  664.            mult*=(-1);  
  665.            break;  
  666.         }  
  667.     }  
  668.     if(bar!=MIN+1)  
  669.     {  
  670.     if(t[bar][col]!=0)  
  671.     {  
  672.        for(size_t i=bar+1;i  
  673.      {  
  674.         t[bar][i]=t[bar][i]/t[bar][col];  
  675.         if(i==bar+1)  
  676.         mult*=t[bar][col];  
  677.      }  
  678.         if(col!=MIN-1)  
  679.         t[bar][col]=1;  
  680.     }  
  681.   
  682.      for(size_t i=bar+1;i  
  683.      {  
  684.          t[i]=t[i]-t[i][col]*t[bar];  
  685.      }  
  686.     }  
  687.       col++;  
  688.       bar++;  
  689.       b_c=bar;  
  690.       if(col>bar);  
  691.       b_c=col;  
  692.     }  
  693. return make_pair(t,mult);  
  694. };  
  695.     pair tt=f(m);  
  696.     return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));  
  697. }  
  698. //Cramer法则解线性方程组:要求自变量矩阵为方阵且相应行列式不为0  
  699. c solve(const matrix&m,const c&v)  
  700. {  
  701.     c temp(rep(0,v.re_dim()));  
  702.     double low=det(m);  
  703.     if(m.re_dim().first!=m.re_dim().second)  
  704.       throw runtime_error(“dim error!\n”);  
  705.     if(low==0)  
  706.       throw runtime_error(“singular!\n”);  
  707.     for(size_t i=0;i  
  708.     {  
  709.         matrix temp0(m);  
  710.         temp0.col_change(i,v);  
  711.         temp[i]=det(temp0)/low;  
  712.     }  
  713.     return temp;  
  714. }  
  715. double tr(const matrix&m)  
  716. {  
  717.     if(m.re_dim().first!=m.re_dim().second)  
  718.         throw runtime_error(“dim error!\n”);  
  719.     return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));  
  720. }  
  721. matrix T(const matrix&m)  
  722. {  
  723.     matrix out(m.re_dim().second,m.re_dim().first,0);  
  724.     for(size_t i=0;i  
  725.         for(size_t j=0;j  
  726.         out[i][j]=m[j][i];  
  727.     return out;  
  728. }  
  729. matrix cbind(initializer_listlst)  
  730. {  
  731.     return T(rbind(lst));  
  732. }  
  733. matrix operator%(const matrix&a,const matrix&b)  
  734. {  
  735.     matrix out(a.re_dim().first,a.re_dim().second,0);  
  736.     for(size_t i=0;i  
  737.         for(size_t j=0;j  
  738.         out[i][j]=inter_product(a[i],b.col(j));  
  739.     return out;  
  740. }  
  741. double F_norm(const matrix&m)  
  742. {  
  743.     return sqrt(tr(T(m)%m));  
  744. }  
  745.    
  746. //a,b为区间 d为步长 e为精度 默认为求整数精度:以1为精度特征值 如所求特征值太少可通过减少d改善  
  747. //只能求实值  
  748. //可以用F_norm作为a、b的界(-F_norm(m),F_norm(m))可确保值在范围中  
  749. //可用于求解但效率低  
  750. c eigen(const matrix&m,const double&a,const double& b,const double &d=0.0001,const double &e=1)  
  751. {  
  752.     if(m.re_dim().first!=m.re_dim().second)  
  753.         throw runtime_error(“dim error!\n”);  
  754.     vectortemp;  
  755.     double pre_beg=100000000;  
  756.     for(double beg=a;beg  
  757.     {  
  758.         matrix mm0(diag(rep(beg,m.re_dim().first))-m);  
  759.         double s0=det(mm0);  
  760.         if(sqrt(s0*s0)  
  761.         {  
  762.             if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)  
  763.             {  
  764.              temp.push_back(beg);  
  765.              pre_beg=beg;  
  766.              cout<<”Done!\n”;  
  767.              cout<<beg<<endl;//eigenvalue  
  768.             }  
  769.         }  
  770.     }  
  771.     return c(temp);  
  772. }  
  773. //重载求eign的函数 d:步长 可保证绝对求出所有特征值(利用循环增大精度)  
  774. c eigen(const matrix&m,const double&d=0.1)  
  775. {  
  776.     double dd=d;  
  777.     double ra=F_norm(m);  
  778.     double trr=tr(m);  
  779.     c temp;  
  780.     double r=0;  
  781.     do  
  782.     {  
  783.         dd*=0.1;  
  784.         temp=eigen(m,-ra,ra,dd);  
  785.         if(temp.re_dim()==m.re_dim().first)  
  786.             break;  
  787.         r=sum(temp)-trr;  
  788.     }  
  789.     while(sqrt((r*r)>1));  
  790.         //与迹的差距(精度)设定为1 绝对值小于此值的特征值被忽略  
  791.         //此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握  
  792.         //只适用于所有特征值为实值的情况  
  793.         //如有复特征值,只输出实数的并循环不能结束  
  794.     return temp;  
  795. }  
  796. matrix M(const matrix&m,size_t k,size_t l)  
  797. {  
  798.     matrix temp(m);  
  799.     set s1;  
  800.     for(size_t i=0;i  
  801.         s1.insert(i);  
  802.     set s2;  
  803.     for(size_t i=0;i  
  804.         s2.insert(i);  
  805.     s1.erase(k);  
  806.     s2.erase(l);  
  807.     vectortemp00;  
  808.     for(auto j:s2)  
  809.         for(auto i:s1)  
  810.         temp00.push_back(temp[i][j]);  
  811.     return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);  
  812. }  
  813. matrix A_accompany(const matrix&m)  
  814. {  
  815.     matrix out(m.re_dim().first,m.re_dim().second,0);  
  816.     for(size_t i=0;i  
  817.         for(size_t j=0;j  
  818.         out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);  
  819.         //代数余子式  
  820.     return T(out);  
  821. }  
  822. //求方阵的逆矩阵(利用伴随矩阵方法)  
  823. matrix solve(const matrix&m)  
  824. {  
  825.     double de=det(m);  
  826.     if(sqrt(de*de)<0.000001)  
  827.      throw runtime_error(“singular!\n”);  
  828.     return A_accompany(m)/de;  
  829. }  
  830. //利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆  
  831. c solve_eq(const matrix&m,const c&v)  
  832. {  
  833.     c out;  
  834.     double de=det(m);  
  835.     if(sqrt(de*de)<0.000001)  
  836.       throw runtime_error(“singular!\n”);  
  837.     return solve(m)%v;  
  838. }  
  839.    
  840.    
  841. 测试程序:  
  842. int main()  
  843. {  
  844. //测试程序:  
  845.     c c0({35.1,12});  
  846.     vectorv0;  
  847.     for(size_t i=0;i<10;i++)  
  848.         v0.push_back(i);  
  849.     c c1(v0);  
  850.     c0.out();  
  851.     c1.out();  
  852.     cout<<c0[1]<<endl;  
  853.     c0[1]=10;  
  854.     cout<<c0[1]<<endl;  
  855.     (c0+c1).out();  
  856.     (-c0).out();  
  857.      (c0-c1).out();  
  858.      rep(1,10).out();  
  859.      (c0+rep(1,10)).out();  
  860.     (c0-rep(1,10)).out();  
  861.     (rep(1,10)+1).out();  
  862.     (rep(1,10)-1).out();  
  863.     (10/rep(1,10)).out();  
  864.     (rep(1,10)/10).out();  
  865.     cout<<((10/rep(1,10))==(rep(1,10)/10))<<endl;  
  866.     matrix m0=matrix(3,4,10);  
  867.     m0[1]=rep(4,4);  
  868.     m0.out();  
  869.     m0.col(2).out();  
  870.     m0.col_change(0,c({1,2,3}));  
  871.     m0.out();  
  872.     matrix m1=(matrix(c({1,5,7,8,9,0,3,4,5}),3,3)+matrix(c({3,5,7,9}),2,2));  
  873.     m1.out();  
  874.     (-m1).out();  
  875.     (m1-matrix(rep(1,10),5,2)).out();  
  876.     (m1+1).out();  
  877.     (m1*0.6).out();  
  878.     (m1-1).out();  
  879.     (1-m1).out();  
  880.   
  881.     (0.6*m1).out();  
  882.     (1/m1).out();  
  883.     (m1/10).out();  
  884.     m1.out();  
  885.     m1.Swap(1,2);  
  886.     m1.out();  
  887.     m1.col_Swap(1,2);  
  888.     m1.out();  
  889.     matrix    m2=m1;  
  890.     m1.Swap(0,1);  
  891.     cout<<(m1==m2)<<endl;  
  892.     cout<<m1.re_dim().first<<” ”<<m1.re_dim().second<<endl;  
  893.     matrix A(c({1,2,1,4,7,1,5,2,1}),3,3);  
  894.     cout<<det((A+diag(rep(10,A.re_dim().first))))<<endl;  
  895.     solve(matrix(c({30,7.3,4.2,0,0,6,8,7,0}),3,3),c({11,12,13})).out();  
  896.     m1.out();  
  897.     cout<<F_norm(m1)<<endl;  
  898.     eigen(matrix(c({50,1,1,9,7,1,4,1,13}),3,3)).out();  
  899.     solve(matrix(c({50,1,1,9,7,1,4,1,13}),3,3),c({1,1,1})).out();  
  900.     solve(matrix(c({50,1,1,9,7,1,4,1,13}),3,3)).out();  
  901.     solve_eq(matrix(c({50,1,1,9,7,1,4,1,13}),3,3),c({1,1,1})).out();  
  902.     m2.out();  
  903.     M(m2,0,0).out();  
  904.     cout<<sum(m2)<<endl;  
  905.     cout<<prod(m2)<<endl;  
  906.     (c({1,1})%matrix(c({5,6,7,8,9,10}),2,3)).out();  
  907.     outer_product(c({1,1}),c({5,6,7})).out();  
  908.     diag(m2).out();  
  909.     diag(c(diag(m2))).out();  
  910.     cout<<length(m2)<<endl;  
  911.     cout<<length(m2[0])<<endl;  
  912.     sort(c({3,2,1})).out();  
  913.     cout<<mean(m2)<<endl;  
  914.     cout<<mean(m2.col(0))<<endl;  
  915.     cout<<var(matrix(c({1,2,3,4,5,6,7,8,9}),3,3))<<endl;  
  916.     connect({c({1,1,1}),c({-1,1,-1})}).out();  
  917.     cout<<STD(matrix(c({1,2,3,4,5,6,7,8,9}),3,3))<<endl;  
  918.     matrix    m3=rbind({c({1,2}),c({3,4,5})});  
  919.     cbind({c({1,2}),c({3,4,5})}).out();  
  920.     cout<<dim(m3).first<<” ”<<dim(m3).second<<endl;  
  921.   
  922.     m3.out();  
  923.     cout<<dim(m3[0])<<endl;  
  924.     cout<<min(m3)<<endl;  
  925.     cout<<min(m3[0])<<endl;  
  926.     cout<<max(m3)<<endl;  
  927.     cout<<max(m3[0])<<endl;  
  928.     c(m3).out();  
  929.    pmin({m3[0],m3[1]}).out();  
  930.    pmax({m3[0],m3[1]}).out();  
  931.    range(m3);  
  932.    cout<<median(m3)<<endl;  
  933.    var(m3,m3).out();  
  934.    cov(m3).out();  
  935.    cor(m3).out();  
  936.    rank(c({0,3,4,4,5,6,7,6,6,3})).out();  
  937.    rev(c({1,5,7,6,9})).out();  
  938.    order(c({0,3,4,4,5,6,7,6,3})).out();  
  939.    m3.out();  
  940.    cumsum(m3[0]).out();  
  941.    cumsum(m3).out();  
  942.    cumprod(m3[0]).out();  
  943.    cumprod(m3).out();  
  944.     cout<<det(matrix(c({50,7,8,9,10,4,3,2,9,4,8,7,1,1,1,10}),4,4))<<endl;  
  945.   
  946.     return 0;  
  947. }  
#includeusing std::cout;using std::endl;using std::cin;using std::istream;using std::ios_base;
#includeusing std::vector;#includeusing std::initializer_list;#includeusing std::for_each;using std::sort;using std::copy;using std::find;#includeusing std::runtime_error;#includeusing std::srand;using std::rand;using std::system;#includeusing std::time;#includeusing std::sqrt;#includeusing std::pair;using std::make_pair;#includeusing std::setw;#includeusing std::multiset;using std::set;#includeusing std::unordered_map;using std::unordered_multimap;#includeusing 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 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);}//标准差templatedouble 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;}templatepair 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[bar];     }    }      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法则解线性方程组:要求自变量矩阵为方阵且相应行列式不为0c 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.

原创粉丝点击