Precompute a divisibility matrix `div` where `div[i][j] = 1` if `nums[i] % nums[j] == 0` or `nums[j] % nums[i] == 0`, else `0`. This helps quickly check valid transitions during DP.

Initialize the DP table with base cases: for each element `i`, set `dp[1 << i][i] = 1` since a single-element permutation is trivially special.

For each mask from `1` to `2^n - 1`, and for each `last` in the mask, iterate over all `next` not in the mask. If `div[last][next]` is true, update `dp[mask | (1 << next)][next] += dp[mask][last]` modulo `10^9 + 7`. The answer is the sum of all `dp[(1 << n) - 1][i]` for all `i`.