C++ math::power2方法代码示例

  void run(alps::ObservableSet& obs) {
    ++mcs;

    // initialize cluster information
    std::fill(fragments.begin(), fragments.end(), fragment_t());

    // cluster generation
    for (double t = r_time(); t < 1; t += r_time()) {
      int s0 = lattice.num_sites() * uniform_01();
      int s1 = lattice.num_sites() * uniform_01();
      if (spins[s0] == spins[s1]) unify(fragments, s0, s1);
    }

    // assign cluster id & accumulate cluster properties
    int nc = 0;
    double mag2 = 0, mag4 = 0;
    BOOST_FOREACH(fragment_t& f, fragments) {
      if (f.is_root()) {
        f.set_id(nc++);
        double w = f.weight();
        mag2 += power2(w);
        mag4 += power4(w);
      }
    }
    BOOST_FOREACH(fragment_t& f, fragments) f.set_id(cluster_id(fragments, f));

    // flip spins
    for (int c = 0; c < nc; ++c) flip[c] = (uniform_01() < 0.5);
    for (int s = 0; s < lattice.num_sites(); ++s)
      if (flip[fragments[s].id()]) spins[s] ^= 1;

    double mu = 0;
    for (int s = 0; s < lattice.num_sites(); ++s) mu += 2 * spins[s] - 1;

    obs["Number of Clusters"] << (double)nc;
    obs["Magnetization (unimproved)"] << mu;
    obs["Magnetization^2 (unimproved)"] << power2(mu);
    obs["Magnetization^4 (unimproved)"] << power4(mu);
    obs["Magnetization^2"] << mag2;
    obs["Magnetization^4"] << (3 * power2(mag2) - 2 * mag4);
  }

//.........这里部分代码省略.........
           itr!= operators.end(); ++itr, itr_new += 2) *itr_new = *itr;
      std::swap(operators, operators_new);
    }

    // initialize cluster information (setup s cluster fragments)
    fragments.resize(nsites);
    std::fill(fragments.begin(), fragments.end(), fragment_t());
    for (unsigned int s = 0; s < nsites; ++s) current[s] = s;

    for (std::vector<local_operator_t>::iterator oi = operators.begin();
         oi != operators.end(); ++oi) {

      // diagonal update
      if (oi->type == identity) {
        unsigned int b = nbonds * random();
        if (spins[left(nbonds, b)] != spins[right(nbonds, b)] &&
            (operators.size() - nop) * random() < lb2) {
          *oi = local_operator_t(b);
          ++nop;
        } else {
          continue;
        }
      } else {
        if (oi->type == diagonal &&
            lb2 * random() < operators.size() - nop + 1) {
          oi->type = identity;
          --nop;
          continue;
        }
      }

      // cluster generation
      unsigned int s0 = left(nbonds, oi->bond);
      unsigned int s1 = right(nbonds, oi->bond);
      oi->lower_loop = unify(fragments, current[s0], current[s1]);
      oi->upper_loop = current[s0] = current[s1] = add(fragments);
      if (oi->type == offdiagonal) {
        spins[s0] ^= 1;
        spins[s1] ^= 1;
      }
    }

    // connect bottom and top cluster fragments
    for (int s = 0; s < nsites; ++s) unify(fragments, s, current[s]);

    //
    // cluster flip
    //

    // assign cluster id & determine if clusters are to be flipped
    int nc = 0;
    BOOST_FOREACH(fragment_t& f, fragments) { if (f.is_root()) f.set_id(nc++); }
    clusters.resize(nc);
    BOOST_FOREACH(fragment_t& f, fragments) { f.set_id(cluster_id(fragments, f)); }
    for (int c = 0; c < nc; ++c) clusters[c] = cluster_t(random() < 0.5);

    // 'flip' operators & do improved measurements
    unsigned int t = 0;
    for (std::vector<local_operator_t>::iterator oi = operators.begin();
         oi != operators.end(); ++oi) {
      if (oi->type == identity) continue;
      int id_l = fragments[oi->lower_loop].id();
      int id_u = fragments[oi->upper_loop].id();
      clusters[id_l].length += 2 * t;
      clusters[id_u].length -= 2 * t;
      if (clusters[id_l].to_flip ^ clusters[id_u].to_flip) oi->flip();
      ++t;
    }

    // flip spins & do improved measurements
    for (unsigned int s = 0; s < nsites; ++s) {
      int id = fragments[s].id();
      clusters[id].size += 1;
      clusters[id].mag += 1 - 2 * spins[s];
      clusters[id].length += nop;
      if (clusters[id].to_flip) spins[s] ^= 1;
    }

    if (mcs < therm) continue;

    //
    // measurements
    //

    // accumurate loop length and magnetization
    double s2 = 0;
    double m2 = 0;
    double l2 = 0;
    for (std::vector<cluster_t>::const_iterator pi = clusters.begin();
         pi != clusters.end(); ++pi) {
      s2 += power2(pi->size);
      m2 += power2(pi->mag);
      l2 += power2(pi->length);
    }

    energy << (0.25 * nbonds - nop / beta) / nsites;
    smag << 0.25 * s2 / nsites;
    usus << 0.25 * beta * m2 / nsites;
    ssus << 0.25 * beta * ((nop ? l2 / nop : 0) + s2) / (nop + 1) / nsites;
  }