LeetCode Daily Problem P3539 Sum of Array Products of Magic Sequences

---

Timeline

Timeline

2025-10-12

init

A combinatorial counting problem that uses DP to simulate binary addition carry rules, advanced dynamic programming,__builtin_popcount

---

Problem:

P3539 Sum of Array Products of Magic Sequences

https://leetcode.cn/problems/find-sum-of-array-product-of-magical-sequences/?envType=daily-question&envId=2025-10-12

Attempt exhaustive enumeration

At first I only thought of exhaustive enumeration, but I didn’t even enumerate completely, frustrating.

#include <algorithm>
#include <stack>
#include <vector>

using std::stack;
using std::vector;

const int MOD = 1000000007; // 10^9 + 7

class Solution {
    public:
        bool magicalSeq(vector<int> &seq, int m, int k, int length)
        {
                // Magic sequence: the sequence length of seq is m.
                // 0 <= seq[i] < nums.length
                // The binary form of 2^seq[0] + 2^seq[1] + ... + 2^seq[m - 1] has m set bits.
                // That is, in {seq[0], seq[1] .. seq[M-1]}, if two numbers are the same,
                // Recursively merge the same numbers into one with the original value +1, until there are no identical numbers.
                int val;
                int set = 0;
                if (seq.size() != m) {
                        return false;
                }

                std::sort(seq.begin(), seq.end());

                if (seq[0] < 0 || seq[m - 1] >= length) {
                        return false;
                }

                stack<int, std::vector<int> > st(seq);

                while (!st.empty()) {
                        val = st.top();
                        st.pop();
                        set++;
                        if (!st.empty() && val == st.top()) {
                                st.top()++;
                                set--;
                        }
                }

                if (set == k) {
                        return true;
                } else {
                        return false;
                }
        }

        int magicalSum(int m, int k, vector<int> &nums)
        {
                int i, j;
                int n = nums.size();
                unsigned long magical_sum = 0;
                unsigned long array_mul;

                vector<vector<int> > result;
                int res_size;

                if (m <= n) {
                        // Generate all subsequences of nums of length m
                        // predict array, similar to SVE instructions
                        vector<int> predict(n, 0);
                        // Start from the smallest lexicographic order
                        std::fill(predict.end() - m, predict.end(), 1);

                        do {
                                vector<int> tmp;
                                for (i = 0; i < n; i++) {
                                        if (predict[i] != 0) {
                                                tmp.push_back(i);
                                        }
                                }

                                result.push_back(tmp);
                        } while (std::next_permutation(predict.begin(), predict.end()));
                }

                res_size = result.size();
                // Determine if it is a magic sequence
                for (i = 0; i < res_size; i++) {
                        if (magicalSeq(result[i], m, k, n)) {
                                // If it is a magic sequence, calculate the array product of all its permutations
                                // But if the numbers are the same, different permutations have the same array product
                                // Therefore directly multiply by the factorial of m
                                array_mul = 1;
                                for (j = 0; j < m; j++) {
                                        array_mul = (array_mul * nums[result[i][j]]) % MOD;
                                }
                                for (j = 1; j <= m; j++) {
                                        array_mul = (array_mul * j) % MOD;
                                }
                                magical_sum = (magical_sum + array_mul) % MOD;
                        }
                }
                return magical_sum;
        }
};

Dynamic programming

The solution to this problem still uses dynamic programming, but it is high-order dynamic programming.

Problem

Given two integersMandK, and an integer arraynums。

An integer sequenceseqif it satisfies the following conditions, is calledmagicSequence:

• seqhas sequence length ofM。
• 0 <= seq[i] < nums.length
• 2^seq[0] + 2^seq[1] + ... + 2^seq[M - 1]ofbinary formhasKaset bits。

This sequence’sarray productis defined asprod(seq) = (nums[seq[0]] * nums[seq[1]] * ... * nums[seq[M - 1]])。

Return all validmagicsequence’sarray productofsum 。

Since the answer may be large, return the result modulo10^9 + 7modulo。

set bitrefers to a bit that is 1 in the binary representation of a number.

What we need to do is:

