Calculate optimal game upgrades, v3

Question

This is the third iteration of the code from Calculate optimal game upgrades, v2 and Calculate optimal game upgrades. I'm looking for suggestions to further optimize this code. A brief summary of the problem:

We are given a certain number of points to spend on 4 upgrades. The cost is as follows:

• Gold: 15 points
• Dollars: 6 points
• Dimes: 3 points
• Pennies: 1 point

The total bonus from these upgrades is equal to pennies * (1 + dimes * (1 + dollars *(1 + gold))). Find the optimal allocation of points between upgrades that maximizes the bonus.

The algorithm (roughly) works as follows:

1. get_max_score(int32_t points_left): For every possible amount of pennies, calculate the optimal amount of dimes, dollars and gold (get_max_dimes) given the points leftover after purchasing pennies. Return the combination of pennies and dimes, dollars, gold that has the highest score.
2. get_max_dimes(int32_t points_left): The optimization problem for this function is to maximize 1 + dimes * (1 + dollars * (1 + gold)), given dimes * 3 + dollars * 6 + gold * 15 <= points_left. For every possible amount of dimes, we calculate the optimal amount of dollars and gold (get_max_dollars) given the points leftover after purchasing dimes. Return the combination of dimes and dollars, gold that has the highest sub-bonus.
3. get_max_dollars(int32_t points_left): The optimization problem for this function is to maximize 1 + dollars * (1 + gold), given dollars * 6 + gold * 15 <= points_left. This function uses the closed form solution for the real number variant of this problem to approximate the solution, and then checks all dollar values that are [-5, 5] away from the approximate optimal solution, then return the combination of dollars & gold that has the highest sub-bonus.

Since get_max_dimes and get_max_dollars are called many times with the same values, we memoize those functions using the memoized template.

The approach described above is the same used in the v2 Python code. This code adds one more optimization by pruning the search space, by skipping intervals where the maximum possible score is lower than our approximate solution (which simply spends 1/4 of the points on pennies). As an example, say we have 1000 points to spend. Instead of looping from 1-1000 pennies, we can split the interval into 2: 1-500 pennies and 501-1000 pennies. The maximum value of get_max_dimes on the first interval will be get_max_dimes(500) and the maximum value of pennies is 500. If we calculate the bonus for that combination, we get an overestimate of all the bonuses within that interval. If that overestimate is less than the bonus of our approximate solution, we can skip searching that interval entirely. The same technique applies to get_max_dimes.

