diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/gop_nmf_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/gop_nmf_impl.h
index 28c375e2da..da47bd5d17 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/gop_nmf_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/gop_nmf_impl.h
@@ -706,7 +706,7 @@ bool RelaxedRescaledNmfL1::IsIntermediateStepOver(Matrix &coordinates,
   double norm_gradient=la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr());
-  if (step/std::pow(2, norm_gradient) < grad_tolerance_) {
+  if (math::Pow<1,2>(norm_gradient)*step < grad_tolerance_) {
     return true;
   }
   return false;
@@ -1269,7 +1269,7 @@ bool RelaxedNmfIsometric::IsOptimizationOver(Matrix &coordinates,
 bool RelaxedNmfIsometric::IsIntermediateStepOver(Matrix &coordinates, 
                                         Matrix &gradient, 
                                         double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   if (norm_gradient*step < grad_tolerance_) {
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/sdp_objectives_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
index b4b1a62d61..5866af2b2b 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
@@ -87,7 +87,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
       //double w=coordinates.get(j, w_i);
       double t1=coordinates.get(j, t1_i);
       double dw=-1.0;
-      double dt1=0.5/std::pow(-2,t1);
+      double dt1=0.5*math::Pow<-1,2>(t1);
       gradient->set(j, w_i, dw);
       gradient->set(j, t1_i, dt1);
     }
@@ -100,7 +100,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
      // double h=coordinates.get(j, h_i);
       double t2=coordinates.get(j, t2_i);
       double dh=-1.0;
-      double dt2=0.5/std::pow(-2,t2);
+      double dt2=0.5*math::Pow<-1,2>(t2);
       gradient->set(j, h_i, dh);
       gradient->set(j, t2_i, dt2);
     }
@@ -134,7 +134,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
     }     
   } 
   // from the SDP cones
-  // determinant=(t1-w*w)*(t2-h*h)-std::pow(2,w*h-v);
+  // determinant=(t1-w*w)*(t2-h*h)-math::Pow<2,1>(w*h-v);
   // dw=-2*w*(t2-h*h)-2*h*(w*h-v);
   // dh=-2*h*(t1-w*w)-2*w*(w*h-v);
   // dt1=t2-h*h;
@@ -156,7 +156,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
       double t2_minus_hh=(t2-h*h);
       double wh_minus_v=(w*h-v);
       double determinant=t1_minus_ww*t2_minus_hh
