C++ Vec::length方法代码示例

// Prime-power factorization
void pp_factorize(vector<long>& factors, long N)
{
  Vec< Pair<long, long> > pf;
  factorize(pf,N); // prime factors, N = \prod_i pf[i].first^{pf[i].second}
  factors.resize(pf.length());
  for (long i=0; i<pf.length(); i++)
    factors[i] = power_long(pf[i].a, pf[i].b); // p_i^e_i
}

void assert_vec_p(const Vec<T> &ref, const Vec<T> &act, int line)
{
  static const double tol = 1e-3;
  ASSERT_EQ(ref.length(), ref.length()) << line;
  for (int n = 0; n < ref.length(); ++n) {
    ASSERT_NEAR(ref(n), act(n), tol) << line;
  }
}

double computeAverageDouble(Vec<double> &result){

    RR average = 0;

    for (int i = 1; i <= result.length() ; i++){
        average += result(i);
    }

    return to_double(average/result.length());
}

void setHyperColumn(const Vec<T>& v, const CubeSlice<T>& s, long pos)
{
   long m = s.getProd(1);
   long n = s.getDim(0);

   assert(pos >= 0 && pos < m);
   if (v.length() < n) n = v.length();

   for (long i = 0; i < n; i++)
      s[pos + i*m] = v[i];
}

void
DEMParticleCreator::addExtraTranslations(const Vec& shift, 
                                         REAL boundaryMargin,
                                         Vec inPlaneDiag,
                                         REAL widthX, REAL widthY, REAL widthZ,
                                         const Face& face,
                                         const IntersectionStatus status, 
                                         std::vector<Vec>& translations) 
{
  if (status.second.first == Face::Location::VERTEX) {
    int vertIndex = status.second.second;
    //std::cout << "is vertex " << vertIndex << "\n";
    if (vertIndex == 0) {
      for (int ii = 1; ii < 4; ++ii) {
        Vec inPlane = face.vertex[ii] - face.vertex[vertIndex];
        auto length = inPlane.length();
        inPlane.normalizeInPlace();
        if (ii == 2) {
          inPlane *= length;
          inPlane += (inPlaneDiag*boundaryMargin);
        } else {
          inPlane *= (length + boundaryMargin);
        }
        Vec outOfPlane = inPlane + shift;
        translations.push_back(inPlane);
        translations.push_back(outOfPlane);
      }
    } else if (vertIndex == 1 || vertIndex == 3) {
      Vec normal = Vec(1,1,1) - inPlaneDiag;
      Vec vec1 = normal * Vec(widthX + boundaryMargin, 
        widthY + boundaryMargin, widthZ + boundaryMargin);
      Vec vec2 = face.vertex[2] - face.vertex[vertIndex];
      auto length = vec2.length();
      vec2.normalizeInPlace();
      vec2 *= (length + boundaryMargin);
      //std::cout << "vec1 = " << vec1 << " vec2 = " << vec2 << "\n";
      translations.push_back(vec1 + vec2);
    }
  } 
  else if (status.second.first == Face::Location::EDGE) {
    int edgeIndex = status.second.second;
    //std::cout << "is edge " << edgeIndex << "\n";
    if (edgeIndex == 0 || edgeIndex == 3) {
      int oppIndex = (edgeIndex+3) % 4;
      Vec inPlane = face.vertex[oppIndex] - face.vertex[edgeIndex];
      auto length = inPlane.length() + boundaryMargin;
      inPlane.normalizeInPlace();
      inPlane *= length;
      Vec outOfPlane = inPlane + shift;
      translations.push_back(inPlane);
      translations.push_back(outOfPlane);
    }
  }
}