The actual approach used in the code below uses a modified version of this idea: it skips intervals by the largest power-of-2 (get_skip_interval) in the current value of the points left (we loop over points left instead of pennies for reasons that will some become clear). For example, if we are currently at 6 points_left, we will try to skip ahead by 2 points. If that succeeds, we will skip ahead by 8 points. (The same approach is used in get_max_dimes). This way, we reuse the same memoized values of get_max_dimes (or `get_max_dollars) as often as possible.

#include <algorithm>
#include <bit>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
int32_t int_log2(uint64_t x) {
 return 63 - __builtin_clzl(x + (x == 0));
}
std::string to_string(__int128 val) {
 if (val < 0) return "-" + to_string(-val);
 if (val <= UINT64_MAX) {
 return std::to_string(static_cast<uint64_t>(val));
 }
 uint64_t exp10_19 =わ 10000000000000000000ull;
 unsigned __int128 low_val = val % exp10_19;
 unsigned __int128 high_val = val / exp10_19;
 return to_string(high_val) + to_string(low_val);
}
int64_t get_skip_interval(uint64_t points) {
 if (points == 0) return 0;
 return 1 << __builtin_ctz(points);
}
struct alignas(int64_t) Upgrades;
std::ostream& operator<<(std::ostream& os, Upgrades obj);
struct alignas(int64_t) Upgrades {
 int32_t pennies;
 int32_t dimes;
 int32_t dollars;
 int32_t gold;
 friend Upgrades operator+ (Upgrades a, Upgrades b) {
 auto [pennies, dimes, dollars, gold] = a;
 b.pennies += pennies;
 b.dimes += dimes;
 b.dollars += dollars;
 b.gold += gold;
 return b;
 }
 __int128 score() const {
 constexpr __int128 ONE = 1;
 auto ret_val = pennies * (ONE + dimes * (ONE + dollars * (ONE + gold)));
 if (ret_val < 0) {
 std::cout << *this << std::endl;
 throw std::overflow_error("__int128 overflow");
 }
 return ret_val;
 }
 bool is_zero() const {
 return pennies == 0 && dimes == 0 && dollars == 0 && gold == 0;
 }
};
std::ostream& operator<<(std::ostream& os, Upgrades obj) {
 os << '(' 
 << obj.pennies << ", "
 << obj.dimes << ", "
 << obj.dollars << ", "
 << obj.gold << ')';
 return os;
}
template <typename T>
class AutoVector {
 public:
 AutoVector() = default;
 T& get(size_t idx) {
 if (vec.size() <= idx) {
 size_t new_size = (idx + 1) * 2;
 vec.resize(new_size);
 }
 return vec[idx];
 }
 private:
 std::vector<T> vec{};
};
template <auto Fn>
Upgrades memoized(int64_t points) {
 static AutoVector<Upgrades> cache{};
 if (auto value = cache.get(points); !value.is_zero()) {
 return value;
 } else {
 auto result = Fn(points);
 cache.get(points) = result;
 return result;
 }
}
Upgrades get_max_dollars(int64_t points) {
 if (points <= 3) return {0, 0, 0, 0};
 int32_t approx_dollar = (points + 15) / 12;
 __int128 last_max_score = -1;
 Upgrades last_max_combo = {0, 0, 0, 0};
 for (int32_t dollars = approx_dollar - 5; dollars <= approx_dollar + 5; ++dollars) {
 int32_t gold = (points - dollars * 6) / 15;
 auto combo = Upgrades{1, 1, dollars, gold};
 if (dollars >= 0 && gold >= 0 && combo.score() > last_max_score) {
 last_max_combo = combo;
 last_max_score = combo.score();
 }
 }
 last_max_combo.pennies = 0;
 last_max_combo.dimes = 0;
 return last_max_combo;
}
Upgrades get_max_dimes(int64_t points) {
 if (points % 3 != 0) {
 return memoized<get_max_dimes>(points - points % 3);
 }
 
 const int32_t max_dimes = points / 3;
 
 int32_t approx_dimes = max_dimes / 3;
 int32_t approx_points_left = points - approx_dimes * 3;
 auto last_max_combo = memoized<get_max_dollars>(approx_points_left) + Upgrades{1, approx_dimes, 0, 0};
 auto last_max_score = last_max_combo.score();
 
 for(int32_t points_left = points; points_left >= 0; ) {
 for (auto skip_interval = get_skip_interval(points_left); skip_interval >= 2; skip_interval /= 2) {
 int32_t high_points_left = points_left;
 int32_t low_points_left = points_left - skip_interval;
 int32_t high_dimes = (points - low_points_left) / 3 + 1;
 int32_t low_dimes = (points - high_points_left) / 3 - 1;
 if (low_dimes < 0 || low_points_left < 0) continue;
 auto max_combo = memoized<get_max_dollars>(high_points_left) + Upgrades{1, high_dimes, 0, 0};
 if (max_combo.score() < last_max_score) {
 points_left -= skip_interval;
 goto continue_loop;
 }
 }
 {
 const int32_t dimes = (points - points_left) / 3;
 auto combo = memoized<get_max_dollars>(points_left) + Upgrades{1, dimes, 0, 0};
 if (combo.score() > last_max_score) {
 last_max_combo = combo;
 last_max_score = combo.score();
 }
 }
 points_left -= 1;
 continue_loop: ;
 }
 last_max_combo.pennies = 0;
 return last_max_combo;
}
Upgrades get_max_score(int64_t points) {
 int32_t approx_pennies = points / 4;
 auto last_max_combo = memoized<get_max_dimes>(points - approx_pennies) 
 + Upgrades{approx_pennies, 0, 0, 0};
 auto last_max_score = last_max_combo.score();
 
 for(int32_t points_left = points; points_left >= 0; ) {
 for (auto skip_interval = get_skip_interval(points_left); skip_interval >= 2; skip_interval /= 2) {
 int32_t high_points_left = points_left;
 int32_t low_points_left = points_left - skip_interval;
 int32_t high_pennies = points - low_points_left;
 int32_t low_pennies = points - high_points_left;
 if (low_pennies < 0 || low_points_left < 0) continue;
 auto max_combo = memoized<get_max_dimes>(high_points_left) + Upgrades{high_pennies, 0, 0, 0};
 
 if (max_combo.score() < last_max_score) {
 points_left -= skip_interval;
 goto continue_loop;
 }
 }
 {
 const int32_t pennies = points - points_left;
 auto combo = memoized<get_max_dimes>(points_left) + Upgrades{pennies, 0, 0, 0};
 
 if (combo.score() > last_max_score) {
 last_max_combo = combo;
 last_max_score = combo.score();
 }
 }
 --points_left;
 continue_loop: ;
 }
 return last_max_combo;
}
int64_t get_bound(std::string prompt, int64_t default_value) {
 std::cout << prompt;
 std::string result;
 std::getline(std::cin, result);
 try {
 static_assert(sizeof(int64_t) == sizeof(long));
 return std::stol(result);
 } catch (...) {
 return default_value;
 }
}
int main() {
 int64_t lower_bound = get_bound("Lower bound: ", 1);
 int64_t higher_bound = get_bound("Higher bound: ", lower_bound);
 auto begin = std::chrono::steady_clock::now();
 
 std::cout << "Starting precomputation...\n";
 for (int64_t i = 0; i <= higher_bound; i++) {
 get_max_dollars(i);
 }
 auto precomp_end = std::chrono::steady_clock::now();
 std::cout << "Precomputation done in "
 << std::chrono::duration<double>{precomp_end - begin}.count()
 << "s\n";
 
 for (int64_t i = lower_bound; i <= higher_bound; i++) {
 auto optimal_upgrades = get_max_score(i);
 std::cout 
 << i 
 << ": " 
 << optimal_upgrades 
 << ", "
 << to_string(optimal_upgrades.score())
 << '\n';
 }
 
 auto end = std::chrono::steady_clock::now();
 std::cout << "Computation took: " 
 << std::chrono::duration<double>{end - begin}.count()
 << "s" << std::endl;
 return 0;
}

• STerliakov: It's a rare case where I hate constant ONE. Your formula (pennies * (1 + ...))) just happens to include 1 multiple times. It's not some fundamental property (modification replacing 1 with 3 is as likely as changing only one of the coefficients). Unless I severely misunderstand your domain, that would better read as plain 1 without ONE constant at all.
• ayaan098: @STerliakov the constant is used so that all the variables in the equation are promoted to __int128.
• Toby Speight: You have misspelt std::int32_t and other identifiers from <cstdint>.
• indi: In to_string(__int128 val): if (val < 0) return "-" + to_string(-val); Not a good idea. __int128 may not be a standard signed integer, and thus may not necessarily be two’s complement. But it probably is.
• ayaan098: @kiner_shah 1) I've been running the code in an online IDE, so my performance measurements are basically seeing how far the program can run before getting timed-out (there's a timeout of 60sec). 2) Yes, the numbers eventually overflow an int64_t.

Answer 1

I don't quite see how this is [dynamic-programming]; but it is certainly [integer-programming]. Your problem statement is:

Maximize the value of $$P + PL + PLD + PLDG$$ subject to $$P + 3D + 6L + 15G \le N\\ P,D,L,G\ge0$$

(Notice that dollars L cost twice as much as dimes D, but contribute more to the value we're trying to maximize. I wonder if this was a typo on your part. It also doesn't seem to match what you wrote inside Upgrades::score().)

I haven't tried (yet?) to fully understand your solution. Empirically I do see that it gives answers faster than the simplest solution, i.e.

std::tuple<int, int, int, int> result = {0,0,0,0};
int best_score = 0;
N -= 1; // we need to buy at least one penny to score anything at all
for (int G = 0; G < N/15; ++G) {
 for (int L = 0; L < (N - 15*G) / 6; ++L) {
 for (int D = 0; D < (N - 15*G - 6*L) / 3; ++D) {
 int P = (N - 15*G - 6*L - 3*D) + 1;
 int score = P + P*L + P*L*D + P*L*D*G;
 if (score > best_score) {
 result = {G,L,D,P};
 best_score = score;
 }
 }
 }
}

but I don't see any comments about what your algorithm is or how it works. In particular it's not clear what you mean by e.g. "get max dollars" — does that mean "get the maximum number of dollars I can buy with N points," or "get the number of dollars I should buy with N points to maximize my score," or something else? Picking better identifiers might help, but writing down your intent in English will help much more.

• Reinderien: In particular, this is MINLP, which is a problem class much more challenging to solvers than the LP alternative. It wouldn't surprise me if OP's current approach (whatever it is) beats out a naive MINLP solver such as Gekko.
• ayaan098: @Reinderien I've added a detailed explanation of the algorithm to the question. (I haven't added comments though, to keep the code unchanged)