int main()
{
    MarchingCubeGridf phi(1,1,1);
    std::ifstream phifile("cube.phi");
    phi.read(phifile);
    Meshf mesh;
    phi.populateMesh(mesh,0);
    std::cout << mesh.vertices.size() << std::endl;
    std::cout << mesh.edges.size() << std::endl;
    std::cout << mesh.faces.size() << std::endl;
    std::ofstream out("output.obj");
    MeshIOUtils::OBJParser<NumericalTraitsXf> parser(&mesh);
    parser.write(out);



}

//**********************************************************************
PHX_EVALUATE_FIELDS(FEInterpolation,cell_data)
{ 

  std::vector<Element_Linear2D>::iterator cell_it = cell_data.begin;

  // Loop over number of cells
  for (std::size_t cell = 0; cell < cell_data.num_cells; ++cell) {
    
    const shards::Array<double,shards::NaturalOrder,QuadPoint,Node>& phi = 
      cell_it->basisFunctions();

    const shards::Array<double,shards::NaturalOrder,QuadPoint,Node,Dim>& 
      grad_phi = cell_it->basisFunctionGradientsRealSpace();

    // Loop over quad points of cell
    for (int qp = 0; qp < num_qp; ++qp) {
      
      val_qp(cell,qp) = 0.0;

      for (int dim = 0; dim < num_dim; ++dim)
	val_grad_qp(cell,qp,dim) = 0.0;

      // Sum nodal contributions to qp
      for (int node = 0; node < num_nodes; ++node) {

	val_qp(cell,qp) += phi(qp,node) * val_node(cell,node);
	
	for (int dim = 0; dim < num_dim; ++dim)
	  val_grad_qp(cell,qp,dim) += 
	    grad_phi(qp,node,dim) * val_node(cell,node);
	
      }
    }
    
    ++cell_it;
 
  }
    
//   std::cout << "FEINterpolation: val_node" << std::endl;
//   val_node.print(std::cout,true);
//   std::cout << "FEINterpolation: val_qp" << std::endl;
//   val_qp.print(std::cout,true);
//   std::cout << "FEINterpolation: val_grad_qp" << std::endl;
//   val_grad_qp.print(std::cout,true);

}

Foam::tmp<Foam::volScalarField>
Foam::diameterModels::IATEsources::randomCoalescence::R() const
{
    tmp<volScalarField> tR
    (
        new volScalarField
        (
            IOobject
            (
                "R",
                iate_.phase().time().timeName(),
                iate_.phase().mesh()
            ),
            iate_.phase().mesh(),
            dimensionedScalar("R", dimless/dimTime, 0)
        )
    );

    volScalarField R = tR();

    scalar Crc = Crc_.value();
    scalar C = C_.value();
    scalar alphaMax = alphaMax_.value();
    volScalarField Ut(this->Ut());
    const volScalarField& alpha = phase();
    const volScalarField& kappai = iate_.kappai();
    scalar cbrtAlphaMax = cbrt(alphaMax);

    forAll(R, celli)
    {
        if (alpha[celli] < alphaMax - SMALL)
        {
            scalar cbrtAlphaMaxMAlpha = cbrtAlphaMax - cbrt(alpha[celli]);

            R[celli] =
                (-12)*phi()*kappai[celli]*alpha[celli]
               *Crc
               *Ut[celli]
               *(1 - exp(-C*cbrt(alpha[celli]*alphaMax)/cbrtAlphaMaxMAlpha))
               /(cbrtAlphaMax*cbrtAlphaMaxMAlpha);
        }
    }

    return tR;
}

FnMap minimizeLength(const FnWord &u, const Basis &basis)
{

  /* Returns an automorphism phi such that phi(u) has minimal
     length. */

  int r = basis.getRank();
  FnMap phi(r);

  if (u == Id) return phi;
  if (u.length() == 1) return phi;

  int old_norm,new_norm;
  FnWord u_min(u),tmp;
  FnMap whAuto(r);

  QList<WhiteheadData> whAutos = whiteheadAutos(basis);
  QListIterator<WhiteheadData> move(whAutos);

  bool reduced_norm = true;

  while (reduced_norm) {

    old_norm = u_min.length();
    new_norm = old_norm;
    move.toFront();

    while (old_norm <= new_norm && move.hasNext()) {
      whAuto = whitehead(move.next(),basis);
      tmp = whAuto(u_min);
      new_norm = tmp.length();
    }
    if (new_norm < old_norm) {
      u_min = tmp;
      phi = whAuto*phi;
    } else
      reduced_norm = false;
  } // end while (reduced_norm)
 
  // now u_min = phi(u)

  return phi;

} 

