diff --git a/fastlib2/mlpack/series_expansion/build.py b/fastlib2/mlpack/series_expansion/build.py
index 9263aa54cb..3e467d9254 100644
--- a/fastlib2/mlpack/series_expansion/build.py
+++ b/fastlib2/mlpack/series_expansion/build.py
@@ -5,9 +5,12 @@ librule(
     headers = ["farfield_expansion.h",
                "farfield_expansion_impl.h",
                "mult_farfield_expansion.h",
+               "mult_farfield_expansion_impl.h",
                "kernel_aux.h",
                "local_expansion.h",
+               "local_expansion_impl.h"
                "mult_local_expansion.h",
+               "mult_local_expansion_impl.h",
                "mult_series_expansion_aux.h",
                "series_expansion_aux.h",
                "mult_series_expansion_aux.h"],
diff --git a/fastlib2/mlpack/series_expansion/local_expansion.h b/fastlib2/mlpack/series_expansion/local_expansion.h
index ea5c6ac1c7..c810015be1 100644
--- a/fastlib2/mlpack/series_expansion/local_expansion.h
+++ b/fastlib2/mlpack/series_expansion/local_expansion.h
@@ -135,435 +135,8 @@ class LocalExpansion {
 
 };
 
-template<typename TKernelAux>
-void LocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data, 
-						  const Vector& weights, 
-						  int begin, int end, 
-						  int order) {
-
-  if(order > order_) {
-    order_ = order;
-  }
-
-  int dim = sea_->get_dimension();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order);
-  
-  // get inverse factorials (precomputed)
-  Vector neg_inv_multiindex_factorials;
-  neg_inv_multiindex_factorials.Alias
-    (sea_->get_neg_inv_multiindex_factorials());
-
-  // declare deritave mapping
-  Matrix derivative_map;
-  derivative_map.Init(dim, order + 1);
-  
-  // some temporary variables
-  Vector arrtmp;
-  arrtmp.Init(total_num_coeffs);
-  Vector x_r_minus_x_Q;
-  x_r_minus_x_Q.Init(dim);
-  
-  // sqrt two times bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-  
-  // for each data point,
-  for(index_t r = begin; r < end; r++) {
-    
-    // calculate x_r - x_Q
-    for(index_t d = 0; d < dim; d++) {
-      x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) / 
-	bandwidth_factor;
-    }
-    
-    // precompute necessary partial derivatives based on coordinate difference
-    ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map);
-    
-    // compute h_{beta}((x_r - x_Q) / sqrt(2h^2))
-    for(index_t j = 0; j < total_num_coeffs; j++) {
-      const ArrayList<int> &mapping = sea_->get_multiindex(j);
-      arrtmp[j] = ka_->ComputePartialDerivative(derivative_map, mapping);
-    }
-
-    for(index_t j = 0; j < total_num_coeffs; j++) {
-      coeffs_[j] += neg_inv_multiindex_factorials[j] * weights[r] * 
-	arrtmp[j];
-    }
-  } // End of looping through each reference point.
-}
-
-template<typename TKernelAux>
-void LocalExpansion<TKernelAux>::PrintDebug(const char *name, 
-					    FILE *stream) const {
-  
-  int dim = sea_->get_dimension();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
-  fprintf(stream, "Local expansion\n");
-  fprintf(stream, "Center: ");
-  
-  for (index_t i = 0; i < center_.length(); i++) {
-    fprintf(stream, "%g ", center_[i]);
-  }
-  fprintf(stream, "\n");
-  
-  fprintf(stream, "f(");
-  for(index_t d = 0; d < dim; d++) {
-    fprintf(stream, "x_q%d", d);
-    if(d < dim - 1)
-      fprintf(stream, ",");
-  }
-  fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
-  
-  for (index_t i = 0; i < total_num_coeffs; i++) {
-    const ArrayList<int> &mapping = sea_->get_multiindex(i);
-    fprintf(stream, "%g", coeffs_[i]);
-    
-    for(index_t d = 0; d < dim; d++) {
-      fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
-    }
-
-    if(i < total_num_coeffs - 1) {
-      fprintf(stream, " + ");
-    }
-  }
-  fprintf(stream, "\n");
-}
-
-template<typename TKernelAux>
-double LocalExpansion<TKernelAux>::EvaluateField(const Matrix& data, 
-						 int row_num) const {
-  
-  // if there are no local expansion here, then return 0
-  if(order_ < 0) {
-    return 0;
-  }
-
-  index_t k, t, tail;
-  
-  // total number of coefficient
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  // number of dimensions
-  int dim = sea_->get_dimension();
-
-  // evaluated sum to be returned
-  double sum = 0;
-  
-  // sqrt two bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-
-  // temporary variable
-  Vector x_Q_to_x_q;
-  x_Q_to_x_q.Init(dim);
-  Vector tmp;
-  tmp.Init(total_num_coeffs);
-  ArrayList<int> heads;
-  heads.Init(dim + 1);
-  
-  // compute (x_q - x_Q) / (sqrt(2h^2))
-  for(index_t i = 0; i < dim; i++) {
-    x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
-  }
-  
-  for(index_t i = 0; i < dim; i++)
-    heads[i] = 0;
-  heads[dim] = MAXINT;
-
-  tmp[0] = 1.0;
-
-  for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
-
-    for(index_t i = 0; i < dim; i++) {
-      int head = heads[i];
-      heads[i] = t;
-
-      for(index_t j = head; j < tail; j++, t++) {
-        tmp[t] = tmp[j] * x_Q_to_x_q[i];
-      }
-    }
-  }
-
-  for(index_t i = 0; i < total_num_coeffs; i++) {
-    sum += coeffs_[i] * tmp[i];
-  }
-
-  return sum;
-}
-
-template<typename TKernelAux>
-double LocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
-  
-  // if there are no local expansion here, then return 0
-  if(order_ < 0) {
-    return 0;
-  }
-
-  index_t k, t, tail;
-  
-  // total number of coefficient
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  // number of dimensions
-  int dim = sea_->get_dimension();
-
-  // evaluated sum to be returned
-  double sum = 0;
-  
-  // sqrt two bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-
-  // temporary variable
-  Vector x_Q_to_x_q;
-  x_Q_to_x_q.Init(dim);
-  Vector tmp;
-  tmp.Init(total_num_coeffs);
-  ArrayList<int> heads;
-  heads.Init(dim + 1);
-  
-  // compute (x_q - x_Q) / (sqrt(2h^2))
-  for(index_t i = 0; i < dim; i++) {
-    x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
-  }
-  
-  for(index_t i = 0; i < dim; i++)
-    heads[i] = 0;
-  heads[dim] = MAXINT;
-
-  tmp[0] = 1.0;
-
-  for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
-
-    for(index_t i = 0; i < dim; i++) {
-      int head = heads[i];
-      heads[i] = t;
-
-      for(index_t j = head; j < tail; j++, t++) {
-        tmp[t] = tmp[j] * x_Q_to_x_q[i];
-      }
-    }
-  }
-
-  for(index_t i = 0; i < total_num_coeffs; i++) {
-    sum += coeffs_[i] * tmp[i];
-  }
-
-  return sum;
-}
-
-template<typename TKernelAux>
-void LocalExpansion<TKernelAux>::Init(const Vector& center,
-				      const TKernelAux &ka) {
-  
-  // copy kernel type, center, and bandwidth squared
-  kernel_ = &(ka.kernel_);
-  center_.Copy(center);
-  order_ = -1;
-  sea_ = &(ka.sea_);
-  ka_ = &ka;
-
-  // initialize coefficient array
-  coeffs_.Init(sea_->get_max_total_num_coeffs());
-  coeffs_.SetZero();
-}
-
-template<typename TKernelAux>
-void LocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
-  
-  // copy kernel type, center, and bandwidth squared
-  kernel_ = &(ka.kernel_);
-  order_ = -1;
-  sea_ = &(ka.sea_);
-  center_.Init(sea_->get_dimension());
-  ka_ = &ka;
-
-  // initialize coefficient array
-  coeffs_.Init(sea_->get_max_total_num_coeffs());
-  coeffs_.SetZero();
-}
-
-template<typename TKernelAux>
-template<typename TBound>
-int LocalExpansion<TKernelAux>::OrderForEvaluating
-(const TBound &far_field_region, 
- const TBound &local_field_region, double min_dist_sqd_regions,
- double max_dist_sqd_regions, double max_error, double *actual_error) const {
-  
-  return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region, 
-				     min_dist_sqd_regions,
-				     max_dist_sqd_regions, max_error, 
-				     actual_error);
-}
-
-template<typename TKernelAux>
-void LocalExpansion<TKernelAux>::TranslateFromFarField
-(const FarFieldExpansion<TKernelAux> &se) {
-
-  Vector pos_arrtmp, neg_arrtmp;
-  Matrix derivative_map;
-  Vector far_center;
-  Vector cent_diff;
-  Vector far_coeffs;
-  int dimension = sea_->get_dimension();
-  int far_order = se.get_order();
-  int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
-  int limit;
-  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
-
-  // get center and coefficients for far field expansion
-  far_center.Alias(*(se.get_center()));
-  far_coeffs.Alias(se.get_coeffs());
-  cent_diff.Init(dimension);
-
-  // if the order of the far field expansion is greater than the
-  // local one we are adding onto, then increase the order.
-  if(far_order > order_) {
-    order_ = far_order;
-  }
-
-  // compute Gaussian derivative
-  limit = 2 * order_ + 1;
-  derivative_map.Init(dimension, limit);
-  pos_arrtmp.Init(total_num_coeffs);
-  neg_arrtmp.Init(total_num_coeffs);
-
-  // compute center difference divided by bw_times_sqrt_two;
-  for(index_t j = 0; j < dimension; j++) {
-    cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
-  }
-
-  // compute required partial derivatives
-  ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
-  ArrayList<int> beta_plus_alpha;
-  beta_plus_alpha.Init(dimension);
-
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-
-    const ArrayList<int> &beta_mapping = sea_->get_multiindex(j);
-    pos_arrtmp[j] = neg_arrtmp[j] = 0;
-
-    for(index_t k = 0; k < total_num_coeffs; k++) {
-
-      const ArrayList<int> &alpha_mapping = sea_->get_multiindex(k);
-      for(index_t d = 0; d < dimension; d++) {
-	beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
-      }
-      double derivative_factor =
-	ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
-      
-      double prod = far_coeffs[k] * derivative_factor;
-
-      if(prod > 0) {
-	pos_arrtmp[j] += prod;
-      }
-      else {
-	neg_arrtmp[j] += prod;
-      }
-    } // end of k-loop
-  } // end of j-loop
-
-  Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-    coeffs_[j] += (pos_arrtmp[j] + neg_arrtmp[j]) * C_k_neg[j];
-  }
-}
-  
-template<typename TKernelAux>
-void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
-  
-  // if there are no local coefficients to translate, return
-  if(order_ < 0) {
-    return;
-  }
-
-  // get the center and the order and the total number of coefficients of 
-  // the expansion we are translating from. Also get coefficients we
-  // are translating
-  Vector new_center;
-  new_center.Alias(*(se.get_center()));
-  int prev_order = se.get_order();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-  const ArrayList < int > *upper_mapping_index = 
-    sea_->get_upper_mapping_index();
-  Vector new_coeffs;
-  new_coeffs.Alias(se.get_coeffs());
-
-  // dimension
-  int dim = sea_->get_dimension();
-
-  // temporary variable
-  ArrayList<int> tmp_storage;
-  tmp_storage.Init(dim);
-
-  // sqrt two times bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-
-  // center difference between the old center and the new one
-  Vector center_diff;
-  center_diff.Init(dim);
-  for(index_t d = 0; d < dim; d++) {
-    center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
-  }
-
-  // set to the new order if the order of the expansion we are translating
-  // from is higher
-  if(prev_order < order_) {
-    se.set_order(order_);
-  }
-
-  // inverse multiindex factorials
-  Vector C_k;
-  C_k.Alias(sea_->get_inv_multiindex_factorials());
-  
-  // do the actual translation
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-
-    const ArrayList<int> &alpha_mapping = sea_->get_multiindex(j);
-    const ArrayList<int> &upper_mappings_for_alpha = upper_mapping_index[j];
-    double pos_coeffs = 0;
-    double neg_coeffs = 0;
-
-    for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
-      
-      if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
-	break;
-      }
-
-      const ArrayList<int> &beta_mapping = 
-	sea_->get_multiindex(upper_mappings_for_alpha[k]);
-      int flag = 0;
-      double diff1 = 1.0;
-
-      for(index_t l = 0; l < dim; l++) {
-	tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
-
-	if(tmp_storage[l] < 0) {
-	  flag = 1;
-	  break;
-	}
-      } // end of looping over dimension
-      
-      if(flag)
-	continue;
-
-      for(index_t l = 0; l < dim; l++) {
-	diff1 *= pow(center_diff[l], tmp_storage[l]);
-      }
-
-      double prod =  coeffs_[upper_mappings_for_alpha[k]] * diff1 *
-	sea_->get_n_multichoose_k_by_pos(upper_mappings_for_alpha[k], j);
-      
-      if(prod > 0) {
-	pos_coeffs += prod;
-      }
-      else {
-	neg_coeffs += prod;
-      }
-
-    } // end of k loop
-
-    new_coeffs[j] += pos_coeffs + neg_coeffs;
-  } // end of j loop
-}
+#define INSIDE_LOCAL_EXPANSION_H
+#include "local_expansion_impl.h"
+#undef INSIDE_LOCAL_EXPANSION_H
 
 #endif