QString TetraGPSEncoder::generate(TGPSReal lat, TGPSReal lon, TGPSReal precMeters){
	const Vec coord(cos(lat)*cos(lon), cos(lat)*sin(lon), sin(lat));
	qDebug()<<"COORD:"<<coord<<coord.length()<<coord.latlon();
	QVector<Tri> faces;
	const TGPSReal ma=(M_PI*2.0)/3.0;
	const Vec north(    +0.0,                                    +1.0,      +0.0);
	const Vec prime(    +0.0,                                    +sin(-ma), +cos(-ma));
	const Vec meridian2(prime.x()*cos(+ma) - prime.z()*sin(+ma), prime.y(), prime.z()*cos(+ma) + prime.x()*sin(+ma));
	const Vec meridian3(prime.x()*cos(-ma) - prime.z()*sin(-ma), prime.y(), prime.z()*cos(-ma) + prime.x()*sin(-ma));
	qDebug()<<"north=    "<<north<<north.length();
	qDebug()<<"prime=    "<<prime<<prime.length();
	qDebug()<<"meridian2="<<meridian2<<meridian2.length();
	qDebug()<<"meridian3="<<meridian3<<meridian3.length();
	faces.push_back(Tri(north,prime,meridian2, 1));
	faces.push_back(Tri(north,prime,meridian3, 2));
	faces.push_back(Tri(north,meridian2,meridian3, 3));
	faces.push_back(Tri(prime,meridian2,meridian3, 4));
	TGPSReal score=1000.0;
	QStringList ret;
	int sub=0;
	Tri best,second;
	while((score*planetRadiusMeters)>precMeters && sub++<9){
		//qDebug()<<"----ROUND "<<sub<< score;
		TGPSReal eps=1000.0;
		for(auto tri:faces) {
			TGPSReal delta=(tri.center-coord).length();
			//qDebug()<<"DELTA: "<<delta<<tri.center<<tri.center.length();
			if(delta<eps){
				if(delta<score){
					score=delta;
				}
				eps=delta;
				best=tri;
				qDebug()<<"NEW BEST: "<<best<<"!";
			}
		}
		QString name=best.name();
		//if("0"==name){			break;		}
		ret<<name;
		faces=best.subdivide();

		if(best!=second){
			second=best;
		}
		else{
			break;
		}
	}
	TGPSReal d=(coord-best.center).length()/best.reach();
	qDebug()<<"BEST: "<<best<<" with "<< d<<" normalized distance from center and " << (score*planetRadiusMeters)<<" meters after "<<sub<<" iterations";
	return ret.join(".");
}

// Get multiple layers of a Benes permutation network. Returns in out[i][j]
// the shift amount to move item j in the i'th layer. Also isID[i]=true if
// the i'th layer is the identity (i.e., contains only 0 shift amounts).
void ColPerm::getBenesShiftAmounts(Vec<Permut>& out, Vec<bool>& isID,
				   const Vec<long>& benesLvls) const
{
  // Go over the columns one by one. For each column extract the columns
  // permutation, prepare a Benes network for it, and then for every layer
  // compute the shift amounts for this columns.

  long n = getDim(dim);     // the permutations are over [0,n-1]

  // Allocate space
  out.SetLength(benesLvls.length());
  isID.SetLength(benesLvls.length());
  for (long k=0; k<benesLvls.length(); k++) {
    out[k].SetLength(getSize());
    isID[k] = true;
  }

  Vec<long> col;
  col.SetLength(n);

  for (long slice_index = 0; slice_index < numSlices(dim); slice_index++) {
    ConstCubeSlice<long> slice(*this, slice_index, dim);
    for (long col_index = 0; col_index < slice.numCols(); col_index++) {
      getHyperColumn(col, slice, col_index);

      GeneralBenesNetwork net(col); // build a Benes network for this column

      // Sanity checks: width of network == n,
      //                and sum of benesLvls entries == # of levels
      assert(net.getSize()==n);
      {long sum=0;
       for (long k=0; k<benesLvls.length(); k++) sum+=benesLvls[k];
       assert(net.getNumLevels()==sum);
      }

      // Compute the layers of the collapased network for this column
      for (long lvl=0,k=0; k<benesLvls.length(); lvl += benesLvls[k], k++) {

	// Returns in col the shift amounts for this layer in the network,
	// restricted to this column. Also returns true if the returned
	// permutation is the idendity, false otherwise.
	bool id = collapseBenesLevels(col, net, lvl, benesLvls[k]);
	isID[k] = isID[k] && id;

        CubeSlice<long> oslice(out[k], getSig());
        CubeSlice<long> osubslice(oslice, slice_index, dim);
        setHyperColumn(col, osubslice, col_index);
      }  // next collapsed layer
    }  // next column
  } // next slice
}

