C++ vd类代码示例

bool gaussian_elimination(vvd &matrix, vd &output) {
  for (int i = 0; i < matrix.size(); i++) {
    int pivot = find_pivot(matrix, i);
    if (fabs(matrix[pivot][i]) < EPS) {
      return false;
    }
    swap_row(matrix, pivot, i);
    double factor = matrix[i][i];
    for (int j = 0; j < matrix[i].size(); j++) {
      matrix[i][j] /= factor;
    }
    for (int j = 0; j < matrix.size(); j++) {
      if (j == i) {
        continue;
      }
      factor = matrix[j][i];
      for (int k = 0; k < matrix[j].size(); k++) {
        matrix[j][k] -= factor*matrix[i][k];
      }
    }
  }
  output.clear();
  for (int i = 0; i < matrix.size(); i++) {
    output.push_back(matrix[i][matrix[i].size() - 1]);
  }
  return true;
}

	double SumSquaredDifferences( vd v, double avg)
	{
		vd::iterator it;
		double sumSq = 0.0;

		for(it = v.begin(); it < v.end(); ++it)
			sumSq += pow(*it - avg,2.0);

		return( sumSq );
	}

	double Sum( vd v)
	{
		vd::iterator it;
		double sum = 0.0;

		for(it = v.begin(); it < v.end(); ++it)
			sum += *it;

		return(sum);
	}

 vd conv(vd a, vd b) {
     int s = max(sz(a),sz(b)), L = get(s), n = 1<<L;
     if (s <= 0) return {};
     
     a.resize(n); a = fwht(a);
     b.resize(n); b = fwht(b);
     
     F0R(i,n) a[i] = a[i]*b[i];
     a = fwht_rev(a);
     return a;
 }

	double average( vd v)
	{
		vd::iterator it;
		double sum = 0.0;

		if ( v.size() == 0) return(0.0);

		for(it = v.begin(); it < v.end(); ++it)
			sum += *it;

		return(sum/v.size());
	}

//a1x1+a2x2+a3=0
//b1x1+b2x2+b3=0
vctr solve(vd &a, vd &b)
{
    assert(a.size()==b.size()&&a.size()==3);
    if(fabs(b[0])>fabs(a[0]))swap(a,b);
    assert(fabs(a[0])>EPS);
    double k=b[0]/a[0];
    rep(i,3)b[i]-=a[i]*k;
    assert(fabs(b[1])>EPS);
    double x2=-b[2]/b[1];
    double x1=-(a[2]+a[1]*x2)/a[0];
    return vctr(x1,x2);
}

void constructPoints(const vector<Edge> &edges, const vd &mesh,
    vector<Point> &points, vvi &edge2point) {
  points.clear();
  for (int x = 2; x <= 8; x += 2) {
    for (int y = 2; y <= 8; y += 2) {
      for (int z = 2; z <= 8; z += 2) {
        D3 p(x, y, z);
        if (isOneOfEndPoint(p, edges)) {
          points.push_back(Point(p));
          for (int i = 0; i < edges.size(); i++) {
            if (edges[i].src == p) {
              addPoint(points, edge2point, i);
            } else if (edges[i].dst == p) {
              addPoint(points, edge2point, i);
            }
          }
        }
      }
    }
  }
  for (int i = 0; i < edges.size(); i++) {
    D3 p = edges[i].src;
    D3 v = edges[i].dst - edges[i].src;
    for (int j = 0; j < mesh.size(); j++) {
      points.push_back(Point(p + mesh[j]*v));
      addPoint(points, edge2point, i);
    }
  }
}

	double stddev( vd v, double avg)
	{
		vd::iterator it;
		double sumSq = SumSquaredDifferences(v, avg);

		return( stddev(sumSq,v.size()) );
	}

