# Find Sum of Array Product of Magical Sequences
**Difficulty:** HARD
[External](https://leetcode.com/problems/find-sum-of-array-product-of-magical-sequences)
Canonical: https://scaleengineer.com/dsa/problems/find-sum-of-array-product-of-magical-sequences
**Patterns:** [Math](https://scaleengineer.com/dsa/patterns/math), [Dynamic Programming](https://scaleengineer.com/dsa/patterns/dynamic-programming), [Bit Manipulation](https://scaleengineer.com/dsa/patterns/bit-manipulation), [Combinatorics](https://scaleengineer.com/dsa/patterns/combinatorics), [Bitmask](https://scaleengineer.com/dsa/patterns/bitmask)
**Data structures:** Array
---
## Problem
You are given two integers, `m` and `k`, and an integer array `nums`.

A sequence of integers `seq` is called **magical** if: 
* `seq` has a size of `m`.
* `0 <= seq[i] < nums.length`
* The **binary representation** of `2seq[0] + 2seq[1] + ... + 2seq[m - 1]` has `k` **set bits**.

The **array product** of this sequence is defined as `prod(seq) = (nums[seq[0]] * nums[seq[1]] * ... * nums[seq[m - 1]])`.

Return the **sum** of the **array products** for all valid **magical** sequences.

Since the answer may be large, return it **modulo** `109 + 7`.

A **set bit** refers to a bit in the binary representation of a number that has a value of 1.

**Example 1:**

**Input:** m = 5, k = 5, nums = \[1,10,100,10000,1000000\]

**Output:** 991600007

**Explanation:**

All permutations of `[0, 1, 2, 3, 4]` are magical sequences, each with an array product of 1013.

**Example 2:**

**Input:** m = 2, k = 2, nums = \[5,4,3,2,1\]

**Output:** 170

**Explanation:**

The magical sequences are `[0, 1]`, `[0, 2]`, `[0, 3]`, `[0, 4]`, `[1, 0]`, `[1, 2]`, `[1, 3]`, `[1, 4]`, `[2, 0]`, `[2, 1]`, `[2, 3]`, `[2, 4]`, `[3, 0]`, `[3, 1]`, `[3, 2]`, `[3, 4]`, `[4, 0]`, `[4, 1]`, `[4, 2]`, and `[4, 3]`.

**Example 3:**

**Input:** m = 1, k = 1, nums = \[28\]

**Output:** 28

**Explanation:**

The only magical sequence is `[0]`.

**Constraints:**

* `1 <= k <= m <= 30`
* `1 <= nums.length <= 50`
* `1 <= nums[i] <= 108`

# Approaches
## Brute-Force Recursion
This approach directly follows the problem definition by generating every possible sequence and checking if it's magical. A sequence `seq` is an array of `m` indices, where each index can be any value from `0` to `nums.length - 1`. There are `(nums.length)^m` such sequences.

We can use a recursive backtracking function to generate all these sequences. For each generated sequence, we verify if it meets the 'magical' criteria. If it does, we compute its array product and add it to a total sum. All calculations involving the sum are performed modulo `10^9 + 7`.
**Time:** O((nums.length)^m * m)
There are `(nums.length)^m` possible sequences. For each sequence, we perform calculations (sum of powers, product) that take O(m) time. · **Space:** O(m)
This is for the recursion stack depth, as the recursion goes `m` levels deep.
**Pros:** Simple to understand and implement.; Directly translates the problem statement into code.
**Cons:** Extremely inefficient due to its exponential time complexity.; Will result in a 'Time Limit Exceeded' (TLE) error for all but the smallest possible inputs.
### Explanation
The algorithm uses a standard recursive template to explore the entire search space of possible sequences. A helper function, say `generate(k, currentSequence)`, is defined, where `k` is the number of elements remaining to be chosen for the sequence.

When `k` is `0`, we have a complete sequence. We then compute `S = 2^seq[0] + ... + 2^seq[m-1]`. A key check is `Long.bitCount(S) == k`. If this holds, the sequence's product `nums[seq[0]] * ... * nums[seq[m-1]]` is calculated and added to the total. The recursion proceeds by iterating through all possible next elements (`0` to `nums.length - 1`) and making a recursive call for `k-1`.

```java
class Solution {
    long totalSum = 0;
    int m, k;
    int[] nums;
    long MOD = 1_000_000_007;

    public int sumOfProducts(int m, int k, int[] nums) {
        this.m = m;
        this.k = k;
        this.nums = nums;
        generate(0, new ArrayList<>());
        return (int) totalSum;
    }

    private void generate(int index, List<Integer> currentSeq) {
        if (index == m) {
            // A full sequence is generated, check if it's magical
            long powerSum = 0;
            for (int val : currentSeq) {
                // This can overflow if val is large, but constraints keep it within long.
                powerSum += 1L << val;
            }

            if (Long.bitCount(powerSum) == k) {
                long product = 1;
                for (int val : currentSeq) {
                    product = (product * nums[val]) % MOD;
                }
                totalSum = (totalSum + product) % MOD;
            }
            return;
        }

        for (int i = 0; i < nums.length; i++) {
            currentSeq.add(i);
            generate(index + 1, currentSeq);
            currentSeq.remove(currentSeq.size() - 1); // Backtrack
        }
    }
}
```
### Algorithm
*   Define a recursive function, say `generate(index, currentSequence)`, to build sequences of length `m`.
*   The `index` parameter tracks the current position in the sequence being built (from `0` to `m-1`).
*   The `currentSequence` is a list or array storing the indices chosen so far.
*   **Base Case:** When `index == m`, a full sequence of length `m` has been generated.
    1.  Calculate the sum of powers of two: `S = 2^seq[0] + 2^seq[1] + ... + 2^seq[m-1]`. Be careful with potential overflows if `m` and `seq[i]` are large, though for the given constraints `m <= 30` and `nums.length <= 50`, `2^49` fits in a 64-bit integer (`long`).
    2.  Count the number of set bits in `S`. In Java, this is `Long.bitCount(S)`.
    3.  If the bit count equals `k`, the sequence is magical. Calculate its array product: `P = nums[seq[0]] * nums[seq[1]] * ... * nums[seq[m-1]]` (modulo `10^9 + 7`).
    4.  Add this product `P` to a running total sum (also modulo `10^9 + 7`).
*   **Recursive Step:** For the current `index`, iterate through all possible values `i` from `0` to `nums.length - 1`.
    1.  Add `i` to the `currentSequence` at the current `index`.
    2.  Make a recursive call: `generate(index + 1, currentSequence)`.
    3.  (If using a mutable list, backtrack by removing the element added in step 1).
*   The initial call to start the process is `generate(0, new empty_sequence)`. The final result is the total sum accumulated.

## Recursion with Count Generation
A more optimized brute-force approach is to recognize that the order of elements in a magical sequence `seq` does not affect the sum of powers of two (`2^a + 2^b = 2^b + 2^a`) or the final array product (`a*b = b*a`). Therefore, many sequences are just permutations of each other and can be grouped.

We can iterate through all unique multisets of `m` indices from `0` to `nums.length - 1`. A multiset can be represented by an array of counts `c`, where `c[i]` is the number of times index `i` appears. The sum of all counts must be `m`.

For each valid count distribution, we check the magical condition. If it holds, we calculate its contribution to the total sum. The number of distinct permutations for a given count distribution `(c_0, ..., c_{n-1})` is given by the multinomial coefficient `m! / (c_0! * c_1! * ... * c_{n-1}!)`. The product for each of these permutations is the same: `nums[0]^c_0 * ... * nums[n-1]^c_{n-1}`. We calculate this total contribution and add it to our answer.
**Time:** O(C(m+n-1, n-1) * n)
Where `n` is `nums.length`. The complexity is dominated by the number of ways to partition `m` into `n` parts, and for each partition, we do O(n) work. · **Space:** O(n)
This is for the recursion stack depth and the `counts` array.
**Pros:** More efficient than the naive brute-force as it groups permutations together.; Introduces the combinatorial insight (multinomial coefficients) that is key to the optimal solution.
**Cons:** The number of ways to choose `m` items from `n` with replacement is given by `C(m+n-1, n-1)`, which is still too large for the given constraints.; Will also result in a 'Time Limit Exceeded' error.
### Explanation
This method uses recursion to find all partitions of the integer `m` into `n` parts, where `n` is `nums.length`. Each part corresponds to the count `c_i`.

```java
// This approach is still too slow and is for conceptual understanding.
// A full implementation would require modular arithmetic helpers for inverse, power, etc.
class Solution {
    long totalSum = 0;
    int m, k, n;
    int[] nums;
    long MOD = 1_000_000_007;
    long[] fact, invFact;

    // Assume modular arithmetic helpers (power, inverse) and precomputed factorials exist.

    public int sumOfProducts(int m, int k, int[] nums) {
        this.m = m;
        this.k = k;
        this.n = nums.length;
        this.nums = nums;
        // Precompute factorials and their modular inverses
        // ... 

        generateCounts(0, m, new int[n]);
        return (int) totalSum;
    }

    private void generateCounts(int numIndex, int remainingM, int[] counts) {
        if (numIndex == n - 1) {
            counts[numIndex] = remainingM;
            
            long powerSum = 0;
            for (int i = 0; i < n; i++) {
                if (counts[i] > 0) {
                    powerSum += (long)counts[i] * (1L << i);
                }
            }

            if (Long.bitCount(powerSum) == k) {
                long productTerm = 1;
                long multinomialCoeff = fact[m];

                for (int i = 0; i < n; i++) {
                    productTerm = (productTerm * power(nums[i], counts[i])) % MOD;
                    multinomialCoeff = (multinomialCoeff * invFact[counts[i]]) % MOD;
                }
                long contribution = (multinomialCoeff * productTerm) % MOD;
                totalSum = (totalSum + contribution) % MOD;
            }
            return;
        }

        for (int c = 0; c <= remainingM; c++) {
            counts[numIndex] = c;
            generateCounts(numIndex + 1, remainingM - c, counts);
        }
    }
    // Helper for modular exponentiation: long power(long base, long exp)
    // Helper for modular inverse: long modInverse(long n)
}
```
### Algorithm
*   Instead of generating sequences, we generate the multiset of chosen indices, represented by their counts `c_0, c_1, ..., c_{n-1}`, where `n = nums.length` and `sum(c_i) = m`.
*   Define a recursive function `generateCounts(numIndex, remainingM, counts)`.
    *   `numIndex`: The index `i` in `nums` for which we are deciding the count `c_i`.
    *   `remainingM`: The number of elements yet to be chosen.
    *   `counts`: An array storing the counts `c_0, ..., c_{numIndex-1}`.
*   **Base Case:** When `numIndex == n` (all counts decided):
    1.  If `remainingM` is not `0`, this path is invalid (it means `sum(c_i) != m`).
    2.  Calculate `S = sum(counts[j] * 2^j for j=0..n-1)`.
    3.  If `Long.bitCount(S) == k`, this multiset is magical.
    4.  Calculate the contribution of all sequences with this multiset of indices: `Term = (m! / (c_0! * ... * c_{n-1}!)) * (nums[0]^c_0 * ... * nums_{n-1}^c_{n-1})`.
    5.  Add `Term` to the total sum (modulo `10^9 + 7`). Precomputing factorials, their inverses, and powers is necessary for this.
*   **Recursive Step:** For the current `numIndex`, iterate through all possible counts `c` from `0` to `remainingM`.
    1.  Set `counts[numIndex] = c`.
    2.  Make a recursive call: `generateCounts(numIndex + 1, remainingM - c, counts)`.
*   The initial call is `generateCounts(0, m, new int[n])`.

## Dynamic Programming
This problem can be solved efficiently using dynamic programming. The structure of the problem, involving choices (counts `c_i`) that build up a sum (`m`) and a bitwise property, is a strong indicator for DP. This approach is analogous to 'Digit DP' but applied to the coefficients `c_i` in the base-2 representation `S = sum(c_i * 2^i)`.
**Time:** O(n * m^3 * k)
We iterate through `n` numbers. For each, we iterate through the DP table states `j` (m+1), `carry` (m+1), `b` (k+1), and for each state, we iterate through `count` (m+1). This gives `n * m * m * k * m`. · **Space:** O(m^2 * k)
We use a DP table of size `(m+1) x (m+1) x (k+1)`. We can optimize this to use two such tables (current and next), but the complexity remains the same.
**Pros:** Provides an efficient solution that passes within the time limits.; Systematically builds the solution by breaking it down into manageable subproblems.; Correctly handles the complex combinatorial and bitwise constraints of the problem.
**Cons:** The DP state and transitions are complex to derive and implement correctly.; Requires careful handling of modular arithmetic, including precomputation of factorials, inverse factorials, and powers.
### Explanation
We define a DP state that captures the necessary information as we iterate through the `nums` array. Let `dp[i][j][carry][b]` be the partial sum of products considering elements `nums[0]` to `nums[i-1]`, having used `j` elements in total, generating a `carry` to the `i`-th bit position, and having `b` set bits in the resulting sum up to bit `i-1`. To save space, we can use a 2D DP table (for `i` and `i-1`) of size `(m+1) x (m+1) x (k+1)`.

```java
class Solution {
    long MOD = 1_000_000_007;

    public int sumOfProducts(int m, int k, int[] nums) {
        int n = nums.length;

        long[] fact = new long[m + 1];
        long[] invFact = new long[m + 1];
        fact[0] = 1;
        invFact[0] = 1;
        for (int i = 1; i <= m; i++) {
            fact[i] = (fact[i - 1] * i) % MOD;
            invFact[i] = power(fact[i], MOD - 2);
        }

        long[][][] dp = new long[m + 1][m + 1][k + 1];
        dp[0][0][0] = 1;

        for (int i = 0; i < n; i++) {
            long[][] powers = new long[m + 1][2];
            for(int j = 0; j <= m; j++) {
                powers[j][0] = power(nums[i], j);
                powers[j][1] = invFact[j];
            }

            long[][][] newDp = new long[m + 1][m + 1][k + 1];
            for (int j = 0; j <= m; j++) { // elements used so far
                for (int carry = 0; carry <= m; carry++) {
                    for (int b = 0; b <= k; b++) {
                        if (dp[j][carry][b] == 0) continue;
                        for (int count = 0; j + count <= m; count++) { // count for nums[i]
                            int newJ = j + count;
                            int val = carry + count;
                            int newCarry = val / 2;
                            int bit = val % 2;
                            int newB = b + bit;

                            if (newB <= k) {
                                long term = (powers[count][0] * powers[count][1]) % MOD;
                                long toAdd = (dp[j][carry][b] * term) % MOD;
                                newDp[newJ][newCarry][newB] = (newDp[newJ][newCarry][newB] + toAdd) % MOD;
                            }
                        }
                    }
                }
            }
            dp = newDp;
        }

        long totalSum = 0;
        for (int carry = 0; carry <= m; carry++) {
            for (int b = 0; b <= k; b++) {
                if (dp[m][carry][b] > 0) {
                    if (b + Integer.bitCount(carry) == k) {
                        totalSum = (totalSum + dp[m][carry][b]) % MOD;
                    }
                }
            }
        }

        return (int) ((totalSum * fact[m]) % MOD);
    }

    private long power(long base, long exp) {
        long res = 1;
        base %= MOD;
        while (exp > 0) {
            if (exp % 2 == 1) res = (res * base) % MOD;
            base = (base * base) % MOD;
            exp /= 2;
        }
        return res;
    }
}
```
### Algorithm
*   The core idea is to build the sum of products iteratively. The expression `sum_{counts} (m! / product(c_i!)) * product(nums[i]^c_i)` can be rewritten as `m! * sum_{counts} product( (nums[i]^c_i) / c_i! )`. We use DP to compute the inner sum.
*   The magical condition `bitCount(sum c_i * 2^i) == k` is handled by processing the sum bit by bit. When we consider `c_i`, we are essentially determining the coefficient of `2^i`. This operation `c_i + carry_from_prev_bit` determines the bit at position `i` and the new carry to position `i+1`. The carry value is bounded by `m`.
*   **DP State:** `dp[j][carry][b]` will store the sum of `product_{l=0..i-1} (nums[l]^c_l / c_l!)` after considering `nums[0]` to `nums[i-1]`, where:
    *   `j`: total number of elements used so far (`sum c_l`).
    *   `carry`: the carry value propagated to bit `i`.
    *   `b`: the number of set bits produced by the sum from bits `0` to `i-1`.
*   **Initialization:** Create a DP table `dp[m+1][m+1][k+1]`. Initialize `dp[0][0][0] = 1`, representing a state before considering any numbers (0 elements used, 0 carry, 0 set bits).
*   **Transition:** Iterate `i` from `0` to `n-1` (for each `nums[i]`). Create a `new_dp` table for the next state.
    *   For each state `(j, carry, b)` in the current `dp` table:
        *   Iterate through `count` (the value of `c_i`) from `0` to `m-j`.
        *   Calculate `new_j = j + count`.
        *   Calculate `val = carry + count`.
        *   `new_carry = val / 2`.
        *   `bit_i = val % 2`.
        *   `new_b = b + bit_i`.
        *   If `new_b <= k`, update the `new_dp` table: `new_dp[new_j][new_carry][new_b] += dp[j][carry][b] * (nums[i]^count / count!)`.
    *   After iterating through all states and counts, replace `dp` with `new_dp`.
*   **Final Calculation:** After the main loop over `i` finishes, we have the final `dp` table. The total number of elements must be `m`. We are interested in states `dp[m][carry][b]`.
    *   Iterate through all final `carry` and `b` values.
    *   If `b + Integer.bitCount(carry) == k`, add `dp[m][carry][b]` to a running sum.
    *   The final answer is this sum multiplied by `m!` (all modulo `10^9 + 7`).