#include <iostream>
#include <cassert>
#include <climits>
#include <cstddef>

// We have previously seen that calculating
// and storing intermediate results in an
// array can significantly reduce the running
// time of a calculation.
//
// One previously slow calculation was finding
// the number of integer partitions into
// distinct positive parts.
//
// For example, the partitions of 9 are
//
//     9
//     8 + 1
//     7 + 2
//     6 + 3
//     6 + 2 + 1
//     5 + 4
//     5 + 3 + 1
//     4 + 3 + 2
//
// Thus, q(9) = 8.
//
// In Abramowitz and Stegun, we find the
// following recurrence:
//
//             1   n
//     q(n) = --- SUM q(n - k)*(
//             n  k=1
//
//         sigma_1(k) - 2*sigma_1(k/2)
//     ).
//
// Here, sigma_1(k) is the sum of the
// positive divisors of k.
//
// For example:
//
//     sigma_1(6) = 1 + 2 + 3 + 6 = 12.
//
// If k is odd, sigma_1(k/2) is taken
// to be zero.
//
// The initial value is
//
//     q(0) = 1.
//
// This represents the single partition
// of zero: the empty partition.
//
// You have already implemented functions
// for calculating sums of powers of divisors.
// These functions are supplied below.
//
// PROBLEM
//
// Implement two functions.
//
// A. Using a working array:
//
//    Allocate an array of capacity n + 1.
//
//    Initialize its first entry to 1,
//    corresponding to q(0).
//
//    Calculate q(1), q(2), ..., q(n)
//    successively, storing each result
//    in the array.
//
//    Instead of making recursive calls,
//    access the previously calculated
//    values in the working array.
//
//    Return q(n) and release all
//    dynamically allocated memory.
//
// B. Returning an array:
//
//    Allocate an array of capacity n + 1.
//
//    Calculate and store q(0), ..., q(n)
//    as described above.
//
//    Return the address of the array.
//
//    The calling function is responsible
//    for releasing the returned array.
//
// For this exercise, assume n <= 30.

// Function declarations
int main();

unsigned long distinct_integer_partitions_memo(
    unsigned long const n
);

unsigned long *distinct_integer_partitions_array(
    unsigned long const n
);

unsigned long *sigma_array(
    unsigned long const n,
    unsigned long const k = 1
);

unsigned long integer_power(
    unsigned long const n,
    unsigned long const m
);

// Function definitions
int main() {
    // The expected values q(0), ..., q(29) are
    //
    //  1   1   1   2   2   3   4   5   6   8
    // 10  12  15  18  22  27  32  38  46  54
    // 64  76  89 104 122 142 165 192 222 256

    unsigned long const N{ 29 };

    for ( unsigned long m{ 0 }; m <= N; ++m ) {
        std::cout << "q(" << m << ") = "
                  << distinct_integer_partitions_memo( m )
                  << std::endl;
    }

    std::cout << std::endl;

    unsigned long *a_q{
        distinct_integer_partitions_array( N )
    };

    assert( a_q != nullptr );

    for ( unsigned long m{ 0 }; m <= N; ++m ) {
        std::cout << "q(" << m << ") = "
                  << a_q[m] << std::endl;

        assert(
            a_q[m]
            == distinct_integer_partitions_memo( m )
        );
    }

    assert( a_q[0] == 1 );
    assert( a_q[9] == 8 );
    assert( a_q[19] == 54 );
    assert( a_q[29] == 256 );

    // The returned array was dynamically
    // allocated and must be released.
    delete[] a_q;
    a_q = nullptr;

    return 0;
}

// Calculate the sum of the kth powers
// of the positive divisors of each
// integer from 1 through n.
//
// The returned array has capacity n + 1.
// Entry zero is initialized to zero.
//
// The caller must release the array.

unsigned long *sigma_array(
    unsigned long const n,
    unsigned long const k
) {
    unsigned long *a_values{
        new unsigned long[n + 1]{}
    };

    for (
        unsigned long divisor{ 1 };
        divisor <= n;
        ++divisor
    ) {
        unsigned long const contribution{
            integer_power( divisor, k )
        };

        // Visit every multiple of divisor.
        for (
            unsigned long multiple{ divisor };
            multiple <= n;
            multiple += divisor
        ) {
            assert(
                contribution
                <= ULONG_MAX - a_values[multiple]
            );

            a_values[multiple] += contribution;
        }
    }

    return a_values;
}