diff --git a/fastlib2/mlpack/series_expansion/local_expansion_impl.h b/fastlib2/mlpack/series_expansion/local_expansion_impl.h
new file mode 100644
index 0000000000..b2e6725758
--- /dev/null
+++ b/fastlib2/mlpack/series_expansion/local_expansion_impl.h
@@ -0,0 +1,441 @@
+#ifndef INSIDE_LOCAL_EXPANSION_H
+#error "This is not a public header file!"
+#endif
+
+#ifndef LOCAL_EXPANSION_IMPL_H
+#define LOCAL_EXPANSION_IMPL_H
+
+
+template<typename TKernelAux>
+void LocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data, 
+						  const Vector& weights, 
+						  int begin, int end, 
+						  int order) {
+
+  if(order > order_) {
+    order_ = order;
+  }
+
+  int dim = sea_->get_dimension();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+  
+  // get inverse factorials (precomputed)
+  Vector neg_inv_multiindex_factorials;
+  neg_inv_multiindex_factorials.Alias
+    (sea_->get_neg_inv_multiindex_factorials());
+
+  // declare deritave mapping
+  Matrix derivative_map;
+  derivative_map.Init(dim, order + 1);
+  
+  // some temporary variables
+  Vector arrtmp;
+  arrtmp.Init(total_num_coeffs);
+  Vector x_r_minus_x_Q;
+  x_r_minus_x_Q.Init(dim);
+  
+  // sqrt two times bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+  
+  // for each data point,
+  for(index_t r = begin; r < end; r++) {
+    
+    // calculate x_r - x_Q
+    for(index_t d = 0; d < dim; d++) {
+      x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) / 
+	bandwidth_factor;
+    }
+    
+    // precompute necessary partial derivatives based on coordinate difference
+    ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map);
+    
+    // compute h_{beta}((x_r - x_Q) / sqrt(2h^2))
+    for(index_t j = 0; j < total_num_coeffs; j++) {
+      const ArrayList<int> &mapping = sea_->get_multiindex(j);
+      arrtmp[j] = ka_->ComputePartialDerivative(derivative_map, mapping);
+    }
+
+    for(index_t j = 0; j < total_num_coeffs; j++) {
+      coeffs_[j] += neg_inv_multiindex_factorials[j] * weights[r] * 
+	arrtmp[j];
+    }
+  } // End of looping through each reference point.
+}
+
+template<typename TKernelAux>
+void LocalExpansion<TKernelAux>::PrintDebug(const char *name, 
+					    FILE *stream) const {
+  
+  int dim = sea_->get_dimension();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
+  fprintf(stream, "Local expansion\n");
+  fprintf(stream, "Center: ");
+  
+  for (index_t i = 0; i < center_.length(); i++) {
+    fprintf(stream, "%g ", center_[i]);
+  }
+  fprintf(stream, "\n");
+  
+  fprintf(stream, "f(");
+  for(index_t d = 0; d < dim; d++) {
+    fprintf(stream, "x_q%d", d);
+    if(d < dim - 1)
+      fprintf(stream, ",");
+  }
+  fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
+  
+  for (index_t i = 0; i < total_num_coeffs; i++) {
+    const ArrayList<int> &mapping = sea_->get_multiindex(i);
+    fprintf(stream, "%g", coeffs_[i]);
+    
+    for(index_t d = 0; d < dim; d++) {
+      fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
+    }
+
+    if(i < total_num_coeffs - 1) {
+      fprintf(stream, " + ");
+    }
+  }
+  fprintf(stream, "\n");
+}
+
+template<typename TKernelAux>
+double LocalExpansion<TKernelAux>::EvaluateField(const Matrix& data, 
+						 int row_num) const {
+  
+  // if there are no local expansion here, then return 0
+  if(order_ < 0) {
+    return 0;
+  }
+
+  index_t k, t, tail;
+  
+  // total number of coefficient
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  // number of dimensions
+  int dim = sea_->get_dimension();
+
+  // evaluated sum to be returned
+  double sum = 0;
+  
+  // sqrt two bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+
+  // temporary variable
+  Vector x_Q_to_x_q;
+  x_Q_to_x_q.Init(dim);
+  Vector tmp;
+  tmp.Init(total_num_coeffs);
+  ArrayList<int> heads;
+  heads.Init(dim + 1);
+  
+  // compute (x_q - x_Q) / (sqrt(2h^2))
+  for(index_t i = 0; i < dim; i++) {
+    x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
+  }
+  
+  for(index_t i = 0; i < dim; i++)
+    heads[i] = 0;
+  heads[dim] = MAXINT;
+
+  tmp[0] = 1.0;
+
+  for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
+
+    for(index_t i = 0; i < dim; i++) {
+      int head = heads[i];
+      heads[i] = t;
+
+      for(index_t j = head; j < tail; j++, t++) {
+        tmp[t] = tmp[j] * x_Q_to_x_q[i];
+      }
+    }
+  }
+
+  for(index_t i = 0; i < total_num_coeffs; i++) {
+    sum += coeffs_[i] * tmp[i];
+  }
+
+  return sum;
+}
+
+template<typename TKernelAux>
+double LocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
+  
+  // if there are no local expansion here, then return 0
+  if(order_ < 0) {
+    return 0;
+  }
+
+  index_t k, t, tail;
+  
+  // total number of coefficient
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  // number of dimensions
+  int dim = sea_->get_dimension();
+
+  // evaluated sum to be returned
+  double sum = 0;
+  
+  // sqrt two bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+
+  // temporary variable
+  Vector x_Q_to_x_q;
+  x_Q_to_x_q.Init(dim);
+  Vector tmp;
+  tmp.Init(total_num_coeffs);
+  ArrayList<int> heads;
+  heads.Init(dim + 1);
+  
+  // compute (x_q - x_Q) / (sqrt(2h^2))
+  for(index_t i = 0; i < dim; i++) {
+    x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
+  }
+  
+  for(index_t i = 0; i < dim; i++)
+    heads[i] = 0;
+  heads[dim] = MAXINT;
+
+  tmp[0] = 1.0;
+
+  for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
+
+    for(index_t i = 0; i < dim; i++) {
+      int head = heads[i];
+      heads[i] = t;
+
+      for(index_t j = head; j < tail; j++, t++) {
+        tmp[t] = tmp[j] * x_Q_to_x_q[i];
+      }
+    }
+  }
+
+  for(index_t i = 0; i < total_num_coeffs; i++) {
+    sum += coeffs_[i] * tmp[i];
+  }
+
+  return sum;
+}
+
+template<typename TKernelAux>
+void LocalExpansion<TKernelAux>::Init(const Vector& center,
+				      const TKernelAux &ka) {
+  
+  // copy kernel type, center, and bandwidth squared
+  kernel_ = &(ka.kernel_);
+  center_.Copy(center);
+  order_ = -1;
+  sea_ = &(ka.sea_);
+  ka_ = &ka;
+
+  // initialize coefficient array
+  coeffs_.Init(sea_->get_max_total_num_coeffs());
+  coeffs_.SetZero();
+}
+
+template<typename TKernelAux>
+void LocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
+  
+  // copy kernel type, center, and bandwidth squared
+  kernel_ = &(ka.kernel_);
+  order_ = -1;
+  sea_ = &(ka.sea_);
+  center_.Init(sea_->get_dimension());
+  ka_ = &ka;
+
+  // initialize coefficient array
+  coeffs_.Init(sea_->get_max_total_num_coeffs());
+  coeffs_.SetZero();
+}
+
+template<typename TKernelAux>
+template<typename TBound>
+int LocalExpansion<TKernelAux>::OrderForEvaluating
+(const TBound &far_field_region, 
+ const TBound &local_field_region, double min_dist_sqd_regions,
+ double max_dist_sqd_regions, double max_error, double *actual_error) const {
+  
+  return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region, 
+				     min_dist_sqd_regions,
+				     max_dist_sqd_regions, max_error, 
+				     actual_error);
+}
+
+template<typename TKernelAux>
+void LocalExpansion<TKernelAux>::TranslateFromFarField
+(const FarFieldExpansion<TKernelAux> &se) {
+
+  Vector pos_arrtmp, neg_arrtmp;
+  Matrix derivative_map;
+  Vector far_center;
+  Vector cent_diff;
+  Vector far_coeffs;
+  int dimension = sea_->get_dimension();
+  int far_order = se.get_order();
+  int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
+  int limit;
+  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
+
+  // get center and coefficients for far field expansion
+  far_center.Alias(*(se.get_center()));
+  far_coeffs.Alias(se.get_coeffs());
+  cent_diff.Init(dimension);
+
+  // if the order of the far field expansion is greater than the
+  // local one we are adding onto, then increase the order.
+  if(far_order > order_) {
+    order_ = far_order;
+  }
+
+  // compute Gaussian derivative
+  limit = 2 * order_ + 1;
+  derivative_map.Init(dimension, limit);
+  pos_arrtmp.Init(total_num_coeffs);
+  neg_arrtmp.Init(total_num_coeffs);
+
+  // compute center difference divided by bw_times_sqrt_two;
+  for(index_t j = 0; j < dimension; j++) {
+    cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
+  }
+
+  // compute required partial derivatives
+  ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
+  ArrayList<int> beta_plus_alpha;
+  beta_plus_alpha.Init(dimension);
+
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+
+    const ArrayList<int> &beta_mapping = sea_->get_multiindex(j);
+    pos_arrtmp[j] = neg_arrtmp[j] = 0;
+
+    for(index_t k = 0; k < total_num_coeffs; k++) {
+
+      const ArrayList<int> &alpha_mapping = sea_->get_multiindex(k);
+      for(index_t d = 0; d < dimension; d++) {
+	beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
+      }
+      double derivative_factor =
+	ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
+      
+      double prod = far_coeffs[k] * derivative_factor;
+
+      if(prod > 0) {
+	pos_arrtmp[j] += prod;
+      }
+      else {
+	neg_arrtmp[j] += prod;
+      }
+    } // end of k-loop
+  } // end of j-loop
+
+  Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+    coeffs_[j] += (pos_arrtmp[j] + neg_arrtmp[j]) * C_k_neg[j];
+  }
+}
+  
+template<typename TKernelAux>
+void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
+  
+  // if there are no local coefficients to translate, return
+  if(order_ < 0) {
+    return;
+  }
+
+  // get the center and the order and the total number of coefficients of 
+  // the expansion we are translating from. Also get coefficients we
+  // are translating
+  Vector new_center;
+  new_center.Alias(*(se.get_center()));
+  int prev_order = se.get_order();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+  const ArrayList < int > *upper_mapping_index = 
+    sea_->get_upper_mapping_index();
+  Vector new_coeffs;
+  new_coeffs.Alias(se.get_coeffs());
+
+  // dimension
+  int dim = sea_->get_dimension();
+
+  // temporary variable
+  ArrayList<int> tmp_storage;
+  tmp_storage.Init(dim);
+
+  // sqrt two times bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+
+  // center difference between the old center and the new one
+  Vector center_diff;
+  center_diff.Init(dim);
+  for(index_t d = 0; d < dim; d++) {
+    center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
+  }
+
+  // set to the new order if the order of the expansion we are translating
+  // from is higher
+  if(prev_order < order_) {
+    se.set_order(order_);
+  }
+
+  // inverse multiindex factorials
+  Vector C_k;
+  C_k.Alias(sea_->get_inv_multiindex_factorials());
+  
+  // do the actual translation
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+
+    const ArrayList<int> &alpha_mapping = sea_->get_multiindex(j);
+    const ArrayList<int> &upper_mappings_for_alpha = upper_mapping_index[j];
+    double pos_coeffs = 0;
+    double neg_coeffs = 0;
+
+    for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
+      
+      if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
+	break;
+      }
+
+      const ArrayList<int> &beta_mapping = 
+	sea_->get_multiindex(upper_mappings_for_alpha[k]);
+      int flag = 0;
+      double diff1 = 1.0;
+
+      for(index_t l = 0; l < dim; l++) {
+	tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
+
+	if(tmp_storage[l] < 0) {
+	  flag = 1;
+	  break;
+	}
+      } // end of looping over dimension
+      
+      if(flag)
+	continue;
+
+      for(index_t l = 0; l < dim; l++) {
+	diff1 *= pow(center_diff[l], tmp_storage[l]);
+      }
+
+      double prod =  coeffs_[upper_mappings_for_alpha[k]] * diff1 *
+	sea_->get_n_multichoose_k_by_pos(upper_mappings_for_alpha[k], j);
+      
+      if(prod > 0) {
+	pos_coeffs += prod;
+      }
+      else {
+	neg_coeffs += prod;
+      }
+
+    } // end of k loop
+
+    new_coeffs[j] += pos_coeffs + neg_coeffs;
+  } // end of j loop
+}
+
+#endif
+
diff --git a/fastlib2/mlpack/series_expansion/mult_farfield_expansion.h b/fastlib2/mlpack/series_expansion/mult_farfield_expansion.h
index 4076817f82..4dd330b332 100644
--- a/fastlib2/mlpack/series_expansion/mult_farfield_expansion.h
+++ b/fastlib2/mlpack/series_expansion/mult_farfield_expansion.h
@@ -207,571 +207,8 @@ class MultFarFieldExpansion {
 
 };
 
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data, 
-							 const Vector& weights,
-							 int begin, int end, 
-							 int order) {
-  
-  int dim = data.n_rows();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order);
-  int max_total_num_coeffs = sea_->get_max_total_num_coeffs();
-  Vector x_r, tmp;
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-
-  // initialize temporary variables
-  x_r.Init(dim);
-  tmp.Init(max_total_num_coeffs);
-  Vector pos_coeffs;
-  Vector neg_coeffs;
-  pos_coeffs.Init(max_total_num_coeffs);
-  pos_coeffs.SetZero();
-  neg_coeffs.Init(max_total_num_coeffs);
-  neg_coeffs.SetZero();
-
-  // set to new order if greater
-  if(order_ < order) {
-    order_ = order;
-  }
-  
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
-
-  // Repeat for each reference point in this reference node.
-  for(index_t r = begin; r < end; r++) {
-
-    // Calculate the coordinate difference between the ref point and the 
-    // centroid.
-    for(index_t i = 0; i < dim; i++) {
-      x_r[i] = (data.get(i, r) - center_[i]) / bandwidth_factor;
-    }
-    
-    tmp.SetZero();
-    tmp[0] = 1.0;
-
-    for(index_t i = 1; i < total_num_coeffs; i++) {
-      
-      int index = traversal_order[i];
-      const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
-
-      // from the direct descendant, recursively compute the multipole moments
-      int direct_ancestor_mapping_pos = 
-	lower_mappings[lower_mappings.size() - 2];
-
-      int position = 0;
-      const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
-      const ArrayList<int> &direct_ancestor_mapping = 
-	sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
-      for(index_t i = 0; i < dim; i++) {
-	if(mapping[i] != direct_ancestor_mapping[i]) {
-	  position = i;
-	  break;
-	}
-      }
-      
-      tmp[index] = tmp[direct_ancestor_mapping_pos] * x_r[position];
-    }
-
-    // Tally up the result in A_k.
-    for(index_t i = 0; i < total_num_coeffs; i++) {
-
-      int index = traversal_order[i];
-      double prod = weights[r] * tmp[index];
-      
-      if(prod > 0) {
-	pos_coeffs[index] += prod;
-      }
-      else {
-	neg_coeffs[index] += prod;
-      }
-    }
-
-  } // End of looping through each reference point
-
-  for(index_t r = 0; r < total_num_coeffs; r++) {
-    int index = traversal_order[r];
-    coeffs_[index] += (pos_coeffs[index] + neg_coeffs[index]) * 
-      sea_->inv_multiindex_factorials_[index];
-  }
-}
-
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::
-ComputeStratifiedLengthSquareDistribution() {
-
-  squared_coeffs_[0] = coeffs_[0] * coeffs_[0];
-  for(index_t i = 1; i < coeffs_.length(); i++) {
-    squared_coeffs_[i] = squared_coeffs_[i - 1] + coeffs_[i] * coeffs_[i];
-  }
-  squared_coeffs_.PrintDebug();
-}
-
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data, 
-						     const Vector& weights, 
-						     int begin, int end, 
-						     int order) {
-
-  // if we already have the order of approximation, then return.
-  if(order_ >= order) {
-    return;
-  }
-
-  // otherwise, recompute from scratch... this could be improved potentially
-  // but I believe it will not squeeze out more performance (as in O(D^p)
-  // expansions).
-  else {
-    order_ = order;
-    
-    coeffs_.SetZero();
-    AccumulateCoeffs(data, weights, begin, end, order);
-  }
-}
-
-template<typename TKernelAux>
-double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Matrix& data, 
-							int row_num, 
-							int order) const {
-  
-  // dimension
-  int dim = sea_->get_dimension();
-
-  // total number of coefficients
-  int total_num_coeffs = sea_->get_total_num_coeffs(order);
-
-  // square root times bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-  
-  // the evaluated sum
-  double pos_multipole_sum = 0;
-  double neg_multipole_sum = 0;
-  double multipole_sum = 0;
-  
-  // computed derivative map
-  Matrix derivative_map;
-  derivative_map.Init(dim, order_ + 1);
-
-  // temporary variable
-  Vector arrtmp;
-  arrtmp.Init(total_num_coeffs);
-
-  // (x_q - x_R) scaled by bandwidth
-  Vector x_q_minus_x_R;
-  x_q_minus_x_R.Init(dim);
-
-  // compute (x_q - x_R) / (sqrt(2h^2))
-  for(index_t d = 0; d < dim; d++) {
-    x_q_minus_x_R[d] = (data.get(d, row_num) - center_[d]) / bandwidth_factor;
-  }
-
-  // compute deriative maps based on coordinate difference.
-  ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
-  
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
-
-  // compute h_{\alpha}((x_q - x_R)/sqrt(2h^2)) ((x_r - x_R)/h)^{\alpha}
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-    
-    int index = traversal_order[j];
-    const ArrayList<int> &mapping = sea_->get_multiindex(index);
-    double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
-    double prod = coeffs_[index] * arrtmp;
-    
-    if(prod > 0) {
-      pos_multipole_sum += prod;
-    }
-    else {
-      neg_multipole_sum += prod;
-    }
-  }
-
-  multipole_sum = pos_multipole_sum + neg_multipole_sum;
-  return multipole_sum;
-}
-
-template<typename TKernelAux>
-double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Vector& x_q, 
-							int order) const {
-  
-  // dimension
-  int dim = sea_->get_dimension();
-
-  // total number of coefficients
-  int total_num_coeffs = sea_->get_total_num_coeffs(order);
-
-  // square root times bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-  
-  // the evaluated sum
-  double pos_multipole_sum = 0;
-  double neg_multipole_sum = 0;
-  double multipole_sum = 0;
-  
-  // computed derivative map
-  Matrix derivative_map;
-  derivative_map.Init(dim, order_ + 1);
-
-  // temporary variable
-  Vector arrtmp;
-  arrtmp.Init(total_num_coeffs);
-
-  // (x_q - x_R) scaled by bandwidth
-  Vector x_q_minus_x_R;
-  x_q_minus_x_R.Init(dim);
-
-  // compute (x_q - x_R) / (sqrt(2h^2))
-  for(index_t d = 0; d < dim; d++) {
-    x_q_minus_x_R[d] = (x_q[d] - center_[d]) / bandwidth_factor;
-  }
-
-  // compute deriative maps based on coordinate difference.
-  ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
-  
-  // get the order of traversal for the given order of approximation
-  ArrayList<int> traversal_order = sea_->traversal_mapping_[order_];
-
-  // compute h_{\alpha}((x_q - x_R)/sqrt(2h^2)) ((x_r - x_R)/h)^{\alpha}
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-    
-    int index = traversal_order[j];
-    ArrayList<int> mapping = sea_->get_multiindex(index);
-    double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
-    double prod = coeffs_[index] * arrtmp;
-    
-    if(prod > 0) {
-      pos_multipole_sum += prod;
-    }
-    else {
-      neg_multipole_sum += prod;
-    }
-  }
-
-  multipole_sum = pos_multipole_sum + neg_multipole_sum;
-  return multipole_sum;
-}
-
-template<typename TKernelAux>
-double MultFarFieldExpansion<TKernelAux>::EvaluateFieldByMonteCarlo
-(const Matrix& data, int row_num, int order, int num_samples) const {
-
-  // I need to implement this...
-  return 0;
-}
-
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::Init(const Vector& center, 
-					     const TKernelAux &ka) {
-  
-  // copy kernel type, center, and bandwidth squared
-  kernel_ = &(ka.kernel_);
-  center_.Copy(center);
-  order_ = -1;
-  sea_ = &(ka.sea_);
-  ka_ = &ka;
-
-  // Initialize coefficient array
-  coeffs_.Init(sea_->get_max_total_num_coeffs());
-  coeffs_.SetZero();
-
-  // Initialize the list of squared coefficients for sampling.
-  squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
-  squared_coeffs_.SetZero();
-}
-
-template<typename TKernelAux>
-  void MultFarFieldExpansion<TKernelAux>::Init(const TKernelAux &ka) {
-  
-  // copy kernel type, center, and bandwidth squared
-  kernel_ = &(ka.kernel_);  
-  order_ = -1;
-  sea_ = &(ka.sea_);
-  center_.Init(sea_->get_dimension());
-  center_.SetZero();
-  ka_ = &ka;
-
-  // Initialize coefficient array.
-  coeffs_.Init(sea_->get_max_total_num_coeffs());
-  coeffs_.SetZero();
-
-  // Initialize the list of squared coefficients for sampling.
-  squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
-  squared_coeffs_.SetZero();
-}
-
-template<typename TKernelAux>
-template<typename TBound>
-int MultFarFieldExpansion<TKernelAux>::OrderForEvaluating
-(const TBound &far_field_region, 
- const TBound &local_field_region, double min_dist_sqd_regions,
- double max_dist_sqd_regions, double max_error, double *actual_error) const {
-  
-  return ka_->OrderForEvaluatingFarField(far_field_region,
-					 local_field_region,
-					 min_dist_sqd_regions, 
-					 max_dist_sqd_regions, max_error,
-					 actual_error);
-}
-
-template<typename TKernelAux>
-template<typename TBound>
-int MultFarFieldExpansion<TKernelAux>::OrderForEvaluatingByMonteCarlo
-(const TBound &far_field_region,
- const TBound &local_field_region, double min_dist_sqd_regions,
- double max_dist_sqd_regions, double max_error, double *actual_error,
- int *num_samples) const {
-
-  int order = ka_->OrderForEvaluatingFarField(far_field_region,
-					      local_field_region,
-					      min_dist_sqd_regions,
-					      max_dist_sqd_regions, max_error,
-					      actual_error);
-  *num_samples = std::max((int) sqrt(coeffs_.length()), order + 1);
-
-  return order;
-}
-
-template<typename TKernelAux>
-template<typename TBound>
-int MultFarFieldExpansion<TKernelAux>::
-OrderForConvertingToLocal(const TBound &far_field_region,
-			  const TBound &local_field_region, 
-			  double min_dist_sqd_regions, 
-			  double max_dist_sqd_regions,
-			  double max_error, 
-			  double *actual_error) const {
-  
-  return ka_->OrderForConvertingFromFarFieldToLocal
-    (far_field_region, local_field_region, min_dist_sqd_regions,
-     max_dist_sqd_regions, max_error, actual_error);
-}
-
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::PrintDebug
-(const char *name, FILE *stream) const {
-    
-  int dim = sea_->get_dimension();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
-  fprintf(stream, "Far field expansion\n");
-  fprintf(stream, "Center: ");
-  
-  for (index_t i = 0; i < center_.length(); i++) {
-    fprintf(stream, "%g ", center_[i]);
-  }
-  fprintf(stream, "\n");
-  
-  fprintf(stream, "f(");
-  for(index_t d = 0; d < dim; d++) {
-    fprintf(stream, "x_q%d", d);
-    if(d < dim - 1)
-      fprintf(stream, ",");
-  }
-  fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
-  
-  for (index_t i = 0; i < total_num_coeffs; i++) {
-    const ArrayList<int> &mapping = sea_->get_multiindex(i);
-    fprintf(stream, "%g ", coeffs_[i]);
-    
-    fprintf(stream, "(-1)^(");
-    for(index_t d = 0; d < dim; d++) {
-      fprintf(stream, "%d", mapping[d]);
-      if(d < dim - 1)
-	fprintf(stream, " + ");
-    }
-    fprintf(stream, ") D^((");
-    for(index_t d = 0; d < dim; d++) {
-      fprintf(stream, "%d", mapping[d]);
-      
-      if(d < dim - 1)
-	fprintf(stream, ",");
-    }
-    fprintf(stream, ")) f(x_q - x_R)");
-    if(i < total_num_coeffs - 1) {
-      fprintf(stream, " + ");
-    }
-  }
-  fprintf(stream, "\n");
-}
-
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::TranslateFromFarField
-(const MultFarFieldExpansion &se) {
-  
-  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
-  int dim = sea_->get_dimension();
-  int order = se.get_order();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order);
-  Vector prev_coeffs;
-  Vector prev_center;
-  const ArrayList < int > *multiindex_mapping = sea_->get_multiindex_mapping();
-  const ArrayList < int > *lower_mapping_index = 
-    sea_->get_lower_mapping_index();
-
-  ArrayList <int> tmp_storage;
-  Vector center_diff;
-  Vector inv_multiindex_factorials;
-
-  center_diff.Init(dim);
-
-  // retrieve coefficients to be translated and helper mappings
-  prev_coeffs.Alias(se.get_coeffs());
-  prev_center.Alias(*(se.get_center()));
-  tmp_storage.Init(sea_->get_dimension());
-  inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
-
-  // no coefficients can be translated
-  if(order == -1) {
-    return;
-  }
-  else {
-    order_ = order;
-  }
-  
-  // compute center difference
-  for(index_t j = 0; j < dim; j++) {
-    center_diff[j] = prev_center[j] - center_[j];
-  }
-
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
-
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-   
-    int index = traversal_order[j];
-    const ArrayList <int> &gamma_mapping = multiindex_mapping[index];
-    const ArrayList <int> &lower_mappings_for_gamma = 
-      lower_mapping_index[index];
-    double pos_coeff = 0;
-    double neg_coeff = 0;
-
-    for(index_t k = 0; k < lower_mappings_for_gamma.size(); k++) {
-
-      const ArrayList <int> &inner_mapping = 
-	multiindex_mapping[lower_mappings_for_gamma[k]];
-
-      int flag = 0;
-      double diff1;
-      
-      // compute gamma minus alpha
-      for(index_t l = 0; l < dim; l++) {
-	tmp_storage[l] = gamma_mapping[l] - inner_mapping[l];
-
-	if(tmp_storage[l] < 0) {
-	  flag = 1;
-	  break;
-	}
-      }
-      
-      if(flag) {
-	continue;
-      }
-      
-      diff1 = 1.0;
-      
-      for(index_t l = 0; l < dim; l++) {
-	diff1 *= pow(center_diff[l] / bandwidth_factor, tmp_storage[l]);
-      }
-
-      double prod = prev_coeffs[lower_mappings_for_gamma[k]] * diff1 * 
-	inv_multiindex_factorials
-	[sea_->ComputeMultiindexPosition(tmp_storage)];
-      
-      if(prod > 0) {
-	pos_coeff += prod;
-      }
-      else {
-	neg_coeff += prod;
-      }
-
-    } // end of k-loop
-    
-    coeffs_[j] += pos_coeff + neg_coeff;
-
-  } // end of j-loop
-}
-
-template<typename TKernelAux>
-void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
-(MultLocalExpansion<TKernelAux> &se, int truncation_order) {
-  
-  Vector pos_arrtmp, neg_arrtmp;
-  Matrix derivative_map;
-  Vector local_center;
-  Vector cent_diff;
-  Vector local_coeffs;
-  int local_order = se.get_order();
-  int dimension = sea_->get_dimension();
-  int total_num_coeffs = sea_->get_total_num_coeffs(truncation_order);
-  int limit;
-  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
-
-  // get center and coefficients for local expansion
-  local_center.Alias(*(se.get_center()));
-  local_coeffs.Alias(se.get_coeffs());
-  cent_diff.Init(dimension);
-
-  // if the order of the far field expansion is greater than the
-  // local one we are adding onto, then increase the order.
-  if(local_order < truncation_order) {
-    se.set_order(truncation_order);
-  }
-
-  // compute Gaussian derivative
-  limit = 2 * truncation_order + 1;
-  derivative_map.Init(dimension, limit);
-  pos_arrtmp.Init(sea_->get_max_total_num_coeffs());
-  neg_arrtmp.Init(sea_->get_max_total_num_coeffs());
-
-  // compute center difference divided by bw_times_sqrt_two;
-  for(index_t j = 0; j < dimension; j++) {
-    cent_diff[j] = (local_center[j] - center_[j]) / bandwidth_factor;
-  }
-
-  // compute required partial derivatives
-  ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
-  ArrayList<int> beta_plus_alpha;
-  beta_plus_alpha.Init(dimension);
-
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = 
-    sea_->traversal_mapping_[truncation_order];
-
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-
-    int index = traversal_order[j];
-    const ArrayList<int> &beta_mapping = sea_->get_multiindex(index);
-    pos_arrtmp[index] = neg_arrtmp[index] = 0;
-
-    for(index_t k = 0; k < total_num_coeffs; k++) {
-
-      int index_k = traversal_order[k];
-
-      const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_k);
-      for(index_t d = 0; d < dimension; d++) {
-	beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
-      }
-      double derivative_factor =
-	ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
-      
-      double prod = coeffs_[index_k] * derivative_factor;
-
-      if(prod > 0) {
-	pos_arrtmp[index] += prod;
-      }
-      else {
-	neg_arrtmp[index] += prod;
-      }
-    } // end of k-loop
-  } // end of j-loop
-
-  Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-    int index = traversal_order[j];
-    local_coeffs[index] += (pos_arrtmp[index] + neg_arrtmp[index]) * 
-      C_k_neg[index];
-  }
-}
+#define INSIDE_MULT_FARFIELD_EXPANSION_H
+#include "mult_farfield_expansion_impl.h"
+#undef INSIDE_MULT_FARFIELD_EXPANSION_H
 
 #endif