• Given an arraynums, lengthn
• Take a sequence of lengthMsequenceseq(elements from nums can be chosen repeatedly)
• Condition:2^seq[0] + ... + 2^seq[M-1]has in binaryKones
• Array product:prod(seq) = nums[seq[0]] * ... * nums[seq[M-1]]
• Goal: sum of products of all magic sequence arrays (mod 1e9+7)

Number of sequence permutations

Suppose we start fromnumsindex0..n-1take numbers:

• Each numberitookr_ipieces
• Total satisfies:r_0 + r_1 + ... + r_{n-1} = M

Number of permutations of length M:

$$\frac{M!}{r_0! \, r_1! \, ... \, r_{n-1}!}$$

Formula for number of permutations of length M

Number of permutations without considering duplicate numbers

• Ifall numbers are distinct, then the number of sequence permutations of length M isM!
• Because M distinct elements can have M! permutations

Consider repeated digits

• If some digits appear repeatedly, for example:

$$seq = [a, a, b, b, b, c]$$

• Swapping repeated digits does not produce a new sequence
• In combinatorics, we need todivide by the factorial of each repeated digit, to eliminate the effect of repeated permutations
• Suppose digit i appearsr_itimes, the total permutation formula is:

$$\text{排列数} = \frac{M!}{r_0! \, r_1! \, \dots \, r_{n-1}!}$$

This is the classicmultiset permutation formula: number of permutations of length M = factorial of total elements / factorial of each repeated element

Array product

$$\prod_{i=0}^{n-1} nums[i]^{r_i} \cdot \frac{M!}{r_0! \, r_1! \, ... \, r_{n-1}!}$$

Array product formula

The array product of a sequence is defined as:

$$\text{prod(seq)} = nums[seq[0]] \cdot nums[seq[1]] \cdot \dots \cdot nums[seq[M-1]]$$

• If the number i appears in the sequencer_itimes, then the product can be written as:

$$\prod_{i=0}^{n-1} nums[i]^{r_i}$$

• Combining the number of permutations, a certain sequence { r0, r_1, r_2, … ,r } the total contribution of the sequence is:

$$\text{总贡献} = \frac{M!}{r_0! \, r_1! \, ... \, r_{n-1}!} \cdot \prod_{i=0}^{n-1} nums[i]^{r_i}$$

• Explanation:
• The array product corresponding to each sequence of length M is∏ nums[i]^r_i
• There are a total ofM! / ∏ r_i!permutations (because permutations of repeated numbers are not considered new sequences)
• Somultiplying them gives the total contribution of that combination

Dynamic programming approach

Use DP to simulate the carry rule of binary addition

• Enumerate all possibler_i(the number of occurrences of each number) is very complex
• Dynamic programming approach:Consider each number in nums one by one

State definition

• f[i][j][p][q]represents:
  • Already considered the firstinumbers (nums[0…i])
  • Total takenj
  • The sum of the lower i bits of the current mask (carry state) isp
  • The number of set bits in the lower i bits isq
  • Value = sum of all corresponding sequences∏ r_t! * nums[t]^r_twhich is the sum of contributions of all these cases

Transition formula

When considering thei+1th number:

• Assume we chooseranums[i+1]
• New status:

$$f[i+1][j+r][\lfloor p/2 \rfloor + r][q + (p \bmod 2)] += f[i][j][p][q] * nums[i+1]^r * r!$$

• p/2andp%2the role of:
  • Simulate the binary of mask from low to high bits
  • p % 2= number of low bits set
  • p / 2= mask after shifting up one bit

Initialization

• Fori=0, that is the first number:

$$f[0][j][j][0] = nums[0]^j * j!$$

• Low 0 position bits = 0

---

Result

• Total set bits of mask = set bits of high mask + low q

if (__builtin_popcount(p) + q == K)
    res += f[n-1][M][p][q] * M!

Code implementation

The code uses many techniques, here we first introduce the following

Fast exponentiation

Core idea

Using the exponent’sbinary representationandproperty of splitting powers：

$$x^y = x^{(b_0·2^0 + b_1·2^1 + ... + b_k·2^k)} = \prod_{i | b_i=1} x^{2^i}$$

That is:

• Takeyas binary;
• Each time square the base, which is equivalent to doubling the exponent;
• If the current binary digit is 1, multiply this term into the result

Example

Suppose we want to calculate:

$$3^{13}$$