-        -std::pow(2,wh_minus_v);
+        -math::Pow<2,1>(wh_minus_v);
       
       //NOTIFY("determinant:%lg", determinant);
       if (determinant<=0.0) {
@@ -196,7 +196,7 @@ void SmallSdpNmf::ComputeObjective(Matrix &coordinates, double *objective) {
     for(index_t j=0; j<new_dim_; j++) {
       double w=coordinates.get(j, w_i);
       double t1=coordinates.get(j, t1_i);
-      *objective+=1/std::pow(2,t1)-w;
+      *objective+=math::Pow<1,2>(t1)-w;
     }
   }
   for(index_t i=0; i<num_of_columns_; i++) {
@@ -205,7 +205,7 @@ void SmallSdpNmf::ComputeObjective(Matrix &coordinates, double *objective) {
     for(index_t j=0; j<new_dim_; j++) {
       double h=coordinates.get(j, h_i);
       double t2=coordinates.get(j, t2_i);
-      *objective+=1/std::pow(2,t2)-h;
+      *objective+=math::Pow<1,2>(t2)-h;
     }
   }
   //Vector temp;
@@ -250,7 +250,7 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
   }
   lagrangian-=log(temp_prod);   
   // from the SDP cones
-  // determinant=(t1-w*w)*(t2-h*h)-std::pow(2,w*h-v);
+  // determinant=(t1-w*w)*(t2-h*h)-math::Pow<2,1>(w*h-v);
   temp_prod=1.0;
  for(index_t i=0; i<values_.size(); i++) {
     index_t w_i=rows_[i];
@@ -275,7 +275,7 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
       }
       double wh_minus_v=(w*h-v);
       double determinant=t1_minus_ww*t2_minus_hh
-        -std::pow(2,wh_minus_v);
+        -math::Pow<2,1>(wh_minus_v);
       if (unlikely(determinant<=0)) {
         return DBL_MAX;
       }
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_classification_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_classification_impl.h
index b3295634bf..f46d162fc0 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_classification_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_classification_impl.h
@@ -81,7 +81,7 @@ void MaxFurthestNeighborsSemiSupervised::Init(fx_module *module, Matrix &labeled
   }
   sum_of_nearest_distances_=0;
   for(index_t i=0; i<nearest_distances_.size(); i++) {
-   sum_of_nearest_distances_+=1.0/std::pow(2,nearest_distances_[i]);
+   sum_of_nearest_distances_+=math::Pow<1,2>(nearest_distances_[i]);
   }
   NOTIFY("Sum of all nearest distances:%lg", sum_of_nearest_distances_);
   fx_result_double(module_, "sum_of_nearest_distances", sum_of_nearest_distances_);
@@ -323,7 +323,7 @@ bool MaxFurthestNeighborsSemiSupervised::IsOptimizationOver(
 
 bool MaxFurthestNeighborsSemiSupervised::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -445,7 +445,7 @@ void MaxFurthestNeighborsSvmSemiSupervised::Init(fx_module *module,
   }
   sum_of_nearest_distances_=0;
   for(index_t i=0; i<nearest_distances_.size(); i++) {
-   sum_of_nearest_distances_+=1.0/std::pow(2,nearest_distances_[i]);
+   sum_of_nearest_distances_+=math::Pow<1,2>(nearest_distances_[i]);
   }
   NOTIFY("Sum of all nearest distances:%lg", sum_of_nearest_distances_);
   fx_result_double(module_, "sum_of_nearest_distances", sum_of_nearest_distances_);
@@ -718,7 +718,7 @@ double MaxFurthestNeighborsSvmSemiSupervised::ComputeLagrangian(Matrix &coordina
       if (sigma1_ * ineq  <= ineq_lagrange_mult_[i*num_of_labeled_+j]) {
         lagrangian+=(-ineq_lagrange_mult_[i*num_of_labeled_+j] + sigma1_/2*ineq)*ineq;
       } else {
-        lagrangian+=-std::pow(2,ineq_lagrange_mult_[i*num_of_labeled_+j])/(2*sigma1_);
+        lagrangian+=-math::Pow<2,1>(ineq_lagrange_mult_[i*num_of_labeled_+j])/(2*sigma1_);
       }
     }   
   }
@@ -788,7 +788,7 @@ bool MaxFurthestNeighborsSvmSemiSupervised::IsOptimizationOver(
 
 bool MaxFurthestNeighborsSvmSemiSupervised::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -912,7 +912,7 @@ void MaxFurthestNeighborsSvmSemiSupervised1::Init(fx_module *module,
   }
   sum_of_nearest_distances_=0;
   for(index_t i=0; i<nearest_distances_.size(); i++) {
-   sum_of_nearest_distances_+=1.0/std::pow(2,nearest_distances_[i]);
+   sum_of_nearest_distances_+=math::Pow<1,2>(nearest_distances_[i]);
   }
   NOTIFY("Sum of all nearest distances:%lg", sum_of_nearest_distances_);
   fx_result_double(module_, "sum_of_nearest_distances", sum_of_nearest_distances_);
@@ -1185,7 +1185,7 @@ double MaxFurthestNeighborsSvmSemiSupervised1::ComputeLagrangian(Matrix &coordin
       if (sigma1_ * ineq  <= ineq_lagrange_mult_[i*num_of_labeled_+j]) {
         lagrangian+=(-ineq_lagrange_mult_[i*num_of_labeled_+j] + sigma1_/2*ineq)*ineq;
       } else {
-        lagrangian+=-std::pow(2,ineq_lagrange_mult_[i*num_of_labeled_+j])/(2*sigma1_);
+        lagrangian+=-math::Pow<2,1>(ineq_lagrange_mult_[i*num_of_labeled_+j])/(2*sigma1_);
       }
     }   
   }
@@ -1255,7 +1255,7 @@ bool MaxFurthestNeighborsSvmSemiSupervised1::IsOptimizationOver(
 
 bool MaxFurthestNeighborsSvmSemiSupervised1::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_objectives_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_objectives_impl.h
index 5f5a20ec5f..26b6d6d12e 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_objectives_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/mvu/mvu_objectives_impl.h
@@ -564,7 +564,7 @@ void MaxFurthestNeighbors::Init(datanode *module, Matrix &data) {
         &nearest_distances_,
         &num_of_nearest_pairs_);
   }
-  sum_of_nearest_distances_=1.0/std::pow(2,la::Dot(nearest_distances_.size(), 
+  sum_of_nearest_distances_=math::Pow<1,2>(la::Dot(nearest_distances_.size(), 
       &nearest_distances_[0], &nearest_distances_[0]));
   fx_format_result(module_, "num_of_constraints", "%i", num_of_nearest_pairs_);
   eq_lagrange_mult_.Init(num_of_nearest_pairs_);
@@ -633,7 +633,7 @@ void MaxFurthestNeighbors::Init(fx_module *module) {
   }
   num_of_points_++;
   num_of_nearest_pairs_=nearest_neighbor_pairs_.size();
-  sum_of_nearest_distances_=1.0/std::pow(2,la::Dot(nearest_distances_.size(), 
+  sum_of_nearest_distances_=math::Pow<1,2>(la::Dot(nearest_distances_.size(), 
       &nearest_distances_[0], &nearest_distances_[0]));
 
   fclose(fp);
@@ -743,7 +743,7 @@ void MaxFurthestNeighbors::ComputeFeasibilityError(Matrix &coordinates, double *
                            -nearest_distances_[i];
     *error+=dist_diff*dist_diff;
   }
-  *error= 100.0 / std::pow(2,*error)/sum_of_nearest_distances_;
+  *error= 100 * math::Pow<1,2>(*error)/sum_of_nearest_distances_;
 }
 
 double MaxFurthestNeighbors::ComputeLagrangian(Matrix &coordinates) {
@@ -825,7 +825,7 @@ bool MaxFurthestNeighbors::IsOptimizationOver(Matrix &coordinates,
 
 bool MaxFurthestNeighbors::IsIntermediateStepOver(Matrix &coordinates, 
       Matrix &gradient, double step) {
-   double norm_gradient=1/std::pow(2,la::Dot(gradient.n_elements(), 
+   double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -894,7 +894,7 @@ void MaxVarianceUtils::EstimateKnns(ArrayList<index_t> &neares_neighbors,
       double probability=0;
       for(index_t j=0; j<k; j++) {
         probability+=exp(-nearest_distances[i*maximum_knns+j]
-            *std::pow(2,nearest_distances[i*maximum_knns+k]-1))/scale_factor;
+            /(2*math::Pow<1,2>(nearest_distances[i*maximum_knns+k])))/scale_factor;
       }
       loocv_score+=log(probability);
       mean_band+=nearest_distances[i*maximum_knns+k];
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/nmf_tree/relaxed_nmf_bound_tightener_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/nmf_tree/relaxed_nmf_bound_tightener_impl.h
index 8033b08651..114ca6e7aa 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/nmf_tree/relaxed_nmf_bound_tightener_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/nmf_tree/relaxed_nmf_bound_tightener_impl.h
@@ -665,7 +665,7 @@ bool RelaxedNmfIsometricBoundTightener::IsOptimizationOver(Matrix &coordinates,
 bool RelaxedNmfIsometricBoundTightener::IsIntermediateStepOver(Matrix &coordinates, 
                                         Matrix &gradient, 
                                         double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   if (norm_gradient*step < grad_tolerance_) {
@@ -1055,7 +1055,7 @@ bool RelaxedNmfIsometricBoxTightener::IsOptimizationOver(Matrix &coordinates,
 bool RelaxedNmfIsometricBoxTightener::IsIntermediateStepOver(Matrix &coordinates, 
                                         Matrix &gradient, 
                                         double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   if (norm_gradient*step < grad_tolerance_) {
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h
index 10781cd631..6f0f4b8b62 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h
@@ -134,7 +134,7 @@ void NonConvexMVU::ComputeLocalOptimumSGD() {
   index_t exterior_points=0;
   for(index_t kk=0; kk<3; kk++) {
   for(it1=0; it1<100; it1++) {  
-    step=1.0/std::pow(2,step0/(it1+1));
+    step=math::Pow<1,2>(step0/(it1+1));
     for(index_t it2=0; it2<max_iterations_; it2++) {
       for(index_t i=0; i<num_of_points_; i++) {
         for(index_t k=0; k<kfns_; k++) {
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/main.cc b/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/main.cc
index 3be02a7b8f..42e25c27ca 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/main.cc
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/main.cc
@@ -116,7 +116,7 @@ void ComputeSpectrum(Matrix &w_mat, Matrix &h_mat, Matrix *spectrum) {
   double total_sum=0;
   for(index_t i=0; i<w_mat.n_rows(); i++) {
     for(index_t j=0; j<w_mat.n_cols(); j++) {
-      spectrum->set(i, 0, spectrum->get(i, 0)+std::pow(2,w_mat.get(i, j)));
+      spectrum->set(i, 0, spectrum->get(i, 0)+math::Pow<2,1>(w_mat.get(i, j)));
     }
   }
   Matrix h_norms;
@@ -128,7 +128,7 @@ void ComputeSpectrum(Matrix &w_mat, Matrix &h_mat, Matrix *spectrum) {
     }
   }
   for(index_t i=0; i<h_norms.n_rows(); i++) {
-    spectrum->set(i, 0, spectrum->get(i, 0)*std::pow(2,h_norms.get(i,0)));
+    spectrum->set(i, 0, spectrum->get(i, 0)/math::Pow<1,2>(h_norms.get(i,0)));
   }
   std::sort(spectrum->ptr(), 
             spectrum->ptr()+spectrum->n_rows(), 
@@ -145,12 +145,12 @@ void ComputeNeighborhoodErrors(ArrayList<std::pair<index_t, index_t> > &neighbor
   for(index_t i=0; i<neighbor_pairs.size(); i++) {
     index_t n1=neighbor_pairs[i].first;
     index_t n2=neighbor_pairs[i].second;
-    *error += std::pow(2,la::DistanceSqEuclidean(new_points.n_rows(), 
+    *error += math::Pow<2,1>(la::DistanceSqEuclidean(new_points.n_rows(), 
                                      new_points.GetColumnPtr(n1),
                                      new_points.GetColumnPtr(n2))
                      -distances[i]);
       total_distances += distances[i]*distances[i];
   }
-  *error=100.0 / std::pow(2,*error/total_distances);
+  *error=100.0 * math::Pow<1,2>(*error/total_distances);
 }
   
diff --git a/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/nmf_objectives_impl.h b/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/nmf_objectives_impl.h
index f0b9cee774..c2f9c74f62 100644
--- a/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/nmf_objectives_impl.h
+++ b/fastlib/branches/fastlib-stl/contrib/nvasil/non_convex_nmf/nmf_objectives_impl.h
@@ -310,7 +310,7 @@ double BigSdpNmfObjectiveMinVarIneq::ComputeLagrangian(Matrix &coordinates) {
 		if (sigma_*diff <= ineq_lagrange_mult_[i]) {
       lagrangian+=(sigma_*diff-ineq_lagrange_mult_[i])*diff;
     } else {
-      lagrangian+=-std::pow(2,ineq_lagrange_mult_[i])/sigma_;
+      lagrangian+=-math::Pow<2,1>(ineq_lagrange_mult_[i])/sigma_;
     }
 	}
 	return lagrangian;
@@ -529,7 +529,7 @@ double BigSdpNmfObjectiveMinVarDiagonalDominance::ComputeLagrangian(Matrix &coor
   	if (sigma2_*diff <= ineq_lagrange_mult_[i]) {
       lagrangian+=(sigma2_*diff/2-ineq_lagrange_mult_[i])*diff;
     } else {
-      lagrangian+=-std::pow(2,ineq_lagrange_mult_[i])/(2*sigma2_);
+      lagrangian+=-math::Pow<2,1>(ineq_lagrange_mult_[i])/(2*sigma2_);
     }
   }
   
@@ -540,7 +540,7 @@ double BigSdpNmfObjectiveMinVarDiagonalDominance::ComputeLagrangian(Matrix &coor
       if (sigma2_*diff < ineq_lagrange_mult1_[i*new_dim_+j]) {
         lagrangian+=(sigma2_*diff-ineq_lagrange_mult1_[i*new_dim_+j])*diff;
       } else {
-        lagrangian+=-std::pow(2,
+        lagrangian+=-math::Pow<2,1>(
             ineq_lagrange_mult1_[i*new_dim_+j])/(2*sigma2_);
       }
     }
@@ -807,8 +807,8 @@ void BigSdpNmfObjectiveMaxVarIsometric::ComputeFeasibilityError(Matrix &coordina
 												coordinates.GetColumnPtr(h))-values_[i];
 		*error+=diff*diff;
   }
-  infeasibility1_=100.0/std::pow(2,*error/v_norm_);
-  printf("dot_infeasibility:%lg%% ", 100.0/std::pow(2,*error/v_norm_)); ;
+  infeasibility1_=math::Pow<1,2>(*error/v_norm_)*100;
+  printf("dot_infeasibility:%lg%% ", math::Pow<1,2>(*error/v_norm_)*100); ;
   // local isometry
   double error2=0;
   for(index_t i=0; i<num_of_nearest_pairs_; i++) {
@@ -824,8 +824,8 @@ void BigSdpNmfObjectiveMaxVarIsometric::ComputeFeasibilityError(Matrix &coordina
                                           -nearest_distances_[i]);
 
   }
-  infeasibility2_=100.0/std::pow(2,error2/sum_all_distances_);
-  printf("dist_infeasibility:%lg %%\n", 100.0/std::pow(2,error2/sum_all_distances_));
+  infeasibility2_=math::Pow<1,2>(error2/sum_all_distances_)*100;
+  printf("dist_infeasibility:%lg %%\n", math::Pow<1,2>(error2/sum_all_distances_)*100);
   printf("sigma_ratio:%lg\n", sigma_ratio_);
 //  sigma_ratio_=infeasibility2_/infeasibility1_;
 }
@@ -939,7 +939,7 @@ bool BigSdpNmfObjectiveMaxVarIsometric::IsOptimizationOver(
 
 bool BigSdpNmfObjectiveMaxVarIsometric::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -1166,7 +1166,7 @@ void BigSdpNmfObjectiveMaxVarIsometricRemoveMean::ComputeFeasibilityError(Matrix
 												coordinates.GetColumnPtr(h))-values_[i];
 		*error+=diff*diff;
   }
-  infeasibility1_=1.0/std::pow(2,*error/v_norm_)*100;
+  infeasibility1_=math::Pow<1,2>(*error/v_norm_)*100;
   printf("dot_infeasibility:%lg%% ", infeasibility1_); ;
   // local isometry
   double error2=0;
@@ -1183,7 +1183,7 @@ void BigSdpNmfObjectiveMaxVarIsometricRemoveMean::ComputeFeasibilityError(Matrix
                                           -nearest_distances_[i]);
 
   }
-  infeasibility2_=1.0/std::pow(2,error2/sum_all_distances_)*100;
+  infeasibility2_=math::Pow<1,2>(error2/sum_all_distances_)*100;
   printf("dist_infeasibility:%lg %%\n", infeasibility2_);
 }
 
@@ -1291,7 +1291,7 @@ bool BigSdpNmfObjectiveMaxVarIsometricRemoveMean::IsOptimizationOver(
 
 bool BigSdpNmfObjectiveMaxVarIsometricRemoveMean::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -1524,7 +1524,7 @@ void BigSdpNmfObjectiveMaxFurthestIsometric::ComputeFeasibilityError(Matrix &coo
 												coordinates.GetColumnPtr(h))-values_[i];
 		*error+=diff*diff;
   }
-  infeasibility1_=100.0/std::pow(2,*error/v_norm_);
+  infeasibility1_=math::Pow<1,2>(*error/v_norm_)*100;
   // local isometry
   double error2=0;
   for(index_t i=0; i<num_of_nearest_pairs_; i++) {
@@ -1540,7 +1540,7 @@ void BigSdpNmfObjectiveMaxFurthestIsometric::ComputeFeasibilityError(Matrix &coo
                                           -nearest_distances_[i]);
 
   }
-  infeasibility2_=100.0/std::pow(2,error2/sum_all_distances_);
+  infeasibility2_=math::Pow<1,2>(error2/sum_all_distances_)*100;
   NOTIFY("dot_infeasibility:%lg%% dist_infeasibility:%lg %%", 
       infeasibility1_, infeasibility2_);
   // NOTIFY("sigma_ratio:%lg\n", sigma_ratio_);
@@ -1656,7 +1656,7 @@ bool BigSdpNmfObjectiveMaxFurthestIsometric::IsOptimizationOver(
 
 bool BigSdpNmfObjectiveMaxFurthestIsometric::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -1856,7 +1856,7 @@ void NmfObjectiveIsometric::ComputeFeasibilityError(Matrix &coordinates,
                                           -nearest_distances_[i]);
 
   }
-  infeasibility2_=100.0/std::pow(2,*error/sum_all_distances_);
+  infeasibility2_=math::Pow<1,2>(*error/sum_all_distances_)*100;
   printf("dist_infeasibility:%lg %%\n", infeasibility2_*100);
 }
 
@@ -1933,7 +1933,7 @@ bool NmfObjectiveIsometric::IsOptimizationOver(
 
 bool NmfObjectiveIsometric::IsIntermediateStepOver(
     Matrix &coordinates, Matrix &gradient, double step) {
-  double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+  double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
diff --git a/fastlib/branches/fastlib-stl/fastlib/fastlib.h b/fastlib/branches/fastlib-stl/fastlib/fastlib.h
index 17e315c540..0289f1c9f6 100644
--- a/fastlib/branches/fastlib-stl/fastlib/fastlib.h
+++ b/fastlib/branches/fastlib-stl/fastlib/fastlib.h
@@ -38,7 +38,6 @@
 #ifndef FASTLIB_FASTLIB_H
 #define FASTLIB_FASTLIB_H
 
-#include <cmath>
 #include "base/base.h"
 #include "la/la.h"
 #include "la/matrix.h"
diff --git a/fastlib/branches/fastlib-stl/fastlib/la/la.h b/fastlib/branches/fastlib-stl/fastlib/la/la.h
index c6018f63e9..cd0cfc34bf 100644
--- a/fastlib/branches/fastlib-stl/fastlib/la/la.h
+++ b/fastlib/branches/fastlib-stl/fastlib/la/la.h
@@ -91,7 +91,7 @@ namespace la {
     double s = 0;
     do {
       double d = *va++ - *vb++;
-      s += std::pow(t_pow,d);
+      s += math::PowAbs<t_pow, 1>(d);
     } while (--length);
     return s;
   }
@@ -107,7 +107,7 @@ namespace la {
   template<int t_pow>
   inline double LMetric(
       index_t length, const double *va, const double *vb) {
-    return 1.0/std::pow(t_pow, RawLMetric<t_pow>(length, va, vb));
+    return math::Pow<1, t_pow>(RawLMetric<t_pow>(length, va, vb));
   }
   /** Finds the trace of the matrix.
    *  Trace(A) is the sum of the diagonal elements
diff --git a/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/lbfgs_impl.h b/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/lbfgs_impl.h
index f225eea7ec..a977d69653 100644
--- a/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/lbfgs_impl.h
+++ b/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/lbfgs_impl.h
@@ -620,7 +620,7 @@ std::string Lbfgs<OptimizedFunction>::ComputeProgress_() {
   double feasibility_error;
   optimized_function_->ComputeFeasibilityError(
       coordinates_, &feasibility_error);
-  double norm_grad = 1.0 / std::pow(2, la::Dot(gradient_.n_elements(), 
+  double norm_grad=math::Pow<1,2>(la::Dot(gradient_.n_elements(), 
       gradient_.ptr(), gradient_.ptr()));
   char buffer[1024];
   sprintf(buffer, "iteration:%i sigma:%lg lagrangian:%lg objective:%lg error:%lg "
diff --git a/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/optimization_utils.h b/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/optimization_utils.h
index 7c3aa13a48..aae2e89ef6 100644
--- a/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/optimization_utils.h
+++ b/fastlib/branches/fastlib-stl/fastlib/optimization/lbfgs/optimization_utils.h
@@ -170,8 +170,8 @@ class OptUtils {
     ones.SetAll(1.0);
     // This part of the sparsity constraint function formula can be
     // precomputed and it is the same for every iteration
-    double precomputed_sparse_factor=-sparse_factor*(1.0/std::pow(2,dimension)-1)+
-          1.0/std::pow(2,dimension);
+    double precomputed_sparse_factor=-sparse_factor*(math::Pow<1,2>(dimension)-1)+
+          math::Pow<1,2>(dimension);
 
     for (index_t i=0; i<data->n_cols(); i++) {
       double *point=data->GetColumnPtr(i);
@@ -199,7 +199,7 @@ class OptUtils {
         double w_norm=la::LengthEuclidean(w_vector);
         double w_times_v = 2*la::Dot(v_vector, w_vector);
         double v_norm_minus_l2=la::LengthEuclidean(v_vector)-l2_norm;
-        double alpha = (-w_times_v+1.0/std::pow(2,w_times_v*w_times_v
+        double alpha = (-w_times_v+math::Pow<1,2>(w_times_v*w_times_v
             -4*w_norm*v_norm_minus_l2))/(2*w_norm);
         la::AddExpert(alpha, w_vector, &v_vector);
         bool all_positive=true;
diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h b/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h
index 0a21f791ac..c4987d7f93 100644
--- a/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h
+++ b/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h
@@ -76,7 +76,7 @@ double DBallBound<TMetric, TPoint>::MinDistance(const Point& point) const {
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MinDistanceSq(const Point& point) const {
-  return std::pow(2,MinDistance(point));
+  return math::Pow<2, 1>(MinDistance(point));
 }
 
 /**
@@ -90,7 +90,7 @@ double DBallBound<TMetric, TPoint>::MinDistance(const DBallBound& other) const {
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MinDistanceSq(const DBallBound& other) const {
-  return std::pow(2,MinDistance(other));
+  return math::Pow<2, 1>(MinDistance(other));
 }
 
 /**
@@ -103,7 +103,7 @@ double DBallBound<TMetric, TPoint>::MaxDistance(const Point& point) const {
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MaxDistanceSq(const Point& point) const {
-  return std::pow(2,MaxDistance(point));
+  return math::Pow<2, 1>(MaxDistance(point));
 }
 
 /**
@@ -116,7 +116,7 @@ double DBallBound<TMetric, TPoint>::MaxDistance(const DBallBound& other) const {
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MaxDistanceSq(const DBallBound& other) const {
-  return std::pow(2, MaxDistance(other));
+  return math::Pow<2, 1>(MaxDistance(other));
 }
 
 /**
@@ -138,8 +138,8 @@ DRange DBallBound<TMetric, TPoint>::RangeDistanceSq(const DBallBound& other) con
   double delta = MidDistance(other.center_);
   double sumradius = radius_ + other.radius_;
   return DRange(
-      std::pow(2,math::ClampNonNegative(delta - sumradius)),
-      std::pow(2,delta + sumradius));
+      math::Pow<2, 1>(math::ClampNonNegative(delta - sumradius)),
+      math::Pow<2, 1>(delta + sumradius));
 }
 
 /**
@@ -161,7 +161,7 @@ double DBallBound<TMetric, TPoint>::MinToMid(const DBallBound& other) const {
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MinToMidSq(const DBallBound& other) const {
-  return std::pow(2, MinToMid(other));
+  return math::Pow<2, 1>(MinToMid(other));
 }
 
 /**
@@ -175,7 +175,7 @@ double DBallBound<TMetric, TPoint>::MinimaxDistance(const DBallBound& other) con
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MinimaxDistanceSq(const DBallBound& other) const {
-  return std::pow(2, MinimaxDistance(other));
+  return math::Pow<2, 1>(MinimaxDistance(other));
 }
 
 /**
@@ -197,7 +197,7 @@ double DBallBound<TMetric, TPoint>::MidDistance(const DBallBound& other) const {
 
 template<typename TMetric, typename TPoint>
 double DBallBound<TMetric, TPoint>::MidDistanceSq(const DBallBound& other) const {
-  return std::pow(2, MidDistance(other));
+  return math::Pow<2, 1>(MidDistance(other));
 }
 
 template<typename TMetric, typename TPoint>
diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h b/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h
index 8df87da255..d0c2c6f62e 100644
--- a/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h
+++ b/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h
@@ -271,6 +271,23 @@ double DHrectBound<t_pow>::MaxDistanceSq(const vec& point) const {
   return pow(sum, 2.0 / (double) t_pow);
 }
 
+/**
+ * Calculates maximum bound-to-point squared distance.
+ */
+/*
+template<int t_pow>
+double DHrectBound<t_pow>::MaxDistanceSq(const double *point) const {
+  double sum = 0;
+
+  for (index_t d = 0; d < dim_; d++) {
+    double v = std::max(point[d] - bounds_[d].lo, bounds_[d].hi - point[d]);
+    sum += math::Pow<t_pow, 1>(v); // v is non-negative
+  }
+
+  return math::Pow<2, t_pow>(sum);
+}
+*/
+
 /**
  * Computes maximum distance.
  */
diff --git a/fastlib/branches/fastlib-stl/mlpack/mvu/mvu_objectives_impl.h b/fastlib/branches/fastlib-stl/mlpack/mvu/mvu_objectives_impl.h
index 9dd6a2f697..28d412e1a5 100644
--- a/fastlib/branches/fastlib-stl/mlpack/mvu/mvu_objectives_impl.h
+++ b/fastlib/branches/fastlib-stl/mlpack/mvu/mvu_objectives_impl.h
@@ -302,7 +302,7 @@ void MaxFurthestNeighbors::Init(datanode *module, Matrix &data) {
         &nearest_distances_,
         &num_of_nearest_pairs_);
   }
-  sum_of_nearest_distances_= 1.0/std::pow(2,la::Dot(nearest_distances_.size(), 
+  sum_of_nearest_distances_=math::Pow<1,2>(la::Dot(nearest_distances_.size(), 
       &nearest_distances_[0], &nearest_distances_[0]));
   fx_format_result(module_, "num_of_constraints", "%i", num_of_nearest_pairs_);
   eq_lagrange_mult_.Init(num_of_nearest_pairs_);
@@ -371,7 +371,7 @@ void MaxFurthestNeighbors::Init(fx_module *module) {
   }
   num_of_points_++;
   num_of_nearest_pairs_=nearest_neighbor_pairs_.size();
-  sum_of_nearest_distances_= 1.0/std::pow(2,la::Dot(nearest_distances_.size(), 
+  sum_of_nearest_distances_=math::Pow<1,2>(la::Dot(nearest_distances_.size(), 
       &nearest_distances_[0], &nearest_distances_[0]));
 
   fclose(fp);
@@ -481,7 +481,7 @@ void MaxFurthestNeighbors::ComputeFeasibilityError(Matrix &coordinates, double *
                            -nearest_distances_[i];
     *error+=dist_diff*dist_diff;
   }
-  *error= 100 / std::pow(2,*error)/sum_of_nearest_distances_;
+  *error= 100 * math::Pow<1,2>(*error)/sum_of_nearest_distances_;
 }
 
 double MaxFurthestNeighbors::ComputeLagrangian(Matrix &coordinates) {
@@ -563,7 +563,7 @@ bool MaxFurthestNeighbors::IsOptimizationOver(Matrix &coordinates,
 
 bool MaxFurthestNeighbors::IsIntermediateStepOver(Matrix &coordinates, 
       Matrix &gradient, double step) {
-   double norm_gradient=1.0/std::pow(2,la::Dot(gradient.n_elements(), 
+   double norm_gradient=math::Pow<1,2>(la::Dot(gradient.n_elements(), 
                                gradient.ptr(), 
                                gradient.ptr()));
   double feasibility_error;
@@ -632,7 +632,7 @@ void MaxVarianceUtils::EstimateKnns(ArrayList<index_t> &neares_neighbors,
       double probability=0;
       for(index_t j=0; j<k; j++) {
         probability+=exp(-nearest_distances[i*maximum_knns+j]
-            *std::pow(2,nearest_distances[i*maximum_knns+k]-1))/scale_factor;
+            /(2*math::Pow<1,2>(nearest_distances[i*maximum_knns+k])))/scale_factor;
       }
       loocv_score+=log(probability);
       mean_band+=nearest_distances[i*maximum_knns+k];
diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/bounds_aux.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/bounds_aux.h
index 78a4696dae..67a0850bc2 100644
--- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/bounds_aux.h
+++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/bounds_aux.h
@@ -32,9 +32,9 @@ class bounds_aux {
 	furthest_point_in_bound1[d] = bound1_range.hi;
 	v = v2;
       }
-      furthest_dsqd += std::pow(std::fabs(t_pow),v); // v is non-negative
+      furthest_dsqd += math::PowAbs<t_pow, 1>(v); // v is non-negative
     }
-    furthest_dsqd = std::pow(2.0/t_pow, furthest_dsqd);
+    furthest_dsqd = math::Pow<2, t_pow>(furthest_dsqd);
   }
 
   template<int t_pow>
@@ -62,9 +62,9 @@ class bounds_aux {
 	furthest_point_in_bound1[d] = bound1_range.hi;
 	v = v2;
       }
-      furthest_dsqd += std::pow(std::fabs(t_pow), v); // v is non-negative
+      furthest_dsqd += math::PowAbs<t_pow, 1>(v); // v is non-negative
     }
-    furthest_dsqd = std::pow(2.0/t_pow,furthest_dsqd);
+    furthest_dsqd = math::Pow<2, t_pow>(furthest_dsqd);
   }
 
   template<int t_pow, typename TVector>
@@ -90,8 +90,8 @@ class bounds_aux {
     furthest_point_in_bound1.CopyValues(bound1.center());
     la::AddExpert(bound1.radius(), unit_vector, &furthest_point_in_bound1);
     
-    furthest_dsqd = std::pow(2.0/t_pow,
-      la::RawLMetric<t_pow>(bound2_centroid.length(), 
+    furthest_dsqd = math::Pow<2, t_pow>
+      (la::RawLMetric<t_pow>(bound2_centroid.length(), 
 			     furthest_point_in_bound1.ptr(), 
 			     bound2_centroid.ptr()));
   }