void Foam::Leonov::correct()
{
    // Velocity gradient tensor
    volTensorField L = fvc::grad(U());

    // Convected derivate term
    volTensorField C = sigma_ & L;

    // Twice the rate of deformation tensor
    volSymmTensorField twoD = twoSymm(L);

    // Two phase transport properties treatment
    volScalarField alpha1f = min(max(alpha(), scalar(0)), scalar(1));

    volScalarField lambda =  alpha1f*lambda1_ + (scalar(1) - alpha1f)*lambda2_;

    volScalarField etaP =  alpha1f*etaP1_ + (scalar(1) - alpha1f)*etaP2_;

     // Stress transport equation
    fvSymmTensorMatrix sigmaEqn
    (
        fvm::ddt(sigma_)
      + fvm::div(phi(), sigma_)
     ==
        twoSymm(C)
      - 1/etaP/2*(symm(sigma_ & sigma_) - Foam::pow((etaP/lambda), 2)*I_)
      + fvm::Sp
        (
            1/etaP/6*
            (
                tr(sigma_)
              - Foam::pow(etaP/lambda,2) * tr(inv(sigma_))
            ),
            sigma_
        )
    );


    sigmaEqn.relax();
    sigmaEqn.solve();

    // Viscoelastic stress
    tau_ = sigma_ - etaP/lambda*I_;
}

void getError(LevelData<double, 1> & a_error,
              double               & a_maxError,
              const double         & a_dx)
{

  BoxLayout layout = a_error.getBoxLayout();

  LevelData<double, 1> phi(layout, s_nghost);
  LevelData<double, 1> lphcalc(layout, 0);
  LevelData<double, 1> lphexac(layout, 0);

  //  cout << "initializing phi to sum_dir(sin 2*pi*xdir)" << endl;
  initialize(phi,lphexac, a_dx);

  //set ghost cells of phi
  phi.exchange();
  //  dumpLevelRDA(&phi);
  Stencil<double> laplace;
  setStencil(laplace, a_dx);

  //apply stencil operator independently on each box
  a_maxError = 0;
  
  for(BLIterator blit(layout); blit != blit.end(); ++blit)
    {
      RectMDArray<double>& phiex =     phi[*blit];
      RectMDArray<double>& lphca = lphcalc[*blit];
      RectMDArray<double>& lphex = lphexac[*blit];
      RectMDArray<double>& error = a_error[*blit];

      //apply is set as an increment so need to set this to zero initially
      lphca.setVal(0.);

      Box bxdst=layout[*blit];
      Stencil<double>::apply(laplace, phiex, lphca, bxdst);
        
      //error = lphicalc -lphiexac
      lphca.copyTo(error);
      //here err holds lphi calc
      double maxbox = forall_max(error, lphex, &errorF, bxdst);

      a_maxError = max(maxbox, a_maxError);
    }
}

///* Computation of Critical Rate in the Jamishidian decomposition, with the newton method to find zero of a function
static double Critical_Rate(ZCMarketData* ZCMarket, double r_initial, double periodicity, double option_maturity, double contract_maturity, double SwaptionFixedRate, double a, double sigma)
{
  double previous,current;
  int nbr_iterations;

  const double precision = 0.000001;

  current = r_initial;
  nbr_iterations = 0;

  do
    {
      nbr_iterations++;
      previous =current;
      current=current-phi(ZCMarket, current, periodicity, option_maturity, contract_maturity, SwaptionFixedRate, a, sigma);
    } while((fabs(previous-current) > precision) && (nbr_iterations <= 50));

  return current;
}