The binary of 13 is:1101₂ = 8 + 4 + 0 + 1

So:

$$3^{13} = 3^{8} × 3^{4} × 3^{1}$$

We can obtain these powers by repeatedly squaring:

Step	Power	Result
3¹	3
3²	9
3⁴	81
3⁸	6561

Therefore:

$$3^{13} = 3⁸ × 3⁴ × 3¹ = 6561 × 81 × 3 = 1594323$$

Corresponding algorithm logic

long long quickPow(long long x, long long y) {
    long long res = 1;
    long long cur = x;  // Current base

    while (y > 0) {
        if (y & 1) {
            res *= cur;  // Current bit is 1 → multiply the current power into the result
        }
        y >>= 1;         // Shift right by one, equivalent to dividing by 2
        cur *= cur;      // Square the base → double the exponent
    }
    return res;
}

Fermat’s Little Theorem

According toFermat’s Little Theorem (Fermat’s Little Theorem)：

$$a^{p-1} ≡ 1 \pmod{p}$$

Therefore:

$$a^{p-2} ≡ a^{-1} \pmod{p}$$

Fermat’s Little Theorem requires that p must be a prime number

That is to say:

The modular inverse of a number = this number raised to the power (mod - 2) modulo mod.

Final code implementation

#include <vector>
using std::vector;

class Solution {
    public:
        long long quickmul(long long x, long long y, long long mod)
        {
                long long res = 1, cur = x % mod;
                while (y) {
                        if (y & 1) {
                                res = res * cur % mod;
                        }
                        y >>= 1;
                        cur = cur * cur % mod;
                }
                return res;
        };

        int magicalSum(int m, int k, vector<int> &nums)
        {
                int n = nums.size();
                const long long mod = 1e9 + 7;

                // Calculate factorials from 1 to m
                vector<long long> fac(m + 1, 1);
                for (int i = 1; i <= m; i++) {
                        fac[i] = fac[i - 1] * i % mod;
                }
                // ifac[i]=(i!)^(−1) mod mod
                vector<long long> ifac(m + 1, 1);
                for (int i = 2; i <= m; i++) {
                        // Use Fermat's little theorem to compute the modular inverse of each i
                        ifac[i] = quickmul(i, mod - 2, mod);
                }
                // Use recurrence relation to convert the inverse of a single integer into factorial inverses
                for (int i = 2; i <= m; i++) {
                        ifac[i] = ifac[i - 1] * ifac[i] % mod;
                }

                // numsPower[i][j] represents nums[i] raised to the power j modulo mod
                // If i appears j times in the sequence, then the product of the sequence is, iterating over all sequences i,
                // multiply cumulatively numsPower[i][j]
                vector numsPower(n, vector<long long>(m + 1, 1));
                for (int i = 0; i < n; i++) {
                        for (int j = 1; j <= m; j++) {
                                numsPower[i][j] = numsPower[i][j - 1] * nums[i] % mod;
                        }
                }

                vector f(n, vector(m + 1, vector(m * 2 + 1, vector<long long>(k + 1, 0))));

                for (int j = 0; j <= m; j++) {
                        f[0][j][j][0] = numsPower[0][j] * ifac[j] % mod;
                }

                for (int i = 0; i + 1 < n; i++) {
                        for (int j = 0; j <= m; j++) {
                                for (int p = 0; p <= m * 2; p++) {
                                        for (int q = 0; q <= k; q++) {
                                                int q2 = p % 2 + q;
                                                if (q2 > k) {
                                                        break;
                                                }
                                                for (int r = 0; r + j <= m; r++) {
                                                        int p2 = p / 2 + r;
                                                        f[i + 1][j + r][p2][q2] +=
                                                                f[i][j][p][q] *
                                                                numsPower[i + 1][r] % mod *
                                                                ifac[r] % mod;
                                                        f[i + 1][j + r][p2][q2] %= mod;
                                                }
                                        }
                                }
                        }
                }
                long long res = 0;
                for (int p = 0; p <= m * 2; p++) {
                        for (int q = 0; q <= k; q++) {
                                if (__builtin_popcount(p) + q == k) {
                                        res = (res + f[n - 1][m][p][q] * fac[m] % mod) % mod;
                                }
                        }
                }
                return res;
        }
};