diff --git a/fastlib2/mlpack/series_expansion/mult_farfield_expansion_impl.h b/fastlib2/mlpack/series_expansion/mult_farfield_expansion_impl.h
new file mode 100644
index 0000000000..d7007a86d3
--- /dev/null
+++ b/fastlib2/mlpack/series_expansion/mult_farfield_expansion_impl.h
@@ -0,0 +1,575 @@
+#ifndef INSIDE_MULT_FARFIELD_EXPANSION_H
+#error "This file is not a public header file!"
+#endif
+
+#ifndef MULT_FARFIELD_EXPANSION_IMPL_H
+#define MULT_FARFIELD_EXPANSION_IMPL_H
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data, 
+							 const Vector& weights,
+							 int begin, int end, 
+							 int order) {
+  
+  int dim = data.n_rows();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+  int max_total_num_coeffs = sea_->get_max_total_num_coeffs();
+  Vector x_r, tmp;
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+
+  // initialize temporary variables
+  x_r.Init(dim);
+  tmp.Init(max_total_num_coeffs);
+  Vector pos_coeffs;
+  Vector neg_coeffs;
+  pos_coeffs.Init(max_total_num_coeffs);
+  pos_coeffs.SetZero();
+  neg_coeffs.Init(max_total_num_coeffs);
+  neg_coeffs.SetZero();
+
+  // set to new order if greater
+  if(order_ < order) {
+    order_ = order;
+  }
+  
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
+
+  // Repeat for each reference point in this reference node.
+  for(index_t r = begin; r < end; r++) {
+
+    // Calculate the coordinate difference between the ref point and the 
+    // centroid.
+    for(index_t i = 0; i < dim; i++) {
+      x_r[i] = (data.get(i, r) - center_[i]) / bandwidth_factor;
+    }
+    
+    tmp.SetZero();
+    tmp[0] = 1.0;
+
+    for(index_t i = 1; i < total_num_coeffs; i++) {
+      
+      int index = traversal_order[i];
+      const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
+
+      // from the direct descendant, recursively compute the multipole moments
+      int direct_ancestor_mapping_pos = 
+	lower_mappings[lower_mappings.size() - 2];
+
+      int position = 0;
+      const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
+      const ArrayList<int> &direct_ancestor_mapping = 
+	sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
+      for(index_t i = 0; i < dim; i++) {
+	if(mapping[i] != direct_ancestor_mapping[i]) {
+	  position = i;
+	  break;
+	}
+      }
+      
+      tmp[index] = tmp[direct_ancestor_mapping_pos] * x_r[position];
+    }
+
+    // Tally up the result in A_k.
+    for(index_t i = 0; i < total_num_coeffs; i++) {
+
+      int index = traversal_order[i];
+      double prod = weights[r] * tmp[index];
+      
+      if(prod > 0) {
+	pos_coeffs[index] += prod;
+      }
+      else {
+	neg_coeffs[index] += prod;
+      }
+    }
+
+  } // End of looping through each reference point
+
+  for(index_t r = 0; r < total_num_coeffs; r++) {
+    int index = traversal_order[r];
+    coeffs_[index] += (pos_coeffs[index] + neg_coeffs[index]) * 
+      sea_->inv_multiindex_factorials_[index];
+  }
+}
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::
+ComputeStratifiedLengthSquareDistribution() {
+
+  squared_coeffs_[0] = coeffs_[0] * coeffs_[0];
+  for(index_t i = 1; i < coeffs_.length(); i++) {
+    squared_coeffs_[i] = squared_coeffs_[i - 1] + coeffs_[i] * coeffs_[i];
+  }
+  squared_coeffs_.PrintDebug();
+}
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data, 
+						     const Vector& weights, 
+						     int begin, int end, 
+						     int order) {
+
+  // if we already have the order of approximation, then return.
+  if(order_ >= order) {
+    return;
+  }
+
+  // otherwise, recompute from scratch... this could be improved potentially
+  // but I believe it will not squeeze out more performance (as in O(D^p)
+  // expansions).
+  else {
+    order_ = order;
+    
+    coeffs_.SetZero();
+    AccumulateCoeffs(data, weights, begin, end, order);
+  }
+}
+
+template<typename TKernelAux>
+double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Matrix& data, 
+							int row_num, 
+							int order) const {
+  
+  // dimension
+  int dim = sea_->get_dimension();
+
+  // total number of coefficients
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+
+  // square root times bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+  
+  // the evaluated sum
+  double pos_multipole_sum = 0;
+  double neg_multipole_sum = 0;
+  double multipole_sum = 0;
+  
+  // computed derivative map
+  Matrix derivative_map;
+  derivative_map.Init(dim, order_ + 1);
+
+  // temporary variable
+  Vector arrtmp;
+  arrtmp.Init(total_num_coeffs);
+
+  // (x_q - x_R) scaled by bandwidth
+  Vector x_q_minus_x_R;
+  x_q_minus_x_R.Init(dim);
+
+  // compute (x_q - x_R) / (sqrt(2h^2))
+  for(index_t d = 0; d < dim; d++) {
+    x_q_minus_x_R[d] = (data.get(d, row_num) - center_[d]) / bandwidth_factor;
+  }
+
+  // compute deriative maps based on coordinate difference.
+  ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
+  
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
+
+  // compute h_{\alpha}((x_q - x_R)/sqrt(2h^2)) ((x_r - x_R)/h)^{\alpha}
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+    
+    int index = traversal_order[j];
+    const ArrayList<int> &mapping = sea_->get_multiindex(index);
+    double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
+    double prod = coeffs_[index] * arrtmp;
+    
+    if(prod > 0) {
+      pos_multipole_sum += prod;
+    }
+    else {
+      neg_multipole_sum += prod;
+    }
+  }
+
+  multipole_sum = pos_multipole_sum + neg_multipole_sum;
+  return multipole_sum;
+}
+
+template<typename TKernelAux>
+double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Vector& x_q, 
+							int order) const {
+  
+  // dimension
+  int dim = sea_->get_dimension();
+
+  // total number of coefficients
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+
+  // square root times bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+  
+  // the evaluated sum
+  double pos_multipole_sum = 0;
+  double neg_multipole_sum = 0;
+  double multipole_sum = 0;
+  
+  // computed derivative map
+  Matrix derivative_map;
+  derivative_map.Init(dim, order_ + 1);
+
+  // temporary variable
+  Vector arrtmp;
+  arrtmp.Init(total_num_coeffs);
+
+  // (x_q - x_R) scaled by bandwidth
+  Vector x_q_minus_x_R;
+  x_q_minus_x_R.Init(dim);
+
+  // compute (x_q - x_R) / (sqrt(2h^2))
+  for(index_t d = 0; d < dim; d++) {
+    x_q_minus_x_R[d] = (x_q[d] - center_[d]) / bandwidth_factor;
+  }
+
+  // compute deriative maps based on coordinate difference.
+  ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
+  
+  // get the order of traversal for the given order of approximation
+  ArrayList<int> traversal_order = sea_->traversal_mapping_[order_];
+
+  // compute h_{\alpha}((x_q - x_R)/sqrt(2h^2)) ((x_r - x_R)/h)^{\alpha}
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+    
+    int index = traversal_order[j];
+    ArrayList<int> mapping = sea_->get_multiindex(index);
+    double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
+    double prod = coeffs_[index] * arrtmp;
+    
+    if(prod > 0) {
+      pos_multipole_sum += prod;
+    }
+    else {
+      neg_multipole_sum += prod;
+    }
+  }
+
+  multipole_sum = pos_multipole_sum + neg_multipole_sum;
+  return multipole_sum;
+}
+
+template<typename TKernelAux>
+double MultFarFieldExpansion<TKernelAux>::EvaluateFieldByMonteCarlo
+(const Matrix& data, int row_num, int order, int num_samples) const {
+
+  // I need to implement this...
+  return 0;
+}
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::Init(const Vector& center, 
+					     const TKernelAux &ka) {
+  
+  // copy kernel type, center, and bandwidth squared
+  kernel_ = &(ka.kernel_);
+  center_.Copy(center);
+  order_ = -1;
+  sea_ = &(ka.sea_);
+  ka_ = &ka;
+
+  // Initialize coefficient array
+  coeffs_.Init(sea_->get_max_total_num_coeffs());
+  coeffs_.SetZero();
+
+  // Initialize the list of squared coefficients for sampling.
+  squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
+  squared_coeffs_.SetZero();
+}
+
+template<typename TKernelAux>
+  void MultFarFieldExpansion<TKernelAux>::Init(const TKernelAux &ka) {
+  
+  // copy kernel type, center, and bandwidth squared
+  kernel_ = &(ka.kernel_);  
+  order_ = -1;
+  sea_ = &(ka.sea_);
+  center_.Init(sea_->get_dimension());
+  center_.SetZero();
+  ka_ = &ka;
+
+  // Initialize coefficient array.
+  coeffs_.Init(sea_->get_max_total_num_coeffs());
+  coeffs_.SetZero();
+
+  // Initialize the list of squared coefficients for sampling.
+  squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
+  squared_coeffs_.SetZero();
+}
+
+template<typename TKernelAux>
+template<typename TBound>
+int MultFarFieldExpansion<TKernelAux>::OrderForEvaluating
+(const TBound &far_field_region, 
+ const TBound &local_field_region, double min_dist_sqd_regions,
+ double max_dist_sqd_regions, double max_error, double *actual_error) const {
+  
+  return ka_->OrderForEvaluatingFarField(far_field_region,
+					 local_field_region,
+					 min_dist_sqd_regions, 
+					 max_dist_sqd_regions, max_error,
+					 actual_error);
+}
+
+template<typename TKernelAux>
+template<typename TBound>
+int MultFarFieldExpansion<TKernelAux>::OrderForEvaluatingByMonteCarlo
+(const TBound &far_field_region,
+ const TBound &local_field_region, double min_dist_sqd_regions,
+ double max_dist_sqd_regions, double max_error, double *actual_error,
+ int *num_samples) const {
+
+  int order = ka_->OrderForEvaluatingFarField(far_field_region,
+					      local_field_region,
+					      min_dist_sqd_regions,
+					      max_dist_sqd_regions, max_error,
+					      actual_error);
+  *num_samples = std::max((int) sqrt(coeffs_.length()), order + 1);
+
+  return order;
+}
+
+template<typename TKernelAux>
+template<typename TBound>
+int MultFarFieldExpansion<TKernelAux>::
+OrderForConvertingToLocal(const TBound &far_field_region,
+			  const TBound &local_field_region, 
+			  double min_dist_sqd_regions, 
+			  double max_dist_sqd_regions,
+			  double max_error, 
+			  double *actual_error) const {
+  
+  return ka_->OrderForConvertingFromFarFieldToLocal
+    (far_field_region, local_field_region, min_dist_sqd_regions,
+     max_dist_sqd_regions, max_error, actual_error);
+}
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::PrintDebug
+(const char *name, FILE *stream) const {
+    
+  int dim = sea_->get_dimension();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
+  fprintf(stream, "Far field expansion\n");
+  fprintf(stream, "Center: ");
+  
+  for (index_t i = 0; i < center_.length(); i++) {
+    fprintf(stream, "%g ", center_[i]);
+  }
+  fprintf(stream, "\n");
+  
+  fprintf(stream, "f(");
+  for(index_t d = 0; d < dim; d++) {
+    fprintf(stream, "x_q%d", d);
+    if(d < dim - 1)
+      fprintf(stream, ",");
+  }
+  fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
+  
+  for (index_t i = 0; i < total_num_coeffs; i++) {
+    const ArrayList<int> &mapping = sea_->get_multiindex(i);
+    fprintf(stream, "%g ", coeffs_[i]);
+    
+    fprintf(stream, "(-1)^(");
+    for(index_t d = 0; d < dim; d++) {
+      fprintf(stream, "%d", mapping[d]);
+      if(d < dim - 1)
+	fprintf(stream, " + ");
+    }
+    fprintf(stream, ") D^((");
+    for(index_t d = 0; d < dim; d++) {
+      fprintf(stream, "%d", mapping[d]);
+      
+      if(d < dim - 1)
+	fprintf(stream, ",");
+    }
+    fprintf(stream, ")) f(x_q - x_R)");
+    if(i < total_num_coeffs - 1) {
+      fprintf(stream, " + ");
+    }
+  }
+  fprintf(stream, "\n");
+}
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::TranslateFromFarField
+(const MultFarFieldExpansion &se) {
+  
+  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
+  int dim = sea_->get_dimension();
+  int order = se.get_order();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+  Vector prev_coeffs;
+  Vector prev_center;
+  const ArrayList < int > *multiindex_mapping = sea_->get_multiindex_mapping();
+  const ArrayList < int > *lower_mapping_index = 
+    sea_->get_lower_mapping_index();
+
+  ArrayList <int> tmp_storage;
+  Vector center_diff;
+  Vector inv_multiindex_factorials;
+
+  center_diff.Init(dim);
+
+  // retrieve coefficients to be translated and helper mappings
+  prev_coeffs.Alias(se.get_coeffs());
+  prev_center.Alias(*(se.get_center()));
+  tmp_storage.Init(sea_->get_dimension());
+  inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
+
+  // no coefficients can be translated
+  if(order == -1) {
+    return;
+  }
+  else {
+    order_ = order;
+  }
+  
+  // compute center difference
+  for(index_t j = 0; j < dim; j++) {
+    center_diff[j] = prev_center[j] - center_[j];
+  }
+
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
+
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+   
+    int index = traversal_order[j];
+    const ArrayList <int> &gamma_mapping = multiindex_mapping[index];
+    const ArrayList <int> &lower_mappings_for_gamma = 
+      lower_mapping_index[index];
+    double pos_coeff = 0;
+    double neg_coeff = 0;
+
+    for(index_t k = 0; k < lower_mappings_for_gamma.size(); k++) {
+
+      const ArrayList <int> &inner_mapping = 
+	multiindex_mapping[lower_mappings_for_gamma[k]];
+
+      int flag = 0;
+      double diff1;
+      
+      // compute gamma minus alpha
+      for(index_t l = 0; l < dim; l++) {
+	tmp_storage[l] = gamma_mapping[l] - inner_mapping[l];
+
+	if(tmp_storage[l] < 0) {
+	  flag = 1;
+	  break;
+	}
+      }
+      
+      if(flag) {
+	continue;
+      }
+      
+      diff1 = 1.0;
+      
+      for(index_t l = 0; l < dim; l++) {
+	diff1 *= pow(center_diff[l] / bandwidth_factor, tmp_storage[l]);
+      }
+
+      double prod = prev_coeffs[lower_mappings_for_gamma[k]] * diff1 * 
+	inv_multiindex_factorials
+	[sea_->ComputeMultiindexPosition(tmp_storage)];
+      
+      if(prod > 0) {
+	pos_coeff += prod;
+      }
+      else {
+	neg_coeff += prod;
+      }
+
+    } // end of k-loop
+    
+    coeffs_[j] += pos_coeff + neg_coeff;
+
+  } // end of j-loop
+}
+
+template<typename TKernelAux>
+void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
+(MultLocalExpansion<TKernelAux> &se, int truncation_order) {
+  
+  Vector pos_arrtmp, neg_arrtmp;
+  Matrix derivative_map;
+  Vector local_center;
+  Vector cent_diff;
+  Vector local_coeffs;
+  int local_order = se.get_order();
+  int dimension = sea_->get_dimension();
+  int total_num_coeffs = sea_->get_total_num_coeffs(truncation_order);
+  int limit;
+  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
+
+  // get center and coefficients for local expansion
+  local_center.Alias(*(se.get_center()));
+  local_coeffs.Alias(se.get_coeffs());
+  cent_diff.Init(dimension);
+
+  // if the order of the far field expansion is greater than the
+  // local one we are adding onto, then increase the order.
+  if(local_order < truncation_order) {
+    se.set_order(truncation_order);
+  }
+
+  // compute Gaussian derivative
+  limit = 2 * truncation_order + 1;
+  derivative_map.Init(dimension, limit);
+  pos_arrtmp.Init(sea_->get_max_total_num_coeffs());
+  neg_arrtmp.Init(sea_->get_max_total_num_coeffs());
+
+  // compute center difference divided by bw_times_sqrt_two;
+  for(index_t j = 0; j < dimension; j++) {
+    cent_diff[j] = (local_center[j] - center_[j]) / bandwidth_factor;
+  }
+
+  // compute required partial derivatives
+  ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
+  ArrayList<int> beta_plus_alpha;
+  beta_plus_alpha.Init(dimension);
+
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = 
+    sea_->traversal_mapping_[truncation_order];
+
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+
+    int index = traversal_order[j];
+    const ArrayList<int> &beta_mapping = sea_->get_multiindex(index);
+    pos_arrtmp[index] = neg_arrtmp[index] = 0;
+
+    for(index_t k = 0; k < total_num_coeffs; k++) {
+
+      int index_k = traversal_order[k];
+
+      const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_k);
+      for(index_t d = 0; d < dimension; d++) {
+	beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
+      }
+      double derivative_factor =
+	ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
+      
+      double prod = coeffs_[index_k] * derivative_factor;
+
+      if(prod > 0) {
+	pos_arrtmp[index] += prod;
+      }
+      else {
+	neg_arrtmp[index] += prod;
+      }
+    } // end of k-loop
+  } // end of j-loop
+
+  Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+    int index = traversal_order[j];
+    local_coeffs[index] += (pos_arrtmp[index] + neg_arrtmp[index]) * 
+      C_k_neg[index];
+  }
+}
+
+#endif
diff --git a/fastlib2/mlpack/series_expansion/mult_local_expansion.h b/fastlib2/mlpack/series_expansion/mult_local_expansion.h
index 9f9d4245d7..78c5be3097 100644
--- a/fastlib2/mlpack/series_expansion/mult_local_expansion.h
+++ b/fastlib2/mlpack/series_expansion/mult_local_expansion.h
@@ -135,471 +135,8 @@ class MultLocalExpansion {
 
 };
 
