#include<cmath>#include<vector>#include<stdexcept>#include<limits>#include<iostream>#include<iomanip>// 计算一组 double 的几何平均数// 行为对齐 MATLAB geomean:// - 任一元素为 0 → 返回 0// - 任一元素为负 → 返回 NaN// - 空输入 → 返回 NaNdoublegeomean_vector(conststd::vector<double>&x){if(x.empty()){returnstd::numeric_limits<double>::quiet_NaN();}doublesum_log=0.0;for(doublev:x){if(v==0.0){return0.0;}if(v<0.0){returnstd::numeric_limits<double>::quiet_NaN();}sum_log+=std::log(v);}returnstd::exp(sum_log/static_cast<double>(x.size()));}// 矩阵按列求几何平均数（行优先存储）// data: rows × cols 的矩阵，按行优先展平// 返回长度为 cols 的 vector，对应每一列的几何平均数std::vector<double>geomean_columns(conststd::vector<double>&data,std::size_t rows,std::size_t cols){if(data.size()!=rows*cols||rows==0||cols==0){throwstd::invalid_argument("维度不匹配或为空");}std::vector<double>result(cols);for(std::size_t j=0;j<cols;++j){doublesum_log=0.0;boolhas_zero=false;boolhas_negative=false;for(std::size_t i=0;i<rows;++i){doublev=data[i*cols+j];if(v==0.0){has_zero=true;break;}if(v<0.0){has_negative=true;break;}sum_log+=std::log(v);}if(has_zero){result[j]=0.0;}elseif(has_negative){result[j]=std::numeric_limits<double>::quiet_NaN();}else{result[j]=std::exp(sum_log/static_cast<double>(rows));}}returnresult;}// 所有元素的几何平均数，对应 geomean(X, 'all')doublegeomean_all(conststd::vector<double>&data){returngeomean_vector(data);}// data: rows × cols 矩阵，列优先存储（data[r + c * rows]）// dim = 1: 沿列方向，返回长度 cols// dim = 2: 沿行方向，返回长度 rowsstd::vector<double>geomean_dim(conststd::vector<double>&data,std::size_t rows,std::size_t cols,intdim){if(data.size()!=rows*cols){throwstd::invalid_argument("数据长度与维度不符");}if(dim==1){std::vector<double>result(cols);for(std::size_t c=0;c<cols;++c){doublesum_log=0.0;boolzero=false,neg=false;for(std::size_t r=0;r<rows;++r){doublev=data[r+c*rows];if(v==0.0){zero=true;break;}if(v<0.0){neg=true;break;}sum_log+=std::log(v);}result[c]=zero?0.0:neg?std::numeric_limits<double>::quiet_NaN():std::exp(sum_log/static_cast<double>(rows));}returnresult;}if(dim==2){std::vector<double>result(rows);for(std::size_t r=0;r<rows;++r){doublesum_log=0.0;boolzero=false,neg=false;for(std::size_t c=0;c<cols;++c){doublev=data[r+c*rows];if(v==0.0){zero=true;break;}if(v<0.0){neg=true;break;}sum_log+=std::log(v);}result[r]=zero?0.0:neg?std::numeric_limits<double>::quiet_NaN():std::exp(sum_log/static_cast<double>(cols));}returnresult;}throwstd::invalid_argument("dim 仅支持 1 或 2");}intmain(){std::cout<<std::setprecision(10);// 向量: [1, 2, 3, 4, 5, 6]std::vector<double>v={1,2,3,4,5,6};std::cout<<"geomean_vector = "<<geomean_vector(v)<<"\n";// 含 0std::vector<double>v0={1,0,3};std::cout<<"含0 = "<<geomean_vector(v0)<<"\n";// 含负值std::vector<double>vn={1,-2,3};doublegn=geomean_vector(vn);std::cout<<"含负 = ";if(std::isnan(gn)){std::cout<<"NaN";}else{std::cout<<gn;}std::cout<<"\n";// 矩阵按列 (2×3, 行优先)std::vector<double>m={1,2,3,4,5,6};autocols=geomean_columns(m,2,3);std::cout<<"矩阵按列 = ";for(std::size_t j=0;j<cols.size();++j){if(j>0){std::cout<<" ";}std::cout<<cols[j];}std::cout<<"\n";return0;}