int main()
{
    long long i,j;
    boolean[0]=1;
    boolean[1]=1;
    for(i=2;i*i<55000;i++)
    {
        if(!boolean[i])
        {
            for(j=i*i;j<55000;j+=i)
            {
                boolean[j]=1;
            }
        }
    }

    int index=0;
    for(i=2;i<55000;i++)
    {
        if(!boolean[i])
        {
            primes[index++]=i;
        }
    }
    coprimes[0]=1;
    long long global=1;
    for(i=1;i<50000;i++)
    {
        coprimes[i]=2*phi(i+1);
    }

    long long input;
    while(scanf("%lld",&input)==1 && input)
    {
        long long output=0;
        for(i=input-1;i>=0;i--)
        {
            output+=coprimes[i];
        }
        printf("%lld\n",output);
    }
    return 0;
}

    Euler(const Dcm<Type> & dcm) : Vector<Type, 3>()
    {
        Type phi_val = Type(atan2(dcm(2,1), dcm(2,2)));
        Type theta_val = Type(asin(-dcm(2,0)));
        Type psi_val = Type(atan2(dcm(1,0), dcm(0,0)));
        Type pi = Type(M_PI);

        if (fabs(theta_val - pi/2) < 1.0e-3) {
            phi_val = Type(0.0);
            psi_val = Type(atan2(dcm(1,2), dcm(0,2)));
        } else if (Type(fabs(theta_val + pi/2)) < Type(1.0e-3)) {
            phi_val = Type(0.0);
            psi_val = Type(atan2(-dcm(1,2), -dcm(0,2)));
        }

        phi() = phi_val;
        theta() = theta_val;
        psi() = psi_val;
    }

void NSForchheimerTerm<EvalT, Traits>::
evaluateFields(typename Traits::EvalData workset)
{
  for (std::size_t cell=0; cell < workset.numCells; ++cell) {
    for (std::size_t qp=0; qp < numQPs; ++qp) {     
      normV(cell,qp) = 0.0; 
      for (std::size_t i=0; i < numDims; ++i) {
          normV(cell,qp) += V(cell,qp,i)*V(cell,qp,i); 
      } 
      if (normV(cell,qp) > 0)
        normV(cell,qp) = std::sqrt(normV(cell,qp));
      else
        normV(cell,qp) = 0.0;
      for (std::size_t i=0; i < numDims; ++i) {
          ForchTerm(cell,qp,i) = phi(cell,qp)*rho(cell,qp)*F(cell,qp)*normV(cell,qp)*V(cell,qp,i)/std::sqrt(K(cell,qp));
      } 
    }
  }
}

  void mk_array_instantiation::instantiate_rule(const rule& r, rule_set & dest)
  {
    //Reset everything
    selects.reset();
    eq_classes.reset();
    cnt = src_manager->get_counter().get_max_rule_var(r)+1;
    done_selects.reset();
    ownership.reset();

    expr_ref_vector phi(m);
    expr_ref_vector preds(m);
    expr_ref new_head = create_head(to_app(r.get_head()));
    unsigned nb_predicates = r.get_uninterpreted_tail_size();
    unsigned tail_size = r.get_tail_size();
    for(unsigned i=0;i<nb_predicates;i++)
    {
      preds.push_back(r.get_tail(i));
    }
    for(unsigned i=nb_predicates;i<tail_size;i++)
    {
      phi.push_back(r.get_tail(i));
    }

    //Retrieve selects
    for(unsigned i=0;i<phi.size();i++)
      retrieve_selects(phi[i].get());

    //Rewrite the predicates
    expr_ref_vector new_tail(m);
    for(unsigned i=0;i<preds.size();i++)
    {
      new_tail.append(instantiate_pred(to_app(preds[i].get())));
    }
    new_tail.append(phi);
    for(obj_map<expr, var*>::iterator it = done_selects.begin(); it!=done_selects.end(); ++it)
    {
      expr_ref tmp(m);
      tmp = &it->get_key();
      new_tail.push_back(m.mk_eq(it->get_value(), tmp));
    }
    proof_ref pr(m);
    src_manager->mk_rule(m.mk_implies(m.mk_and(new_tail.size(), new_tail.c_ptr()), new_head), pr, dest, r.name());
  }