void vecRed(Vec<ZZ>& out, const Vec<ZZ>& in, long q, bool abs)
{
  out.SetLength(in.length());  // allocate space if needed

  for (long i=0; i<in.length(); i++) {
    long c = in[i]%q;
    if (abs)       { if (c<0) c+=q; }
    else if (q==2) { if (in[i]<0) c = -c; }
    else { 
      if (c >= q/2)        c -= q;
      else if (c < -(q/2)) c += q;
    }
    out[i] = c;
  }
}

// powVec[d] = m_d = p_d^{e_d}
// computes multiEvalPoints[d] as a vector of length phi(m_d)
//   containing base^{(m/m_d) j} for j in Z_{m_d}^*
void computeMultiEvalPoints(Vec< Vec<zz_p> >& multiEvalPoints,
                            const zz_p& base,
                            long m,
                            const Vec<long>& powVec,
                            const Vec<long>& phiVec)
{
   long k = powVec.length();

   multiEvalPoints.SetLength(k);

   for (long d = 0; d < k; d++) {
      long m_d = powVec[d];
      long phi_d = phiVec[d];
      long count = 0;

      zz_p pow = conv<zz_p>(1);
      zz_p mult = power(base, m/m_d);

      multiEvalPoints[d].SetLength(phi_d);

      for (long j = 0; j < m_d; j++) {
         if (GCD(j, m_d) == 1) {
            multiEvalPoints[d][count] = pow;
            count++;
         }
         pow = pow * mult;
      }
   }
}

// factors[d] = (p_d, e_d)
// computes powVec[d] = p_d^{e_d}
void computePowVec(Vec<long>& powVec, 
                   const Vec< Pair<long, long> >& factors)
{
  long k = factors.length();
  powVec.SetLength(k);
  for (long d = 0; d < k; d++)
    powVec[d] = computePow(factors[d]);
}

// returns \prod_d vec[d]
long computeProd(const Vec<long>& vec)
{
  long prod = 1;
  long k = vec.length();
  for (long d = 0; d < k; d++)
    prod = prod * vec[d];
  return prod;
}

void computePowVec(Vec<long>& powVec, 
                   const Vec< Pair<long, long> >& factors)
{
  long k = factors.length();
  powVec.SetLength(k);
  for (long i = 0; i < k; i++)
    powVec[i] = computePow(factors[i]);
}

// computeProd(vec) returns the product of the entries of vec
long computeProd(const Vec<long>& vec)
{
  long prod = 1;
  long k = vec.length();
  for (long i = 0; i < k; i++)
    prod = prod * vec[i];
  return prod;
}

// powVec[d] = p_d^{e_d}, m = \prod_d p_d^{e_d}
// computes divVec[d] = m/p_d^{e_d}
void computeDivVec(Vec<long>& divVec, long m,
                   const Vec<long>& powVec)
{
  long k = powVec.length();
  divVec.SetLength(k);

  for (long d = 0; d < k; d++)
    divVec[d] = m/powVec[d];
}

void computeDivVec(Vec<long>& divVec, long m,
                   const Vec<long>& powVec)
{
  long k = powVec.length();
  divVec.SetLength(k);

  for (long i = 0; i < k; i++)
    divVec[i] = m/powVec[i];
}