// Calculate n raised to the power m
// using exponentiation by squaring.
//
// For example, since
//
//     42 = 32 + 8 + 2,
//
// we have
//
//     n^42 = n^32 * n^8 * n^2.
//
// Repeated squaring calculates these
// powers efficiently.
//
// Assume the result and intermediate
// calculations fit in unsigned long.

unsigned long integer_power(
    unsigned long const n,
    unsigned long const m
) {
    unsigned long result{ 1 };
    unsigned long base{ n };
    unsigned long exponent{ m };

    while ( exponent != 0 ) {
        if ( exponent % 2 == 1 ) {
            assert(
                base == 0
                || result <= ULONG_MAX/base
            );

            result *= base;
        }

        exponent /= 2;

        if ( exponent != 0 ) {
            assert(
                base == 0
                || base <= ULONG_MAX/base
            );

            base *= base;
        }
    }

    return result;
}

unsigned long distinct_integer_partitions_memo(
    unsigned long const n
) {
    assert( n <= 30 );

    unsigned long *a_sigmas{
        sigma_array( n )
    };

    // Your implementation goes here.
    // Remember to release a_sigmas.

    delete[] a_sigmas;
    return 0;
}

unsigned long *distinct_integer_partitions_array(
    unsigned long const n
) {
    assert( n <= 30 );

    unsigned long *a_sigmas{
        sigma_array( n )
    };

    // Your implementation goes here.
    // Remember to release a_sigmas.

    delete[] a_sigmas;
    return nullptr;
}

























// First hint:
//
// By definition,
//
//     q(0) = 1.
//
// The recurrence involves sigma_1(k/2).
//
// However, k/2 is not always an integer.
//
// What should this term contribute
// when k is odd?
//
// Also, notice that q(n) depends only
// on values q(0), ..., q(n - 1).
//
// Therefore, all the required values
// will already be available when
// calculating the next array entry.

























// Second hint:
//
// When k is odd, take
//
//     sigma_1(k/2) = 0.
//
// Thus, the coefficient is
//
//     sigma_1(k)
//
// for odd k, and
//
//     sigma_1(k) - 2*sigma_1(k/2)
//
// for even k.
//
// For each array entry q(k), calculate
//
//     sum = 0.
//
// Then, for i = 1, ..., k, add
//
//     q(k - i)*sigma_1(i)
//
// if i is odd, or
//
//     q(k - i)*(
//         sigma_1(i) - 2*sigma_1(i/2)
//     )
//
// if i is even.
//
// Finally, divide the accumulated
// sum by k to obtain q(k).
//
// You can access every required
// partition number directly from
// the working array.

























// Here is one solution:

/*
unsigned long distinct_integer_partitions_memo(
    unsigned long const n
) {
    assert( n <= 30 );

    unsigned long *a_sigmas{
        sigma_array( n )
    };

    unsigned long *a_values{
        new unsigned long[n + 1]{}
    };

    a_values[0] = 1;

    for ( unsigned long k{ 1 }; k <= n; ++k ) {
        unsigned long sum{ 0 };

        for ( unsigned long i{ 1 }; i <= k; ++i ) {
            unsigned long coefficient{ a_sigmas[i] };

            if ( i % 2 == 0 ) {
                assert(
                    2*a_sigmas[i/2] <= coefficient
                );

                coefficient -= 2*a_sigmas[i/2];
            }

            sum += coefficient*a_values[k - i];
        }

        assert( sum % k == 0 );
        a_values[k] = sum/k;
    }

    unsigned long const result{ a_values[n] };

    delete[] a_values;
    delete[] a_sigmas;

    return result;
}

unsigned long *distinct_integer_partitions_array(
    unsigned long const n
) {
    assert( n <= 30 );

    unsigned long *a_sigmas{
        sigma_array( n )
    };

    unsigned long *a_values{
        new unsigned long[n + 1]{}
    };

    a_values[0] = 1;

    for ( unsigned long k{ 1 }; k <= n; ++k ) {
        unsigned long sum{ 0 };

        for ( unsigned long i{ 1 }; i <= k; ++i ) {
            unsigned long coefficient{ a_sigmas[i] };

            if ( i % 2 == 0 ) {
                assert(
                    2*a_sigmas[i/2] <= coefficient
                );

                coefficient -= 2*a_sigmas[i/2];
            }

            sum += coefficient*a_values[k - i];
        }

        assert( sum % k == 0 );
        a_values[k] = sum/k;
    }

    delete[] a_sigmas;

    return a_values;
}
*/