void Feta_PTT::correct()
{
    // Velocity gradient tensor
    volTensorField L = fvc::grad( U() );

    // Convected derivate term
    volTensorField C = tau_ & L;

    // Twice the rate of deformation tensor
    volSymmTensorField twoD = twoSymm( L );

    // Two phase transport properties treatment
    alpha1f_ = min(max(alpha(), scalar(0)), scalar(1));

    // etaP effective
    etaPEff_ =  (etaP1_/( Foam::pow(scalar(1) + A1_*Foam::pow(0.5*( Foam::sqr(tr(tau_)) - tr(tau_ & tau_) ) * Foam::sqr(lambda1_) / Foam::sqr(etaP1_), a1_), b1_) ) )*alpha1f_ + (etaP2_/( Foam::pow(scalar(1) + A2_*Foam::pow(0.5*( Foam::sqr(tr(tau_)) - tr(tau_ & tau_) ) * Foam::sqr(lambda2_) / Foam::sqr(etaP2_), a2_), b2_) ) )*(scalar(1) - alpha1f_);

    // lambda effective
    lambdaEff_ =  (lambda1_ / (scalar(1)  + epsilon1_*lambda1_*tr(tau_) / etaP1_) )*alpha1f_ + (lambda2_ / (scalar(1)  + epsilon2_*lambda2_*tr(tau_) / etaP2_) )*(scalar(1) - alpha1f_);

    volScalarField epsilon =  alpha1f_*epsilon1_ + (scalar(1) - alpha1f_)*epsilon2_;

    volScalarField zeta =  alpha1f_*zeta1_ + (scalar(1) - alpha1f_)*zeta2_;

     // Stress transport equation
    tmp<fvSymmTensorMatrix> tauEqn
    (
          fvm::ddt(lambdaEff_, tau_)
        + lambdaEff_ * fvm::div(phi(), tau_)
        ==
          etaPEff_ * twoD
        + lambdaEff_ * twoSymm( C )
        - zeta * lambdaEff_ /2* symm( (tau_ & twoD) + (twoD & tau_) )
	  //add "twoSymm"(...), by Martin R. Du, [email protected]
	  //replace "*twoSymm" with "*symm", Martin R. Du, 2015/08/08
        - fvm::Sp( epsilon * lambdaEff_ / etaPEff_ * tr(tau_) + 1, tau_ )
    );

    tauEqn().relax();
    solve(tauEqn);

}

/***********************************************************************//**
 * @brief Sum effective area multiplied by livetime over zenith angles
 *
 * @param[in] dir True sky direction.
 * @param[in] energy True photon energy.
 * @param[in] offset Offset from true direction (deg).
 * @param[in] psf Point spread function.
 * @param[in] aeff Effective area.
 *
 * Computes
 * \f[\sum_{\cos \theta} T_{\rm live}(\cos \theta) 
 *    PSF(\log E, \delta, \cos \theta) A_{\rm eff}(\cos \theta, \phi)\f]
 * where
 * \f$T_{\rm live}(\cos \theta)\f$ is the livetime as a function of the
 * cosine of the zenith angle, and
 * \f$PSF(\log E, \delta, \cos \theta)\f$ is the point spread function that
 * depends on
 * the log10 of the energy (in MeV),
 * the offset angle from the true direction (in degrees), and
 * the cosine of the zenith angle, and
 * \f$A_{\rm eff}(\cos \theta, \phi)\f$ is the effective area that depends
 * on the cosine of the zenith angle and (optionally) of the azimuth angle.
 ***************************************************************************/
