/* powi_table generator.
 * Generates a look-up table for generating efficient integer
 * exponentiations by a positive constant integer exponent.
 *
 * Copyright (c) 2016 Joel Yliluoma (http://iki.fi/bisqwit/)
 * License: CC-BY-SA 4.0
 *
 * See section 4.6.3, "Evaluation of Powers" of
 * Donald E. Knuth, "Seminumerical Algorithms", Vol. 2,
 * "The Art of Computer Programming", 3rd Edition, 1998.
 */

#include <utility>
#include <cstdio>
#include <climits>

// std::bitset<> does not have intersect or union methods, so we create our own bitset.
template<std::size_t Size>
class Set
{
    using elem_t = unsigned long;
    static constexpr unsigned bp       = CHAR_BIT * sizeof(elem_t);
    static constexpr std::size_t bound = (Size + bp-1) / bp;
    elem_t data[bound] {};
public:
    void set(std::size_t n) { data[elem(n)] |= mask(n); }
    Set<Size> combine(const Set<Size>& b) const
    {
        if(&b == this) return *this;
        Set<Size> result;
        #pragma omp simd
        for(std::size_t n=0; n<bound; ++n)
            result.data[n] = data[n] | b.data[n];
        return result;
    }
    std::size_t count() const
    {
        std::size_t result=0;
        #pragma omp simd reduction(+:result)
        for(std::size_t n=0; n<bound; ++n)
            result += __builtin_popcountl(data[n]); // GCC, ICC or Clang only
        return result;
    }
private:
    std::size_t elem(std::size_t n) const { return n/bp; };
    elem_t      mask(std::size_t n) const { return elem_t(1) << (n%bp); }
};

constexpr unsigned MaxPow = 4096;
static Set<MaxPow> results[MaxPow];

int main()
{
    /* Each bit in the bitset represents one multiplication.
     * For each powi, we subdivide it into two components that are then
     * multiplied together. The subdivision is chosen such that the number
     * of total multiplications remains minimal. This is done by using
     * the set-intersection method (combine()), which is the sum of
     * number of bits set in both sides of the multiplication, except
     * if the same bit is set on both sides, it is counted only once.
     */
    results[1].set(1);
    unsigned total_multiplications = 0;
    for(unsigned pow=0; pow<MaxPow; ++pow)
    {
        typedef std::pair<unsigned,unsigned> p;
        #pragma omp declare reduction(better:p: omp_in < omp_out ? omp_out=omp_in : omp_out) \
                            initializer(omp_priv = omp_orig)
        p best(~0u, pow);
        #pragma omp parallel for simd reduction(better:best)
        for(unsigned a=pow/2; a>=1; --a)
        {
            unsigned count = results[a].combine(results[pow-a]).count();
            if(count < best.first) best = std::make_pair(count, pow-a);
        }
        if(pow >= 2)
        {
            results[pow] = results[best.second].combine(results[pow-best.second]);
            results[pow].set(pow);
            total_multiplications += best.first;
        }
        std::printf("%4d,", best.second);
        if(pow%16==15) std::printf(" /*%4d -%4d */\n", pow-15,pow);
    }
    std::printf("/* Total multiplications: %u */\n", total_multiplications);
}