-template<typename TKernelAux>
-void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data, 
-						      const Vector& weights, 
-						      int begin, int end, 
-						      int order) {
-
-  if(order > order_) {
-    order_ = order;
-  }
-
-  int dim = sea_->get_dimension();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order);
-  
-  // get inverse factorials (precomputed)
-  Vector neg_inv_multiindex_factorials;
-  neg_inv_multiindex_factorials.Alias
-    (sea_->get_neg_inv_multiindex_factorials());
-
-  // declare deritave mapping
-  Matrix derivative_map;
-  derivative_map.Init(dim, order + 1);
-  
-  // some temporary variables
-  Vector x_r_minus_x_Q;
-  x_r_minus_x_Q.Init(dim);
-  
-  // sqrt two times bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-  
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
-
-  // for each data point,
-  for(index_t r = begin; r < end; r++) {
-    
-    // calculate x_r - x_Q
-    for(index_t d = 0; d < dim; d++) {
-      x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) / 
-	bandwidth_factor;
-    }
-    
-    // precompute necessary partial derivatives based on coordinate difference
-    ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map);
-    
-    // compute h_{beta}((x_r - x_Q) / sqrt(2h^2))
-    for(index_t j = 0; j < total_num_coeffs; j++) {
-      int index = traversal_order[j];
-      const ArrayList<int> &mapping = sea_->get_multiindex(index);
-      double partial_derivative = 
-	ka_->ComputePartialDerivative(derivative_map, mapping);
-      coeffs_[index] += neg_inv_multiindex_factorials[index] * weights[r] * 
-	partial_derivative;
-    }
-  } // End of looping through each reference point.
-}
-
-template<typename TKernelAux>
-void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name, 
-						FILE *stream) const {
-  
-    
-  int dim = sea_->get_dimension();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
-  fprintf(stream, "Local expansion\n");
-  fprintf(stream, "Center: ");
-  
-  for (index_t i = 0; i < center_.length(); i++) {
-    fprintf(stream, "%g ", center_[i]);
-  }
-  fprintf(stream, "\n");
-  
-  fprintf(stream, "f(");
-  for(index_t d = 0; d < dim; d++) {
-    fprintf(stream, "x_q%d", d);
-    if(d < dim - 1)
-      fprintf(stream, ",");
-  }
-  fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
-  
-  for (index_t i = 0; i < total_num_coeffs; i++) {
-    ArrayList<int> mapping = sea_->get_multiindex(i);
-    fprintf(stream, "%g", coeffs_[i]);
-    
-    for(index_t d = 0; d < dim; d++) {
-      fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
-    }
-
-    if(i < total_num_coeffs - 1) {
-      fprintf(stream, " + ");
-    }
-  }
-  fprintf(stream, "\n");
-}
-
-template<typename TKernelAux>
-double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data, 
-						     int row_num) const {
-    
-  // if there are no local coefficients, then return 0
-  if(order_ < 0) {
-    return 0;
-  }
-
-  // total number of coefficient
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  // number of dimensions
-  int dim = sea_->get_dimension();
-
-  // evaluated sum to be returned
-  double sum = 0;
-  
-  // sqrt two bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-
-  // temporary variable
-  Vector x_Q_to_x_q;
-  x_Q_to_x_q.Init(dim);
-  Vector tmp;
-  tmp.Init(sea_->get_max_total_num_coeffs());
-  ArrayList<int> heads;
-  heads.Init(dim + 1);
-  
-  // compute (x_q - x_Q) / (sqrt(2h^2))
-  for(index_t i = 0; i < dim; i++) {
-    x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
-  }
-  
-  for(index_t i = 0; i < dim; i++)
-    heads[i] = 0;
-  heads[dim] = MAXINT;
-
-  tmp[0] = 1.0;
-
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
-
-  for(index_t i = 1; i < total_num_coeffs; i++) {
-    
-    int index = traversal_order[i];
-    const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
-    
-    // from the direct descendant, recursively compute the multipole moments
-    int direct_ancestor_mapping_pos = 
-      lower_mappings[lower_mappings.size() - 2];
-    int position = 0;
-    const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
-    const ArrayList<int> &direct_ancestor_mapping = 
-      sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
-    for(index_t i = 0; i < dim; i++) {
-      if(mapping[i] != direct_ancestor_mapping[i]) {
-	position = i;
-	break;
-      }
-    }
-    tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
-  }
-
-  for(index_t i = 0; i < total_num_coeffs; i++) {
-    int index = traversal_order[i];
-    sum += coeffs_[index] * tmp[index];
-  }
-
-  return sum;
-}
-
-template<typename TKernelAux>
-double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
-    
-  // if there are no local coefficients, then return 0
-  if(order_ < 0) {
-    return 0;
-  }
-
-  // total number of coefficient
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-
-  // number of dimensions
-  int dim = sea_->get_dimension();
-
-  // evaluated sum to be returned
-  double sum = 0;
-  
-  // sqrt two bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_.bandwidth_sq());
-
-  // temporary variable
-  Vector x_Q_to_x_q;
-  x_Q_to_x_q.Init(dim);
-  Vector tmp;
-  tmp.Init(sea_->get_max_total_num_coeffs());
-  ArrayList<int> heads;
-  heads.Init(dim + 1);
-  
-  // compute (x_q - x_Q) / (sqrt(2h^2))
-  for(index_t i = 0; i < dim; i++) {
-    x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
-  }
-  
-  for(index_t i = 0; i < dim; i++)
-    heads[i] = 0;
-  heads[dim] = MAXINT;
-
-  tmp[0] = 1.0;
-
-  // get the order of traversal for the given order of approximation
-  ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
-
-  for(index_t i = 1; i < total_num_coeffs; i++) {
-    
-    int index = traversal_order[i];
-    ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
-    
-    // from the direct descendant, recursively compute the multipole moments
-    int direct_ancestor_mapping_pos = 
-      lower_mappings[lower_mappings.size() - 2];
-    int position = 0;
-    const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
-    const ArrayList<int> &direct_ancestor_mapping = 
-      sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
-    for(index_t i = 0; i < dim; i++) {
-      if(mapping[i] != direct_ancestor_mapping[i]) {
-	position = i;
-	break;
-      }
-    }
-    tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
-  }
-
-  for(index_t i = 0; i < total_num_coeffs; i++) {
-    int index = traversal_order[i];
-    sum += coeffs_[index] * tmp[index];
-  }
-
-  return sum;
-}
-
-template<typename TKernelAux>
-void MultLocalExpansion<TKernelAux>::Init(const Vector& center, 
-					  const TKernelAux &ka) {
-  
-  // copy kernel type, center, and bandwidth squared
-  kernel_ = &(ka.kernel_);
-  center_.Copy(center);
-  order_ = -1;
-  sea_ = &(ka.sea_);
-  ka_ = &ka;
-
-  // initialize coefficient array
-  coeffs_.Init(sea_->get_max_total_num_coeffs());
-  coeffs_.SetZero();
-}
-
-template<typename TKernelAux>
-void MultLocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
-  
-  // copy kernel type, center, and bandwidth squared
-  kernel_ = &(ka.kernel_);
-  sea_ = &(ka.sea_);
-  center_.Init(sea_->get_dimension());
-  order_ = -1;
-  ka_ = &ka;
-
-  // initialize coefficient array
-  coeffs_.Init(sea_->get_max_total_num_coeffs());
-  coeffs_.SetZero();
-}
-
-template<typename TKernelAux>
-template<typename TBound>
-int MultLocalExpansion<TKernelAux>::OrderForEvaluating
-(const TBound &far_field_region, 
- const TBound &local_field_region, double min_dist_sqd_regions,
- double max_dist_sqd_regions, double max_error, double *actual_error) const {
-  
-  return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region, 
-				     min_dist_sqd_regions,
-				     max_dist_sqd_regions, max_error, 
-				     actual_error);
-}
-
-template<typename TKernelAux>
-void MultLocalExpansion<TKernelAux>::TranslateFromFarField
-(const MultFarFieldExpansion<TKernelAux> &se) {
-
-  Vector pos_arrtmp, neg_arrtmp;
-  Matrix derivative_map;
-  Vector far_center;
-  Vector cent_diff;
-  Vector far_coeffs;
-  int dimension = sea_->get_dimension();
-  int far_order = se.get_order();
-  int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
-  int limit;
-  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
-
-  // get center and coefficients for far field expansion
-  far_center.Alias(*(se.get_center()));
-  far_coeffs.Alias(se.get_coeffs());
-  cent_diff.Init(dimension);
-
-  // if the order of the far field expansion is greater than the
-  // local one we are adding onto, then increase the order.
-  if(far_order > order_) {
-    order_ = far_order;
-  }
-
-  // compute Gaussian derivative
-  limit = 2 * order_ + 1;
-  derivative_map.Init(dimension, limit);
-  pos_arrtmp.Init(total_num_coeffs);
-  neg_arrtmp.Init(total_num_coeffs);
-
-  // compute center difference divided by bw_times_sqrt_two;
-  for(index_t j = 0; j < dimension; j++) {
-    cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
-  }
-
-  // compute required partial derivatives
-  ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
-  ArrayList<int> beta_plus_alpha;
-  beta_plus_alpha.Init(dimension);
-
-  // get the order of traversal for the given order of approximation
-  ArrayList<int> &traversal_order = sea_->traversal_mapping_[far_order];
-
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-
-    int index_j = traversal_order[j];
-    ArrayList<int> beta_mapping = sea_->get_multiindex(index_j);
-    pos_arrtmp[index_j] = neg_arrtmp[index_j] = 0;
-
-    for(index_t k = 0; k < total_num_coeffs; k++) {
-
-      int index_k = traversal_order[k];
-      ArrayList<int> alpha_mapping = sea_->get_multiindex(index_k);
-      for(index_t d = 0; d < dimension; d++) {
-	beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
-      }
-      double derivative_factor =
-	ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
-      
-      double prod = far_coeffs[index_k] * derivative_factor;
-
-      if(prod > 0) {
-	pos_arrtmp[index_j] += prod;
-      }
-      else {
-	neg_arrtmp[index_j] += prod;
-      }
-    } // end of k-loop
-  } // end of j-loop
-
-  Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-    int index_j = traversal_order[j];
-    coeffs_[index_j] += (pos_arrtmp[index_j] + neg_arrtmp[index_j]) * 
-      C_k_neg[index_j];
-  }
-}
-  
-template<typename TKernelAux>
-void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
-
-  // if no local coefficients have formed, then nothing to translate
-  if(order_ < 0) {
-    return;
-  }
-
-  // get the center and the order and the total number of coefficients of 
-  // the expansion we are translating from. Also get coefficients we
-  // are translating
-  Vector new_center;
-  new_center.Alias(*(se.get_center()));
-  int prev_order = se.get_order();
-  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
-  const ArrayList < int > *upper_mapping_index = 
-    sea_->get_upper_mapping_index();
-  Vector new_coeffs;
-  new_coeffs.Alias(se.get_coeffs());
-
-  // dimension
-  int dim = sea_->get_dimension();
-
-  // temporary variable
-  ArrayList<int> tmp_storage;
-  tmp_storage.Init(dim);
-
-  // sqrt two times bandwidth
-  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
-
-  // center difference between the old center and the new one
-  Vector center_diff;
-  center_diff.Init(dim);
-  for(index_t d = 0; d < dim; d++) {
-    center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
-  }
-
-  // set to the new order if the order of the expansion we are translating
-  // from is higher
-  if(prev_order < order_) {
-    se.set_order(order_);
-  }
-
-  // inverse multiindex factorials
-  Vector C_k;
-  C_k.Alias(sea_->get_inv_multiindex_factorials());
-
-  // get the order of traversal for the given order of approximation
-  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
-
-  // do the actual translation
-  for(index_t j = 0; j < total_num_coeffs; j++) {
-
-    int index_j = traversal_order[j];
-    const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_j);
-    const ArrayList <int> &upper_mappings_for_alpha = 
-      upper_mapping_index[index_j];
-    double pos_coeffs = 0;
-    double neg_coeffs = 0;
-
-    for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
-
-      if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
-	break;
-      }
-
-      const ArrayList<int> &beta_mapping = 
-	sea_->get_multiindex(upper_mappings_for_alpha[k]);
-      int flag = 0;
-      double diff1 = 1.0;
-
-      for(index_t l = 0; l < dim; l++) {
-	tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
-
-	if(tmp_storage[l] < 0) {
-	  flag = 1;
-	  break;
-	}
-      } // end of looping over dimension
-      
-      if(flag)
-	continue;
-
-      for(index_t l = 0; l < dim; l++) {
-	diff1 *= pow(center_diff[l], tmp_storage[l]);
-      }
-
-      double prod = coeffs_[upper_mappings_for_alpha[k]] * diff1 *
-	sea_->get_n_multichoose_k_by_pos
-	(upper_mappings_for_alpha[k], index_j);
-      
-      if(prod > 0) {
-	pos_coeffs += prod;
-      }
-      else {
-	neg_coeffs += prod;
-      }
-
-    } // end of k loop
-
-    new_coeffs[index_j] += pos_coeffs + neg_coeffs;
-  } // end of j loop
-}
+#define INSIDE_MULT_LOCAL_EXPANSION_H
+#include "mult_local_expansion_impl.h"
+#undef INSIDE_MULT_LOCAL_EXPANSION_H
 
 #endif