double GLATLtCubeMap::operator()(const GSkyDir& dir, const GEnergy& energy,
                                 const double& offset, const GLATPsf& psf,
                                 const GLATAeff& aeff) const
{
    // Get map index
    int pixel = m_map.dir2pix(dir);

    // Initialise sum
    double sum = 0.0;

    // Circumvent const correctness
    GLATPsf*  fpsf  = const_cast<GLATPsf*>(&psf);
    GLATAeff* faeff = const_cast<GLATAeff*>(&aeff);

    // Get log10 of energy
    double logE = energy.log10MeV();

    // If livetime cube and response have phi dependence then sum over
    // zenith and azimuth. Note that the map index starts with m_num_ctheta
    // as the first m_num_ctheta maps correspond to an evaluation without
    // any phi-dependence.
    if (has_phi() && aeff.has_phi()) {
        for (int iphi = 0, i = m_num_ctheta; iphi < m_num_phi; ++iphi) {
            double p = phi(iphi);
            for (int itheta = 0; itheta < m_num_ctheta; ++itheta, ++i) {
                sum += m_map(pixel, i) * (*faeff)(logE, costheta(i), p) *
                       (*fpsf)(offset, logE, costheta(i));
            }
        }
    }

    // ... otherwise sum only over zenith angle
    else {
        for (int i = 0; i < m_num_ctheta; ++i) {
            sum += m_map(pixel, i) * (*faeff)(logE, costheta(i)) *
                   (*fpsf)(offset, logE, costheta(i));
        }
    }

    // Return sum
    return sum;
}

void FENE_P::correct()
{
    // Velocity gradient tensor
    volTensorField L = fvc::grad(U());

    // Convected derivate term
    volTensorField C = tau_ & L;

    // Twice the rate of deformation tensor
    volSymmTensorField twoD = twoSymm(L);

    // Two phase transport properties treatment
    volScalarField alpha1f = 
		min(max(alpha(),scalar(0)),scalar(1));

    volScalarField lambda =  
		alpha1f*lambda1_ + (scalar(1) - alpha1f)*lambda2_;

    volScalarField etaP =  
		alpha1f*etaP1_ + (scalar(1) - alpha1f)*etaP2_;

    volScalarField Lquad =  
		alpha1f*Lquad1_ + (scalar(1) - alpha1f)*Lquad2_;

     // Stress transport equation
    tmp<fvSymmTensorMatrix> tauEqn
    (
		fvm::ddt(lambda, tau_)
	  + lambda * fvm::div(phi(), tau_)
	 ==
		( 1 / (1 - 3/Lquad) ) * etaP * twoD
	  + lambda * twoSymm( C )
	  - fvm::Sp
		( 
			1  + (3/(1 - 3/Lquad) + lambda*tr(tau_)/etaP )/Lquad, 
			tau_ 
		)
    );

    tauEqn().relax();
    solve(tauEqn);
}

