Find and output the given number’s shortest square sum, i.e. the shortest sequence of numbers whose squares summed together give back the input number

An example case: given \( N = 12 \), the SSS (short for shortest squares sum) is

\[SSS = \{ 2, 2, 2 \}\]

since

\[2^2 + 2^2 + 2^2 = N\]

This is a minimum coin change problem in disguise and can be therefore formulated as follows: let \( D(N, S) \) be the minimum number of coins needed among \( N \) ones to get \( S \) by squaring and adding them up. The recursion can be expressed as

\[D(n,s) = min\{ D(n, s - n^2) + 1, \ D(n-1,s) \}\]

and base cases

\[\begin{align} D(n, 0) & = 0 \\ D(0, s) & = \infty \end{align}\]

Code to solve the problem in a dynamic programming fashion and backtrack the solution path follows

// -------------------------------------------------------------------

// Find and output the given number's shortest square sum.

// e.g. for sum == 12, the shortest square sum is made of 3 values only

// i.e. 2^2+2^2+2^2 (not 3^2+1^2+1^2+1^2) so the output is: {2 2 2}  

// -------------------------------------------------------------------


#include <iostream>
#include <vector>
#include <limits>
#include <tuple>
#include <cmath>
#include <algorithm>

// A predecessor node. Identifies a D(n,s) cell and stores

// whether the n-th coin was used in a solution path

struct PredecessorNode {
  PredecessorNode(bool b, int n, int s) :
    taken(b),
    predecessorN(n),
    predecessorS(s)
  {}
  bool taken;
  int predecessorN;
  int predecessorS;
};

int pow(int x, int y) { // Exponentiation by squaring

  if (y == 0)
    return 1;
  else
    return pow(x*x, y >> 1) * ((y % 2 != 0) ? x : 1);
}

bool findSSS(int sum, std::vector<int>& res) {

  // By taking the N squares of only the numbers that will NOT be greater than the sum,

  // i.e. [1;ceil(sqrt(sum))], it's a min-coin-change problem. 

  // This can be solved via dynamic programming. Let D(n, s) be the MINIMUM number

  // of coins needed among n ones we have to make sum s. The recursion is

  //

  //  D(n,s) = min{ D(n, s - i^2) + 1 , D(n-1, s) }

  //                  ^             ^   ^ we don't use the coin here

  //                  |             |

  //                  |             |one coin more

  //                  |

  //                  - not n-1 since once used a coin can be reused

  // 

  // Base cases, or borders to the giant 2D matrix, are

  // D(0, s) = infinity (i.e. can't reach s with no coins, so we provide an infinity

  //                     value which is not desirable)

  // D(n, 0) = 0 (no coins needed to have a 0 sum)


  //               s

  //   -------------------------

  //   |

  // N |

  //   |  

  const int INF = std::numeric_limits<int>::max();
  const int N = static_cast<int>(ceil(sqrt(sum)));
  std::vector<std::vector<int>> dp(N + 1 /* To have it 1-index based */, 
                                   std::vector<int>(sum + 1 /* ditto */, INF));
  std::vector<std::vector<PredecessorNode>> predecessors(N + 1, 
              std::vector<PredecessorNode>(sum + 1, PredecessorNode{ false, 0, -1 }));

  // Base cases

  for (auto& v : dp)
    v[0] = 0;
  // D(0, s) already in place due to initialization



  // D(i,j)

  for (int i = 1; i <= N; ++i) {
    for (int j = 1; j <= sum; ++j) {

      int squared_element = static_cast<int>(pow(i, 2)); // This could be precomputed

      int previous_N = j - squared_element;

      if (previous_N < 0) {
        // This value cannot be used, set the predecessor to D(n-1, s)

        int rv = dp[i - 1][j];
        predecessors[i][j] = PredecessorNode{ false, i - 1, j };
        dp[i][j] = rv;
        continue;
      }

      int lv = dp[i][previous_N] + 1;
      int rv = dp[i - 1][j];

      if (lv < rv) {
        predecessors[i][j] = PredecessorNode{ true, i, previous_N };
        dp[i][j] = lv;
      }
      else {
        predecessors[i][j] = PredecessorNode{ false, i - 1, j };
        dp[i][j] = rv;
      }

    } // j

  } // i


  // After everything has been calculated, check if D(n,s) is not INF (i.e. there's

  // a min sequence of squares)

  if (dp[N][sum] >= INF)
    return false;

  // Reconstruct the shortest squares sequence

  int i = N;
  int j = sum;
  std::vector<int> result;
  while (true) {

    if (predecessors[i][j].taken == true)
      result.push_back(i);

    std::tie(i, j) = std::make_tuple(predecessors[i][j].predecessorN, 
                                     predecessors[i][j].predecessorS);

    if (i <= 0 || j < 0)
      break;
  }

  std::reverse(result.begin(), result.end());
  res = result;
  return true;
}


int main() {

  int sum = 45;
  std::vector<int> result;
  bool found = findSSS(sum, result); // {3, 6}


  if (found) {
    for (auto& v : result)
      std::cout << v << " ";
  }
  std::cout << std::endl;

  return 0;
}

The solution is \( O(N^{\frac{3}{2}}) \) since it features the same square optimization seen in the primes generation article.