diff --git a/fastlib2/mlpack/series_expansion/mult_local_expansion_impl.h b/fastlib2/mlpack/series_expansion/mult_local_expansion_impl.h
new file mode 100644
index 0000000000..f8784b477b
--- /dev/null
+++ b/fastlib2/mlpack/series_expansion/mult_local_expansion_impl.h
@@ -0,0 +1,475 @@
+#ifndef INSIDE_MULT_LOCAL_EXPANSION_H
+#error "This is not a public header file!"
+#endif
+
+#ifndef MULT_LOCAL_EXPANSION_IMPL_H
+#define MULT_LOCAL_EXPANSION_IMPL_H
+
+template<typename TKernelAux>
+void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data, 
+						      const Vector& weights, 
+						      int begin, int end, 
+						      int order) {
+
+  if(order > order_) {
+    order_ = order;
+  }
+
+  int dim = sea_->get_dimension();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+  
+  // get inverse factorials (precomputed)
+  Vector neg_inv_multiindex_factorials;
+  neg_inv_multiindex_factorials.Alias
+    (sea_->get_neg_inv_multiindex_factorials());
+
+  // declare deritave mapping
+  Matrix derivative_map;
+  derivative_map.Init(dim, order + 1);
+  
+  // some temporary variables
+  Vector x_r_minus_x_Q;
+  x_r_minus_x_Q.Init(dim);
+  
+  // sqrt two times bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+  
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
+
+  // for each data point,
+  for(index_t r = begin; r < end; r++) {
+    
+    // calculate x_r - x_Q
+    for(index_t d = 0; d < dim; d++) {
+      x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) / 
+	bandwidth_factor;
+    }
+    
+    // precompute necessary partial derivatives based on coordinate difference
+    ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map);
+    
+    // compute h_{beta}((x_r - x_Q) / sqrt(2h^2))
+    for(index_t j = 0; j < total_num_coeffs; j++) {
+      int index = traversal_order[j];
+      const ArrayList<int> &mapping = sea_->get_multiindex(index);
+      double partial_derivative = 
+	ka_->ComputePartialDerivative(derivative_map, mapping);
+      coeffs_[index] += neg_inv_multiindex_factorials[index] * weights[r] * 
+	partial_derivative;
+    }
+  } // End of looping through each reference point.
+}
+
+template<typename TKernelAux>
+void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name, 
+						FILE *stream) const {
+  
+    
+  int dim = sea_->get_dimension();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
+  fprintf(stream, "Local expansion\n");
+  fprintf(stream, "Center: ");
+  
+  for (index_t i = 0; i < center_.length(); i++) {
+    fprintf(stream, "%g ", center_[i]);
+  }
+  fprintf(stream, "\n");
+  
+  fprintf(stream, "f(");
+  for(index_t d = 0; d < dim; d++) {
+    fprintf(stream, "x_q%d", d);
+    if(d < dim - 1)
+      fprintf(stream, ",");
+  }
+  fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
+  
+  for (index_t i = 0; i < total_num_coeffs; i++) {
+    ArrayList<int> mapping = sea_->get_multiindex(i);
+    fprintf(stream, "%g", coeffs_[i]);
+    
+    for(index_t d = 0; d < dim; d++) {
+      fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
+    }
+
+    if(i < total_num_coeffs - 1) {
+      fprintf(stream, " + ");
+    }
+  }
+  fprintf(stream, "\n");
+}
+
+template<typename TKernelAux>
+double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data, 
+						     int row_num) const {
+    
+  // if there are no local coefficients, then return 0
+  if(order_ < 0) {
+    return 0;
+  }
+
+  // total number of coefficient
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  // number of dimensions
+  int dim = sea_->get_dimension();
+
+  // evaluated sum to be returned
+  double sum = 0;
+  
+  // sqrt two bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+
+  // temporary variable
+  Vector x_Q_to_x_q;
+  x_Q_to_x_q.Init(dim);
+  Vector tmp;
+  tmp.Init(sea_->get_max_total_num_coeffs());
+  ArrayList<int> heads;
+  heads.Init(dim + 1);
+  
+  // compute (x_q - x_Q) / (sqrt(2h^2))
+  for(index_t i = 0; i < dim; i++) {
+    x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
+  }
+  
+  for(index_t i = 0; i < dim; i++)
+    heads[i] = 0;
+  heads[dim] = MAXINT;
+
+  tmp[0] = 1.0;
+
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
+
+  for(index_t i = 1; i < total_num_coeffs; i++) {
+    
+    int index = traversal_order[i];
+    const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
+    
+    // from the direct descendant, recursively compute the multipole moments
+    int direct_ancestor_mapping_pos = 
+      lower_mappings[lower_mappings.size() - 2];
+    int position = 0;
+    const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
+    const ArrayList<int> &direct_ancestor_mapping = 
+      sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
+    for(index_t i = 0; i < dim; i++) {
+      if(mapping[i] != direct_ancestor_mapping[i]) {
+	position = i;
+	break;
+      }
+    }
+    tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
+  }
+
+  for(index_t i = 0; i < total_num_coeffs; i++) {
+    int index = traversal_order[i];
+    sum += coeffs_[index] * tmp[index];
+  }
+
+  return sum;
+}
+
+template<typename TKernelAux>
+double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
+    
+  // if there are no local coefficients, then return 0
+  if(order_ < 0) {
+    return 0;
+  }
+
+  // total number of coefficient
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+
+  // number of dimensions
+  int dim = sea_->get_dimension();
+
+  // evaluated sum to be returned
+  double sum = 0;
+  
+  // sqrt two bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_.bandwidth_sq());
+
+  // temporary variable
+  Vector x_Q_to_x_q;
+  x_Q_to_x_q.Init(dim);
+  Vector tmp;
+  tmp.Init(sea_->get_max_total_num_coeffs());
+  ArrayList<int> heads;
+  heads.Init(dim + 1);
+  
+  // compute (x_q - x_Q) / (sqrt(2h^2))
+  for(index_t i = 0; i < dim; i++) {
+    x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
+  }
+  
+  for(index_t i = 0; i < dim; i++)
+    heads[i] = 0;
+  heads[dim] = MAXINT;
+
+  tmp[0] = 1.0;
+
+  // get the order of traversal for the given order of approximation
+  ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
+
+  for(index_t i = 1; i < total_num_coeffs; i++) {
+    
+    int index = traversal_order[i];
+    ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
+    
+    // from the direct descendant, recursively compute the multipole moments
+    int direct_ancestor_mapping_pos = 
+      lower_mappings[lower_mappings.size() - 2];
+    int position = 0;
+    const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
+    const ArrayList<int> &direct_ancestor_mapping = 
+      sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
+    for(index_t i = 0; i < dim; i++) {
+      if(mapping[i] != direct_ancestor_mapping[i]) {
+	position = i;
+	break;
+      }
+    }
+    tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
+  }
+
+  for(index_t i = 0; i < total_num_coeffs; i++) {
+    int index = traversal_order[i];
+    sum += coeffs_[index] * tmp[index];
+  }
+
+  return sum;
+}
+
+template<typename TKernelAux>
+void MultLocalExpansion<TKernelAux>::Init(const Vector& center, 
+					  const TKernelAux &ka) {
+  
+  // copy kernel type, center, and bandwidth squared
+  kernel_ = &(ka.kernel_);
+  center_.Copy(center);
+  order_ = -1;
+  sea_ = &(ka.sea_);
+  ka_ = &ka;
+
+  // initialize coefficient array
+  coeffs_.Init(sea_->get_max_total_num_coeffs());
+  coeffs_.SetZero();
+}
+
+template<typename TKernelAux>
+void MultLocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
+  
+  // copy kernel type, center, and bandwidth squared
+  kernel_ = &(ka.kernel_);
+  sea_ = &(ka.sea_);
+  center_.Init(sea_->get_dimension());
+  order_ = -1;
+  ka_ = &ka;
+
+  // initialize coefficient array
+  coeffs_.Init(sea_->get_max_total_num_coeffs());
+  coeffs_.SetZero();
+}
+
+template<typename TKernelAux>
+template<typename TBound>
+int MultLocalExpansion<TKernelAux>::OrderForEvaluating
+(const TBound &far_field_region, 
+ const TBound &local_field_region, double min_dist_sqd_regions,
+ double max_dist_sqd_regions, double max_error, double *actual_error) const {
+  
+  return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region, 
+				     min_dist_sqd_regions,
+				     max_dist_sqd_regions, max_error, 
+				     actual_error);
+}
+
+template<typename TKernelAux>
+void MultLocalExpansion<TKernelAux>::TranslateFromFarField
+(const MultFarFieldExpansion<TKernelAux> &se) {
+
+  Vector pos_arrtmp, neg_arrtmp;
+  Matrix derivative_map;
+  Vector far_center;
+  Vector cent_diff;
+  Vector far_coeffs;
+  int dimension = sea_->get_dimension();
+  int far_order = se.get_order();
+  int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
+  int limit;
+  double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
+
+  // get center and coefficients for far field expansion
+  far_center.Alias(*(se.get_center()));
+  far_coeffs.Alias(se.get_coeffs());
+  cent_diff.Init(dimension);
+
+  // if the order of the far field expansion is greater than the
+  // local one we are adding onto, then increase the order.
+  if(far_order > order_) {
+    order_ = far_order;
+  }
+
+  // compute Gaussian derivative
+  limit = 2 * order_ + 1;
+  derivative_map.Init(dimension, limit);
+  pos_arrtmp.Init(total_num_coeffs);
+  neg_arrtmp.Init(total_num_coeffs);
+
+  // compute center difference divided by bw_times_sqrt_two;
+  for(index_t j = 0; j < dimension; j++) {
+    cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
+  }
+
+  // compute required partial derivatives
+  ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
+  ArrayList<int> beta_plus_alpha;
+  beta_plus_alpha.Init(dimension);
+
+  // get the order of traversal for the given order of approximation
+  ArrayList<int> &traversal_order = sea_->traversal_mapping_[far_order];
+
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+
+    int index_j = traversal_order[j];
+    ArrayList<int> beta_mapping = sea_->get_multiindex(index_j);
+    pos_arrtmp[index_j] = neg_arrtmp[index_j] = 0;
+
+    for(index_t k = 0; k < total_num_coeffs; k++) {
+
+      int index_k = traversal_order[k];
+      ArrayList<int> alpha_mapping = sea_->get_multiindex(index_k);
+      for(index_t d = 0; d < dimension; d++) {
+	beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
+      }
+      double derivative_factor =
+	ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
+      
+      double prod = far_coeffs[index_k] * derivative_factor;
+
+      if(prod > 0) {
+	pos_arrtmp[index_j] += prod;
+      }
+      else {
+	neg_arrtmp[index_j] += prod;
+      }
+    } // end of k-loop
+  } // end of j-loop
+
+  Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+    int index_j = traversal_order[j];
+    coeffs_[index_j] += (pos_arrtmp[index_j] + neg_arrtmp[index_j]) * 
+      C_k_neg[index_j];
+  }
+}
+  
+template<typename TKernelAux>
+void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
+
+  // if no local coefficients have formed, then nothing to translate
+  if(order_ < 0) {
+    return;
+  }
+
+  // get the center and the order and the total number of coefficients of 
+  // the expansion we are translating from. Also get coefficients we
+  // are translating
+  Vector new_center;
+  new_center.Alias(*(se.get_center()));
+  int prev_order = se.get_order();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order_);
+  const ArrayList < int > *upper_mapping_index = 
+    sea_->get_upper_mapping_index();
+  Vector new_coeffs;
+  new_coeffs.Alias(se.get_coeffs());
+
+  // dimension
+  int dim = sea_->get_dimension();
+
+  // temporary variable
+  ArrayList<int> tmp_storage;
+  tmp_storage.Init(dim);
+
+  // sqrt two times bandwidth
+  double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
+
+  // center difference between the old center and the new one
+  Vector center_diff;
+  center_diff.Init(dim);
+  for(index_t d = 0; d < dim; d++) {
+    center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
+  }
+
+  // set to the new order if the order of the expansion we are translating
+  // from is higher
+  if(prev_order < order_) {
+    se.set_order(order_);
+  }
+
+  // inverse multiindex factorials
+  Vector C_k;
+  C_k.Alias(sea_->get_inv_multiindex_factorials());
+
+  // get the order of traversal for the given order of approximation
+  const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
+
+  // do the actual translation
+  for(index_t j = 0; j < total_num_coeffs; j++) {
+
+    int index_j = traversal_order[j];
+    const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_j);
+    const ArrayList <int> &upper_mappings_for_alpha = 
+      upper_mapping_index[index_j];
+    double pos_coeffs = 0;
+    double neg_coeffs = 0;
+
+    for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
+
+      if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
+	break;
+      }
+
+      const ArrayList<int> &beta_mapping = 
+	sea_->get_multiindex(upper_mappings_for_alpha[k]);
+      int flag = 0;
+      double diff1 = 1.0;
+
+      for(index_t l = 0; l < dim; l++) {
+	tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
+
+	if(tmp_storage[l] < 0) {
+	  flag = 1;
+	  break;
+	}
+      } // end of looping over dimension
+      
+      if(flag)
+	continue;
+
+      for(index_t l = 0; l < dim; l++) {
+	diff1 *= pow(center_diff[l], tmp_storage[l]);
+      }
+
+      double prod = coeffs_[upper_mappings_for_alpha[k]] * diff1 *
+	sea_->get_n_multichoose_k_by_pos
+	(upper_mappings_for_alpha[k], index_j);
+      
+      if(prod > 0) {
+	pos_coeffs += prod;
+      }
+      else {
+	neg_coeffs += prod;
+      }
+
+    } // end of k loop
+
+    new_coeffs[index_j] += pos_coeffs + neg_coeffs;
+  } // end of j loop
+}
+
+#endif