#include <cmath>

    #include <algorithm>
    #include <iostream>
    #include <vector>

    // The resulting vector contains at index i the sum of 2^-j for j in [i+1, N]
    // and is padded with one 0 to get the same length as `v`
    static std::vector<double> partialSums(std::vector<double> const& v) {
        std::vector<double> result;

        // When summing doubles, we need to start with the smaller ones
        // because of the precision of representations...

        double sum = 0;
        for (std::vector<double>::const_reverse_iterator it = v.rbegin(), end = v.rend(); it != end; ++it)
        {
            sum += *it;
            result.push_back(sum);
        }

        result.pop_back(); // there is a +1 offset in the indexes of the result

        std::reverse(result.begin(), result.end());

        result.push_back(0); // pad the vector to have the same length as `v`

        return result;   
    }

    // The resulting vector contains the indexes elected
    static std::vector<size_t> getIndexesImpl(std::vector<double> const& v,
                                              std::vector<double> const& ps,
                                              double r)
    {
      std::vector<size_t> indexes;

      for (size_t i = 0, max = v.size(); i != max; ++i) {
          if (r >= v[i]) {
              std::cout << ".." << r << ">=" << v[i] << "\n";
              r -= v[i];
              indexes.push_back(i);
              continue;
          }

          // We favor the closest to 0 in case of equality
          if (std::fabs(r - v[i]) < std::fabs(r - ps[i])) {
              std::cout << ".." << v[i] - r << " < " << r - ps[i] << "\n";
              indexes.push_back(i);
              return indexes;
          }
      }

      return indexes;
    }

    std::vector<size_t> getIndexes(std::vector<double>& v, double r) {
        std::vector<double> const ps = partialSums(v);
        return getIndexesImpl(v, ps, r);
    }

    void printIndexes(std::vector<double>& v, double r) {
        std::vector<size_t> indexes = getIndexes(v, r);
        std::cout << r << ":\n";
        for (size_t i = 0, max = indexes.size(); i != max; ++i) {
            std::cout << "   " << indexes[i] << ": " << v[indexes[i]] << "\n";
        }

        std::reverse(indexes.begin(), indexes.end());
        double sum = 0;
        for (size_t i = 0, max = indexes.size(); i != max; ++i) {
            sum += v[indexes[i]];
        }
        std::cout << "=> " << sum << "\n";
    }

    int main() {
        std::vector<double> v;
        v.push_back(1.0/2);
        v.push_back(1.0/4);
        v.push_back(1.0/8);
        v.push_back(1.0/16);
        v.push_back(1.0/32);

        printIndexes(v, 0.3);
        printIndexes(v, 0.25665238);
    }