# Find the Number of Subsequences With Equal GCD
**Difficulty:** HARD
[External](https://leetcode.com/problems/find-the-number-of-subsequences-with-equal-gcd)
Canonical: https://scaleengineer.com/dsa/problems/find-the-number-of-subsequences-with-equal-gcd
**Patterns:** [Math](https://scaleengineer.com/dsa/patterns/math), [Dynamic Programming](https://scaleengineer.com/dsa/patterns/dynamic-programming), [Number Theory](https://scaleengineer.com/dsa/patterns/number-theory)
**Data structures:** Array
**Companies:** [Infosys](https://scaleengineer.com/companies/infosys)
---
## Problem
You are given an integer array `nums`.

Your task is to find the number of pairs of **non-empty** subsequences `(seq1, seq2)` of `nums` that satisfy the following conditions:

* The subsequences `seq1` and `seq2` are **disjoint**, meaning **no index** of `nums` is common between them.
* The GCD of the elements of `seq1` is equal to the GCD of the elements of `seq2`.

Return the total number of such pairs.

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

**Example 1:**

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

**Output:** 10

**Explanation:**

The subsequence pairs which have the GCD of their elements equal to 1 are:

* `([**1**, 2, 3, 4], [1, **2**, **3**, 4])`
* `([**1**, 2, 3, 4], [1, **2**, **3**, **4**])`
* `([**1**, 2, 3, 4], [1, 2, **3**, **4**])`
* `([**1**, **2**, 3, 4], [1, 2, **3**, **4**])`
* `([**1**, 2, 3, **4**], [1, **2**, **3**, 4])`
* `([1, **2**, **3**, 4], [**1**, 2, 3, 4])`
* `([1, **2**, **3**, 4], [**1**, 2, 3, **4**])`
* `([1, **2**, **3**, **4**], [**1**, 2, 3, 4])`
* `([1, 2, **3**, **4**], [**1**, 2, 3, 4])`
* `([1, 2, **3**, **4**], [**1**, **2**, 3, 4])`

**Example 2:**

**Input:** nums = \[10,20,30\]

**Output:** 2

**Explanation:**

The subsequence pairs which have the GCD of their elements equal to 10 are:

* `([**10**, 20, 30], [10, **20**, **30**])`
* `([10, **20**, **30**], [**10**, 20, 30])`

**Example 3:**

**Input:** nums = \[1,1,1,1\]

**Output:** 50

**Constraints:**

* `1 <= nums.length <= 200`
* `1 <= nums[i] <= 200`

# Approaches
## Dynamic Programming
This approach uses dynamic programming to solve the problem. We build the solution by iterating through each number of the input array `nums` and making a decision for each: either include it in the first subsequence (`seq1`), the second subsequence (`seq2`), or neither. The state of our DP table, `dp[g1][g2]`, keeps track of the number of ways to form two disjoint subsequences with GCDs equal to `g1` and `g2` respectively, using the elements processed so far.
**Time:** O(N * M^2 * log(M)), where N is the length of `nums` and M is the maximum value in `nums`. We iterate through N numbers, and for each, we iterate through M*M states, with each transition taking O(log M) for the GCD calculation. Given the constraints, this is on the edge of passing. · **Space:** O(M^2), where M is the maximum value in `nums`. We need a 2D array of size `(M+1) x (M+1)` for the DP table.
**Pros:** Conceptually straightforward to understand as it directly models the choices for each element.; Guaranteed to find the correct solution.
**Cons:** The time complexity is high due to the three nested loops, making it potentially slow for larger constraints.; The space complexity is quadratic with respect to the maximum value in the input array, which can be memory-intensive.
### Explanation
The DP state `dp[g1][g2]` represents the number of ways to partition the elements considered so far into three groups (for `seq1`, for `seq2`, unused) such that the GCD of elements chosen for `seq1` is `g1` and for `seq2` is `g2`. We initialize `dp[0][0] = 1` to represent two empty sets before we start. Then, for each number `x` from `nums`, we update the DP table. For every existing state `dp[g1][g2]`, we can transition to new states by either not using `x`, adding `x` to the first set, or adding `x` to the second set. After processing all numbers, the answer is the sum of `dp[g][g]` for all `g > 0`.

```java
class Solution {
    private int gcd(int a, int b) {
        while (b != 0) {
            int temp = b;
            b = a % b;
            a = temp;
        }
        return a;
    }

    public int numberOfPairs(int[] nums) {
        int maxVal = 0;
        for (int num : nums) {
            maxVal = Math.max(maxVal, num);
        }
        
        long[][] dp = new long[maxVal + 1][maxVal + 1];
        dp[0][0] = 1;
        long MOD = 1_000_000_007;

        for (int x : nums) {
            long[][] nextDp = new long[maxVal + 1][maxVal + 1];
            for (int g1 = 0; g1 <= maxVal; g1++) {
                for (int g2 = 0; g2 <= maxVal; g2++) {
                    if (dp[g1][g2] == 0) continue;

                    // Case 1: x is not used
                    nextDp[g1][g2] = (nextDp[g1][g2] + dp[g1][g2]) % MOD;

                    // Case 2: x is added to seq1
                    int newG1 = (g1 == 0) ? x : gcd(g1, x);
                    nextDp[newG1][g2] = (nextDp[newG1][g2] + dp[g1][g2]) % MOD;

                    // Case 3: x is added to seq2
                    int newG2 = (g2 == 0) ? x : gcd(g2, x);
                    nextDp[g1][newG2] = (nextDp[g1][newG2] + dp[g1][g2]) % MOD;
                }
            }
            dp = nextDp;
        }

        long ans = 0;
        for (int g = 1; g <= maxVal; g++) {
            ans = (ans + dp[g][g]) % MOD;
        }
        return (int) ans;
    }
}
```
### Algorithm
1.  Determine the maximum value `max_val` in the input array `nums`.
2.  Initialize a 2D DP table `dp[max_val + 1][max_val + 1]` with all values set to zero. `dp[g1][g2]` will store the number of ways to partition the elements processed so far into two disjoint subsequences with GCDs `g1` and `g2`.
3.  Set `dp[0][0] = 1`. This represents the initial state with two empty subsequences. A GCD of 0 is used to denote an empty subsequence.
4.  Iterate through each number `x` in the `nums` array:
    a. Create a new DP table `next_dp` of the same size, initialized to zeros.
    b. Iterate through all states `(g1, g2)` from `0` to `max_val`.
    c. If `dp[g1][g2]` is greater than zero, it means this state is reachable. We consider three choices for the current number `x`:
        i.  **Don't use `x`**: Add `dp[g1][g2]` to `next_dp[g1][g2]`.
        ii. **Add `x` to the first subsequence**: Calculate the new GCD `new_g1 = gcd(g1, x)` (if `g1` is 0, `new_g1` is `x`). Add `dp[g1][g2]` to `next_dp[new_g1][g2]`.
        iii. **Add `x` to the second subsequence**: Calculate the new GCD `new_g2 = gcd(g2, x)` (if `g2` is 0, `new_g2` is `x`). Add `dp[g1][g2]` to `next_dp[g1][new_g2]`.
    d. After processing all states for the current `x`, replace `dp` with `next_dp`.
5.  After iterating through all numbers in `nums`, the final `dp` table is ready.
6.  The total number of valid pairs is the sum of `dp[g][g]` for all `g` from 1 to `max_val`, as we are looking for pairs where the GCDs are equal and non-zero.

## Inclusion-Exclusion Principle with Combinatorics
This highly efficient approach uses a combination of combinatorics and the principle of inclusion-exclusion. The main idea is to change the problem from finding pairs with an exact GCD `g` to counting pairs with a GCD that is a multiple of `g`. The latter is easier to count. Then, using the inclusion-exclusion principle, we can subtract the overcounted cases to find the exact counts.
**Time:** O(M * log(M) + N), where N is the length of `nums` and M is the maximum value. The dominant parts are pre-calculating multiples and the inclusion-exclusion loop, both of which are efficient. · **Space:** O(M), where M is the maximum value in `nums`. We use several arrays of size `M+1`.
**Pros:** Very efficient time complexity, making it suitable for larger constraints.; Low space complexity.
**Cons:** The logic involving combinatorics and the inclusion-exclusion principle can be more complex to understand and implement correctly compared to a direct DP approach.
### Explanation
First, we precompute the number of elements in `nums` that are multiples of each integer `g` up to `max(nums)`. Let this be `multiples[g]`. 

For a given `g`, the number of pairs of non-empty, disjoint subsequences `(seq1, seq2)` where both `gcd(seq1)` and `gcd(seq2)` are multiples of `g` can be calculated. Such subsequences must be formed from the `n_g = multiples[g]` elements. For each of these `n_g` elements, we have 3 choices (in `seq1`, in `seq2`, or neither), leading to `3^{n_g}` pairs of disjoint subsequences (including empty ones). Using combinatorics, the number of pairs where both are non-empty is `(3^{n_g} - 2 * 2^{n_g} + 1)`. Let's call this `f[g]`. 

`f[g]` counts pairs where the common GCD is `g`, `2g`, `3g`, etc. Let `ans[g]` be the count for exactly `g`. Then `f[g] = ans[g] + ans[2g] + ...`. We can solve for `ans[g]` by iterating `g` from `max_val` down to 1 and subtracting the counts for multiples: `ans[g] = f[g] - sum(ans[k*g])` for `k > 1`. The total answer is the sum of all `ans[g]`.

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

    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;
    }

    public int numberOfPairs(int[] nums) {
        int maxVal = 0;
        for (int num : nums) {
            maxVal = Math.max(maxVal, num);
        }

        int[] counts = new int[maxVal + 1];
        for (int num : nums) {
            counts[num]++;
        }

        int[] multiples = new int[maxVal + 1];
        for (int g = 1; g <= maxVal; g++) {
            for (int k = g; k <= maxVal; k += g) {
                multiples[g] += counts[k];
            }
        }

        long[] f = new long[maxVal + 1];
        for (int g = 1; g <= maxVal; g++) {
            int n_g = multiples[g];
            if (n_g == 0) continue;
            long term1 = power(3, n_g);
            long term2 = (2 * power(2, n_g)) % MOD;
            f[g] = (term1 - term2 + 1 + MOD) % MOD;
        }

        long[] ans = new long[maxVal + 1];
        for (int g = maxVal; g >= 1; g--) {
            ans[g] = f[g];
            for (int k = 2 * g; k <= maxVal; k += g) {
                ans[g] = (ans[g] - ans[k] + MOD) % MOD;
            }
        }

        long totalPairs = 0;
        for (int g = 1; g <= maxVal; g++) {
            totalPairs = (totalPairs + ans[g]) % MOD;
        }

        return (int) totalPairs;
    }
}
```
### Algorithm
1.  **Preprocessing**:
    a. Find `max_val`, the maximum element in `nums`.
    b. Count the frequency of each number in `nums` into a `counts` array.
    c. Compute `multiples[g]` for each `g` from 1 to `max_val`. `multiples[g]` is the count of numbers in `nums` that are divisible by `g`. This can be done in `O(max_val * log(max_val))` by iterating through multiples.
2.  **Count Superset Pairs**:
    a. Let `f[g]` be the number of pairs of non-empty, disjoint subsequences `(seq1, seq2)` where `gcd(seq1)` and `gcd(seq2)` are both multiples of `g`.
    b. These subsequences can only be formed from the `n_g = multiples[g]` numbers in `nums` divisible by `g`.
    c. For each of these `n_g` numbers, there are 3 choices: put it in `seq1`, `seq2`, or neither. This gives `3^{n_g}` ways to form two disjoint subsequences. To ensure both are non-empty, we use PIE: `f[g] = (3^{n_g} - 2 * 2^{n_g} + 1) % MOD`.
3.  **Inclusion-Exclusion for Exact Counts**:
    a. Let `ans[g]` be the number of pairs where the common GCD is exactly `g`.
    b. The count `f[g]` includes pairs with common GCD `g, 2g, 3g, ...`. So, `f[g] = ans[g] + ans[2g] + ans[3g] + ...`.
    c. We can find `ans[g]` by iterating `g` from `max_val` down to 1: `ans[g] = f[g] - (ans[2g] + ans[3g] + ...)`. The values `ans[k*g]` for `k>1` are already computed due to the downward iteration.
4.  **Final Result**:
    a. The total number of pairs is the sum of `ans[g]` for all `g` from 1 to `max_val`, as the sets of pairs for different common GCDs are disjoint.