#include <functional>
#include <vector>
#include <set>
#include <map>
#include <utility>
#include <iostream>

namespace jcui {
namespace algorithm {

// Return first location where comp(a, b) is not true.
template<class ForwardIterator1, class ForwardIterator2, class Compare>
ForwardIterator1 lower_bound(ForwardIterator1 begin1,
			     ForwardIterator1 end1,
			     ForwardIterator2 begin2,
			     ForwardIterator2 end2,
			     Compare comp) {
  int len = end1 - begin1;
  while (len > 0) {
    int half = len >> 1;
    ForwardIterator1 m1 = begin1 + half;
    ForwardIterator2 m2 = begin2 + half;
    if (comp(*m1, *m2)) {
      begin1 = m1 + 1;
      begin2 = m2 + 1;
      len -= (half + 1);
    } else {
      len = half;
    }
  }
  return begin1;
}

// Return min_x max(inc[x], dec[x]) where inc and dec are non-decreasing and non-increasing   
template<class T, class ForwardIterator1, class ForwardIterator2> 
T min_max(ForwardIterator1 inc_begin,
	  ForwardIterator1 inc_end,
	  ForwardIterator2 dec_begin,
	  ForwardIterator2 dec_end,
	  T default_value) {
  int len = inc_end - inc_begin;
  if (len == 0) return default_value;
  ForwardIterator1 it = lower_bound(inc_begin, inc_end, dec_begin, dec_end, std::less<T>());
  if (it == inc_begin) return *inc_begin;
  if (it == inc_end) return *(dec_begin + len - 1);
  return min(*it, *(dec_begin + (it - inc_begin - 1)));
}

template<class T>
T pow(const T& x, int k) {
  if (k <= 1) return x;
  T y = pow(x, k / 2);
  y = y * y;
  if (k % 2) y = y * x;
  return y;
}
  
template<class T>
class Mat {
public:
  Mat(int H, int W) {
    v.resize(H);
    for (int i = 0; i < H; ++i) v[i].resize(W);
  }
  int H() const { return v.size(); }
  int W() const { return v.size() ? v[0].size() : 0; }  
  void set(int i, int j, T val) { if (valid(i, j)) v[i][j] = val; }
  T get(int i, int j) const { return valid(i, j) ? v[i][j] : 0; }
  void normalize_by_row_sum() {
    std::vector<T> s = Mat<T>::row_sum(*this);
    for (int i = 0; i < H(); ++i) {
      for (int j = 0; j < W(); ++j) {
	v[i][j] /= s[j];
      }
    }
  }
  static std::vector<T> row_sum(const Mat& A) {
    std::vector<T> y(A.W(), 0);
    for (int j = 0; j < A.W(); ++j) {
      for (int i = 0; i < A.H(); ++i) {
	y[j] += A.v[i][j];
      }
    }
    return y;
  }
  Mat operator*(const Mat& y) const {
    Mat p(H(), y.W());
    for (int i = 0; i < H(); ++i) {
      for (int j = 0; j < y.W(); ++j) {
	T s = 0;
	for (int k = 0; k < W(); ++k) {
	  s += v[i][k] * y.v[k][j];
	}
	p.v[i][j] = s;
      }      
    }
    return p;
  }
  std::vector<T> operator*(const std::vector<T>& x) const {
    std::vector<T> y(H(), 0);
    for (int i = 0; i < H(); ++i) {
      for (int j = 0; j < W(); ++j) {
	y[i] += v[i][j] * x[j];
      }
    }
    return y;
  }
  inline bool valid(int i, int j) const { return i >= 0 && i < H() && j >= 0 && j < W(); }
  void print() const {
    for (int i = 0; i < H(); ++i) {
      for (int j = 0; j < W(); ++j) {
	std::cout << v[i][j] << " ";
      }
      std::cout << std::endl;
    }
  }
  std::vector<std::vector<T> > v;
};

template<class T>
class SparseMat {
public:
  typedef typename std::map<std::pair<int, int>, T>::const_iterator MapConstIter;
  typedef std::set<int>::iterator SetIter;