void getDualProblem(const vvd& a, vd& b, vd& c, vvd& dual_a, vd& dual_b, vd& dual_c)
{
    for (auto x : b)
    {
        dual_c.push_back(-x);
    }
    for (auto x : c)
    {
        dual_b.push_back(-x);
    }
    dual_a.assign(a[0].size(), vd());
    for (size_t i = 0; i < a.size(); ++i)
    {
        for (size_t j = 0; j < a[i].size(); ++j)
        {
            dual_a[j].push_back(-a[i][j]);
        }
    }
}

int main() {
    ios::sync_with_stdio(0);
    while(cin >> n >> p && n) {
        po.assign(n + 1, 1);
        for(int i = 1; i <= n; i++) po[i] = po[i - 1] * p;
        memo.assign(n + 1, vd(n + 1, -1));
        cout << fixed << setprecision(10) << solve(0, n) << endl;
    }


    return 0;
}

void dfs(vvid &adj, int fr, int to, vd path, vb &vis){
  vis[fr] = 1;
  if (fr == to){
    for (double e : path)
      res = max(res, e);
    return;
  }
  for (int i=0;i<adj[fr].size();i++){
    if (!vis[adj[fr][i].first]){
      path.push_back(adj[fr][i].second);
      dfs(adj, adj[fr][i].first, to, path, vis);
    }
  }
  return;
}

void printTask(const vvd& a, const vd& b, const vd& c, const std::string& filename = "task.txt")
{
    size_t n = c.size();
    size_t m = a.size();
    std::ofstream out(filename);
    for (size_t i = 0; i < n; ++i)
    {
        auto& t = c[i];
        if (std::abs(t) < EPS)
        {
            continue;
        }
        if (t > 0)
        {
            out << "+\t" << t;
        }
        if (t < 0)
        {
            out << "-\t" << -t;
        }
        out << " * x" << i + 1 << "\t";
    }
    out << "-----> MAX\n";
    for (size_t i = 0; i < m; ++i)
    {
        out << "x" << n + i + 1 << "\t=\t" << b[i] << "\t";
        for (size_t j = 0; j < n; ++j)
        {
            auto t = -a[i][j];
            if (std::abs(t) < EPS)
            {
                continue;
            }
            if (t > 0)
            {
                out << "+\t" << t;
            }
            if (t < 0)
            {
                out << "-\t" << -t;
            }
            out << " * x" << j + 1 << "\t";
        }
        out << endl;
    }
}

double pcet::operator()(const vd& v) {
  double xn1 = v[1],pn1 = v[2],xn2 = v[3],pn2 = v[4],qr = v[5],pr = v[6];
  double n1=e.n(xn1,pn1),n2=e.n(xn2,pn2);
  double h_bath=0.0;
  for (size_t i = 7; i < v.size(); ++(++i)) {
    size_t ib=(i-7)/2;
    double mw2=bp.m*pow(bw[ib],2);
    h_bath+= v[i+1]*v[i+1]*0.5/bp.m
      + mw2*0.5*pow((v[i]+n1*bc[ib]/mw2),2);
  }
  double h11 = sp.eps_r + sp.w_r*sp.w_r*sp.m*pow((qr-sp.qr_initial),2)*0.5;
  double h22 = sp.eps_k + sp.w_k*sp.w_k*sp.m*pow((qr-sp.qr_initial),2)*0.5;
  double h12 = sp.delta;
  double h_avg= (h11+h22)*0.5;
  double h_delta= h11-h22;
  return pr*pr*0.5/sp.m+
    h_avg+
    h_delta*(n1-n2)*0.5+
    sp.delta*(xn1*xn2+pn1*pn2) + h_bath;
    //2*delta*sqrt((n1+e.n_shift)*(n2+e.n_shift))*(xn1*xn2+pn1*pn2); ///> TODO
}

void print(const vd &a){rep(i,a.size())printf("%.2lf%c",a[i],i==a.size()-1?'\n':' ');}

vi bitstream(vd in)
{
	vi out(in.size());
	for (int n=0; n<in.size(); n++) {out[n]=(copysign(1,in[n])+1)/2;}
	return out;
}