On the Bench (Codeforces 840C)
Tags: DP, combinatorics

A year ago on the bench in public park Leha found an array of n numbers. Leha believes that permutation p is right if for all 1 ≤ i < n condition, that api·api + 1 is not perfect square, holds. Leha wants to find number of right permutations modulo 109 + 7.

Input: First line of input data contains single integer n (1 ≤ n ≤ 300) — length of the array.Next line contains n integers a1, a2, ... , an (1 ≤ ai ≤ 109) — found array.

Output: Output single integer — number of right permutations modulo 109 + 7.

Note: For first example:[1, 2, 4] — right permutation, because 2 and 8 are not perfect squares.[1, 4, 2] — wrong permutation, because 4 is square of 2.[2, 1, 4] — wrong permutation, because 4 is square of 2.[2, 4, 1] — wrong permutation, because 4 is square of 2.[4, 1, 2] — wrong permutation, because 4 is square of 2.[4, 2, 1] — right permutation, because 8 and 2 are not perfect squares.

Sample Inputs:
3
1 2 4

7
5 2 4 2 4 1 1

Sample Outputs:
2

144

Solution (Python):
from collections import defaultdict
import sys
readline = sys.stdin.readline

def Extended_Euclid(n, m):
    stack = []
    while m:
        stack.append((n, m))
        (n, m) = (m, n % m)
    if n >= 0:
        (x, y) = (1, 0)
    else:
        (x, y) = (-1, 0)
    for i in range(len(stack) - 1, -1, -1):
        (n, m) = stack[i]
        (x, y) = (y, x - n // m * y)
    return (x, y)

class MOD:

    def __init__(self, p, e=None):
        self.p = p
        self.e = e
        if self.e == None:
            self.mod = self.p
        else:
            self.mod = self.p ** self.e

    def Pow(self, a, n):
        a %= self.mod
        if n >= 0:
            return pow(a, n, self.mod)
        else:
            x = Extended_Euclid(a, self.mod)[0]
            return pow(x, -n, self.mod)

    def Build_Fact(self, N):
        assert N >= 0
        self.factorial = [1]
        if self.e == None:
            for i in range(1, N + 1):
                self.factorial.append(self.factorial[-1] * i % self.mod)
        else:
            self.cnt = [0] * (N + 1)
            for i in range(1, N + 1):
                self.cnt[i] = self.cnt[i - 1]
                ii = i
                while ii % self.p == 0:
                    ii //= self.p
                    self.cnt[i] += 1
                self.factorial.append(self.factorial[-1] * ii % self.mod)
        self.factorial_inve = [None] * (N + 1)
        self.factorial_inve[-1] = self.Pow(self.factorial[-1], -1)
        for i in range(N - 1, -1, -1):
            ii = i + 1
            while ii % self.p == 0:
                ii //= self.p
            self.factorial_inve[i] = self.factorial_inve[i + 1] * ii % self.mod

    def Fact(self, N):
        if N < 0:
            return 0
        retu = self.factorial[N]
        if self.e != None and self.cnt[N]:
            retu *= pow(self.p, self.cnt[N], self.mod) % self.mod
            retu %= self.mod
        return retu

    def Fact_Inve(self, N):
        if self.e != None and self.cnt[N]:
            return None
        return self.factorial_inve[N]

    def Comb(self, N, K, divisible_count=False):
        if K < 0 or K > N:
            return 0
        retu = self.factorial[N] * self.factorial_inve[K] % self.mod * self.factorial_inve[N - K] % self.mod
        if self.e != None:
            cnt = self.cnt[N] - self.cnt[N - K] - self.cnt[K]
            if divisible_count:
                return (retu, cnt)
            else:
                retu *= pow(self.p, cnt, self.mod)
                retu %= self.mod
        return retu

class Prime:

    def __init__(self, N):
        assert N <= 10 ** 8
        self.smallest_prime_factor = [None] * (N + 1)
        for i in range(2, N + 1, 2):
            self.smallest_prime_factor[i] = 2
        n = int(N ** 0.5) + 1
        for p in range(3, n, 2):
            if self.smallest_prime_factor[p] == None:
                self.smallest_prime_factor[p] = p
                for i in range(p ** 2, N + 1, 2 * p):
                    if self.smallest_prime_factor[i] == None:
                        self.smallest_prime_factor[i] = p
        for p in range(n, N + 1):
            if self.smallest_prime_factor[p] == None:
                self.smallest_prime_factor[p] = p
        self.primes = [p for p in range(N + 1) if p == self.smallest_prime_factor[p]]

    def Factorize(self, N):
        assert N >= 1
        factors = defaultdict(int)
        if N <= len(self.smallest_prime_factor) - 1:
            while N != 1:
                factors[self.smallest_prime_factor[N]] += 1
                N //= self.smallest_prime_factor[N]
        else:
            for p in self.primes:
                while N % p == 0:
                    N //= p
                    factors[p] += 1
                if N < p * p:
                    if N != 1:
                        factors[N] += 1
                    break
                if N <= len(self.smallest_prime_factor) - 1:
                    while N != 1:
                        factors[self.smallest_prime_factor[N]] += 1
                        N //= self.smallest_prime_factor[N]
                    break
            else:
                if N != 1:
                    factors[N] += 1
        return factors

    def Divisors(self, N):
        assert N > 0
        divisors = [1]
        for (p, e) in self.Factorize(N).items():
            A = [1]
            for _ in range(e):
                A.append(A[-1] * p)
            divisors = [i * j for i in divisors for j in A]
        return divisors

    def Is_Prime(self, N):
        return N == self.smallest_prime_factor[N]

    def Totient(self, N):
        for p in self.Factorize(N).keys():
            N *= p - 1
            N //= p
        return N

    def Mebius(self, N):
        fact = self.Factorize(N)
        for e in fact.values():
            if e >= 2:
                return 0
        else:
            if len(fact) % 2 == 0:
                return 1
            else:
                return -1
N = int(readline())
mod = 10 ** 9 + 7
MD = MOD(mod)
MD.Build_Fact(N + 1)
dct = defaultdict(int)
P = Prime(10 ** 5)
for a in map(int, readline().split()):
    for p in P.primes:
        while a % p ** 2 == 0:
            a //= p ** 2
    dct[a] += 1
dp = [0] * N
dp[0] = 1
s = 0
ans = 1
for c in dct.values():
    prev = dp
    dp = [0] * N
    ans *= MD.Fact(c)
    ans %= mod
    for i in range(max(s, 1)):
        for j in range(i + 1):
            for k in range(c + 1):
                if i - j + c - k < N:
                    dp[i - j + c - k] += prev[i] * MD.Comb(c - 1, k - 1) % mod * MD.Comb(s + 1 - i, k - j) % mod * MD.Comb(i, j) % mod
    for i in range(N):
        dp[i] %= mod
    s += c
ans *= dp[0]
ans %= mod
print(ans)

Solution (Python):
from collections import defaultdict
import sys
readline = sys.stdin.readline
from math import gcd as GCD

def Extended_Euclid(n, m):
    stack = []
    while m:
        stack.append((n, m))
        (n, m) = (m, n % m)
    if n >= 0:
        (x, y) = (1, 0)
    else:
        (x, y) = (-1, 0)
    for i in range(len(stack) - 1, -1, -1):
        (n, m) = stack[i]
        (x, y) = (y, x - n // m * y)
    return (x, y)

class MOD:

    def __init__(self, p, e=None):
        self.p = p
        self.e = e
        if self.e == None:
            self.mod = self.p
        else:
            self.mod = self.p ** self.e

    def Pow(self, a, n):
        a %= self.mod
        if n >= 0:
            return pow(a, n, self.mod)
        else:
            assert GCD(a, self.mod) == 1
            x = Extended_Euclid(a, self.mod)[0]
            return pow(x, -n, self.mod)

    def Build_Fact(self, N):
        assert N >= 0
        self.factorial = [1]
        if self.e == None:
            for i in range(1, N + 1):
                self.factorial.append(self.factorial[-1] * i % self.mod)
        else:
            self.cnt = [0] * (N + 1)
            for i in range(1, N + 1):
                self.cnt[i] = self.cnt[i - 1]
                ii = i
                while ii % self.p == 0:
                    ii //= self.p
                    self.cnt[i] += 1
                self.factorial.append(self.factorial[-1] * ii % self.mod)
        self.factorial_inve = [None] * (N + 1)
        self.factorial_inve[-1] = self.Pow(self.factorial[-1], -1)
        for i in range(N - 1, -1, -1):
            ii = i + 1
            while ii % self.p == 0:
                ii //= self.p
            self.factorial_inve[i] = self.factorial_inve[i + 1] * ii % self.mod

    def Fact(self, N):
        if N < 0:
            return 0
        retu = self.factorial[N]
        if self.e != None and self.cnt[N]:
            retu *= pow(self.p, self.cnt[N], self.mod) % self.mod
            retu %= self.mod
        return retu

    def Fact_Inve(self, N):
        if self.e != None and self.cnt[N]:
            return None
        return self.factorial_inve[N]

    def Comb(self, N, K, divisible_count=False):
        if K < 0 or K > N:
            return 0
        retu = self.factorial[N] * self.factorial_inve[K] % self.mod * self.factorial_inve[N - K] % self.mod
        if self.e != None:
            cnt = self.cnt[N] - self.cnt[N - K] - self.cnt[K]
            if divisible_count:
                return (retu, cnt)
            else:
                retu *= pow(self.p, cnt, self.mod)
                retu %= self.mod
        return retu

class Prime:

    def __init__(self, N):
        assert N <= 10 ** 8
        self.smallest_prime_factor = [None] * (N + 1)
        for i in range(2, N + 1, 2):
            self.smallest_prime_factor[i] = 2
        n = int(N ** 0.5) + 1
        for p in range(3, n, 2):
            if self.smallest_prime_factor[p] == None:
                self.smallest_prime_factor[p] = p
                for i in range(p ** 2, N + 1, 2 * p):
                    if self.smallest_prime_factor[i] == None:
                        self.smallest_prime_factor[i] = p
        for p in range(n, N + 1):
            if self.smallest_prime_factor[p] == None:
                self.smallest_prime_factor[p] = p
        self.primes = [p for p in range(N + 1) if p == self.smallest_prime_factor[p]]

    def Factorize(self, N):
        assert N >= 1
        factors = defaultdict(int)
        if N <= len(self.smallest_prime_factor) - 1:
            while N != 1:
                factors[self.smallest_prime_factor[N]] += 1
                N //= self.smallest_prime_factor[N]
        else:
            for p in self.primes:
                while N % p == 0:
                    N //= p
                    factors[p] += 1
                if N < p * p:
                    if N != 1:
                        factors[N] += 1
                    break
                if N <= len(self.smallest_prime_factor) - 1:
                    while N != 1:
                        factors[self.smallest_prime_factor[N]] += 1
                        N //= self.smallest_prime_factor[N]
                    break
            else:
                if N != 1:
                    factors[N] += 1
        return factors

    def Divisors(self, N):
        assert N > 0
        divisors = [1]
        for (p, e) in self.Factorize(N).items():
            A = [1]
            for _ in range(e):
                A.append(A[-1] * p)
            divisors = [i * j for i in divisors for j in A]
        return divisors

    def Is_Prime(self, N):
        return N == self.smallest_prime_factor[N]

    def Totient(self, N):
        for p in self.Factorize(N).keys():
            N *= p - 1
            N //= p
        return N

    def Mebius(self, N):
        fact = self.Factorize(N)
        for e in fact.values():
            if e >= 2:
                return 0
        else:
            if len(fact) % 2 == 0:
                return 1
            else:
                return -1
N = int(readline())
mod = 10 ** 9 + 7
MD = MOD(mod)
MD.Build_Fact(N + 1)
dct = defaultdict(int)
P = Prime(10 ** 5)
for a in map(int, readline().split()):
    for p in P.primes:
        while a % p ** 2 == 0:
            a //= p ** 2
    dct[a] += 1
dp = [0] * N
dp[0] = 1
s = 0
ans = 1
for c in dct.values():
    prev = dp
    dp = [0] * N
    ans *= MD.Fact(c)
    ans %= mod
    for i in range(max(s, 1)):
        for j in range(i + 1):
            for k in range(c + 1):
                if i - j + c - k < N:
                    dp[i - j + c - k] += prev[i] * MD.Comb(c - 1, k - 1) % mod * MD.Comb(s + 1 - i, k - j) % mod * MD.Comb(i, j) % mod
    for i in range(N):
        dp[i] %= mod
    s += c
ans *= dp[0]
ans %= mod
print(ans)

Solution (Python):
from collections import defaultdict
import sys
readline = sys.stdin.readline
from math import gcd as GCD

def Extended_Euclid(n, m):
    stack = []
    while m:
        stack.append((n, m))
        (n, m) = (m, n % m)
    if n >= 0:
        (x, y) = (1, 0)
    else:
        (x, y) = (-1, 0)
    for i in range(len(stack) - 1, -1, -1):
        (n, m) = stack[i]
        (x, y) = (y, x - n // m * y)
    return (x, y)

class MOD:

    def __init__(self, p, e=None):
        self.p = p
        self.e = e
        if self.e == None:
            self.mod = self.p
        else:
            self.mod = self.p ** self.e

    def Pow(self, a, n):
        a %= self.mod
        if n >= 0:
            return pow(a, n, self.mod)
        else:
            assert GCD(a, self.mod) == 1
            x = Extended_Euclid(a, self.mod)[0]
            return pow(x, -n, self.mod)

    def Build_Fact(self, N):
        assert N >= 0
        self.factorial = [1]
        if self.e == None:
            for i in range(1, N + 1):
                self.factorial.append(self.factorial[-1] * i % self.mod)
        else:
            self.cnt = [0] * (N + 1)
            for i in range(1, N + 1):
                self.cnt[i] = self.cnt[i - 1]
                ii = i
                while ii % self.p == 0:
                    ii //= self.p
                    self.cnt[i] += 1
                self.factorial.append(self.factorial[-1] * ii % self.mod)
        self.factorial_inve = [None] * (N + 1)
        self.factorial_inve[-1] = self.Pow(self.factorial[-1], -1)
        for i in range(N - 1, -1, -1):
            ii = i + 1
            while ii % self.p == 0:
                ii //= self.p
            self.factorial_inve[i] = self.factorial_inve[i + 1] * ii % self.mod

    def Fact(self, N):
        if N < 0:
            return 0
        retu = self.factorial[N]
        if self.e != None and self.cnt[N]:
            retu *= pow(self.p, self.cnt[N], self.mod) % self.mod
            retu %= self.mod
        return retu

    def Fact_Inve(self, N):
        if self.e != None and self.cnt[N]:
            return None
        return self.factorial_inve[N]

    def Comb(self, N, K, divisible_count=False):
        if K < 0 or K > N:
            return 0
        retu = self.factorial[N] * self.factorial_inve[K] % self.mod * self.factorial_inve[N - K] % self.mod
        if self.e != None:
            cnt = self.cnt[N] - self.cnt[N - K] - self.cnt[K]
            if divisible_count:
                return (retu, cnt)
            else:
                retu *= pow(self.p, cnt, self.mod)
                retu %= self.mod
        return retu

class Prime:

    def __init__(self, N):
        assert N <= 10 ** 8
        self.smallest_prime_factor = [None] * (N + 1)
        for i in range(2, N + 1, 2):
            self.smallest_prime_factor[i] = 2
        n = int(N ** 0.5) + 1
        for p in range(3, n, 2):
            if self.smallest_prime_factor[p] == None:
                self.smallest_prime_factor[p] = p
                for i in range(p ** 2, N + 1, 2 * p):
                    if self.smallest_prime_factor[i] == None:
                        self.smallest_prime_factor[i] = p
        for p in range(n, N + 1):
            if self.smallest_prime_factor[p] == None:
                self.smallest_prime_factor[p] = p
        self.primes = [p for p in range(N + 1) if p == self.smallest_prime_factor[p]]

    def Factorize(self, N):
        assert N >= 1
        factors = defaultdict(int)
        if N <= len(self.smallest_prime_factor) - 1:
            while N != 1:
                factors[self.smallest_prime_factor[N]] += 1
                N //= self.smallest_prime_factor[N]
        else:
            for p in self.primes:
                while N % p == 0:
                    N //= p
                    factors[p] += 1
                if N < p * p:
                    if N != 1:
                        factors[N] += 1
                    break
                if N <= len(self.smallest_prime_factor) - 1:
                    while N != 1:
                        factors[self.smallest_prime_factor[N]] += 1
                        N //= self.smallest_prime_factor[N]
                    break
            else:
                if N != 1:
                    factors[N] += 1
        return factors

    def Divisors(self, N):
        assert N > 0
        divisors = [1]
        for (p, e) in self.Factorize(N).items():
            A = [1]
            for _ in range(e):
                A.append(A[-1] * p)
            divisors = [i * j for i in divisors for j in A]
        return divisors

    def Is_Prime(self, N):
        return N == self.smallest_prime_factor[N]

    def Totient(self, N):
        for p in self.Factorize(N).keys():
            N *= p - 1
            N //= p
        return N

    def Mebius(self, N):
        fact = self.Factorize(N)
        for e in fact.values():
            if e >= 2:
                return 0
        else:
            if len(fact) % 2 == 0:
                return 1
            else:
                return -1
N = int(readline())
mod = 10 ** 9 + 7
MD = MOD(mod)
MD.Build_Fact(N + 1)
dct = defaultdict(int)
P = Prime(10 ** 5)
for a in map(int, readline().split()):
    for p in P.primes:
        while a % p ** 2 == 0:
            a //= p ** 2
    dct[a] += 1
dp = [0] * N
dp[0] = 1
s = 0
ans = 1
for c in dct.values():
    prev = dp
    dp = [0] * N
    ans *= MD.Fact(c)
    ans %= mod
    for i in range(max(s, 1)):
        for j in range(i + 1):
            for k in range(c + 1):
                if i - j + c - k < N:
                    dp[i - j + c - k] += prev[i] * MD.Comb(c - 1, k - 1) % mod * MD.Comb(s + 1 - i, k - j) % mod * MD.Comb(i, j) % mod
                    dp[i - j + c - k] %= mod
    s += c
ans *= dp[0]
ans %= mod
print(ans)