  SparseMat(int H, int W) {
    row_ids.resize(W);
    col_ids.resize(H);
  }
  int H() const { return col_ids.size(); }
  int W() const { return row_ids.size(); }
  void set(int i, int j, T val) {
    if (i < 0 || i >= H() || j < 0 || j >= W()) return;
    v[std::make_pair(i, j)] = val;
    rows.insert(i);
    cols.insert(j);
    row_ids[j].insert(i);
    col_ids[i].insert(j);
  }
  T get(int i, int j) const {
    MapConstIter it = v.find(std::make_pair(i, j));
    if (it == v.end()) return 0;
    return it->second;
  }
  void normalize_by_row_sum() {
    std::vector<T> s = SparseMat<T>::row_sum(*this);
    for (SetIter col = cols.begin(); col != cols.end(); ++col) {
      for (SetIter row = row_ids[*col].begin(); row != row_ids[*col].end(); ++row) {
	v[std::make_pair(*row, *col)] /= s[*col];
      }
    }
  }
  static std::vector<T> row_sum(const SparseMat& A) {
    std::vector<T> y(A.W(), 0);
    for (SetIter col = A.cols.begin(); col != A.cols.end(); ++col) {
      for (SetIter row = A.row_ids[*col].begin(); row != A.row_ids[*col].end(); ++row) {
	y[*col] += A.v.find(std::make_pair(*row, *col))->second;
      }
    }
    return y;
  }
  std::vector<T> operator*(const std::vector<T>& x) const {
    std::vector<T> y(H(), 0);
    for (SetIter r = rows.begin(); r != rows.end(); ++r) {
      for (SetIter k = col_ids[*r].begin(); k != col_ids[*r].end(); ++k) {
	y[*r] += v.find(std::make_pair(*r, *k))->second * x[*k];
      }
    }
    return y;
  }
  SparseMat operator*(const SparseMat& y) const {
    SparseMat p(H(), y.W());
    std::vector<int> intersect(W());
    typedef std::vector<int>::iterator VecIter;
    for (SetIter r = rows.begin(); r != rows.end(); ++r) {
      for (SetIter c = y.cols.begin(); c != y.cols.end(); ++c) {
	T s = 0;
	VecIter intersect_end = set_intersection(col_ids[*r].begin(),
						 col_ids[*r].end(),
						 y.row_ids[*c].begin(), 
						 y.row_ids[*c].end(), 
						 intersect.begin());
	for (VecIter k = intersect.begin(); k != intersect_end; ++k) {
	  s += v.find(std::make_pair(*r, *k))->second * y.v.find(std::make_pair(*k, *c))->second;
	}
	if (s != 0) p.set(*r, *c, s);
      }
    }
    return p;
  }
  void print() const {
    for (MapConstIter i = v.begin(); i != v.end(); ++i) {
      std::cout << "(" << i->first.first << ", " << i->first.second << ") = " << i->second << std::endl;
    }
  }
  std::map<std::pair<int, int>, T> v;
  std::set<int> rows, cols;
  std::vector<std::set<int> > row_ids, col_ids;
};

template<class T>
class SparseMatCSR {
public:
  SparseMatCSR(const SparseMat<T>& m) {
    v.resize(m.v.size());
    col.resize(m.v.size());
    cum_n.resize(1 + m.H());
    if (v.size() == 0) return;

    cum_n[0] = 0;
    int i = 0;
    for (int r = 0; r < H(); ++r) {
      for (std::set<int>::iterator k = m.col_ids[r].begin(); k != m.col_ids[r].end(); ++k) {
	v[i] = m.v.find(std::make_pair(r, *k))->second;
	col[i] = *k;
	i ++;
      }
      cum_n[r + 1] = i;
    }
  }
  int H() const { return cum_n.size() - 1; }
  std::vector<T> operator*(const std::vector<T>& x) const {
    std::vector<T> y(H(), 0);
    for (int i = 0; i < H(); ++i) {
      for (int ind = cum_n[i]; ind < cum_n[i + 1]; ++ind) {
	y[i] += x[col[ind]] * v[ind];
      }
    }
    return y;
  }
  /* Example: For m = [0, 0, 1, 0; 2, 1, 0, 0; 0, 0, 1, 0],
     v = [1 2 1 1]
     col = [2 0 1 2]
     cum_n = [0 1 3 4] */
  std::vector<T> v;
  std::vector<int> col;
  std::vector<int> cum_n;
};

template<class T>
class RingBuffer {
public:
  RingBuffer(int N) : cur(0), size(0) { v.resize(N); } // N > 0
  void push_back(const T& val) {
    v[cur] = val;
    cur ++;
    if (cur >= v.size()) cur = 0;
    if (size < v.size()) size ++;
  }
  // Get value from last_write_position + offset
  T get(int offset, T default_val = 0) const {
    if (offset > 0 || offset <= -size) return default_val; // Cannot read future info, or info that got flushed
    int index = cur - 1 + offset;
    if (index < 0) index += size;
    return v[index];
  }     
  int N() const { return size; }
private:
  int cur;
  int size;
  std::vector<T> v;
};

// Split integer N using non-repeating numbers in [1, k]
std::vector<std::vector<int> > split_no_repeat(int N, int k) {
  std::vector<std::vector<int> > res;
  if (N <= 0 || N > k * (k + 1) / 2) return res;

  if (N == k) res.push_back(std::vector<int>(1, k));

  std::vector<std::vector<int> > partial = split_no_repeat(N - k, k - 1);
  for (int i = 0; i < partial.size(); ++i) {
    partial[i].push_back(k);
    res.push_back(partial[i]);
  }
  
  std::vector<std::vector<int> > rest = split_no_repeat(N, k - 1);
  res.insert(res.end(), rest.begin(), rest.end());

  return res;
}

} // namespace algorithm
} // namespace jcui