void ContinuumMechanicsNeumannBcsAssembler<DIM>::AssembleOnBoundaryElement(BoundaryElement<DIM-1,DIM>& rBoundaryElement,
                                                                           c_vector<double,STENCIL_SIZE>& rBelem,
                                                                           unsigned boundaryConditionIndex)
{
    rBelem.clear();

    c_vector<double, DIM> weighted_direction;
    double jacobian_determinant;
    mpMesh->GetWeightedDirectionForBoundaryElement(rBoundaryElement.GetIndex(), weighted_direction, jacobian_determinant);

    c_vector<double,NUM_NODES_PER_ELEMENT> phi;

    for (unsigned quad_index=0; quad_index<mpQuadRule->GetNumQuadPoints(); quad_index++)
    {
        double wJ = jacobian_determinant * mpQuadRule->GetWeight(quad_index);
        const ChastePoint<DIM-1>& quad_point = mpQuadRule->rGetQuadPoint(quad_index);
        QuadraticBasisFunction<DIM-1>::ComputeBasisFunctions(quad_point, phi);

        c_vector<double,DIM> traction = zero_vector<double>(DIM);
        switch (mpProblemDefinition->GetTractionBoundaryConditionType())
        {
            case ELEMENTWISE_TRACTION:
            {
                traction = mpProblemDefinition->rGetElementwiseTractions()[boundaryConditionIndex];
                break;
            }
            default:
                // Functional traction not implemented yet..
                NEVER_REACHED;
        }

        for (unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
        {
            unsigned spatial_dim = index%DIM;
            unsigned node_index = (index-spatial_dim)/DIM;

            assert(node_index < NUM_NODES_PER_ELEMENT);

            rBelem(index) += traction(spatial_dim) * phi(node_index) * wJ;
        }
    }
}

/******************************************************************
* It computes a sum of a series of values, as needed by formula   *
* 11.25 of "Inside Solid State Drives" book. It helps performing  *
* Horizontal step of Log-Domain Belief Propagation algorithm.     *
******************************************************************/
double log_bp_error_corrector::phi(double** R, int i, int j)
{
	double result;

	result = 0;

	std::list<edge> edges = chk_nodes[j].get_edges();
	std::list<edge>::iterator it;

	for(it = edges.begin(); it != edges.end(); it++){

		int l = (*it).get_dest().get_idx();

		if(l != i){
			result += phi(abs(R[l][j]));
		}
	}

	return result;
}

void LevelSetViewerWidget::openFile(const QString & filename)
{
    m_ready = false;
    qWarning() << "Loading file: "<< filename;
    if(m_obj){delete m_obj;}
    std::ifstream in(filename.toStdString().c_str());
    if(m_grid){delete m_grid;}
    m_grid = new Gridf();
    m_grid->read(in);
    std::cout << m_grid->NI() << " " << m_grid->NJ() << " " << m_grid->NK() << std::endl;
    MarchingSquaresGridf phi(*m_grid); 

    m_obj = new Mesh2f();
    phi.populateMesh(*m_obj,0);
    MeshUtils::normalize(*m_obj);
    std::cout << m_obj->vertices.size() << std::endl;
    std::cout << m_obj->faces.size() << std::endl;
    sendBuffers();

}

inline Polynomial<eT>
polyfit(const eT* x, const eT* y, const uword length, const uword order)
{
    const uword size = order + 1;
    Col<eT> Q(size);   //parameter to be estimated.
    Q.fill(eT(0));        //initialized with zeros.
    Col<eT> phi(size); //
    Mat<eT> P = eye(size, size) * 1e6;
    Col<eT> K;
    Mat<eT> I = eye(size, size);

    for(uword i = 0; i < length; ++i)
    {
        detail::gen_phi(phi, x[i]);
        K = P * phi / (1 + as_scalar(phi.st() * P * phi)); //K(t) = P(t-1)*phi(t)*inv(I+phi'*P(t-1)*phi).
        P = (I - K * phi.st()) * P; //P(t) = (I-K(t)*phi')*P(t-1).
        Q = Q + K * (y[i] - as_scalar(phi.st() * Q)); //Q(t) = Q(t-1) + K(t)*(y(t)-phi'Q(t-1)).
    }
    return std::move(Polynomial<eT>(Q));
}

double PDPS::logpri()const{
  const Mat & Nu(m_->Nu());
  uint d = nrow(Nu);
  Vec phi(d);
  double ans=0;
  for(uint i=0; i<d; ++i){
    double a = sum(Nu.row(i));
    if(a <=0) return BOOM::negative_infinity();
    phi = Nu.row(i);
    for(uint j=0; j<d; ++j) {
      if(phi[j] < min_nu_) {
        return BOOM::negative_infinity();
      }
    }
    phi/=a;
    ans += alpha_row_prior_[i]->logp(a);
    ans += phi_row_prior_[i]->logp(phi);
    ans -= (d-1) * log(a);
  }
  return ans;
}

// for co-ordinates at box nn, ll, evaluates the value using
// coefficients at box n, l
// this function assumes that x is Vector.
double* FunctionAST::evaluate_at_box(int k, double* coeff, int n, int l, int nn, int ll, double* x)
{
    int         i;
    double*     p;
    double      coordinate;
    double      temp;
    double*     returnVector;


    returnVector    = Vector_initialize(NULL, 0, this->func->k+ADD_K);
    for (i=0; i<this->func->k; i++)
    {
        coordinate          = ((x[i] + ll) * (pow(2.0, n - nn))) - l;
        p                   = Vector_initialize(phi(coordinate, this->func->k), 0, this->func->k);
        temp                = Vector_inner(coeff, p, this->func->k);
        returnVector[i]     = temp * sqrt(pow(2.0, n));

        free(p);
    }
    return returnVector;
}