Problem:
 Score : 600 points    Problem Statement  Snuke has a blackboard and a set S consisting of N integers.
The i -th element in S is S_i .  He wrote an integer X on the blackboard, then performed the following operation N times:   Choose one element from S and remove it.  Let x be the number written on the blackboard now, and y be the integer removed from S . Replace the number on the blackboard with x \bmod {y} .   There are N! possible orders in which the elements are removed from S .
For each of them, find the number that would be written on the blackboard after the N operations, and compute the sum of all those N! numbers modulo 10^{9}+7 .      Constraints  All values in input are integers.  1 \leq N \leq 200  1 \leq S_i, X \leq 10^{5}  S_i are pairwise distinct.         Input Input is given from Standard Input in the following format:  N  X  S_1  S_2  \ldots  S_{N}       Output Print the answer.        Sample Input 1 2 19
3 7      Sample Output 1 3   There are two possible orders in which we remove the numbers from S .  If we remove 3 and 7 in this order, the number on the blackboard changes as follows: 19 \rightarrow 1 \rightarrow 1 .  If we remove 7 and 3 in this order, the number on the blackboard changes as follows: 19 \rightarrow 5 \rightarrow 2 .  The output should be the sum of these: 3 .        Sample Input 2 5 82
22 11 6 5 13      Sample Output 2 288       Sample Input 3 10 100000
50000 50001 50002 50003 50004 50005 50006 50007 50008 50009      Sample Output 3 279669259   Be sure to compute the sum modulo 10^{9}+7 .

Solution (Python):
N,X = map(int,input().split())
res = 0
mod = 10** 9 + 7
S = list(map(int,input().split()))
S.sort()
dp = [[0 for i in range(10**5+1)] for j in range(201)]

dp[0][X] = X
for i in range(10**5+1):
    dp[1][i] = i%S[0]
for i in range(2,N+1):
    for j in range(10**5+1):
        dp[i][j] = (dp[i-1][j%S[i-1]] + (i-1)* dp[i-1][j] ) %mod

ans = dp[N][X]
print(ans)

Solution (Python):
def examA():
    ABC =LI(); ABC.sort()
    if ABC[0]==ABC[-1]:
        print("Yes")
    else:
        print("No")
    return

def examB():
    N = I()
    S = SI()
    red = S.count("R")
    if red>N//2:
        print("Yes")
    else:
        print("No")
    return

def examC():
    N,Q = LI()
    S = SI()
    T = [LSI()for _ in range(Q)]
    l = -1; r = N
    while(r-l>1):
        cur = (r+l)//2
        now = copy.deepcopy(cur)
        flag = False
        for t,d in T:
            if S[now]==t:
                if d=="R":
                    now += 1
                else:
                    now -= 1
            if now==N:
                flag = True
                break
            if now==-1:
                break
        if flag:
            r = cur
        else:
            l = cur
    R = r
    l = -1; r = N
    while(r-l>1):
        cur = (r+l)//2
        now = copy.deepcopy(cur)
        flag = False
        for t,d in T:
            if S[now]==t:
                if d=="R":
                    now += 1
                else:
                    now -= 1
            if now==-1:
                flag = True
                break
            if now==N:
                break
        if flag:
            l = cur
        else:
            r = cur
    L = l
    ans = max(0,R-L-1)
    print(ans)
    return

def examD():
    N, X = LI()
    S = LI()
    S.sort(reverse=True)
    dp = [[0]*(X+1)for _ in range(N+1)]
    dp[0][X] = 1
    for i in range(N):
        s = S[i]
        for j in range(X+1):
            dp[i+1][j%s] += dp[i][j]
            dp[i+1][j%s] %= mod
        for j in range(X+1):
            dp[i+1][j] += dp[i][j]*(N-i-1)
            dp[i+1][j] %= mod
    
    ans = 0
    for i,d in enumerate(dp[-1]):
        ans += d*i
        ans %= mod
    print(ans)
    return

def examE():
    ans = 0
    print(ans)
    return

def examF():
    ans = 0
    print(ans)
    return

import sys,copy,bisect,itertools,heapq,math
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
alphabet = [chr(ord('a') + i) for i in range(26)]

if __name__ == '__main__':
    examD()

"""

"""

Solution (Python):
n,x = map(int,input().split())
s = list(map(int,input().split()))
mod = 10**9 + 7

s.sort(reverse=True)
dp = [[0] * (10**5+1) for _ in range(n+1)]

dp[0][x] = 1
for i in range(n):
    for j in range(10**5+1):
        dp[i+1][j] += dp[i][j] * (n-1-i)
        dp[i+1][j%s[i]] += dp[i][j]
        dp[i+1][j] %= mod
        dp[i+1][j%s[i]] %= mod

ans = 0
for i in range(10**5+1):
    ans += i*dp[-1][i]
    ans %= mod

print(ans)

Solution (Python):
import sys
input = sys.stdin.readline

N, X = map(int, input().split())
S = list(map(int, input().split()))
mod = 10**9 + 7
S.sort(reverse=True)
dp = [0] * (X+1)
dp[X] = 1
for i in range(N):
    ndp = dp[:]
    for j in range(X+1):
        
        ndp[j] = ndp[j] * (N - i - 1) % mod
        
        ndp[j % S[i]] = (ndp[j % S[i]] + dp[j]) % mod
    dp = ndp
ans = 0
for j in range(X+1):
    ans = (ans + dp[j] * j) % mod
print(ans)

Solution (Python):
N,X=map(int,input().split())
S=list(map(int,input().split()))
S.sort(reverse=True)
DP=[[0]*(X+1) for i in range(N+1)]
mod=10**9+7
inv=[0,1]
for i in range(2,N+10):
  inv.append((-inv[mod%i]*(mod//i))%mod)
DP[0][X]=1
for i in range(N):
  DP[0][X]=DP[0][X]*(i+1)%mod
for i in range(N):
  for j in range(X+1):
    DP[i+1][j]=(DP[i+1][j]+DP[i][j]*(1-inv[N-i]))%mod
    DP[i+1][j%S[i]]=(DP[i+1][j%S[i]]+DP[i][j]*inv[N-i])%mod
P=0
for i in range(X+1):
  P=(P+DP[N][i]*i)%mod
print(P)

Solution (Python):
import sys
input = sys.stdin.readline
mod = 10**9+7

def main():
    N, X = map(int, input().split())
    A = list(map(int, input().split()))
    A.sort()
    L = max(A[-1], X)

    dp = [[0]*N for _ in range(L+1)]
    a = A[0]
    for x in range(L+1):
        dp[x][0] = x%a

    dp1 = [0 for _ in range(L+1)] 
    for n in range(1, N):
        a = A[n]
        for x in range(L+1):
            dp1[x] = ((n-1) * dp1[x] + dp[x][n-1]) % mod
            dp[x][n] = dp1[x%a]

    ans = 0
    t = 1
    for n in reversed(range(N)):
        ans = (ans + dp[X][n]*t) % mod
        t = t * n % mod

    print(ans)


if __name__ == "__main__":
    main()

Solution (Python):
mod = 10**9+7
n, X = map(int, input().split())
s = [int(x) for x in input().split()]
s = sorted(s, reverse=True)

dp = [[0 for _ in range(X+1)] for _ in range(n)]
dp[-1][X] = 1

for i in range(n):
    for x in range(X+1):
        dp[i][x % s[i]] += dp[i-1][x]
        dp[i][x % s[i]] %= mod
        dp[i][x] += dp[i-1][x]*(n-i-1)
        dp[i][x] %= mod

dp[-1][X] -= 1
ans = 0

for x in range(X+1):
    ans += dp[n-1][x]*x
    ans %= mod

print(ans)

Solution (Python):
mod = 10**9+7
n, X = map(int, input().split())
s = [int(x) for x in input().split()]
s = sorted(s)

dp = [[0 for _ in range(X+1)] for _ in range(n+1)]

for i in range(n+1):
    if i == 0:
        for x in range(X+1):
            dp[i][x] = x
    else:
        for x in range(X+1):
            dp[i][x] = dp[i-1][x % s[i-1]] + dp[i-1][x]*(i-1)
            dp[i][x] %= mod

ans = dp[n][X]

print(ans)

Solution (Python):
mod = 10**9+7
n, X = map(int, input().split())
s = [int(x) for x in input().split()]
s = sorted(s)

dp = [[0 for _ in range(X+1)] for _ in range(n+1)]
inv = [pow(i, mod-2, mod) for i in range(n+1)]

for i in range(n+1):
    if i == 0:
        for x in range(X+1):
            dp[i][x] = x
    else:
        for x in range(X+1):
            dp[i][x] = dp[i-1][x % s[i-1]] * inv[i] + dp[i-1][x]*(1-inv[i])
            dp[i][x] %= mod

ans = dp[n][X]

for i in range(n):
    ans *= i+1
    ans %= mod

print(ans)

Solution (Python):
from collections import defaultdict
mod = 10 ** 9 + 7

N, X = map(int, input().split())
S = list(map(int, input().split()))
S.sort()

dp = defaultdict(int)
dp[X] = 1

for i in reversed(range(N)):
    a = S[i]
    newDP = defaultdict(int)
    for key, val in dp.items():
        newDP[key] += val * i
        newDP[key] %= mod
        newDP[key % a] += val
        newDP[key % a] %= mod
    dp = newDP

ans = 0
for k, v in dp.items():
    ans = (ans + k * v) % mod
print(ans)

Solution (Python):
import sys,bisect,string,math,time,functools,random,fractions
from heapq import heappush,heappop,heapify
from collections import deque,defaultdict,Counter
from itertools import permutations,combinations,groupby
def Golf():n,*t=map(int,open(0).read().split())
def I():return int(input())
def S_():return input()
def IS():return input().split()
def LS():return [i for i in input().split()]
def LI():return [int(i) for i in input().split()]
def LI_():return [int(i)-1 for i in input().split()]
def NI(n):return [int(input()) for i in range(n)]
def NI_(n):return [int(input())-1 for i in range(n)]
def StoLI():return [ord(i)-97 for i in input()]
def ItoS(n):return chr(n+97)
def LtoS(ls):return ''.join([chr(i+97) for i in ls])
def GI(V,E,ls=None,Directed=False,index=1):
    org_inp=[];g=[[] for i in range(V)]
    FromStdin=True if ls==None else False
    for i in range(E):
        if FromStdin:
            inp=LI()
            org_inp.append(inp)
        else:
            inp=ls[i]
        if len(inp)==2:
            a,b=inp;c=1
        else:
            a,b,c=inp
        if index==1:a-=1;b-=1
        aa=(a,c);bb=(b,c);g[a].append(bb)
        if not Directed:g[b].append(aa)
    return g,org_inp
def GGI(h,w,search=None,replacement_of_found='.',mp_def={'

    mp=[boundary]*(w+2);found={}
    for i in range(h):
        s=input()
        for char in search:
            if char in s:
                found[char]=((i+1)*(w+2)+s.index(char)+1)
                mp_def[char]=mp_def[replacement_of_found]
        mp+=[boundary]+[mp_def[j] for j in s]+[boundary]
    mp+=[boundary]*(w+2)
    return h+2,w+2,mp,found
def TI(n):return GI(n,n-1)
def bit_combination(n,base=2):
    rt=[]
    for tb in range(base**n):s=[tb//(base**bt)%base for bt in range(n)];rt+=[s]
    return rt
def gcd(x,y):
    if y==0:return x
    if x%y==0:return y
    while x%y!=0:x,y=y,x%y
    return y
def show(*inp,end='\n'):
    if show_flg:print(*inp,end=end)
 
YN=['YES','NO'];Yn=['Yes','No']
mo=10**9+7
inf=float('inf')
FourNb=[(-1,0),(1,0),(0,1),(0,-1)];EightNb=[(-1,0),(1,0),(0,1),(0,-1),(1,1),(-1,-1),(1,-1),(-1,1)];compas=dict(zip('WENS',FourNb));cursol=dict(zip('LRUD',FourNb))
l_alp=string.ascii_lowercase

read=sys.stdin.buffer.read
readline=sys.stdin.buffer.readline
input=lambda: sys.stdin.readline().rstrip()

class Tree:
    def __init__(self,inp_size=None,ls=None,init=True,index=0):
        self.LCA_init_stat=False
        self.ETtable=[]
        if init:
            if ls==None:
                self.stdin(inp_size,index=index)
            else:
                self.size=len(ls)+1
                self.edges,_=GI(self.size,self.size-1,ls,index=index)
        return
 
    def stdin(self,inp_size=None,index=1):
        if inp_size==None:
            self.size=int(input())
        else:
            self.size=inp_size
        self.edges,_=GI(self.size,self.size-1,index=index)
        return
    
    def listin(self,ls,index=0):
        self.size=len(ls)+1
        self.edges,_=GI(self.size,self.size-1,ls,index=index)
        return
 
    def __str__(self):
        return  str(self.edges)
 
    def dfs(self,x,func=lambda prv,nx,dist:prv+dist,root_v=0):
        q=deque()
        q.append(x)
        v=[-1]*self.size
        v[x]=root_v
        while q:
            c=q.pop()
            for nb,d in self.edges[c]:
                if v[nb]==-1:
                    q.append(nb)
                    v[nb]=func(v[c],nb,d)
        return v
 
    def EulerTour(self,x):
        q=deque()
        q.append(x)
        self.depth=[None]*self.size
        self.depth[x]=0
        self.ETtable=[]
        self.ETdepth=[]
        self.ETin=[-1]*self.size
        self.ETout=[-1]*self.size
        cnt=0
        while q:
            c=q.pop()
            if c<0:
                ce=~c
            else:
                ce=c
                for nb,d in self.edges[ce]:
                    if self.depth[nb]==None:
                        q.append(~ce)
                        q.append(nb)
                        self.depth[nb]=self.depth[ce]+1
            self.ETtable.append(ce)
            self.ETdepth.append(self.depth[ce])
            if self.ETin[ce]==-1:
                self.ETin[ce]=cnt
            else:
                self.ETout[ce]=cnt
            cnt+=1
        return
    
    def LCA_init(self,root):
        self.EulerTour(root)
        self.st=SparseTable(self.ETdepth,init_func=min,init_idl=inf)
        self.LCA_init_stat=True
        return
    
    def LCA(self,root,x,y):
        if self.LCA_init_stat==False:
            self.LCA_init(root)
        xin,xout=self.ETin[x],self.ETout[x]
        yin,yout=self.ETin[y],self.ETout[y]
        a=min(xin,yin)
        b=max(xout,yout,xin,yin)
        id_of_min_dep_in_et=self.st.query_id(a,b+1)
        return self.ETtable[id_of_min_dep_in_et]
 
class SparseTable: 
    def __init__(self,ls,init_func=min,init_idl=float('inf')):
        self.func=init_func
        self.idl=init_idl
        self.size=len(ls)
        self.N0=self.size.bit_length()
        self.table=[ls[:]]
        self.index=[list(range(self.size))]
        self.lg=[0]*(self.size+1)
        
        for i in range(2,self.size+1):
            self.lg[i]=self.lg[i>>1]+1
 
        for i in range(self.N0):
            tmp=[self.func(self.table[i][j],self.table[i][min(j+(1<<i),self.size-1)]) for j in range(self.size)]
            tmp_id=[self.index[i][j] if self.table[i][j]==self.func(self.table[i][j],self.table[i][min(j+(1<<i),self.size-1)]) else self.index[i][min(j+(1<<i),self.size-1)] for j in range(self.size)]
            self.table+=[tmp]
            self.index+=[tmp_id]
    
    
    def query(self,l,r):
        if r>self.size:r=self.size
        
        N=self.lg[r-l]
        return self.func(self.table[N][l],self.table[N][max(0,r-(1<<N))])
    
    
    def query_id(self,l,r):
        if r>self.size:r=self.size
        
        N=self.lg[r-l]
        a,b=self.index[N][l],self.index[N][max(0,r-(1<<N))]
        if self.table[0][a]==self.func(self.table[N][l],self.table[N][max(0,r-(1<<N))]):
            b=a
        return b
    
    def __str__(self):
        return str(self.table[0])
 
    def print(self):
        for i in self.table:
            print(*i)

show_flg=False
show_flg=True
ans=0


n,x=LI()
s=LI()
s.sort(reverse=True)
X=x

dp=[[0]*(X+1) for i in range(n+1)]
dp[0][X]=1

for i in range(n):
    for j in range(X+1):
        dp[i+1][j%s[i]]+=dp[i][j]%mo
        dp[i+1][j]+=dp[i][j]*(n-i-1)%mo

ans=sum([(i)*dp[n][i] for i in range(X+1)])
print(ans%mo)

Solution (Python):
N, X = map(int, input().split())
S = list(map(int, input().split()))
S.sort(reverse=True)
dp = [[0]*(X+1) for _ in range(N+1)]
dp[0][X] = 1
MOD = 10**9+7

for i in range(N):
    for j in range(X+1):
        dp[i+1][j%S[i]] += dp[i][j]
        dp[i+1][j%S[i]] %= MOD
        dp[i+1][j] += (N-i-1)*dp[i][j]
        dp[i+1][j] %= MOD

ans = 0

for i in range(X+1):
    ans += dp[N][i]*i
    ans %= MOD

print(ans)

Solution (Python):
N,X=map(int,input().split())
S=list(map(int,input().split()))
S.sort()
mod=10**9+7
dp=[[0 for i in range(10**5+1)] for j in range(N)]

for i in range(10**5+1):
    dp[0][i]=i%S[0]

for i in range(1,N):
    for j in range(10**5+1):
        dp[i][j]=(i*dp[i-1][j]+dp[i-1][j%S[i]])%mod

print(dp[N-1][X])

Solution (Python):
from bisect import bisect_right
import sys

sys.setrecursionlimit(10 ** 6)
int1 = lambda x: int(x) - 1
p2D = lambda x: print(*x, sep="\n")
def II(): return int(sys.stdin.readline())
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(rows_number): return [LI() for _ in range(rows_number)]
def SI(): return sys.stdin.readline()[:-1]

class mint:
    def __init__(self, x):
        self.__x = x % md

    def __repr__(self):
        return str(self.__x)

    def __neg__(self):
        return mint(-self.__x)

    def __add__(self, other):
        if isinstance(other, mint): other = other.__x
        return mint(self.__x + other)

    def __sub__(self, other):
        if isinstance(other, mint): other = other.__x
        return mint(self.__x - other)

    def __rsub__(self, other):
        return mint(other - self.__x)

    def __mul__(self, other):
        if isinstance(other, mint): other = other.__x
        return mint(self.__x * other)

    __radd__ = __add__
    __rmul__ = __mul__

    def __truediv__(self, other):
        if isinstance(other, mint): other = other.__x
        return mint(self.__x * pow(other, md - 2, md))

    def __rtruediv__(self, other):
        return mint(other * pow(self.__x, md - 2, md))

    def __pow__(self, power, modulo=None):
        return mint(pow(self.__x, power, md))

md = 10**9+7

def main():
    n,x=MI()
    ss=LI()
    ss.sort()
    dp=[0]*(x+1)
    dp[x]=mint(1)
    for i in range(x,0,-1):
        pre=dp[i]
        if pre==0:continue
        mxj=bisect_right(ss,i)
        for j in range(mxj):
            dp[i%ss[j]]+=pre/mxj

    ev=sum(dp[i]*i for i in range(ss[0]))
    fac=mint(1)
    for i in range(1,n+1):fac*=i
    print(ev*fac)

main()

Solution (Python):
N,X=map(int,input().split())
S=list(map(int,input().split()))
pr=10**9+7
b={X:1}
S.sort(reverse=1)
def ir(k,v):
    global d
    if k in d:
        d[k]=(d[k]+v)%pr
    else:
        d[k]=v
for i in range(N):
    t=S[i]
    n=N-i-1
    d={}
    for k,v in b.items():
        ir(k%t,v)
        ir(k,v*n)
    b=d
a=0
for k,v in b.items():
    a=(a+k*v)%pr
print(a)

Solution (Python):
n, x = map(int, input().split())
s = list(map(int, input().split()))
MOD = 10 ** 9 + 7

s = sorted(s, reverse=True)
max_ = max(s[0], x)
dp = [[0] * (max_ + 1) for i in range(n + 1)]
dp[0][x] = 1

for i in range(n):
    for val in range(max_ + 1):
        dp[i + 1][val % s[i]] += dp[i][val]
        dp[i + 1][val % s[i]] %= MOD
        dp[i + 1][val] += dp[i][val] * (n - i - 1)
        dp[i + 1][val] %= MOD

ans = 0
for i, num in enumerate(dp[-1]):
    ans += num * i
    ans %= MOD
print(ans)

Solution (Python):

from collections import defaultdict, deque
from heapq import heappush, heappop
from itertools import permutations, accumulate
import sys
import math
import bisect
def LI(): return [int(x) for x in sys.stdin.buffer.readline().split()]
def I(): return int(sys.stdin.buffer.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
def IR(n):
    return [I() for i in range(n)]
def LIR(n):
    return [LI() for i in range(n)]
def SR(n):
    return [S() for i in range(n)]
def LSR(n):
    return [LS() for i in range(n)]

sys.setrecursionlimit(1000000)
mod = 1000000007

def solve():
    n,x = LI()
    s = LI()
    s.sort()
    dp = [[0]*(n+1) for i in range(x+1)]
    for i in range(n):
        si = s[i]
        ni = i+1
        for j in range(x+1):
            dp[j][0] = j
            nd = (dp[j%si][i]+i*dp[j][i])%mod
            dp[j][ni] += nd
    print(dp[x][n]%mod)
    return


if __name__ == "__main__":
    solve()

Solution (Python):
import sys
input = sys.stdin.readline
N, X = map(int, input().split())
a = list(map(int, input().split()))
mod = 10 ** 9 + 7
a.sort()
dp = [[0] * (X + 1) for _ in range(N + 1)]
for j in range(X + 1): dp[0][j] = j
for i in range(N):
  for j in range(X + 1):
    dp[i + 1][j] = dp[i][j % a[i]] + dp[i][j] * i
    dp[i + 1][j] %= mod
print(dp[-1][X])

Solution (Python):
n,x=map(int,input().split());s,f,M=sorted(list(map(int,input().split())))[::-1],[1,1],10**9+7
for i in range(2,n):f+=[f[-1]*i%M]
D={x:f[n-1]}
for i, k in enumerate(s):
 m,inv=dict(),pow(n-i-1,M-2,M)
 for j,a in D.items():m[j%k]=(m.get(j%k,0)+a)
 for j,a in m.items():D[j]=(D.get(j,0)+a*inv)%M
print(sum([i*j for i, j in m.items()])%M)

Solution (Python):
def main():
    mod = 10**9+7
    fact = [1, 1]
    for i in range(2, 300):
        fact.append(fact[-1]*i % mod)

    n, x = map(int, input().split())
    s = sorted(list(map(int, input().split())))[::-1]
    next_dict = {x: fact[n-1]}
    for i in range(n):
        memo = dict()
        inv = pow(n-i-1, mod-2, mod)
        ss = s[i]
        for j, a in next_dict.items():
            memo[j % ss] = (memo.get(j % ss, 0)+a) % mod
        for j, a in memo.items():
            next_dict[j] = (next_dict.get(j, 0)+a*inv) % mod
    ans = 0
    for i, j in memo.items():
        ans = (ans+i*j) % mod
    print(ans)


main()

Solution (Python):
def main():
    mod = 10**9+7
    fact = [1, 1]
    for i in range(2, 300):
        fact.append(fact[-1]*i % mod)

    n, x = map(int, input().split())
    s = sorted(list(map(int, input().split())))[::-1]
    next_dict = {x: fact[n-1]}
    for i in range(n):
        memo = dict()
        inv = pow(n-i-1, mod-2, mod)
        ss = s[i]
        for j, a in next_dict.items():
            memo[j % ss] = (memo.get(j % ss, 0)+a) % mod
        for j, a in memo.items():
            next_dict[j] = (next_dict.get(j, 0)+a*inv) % mod
    ans = 0
    for i, j in memo.items():
        ans = (ans+i*j) % mod
    print(ans)


main()

Solution (Python):
from collections import defaultdict
n,x = map(int,input().split())
a = list(map(int,input().split()))
a.sort(reverse=True)
mod = 10**9+7
dp = [defaultdict(int) for i in range(n+1)]
dp[0][x] = 1
for i in range(1,n+1):
  for j in dp[i-1].keys():
    k = j%a[i-1]
    dp[i][j] = (dp[i][j]+dp[i-1][j]*(n-i))%mod
    dp[i][k] = (dp[i][k]+dp[i-1][j])%mod
ans = 0
for i,x in dp[n].items():
  ans = (ans+i*x)%mod
print(ans)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
from bisect import bisect_right, bisect_left
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, gamma, log
from operator import mul
from functools import reduce


sys.setrecursionlimit(2147483647)
INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 10 ** 9 + 7


n, x = LI()
S = sorted(LI())
dp = [[0] * (x + 1) for _ in range(n + 1)]
for i in range(x + 1):
    dp[0][i] = i



for j in range(n):
    for k in range(x + 1):
        dp[j + 1][k] = (dp[j][k % S[j]] + dp[j][k] * j) % mod


print(dp[n][x])

Solution (Python):

"""
Writer: SPD_9X2
https://atcoder.jp/contests/exawizards2019/tasks/exawizards2019_d

結局はS中で最小の物以下になる
また、それ以下になったらもう変化しない

通り数を数えなくてはダメ
→0以上min(S)未満の整数kに関して、最後にkになっているのは何通りあるかを数える？

大きい方からdp

具体的には、小さい方で割ったのちに大きい方で割るのは意味がないので大きい方から見ていく
dp[i][n] = 大きい方からi番目までの中で黒板に書かれてる数字がnとなる順序数

推移は、普通に割ったあまりのとこに推移させる & 自分が最初に来た後　
O(max(X,S)*N)

→解けて無いなぁ(途中に来た場合をどうやって数え上げるんですか…)
計算量にNが2回かかってしまう…

影響のあるやつだけ考える？
大きい方から、変化のあるやつだけ取っていく。min(S)以下になったら終了。変化のない奴の位置の組み合わせは…？

N枚のカードをN個あるスロットルに入れ、右側から黒板に向けて射出していくことを考える
dp[i][j] = i枚目までスロットルに入れた時、今射出して黒板の数字がjになる置き方の組み合わせ数
とすれば、Sより大きいjに関しては割ればよいし、そうでなければそれ以降の好きなところにおいてよいので…？
→いやなんか違うな

=========解説を見た=========
https://drken1215.hatenablog.com/entry/2019/04/01/201600
dp[i][x]=
スタート時点での整数がxであるときに、最初のi個についての順列すべてを考えたときに得られる値の総和

うーん総和を直接dpしちゃうのかぁ…
小さい方から見た時に、最初以外に入れても変化しないので
dp[i+1][x] = i * dp[i][x] + dp[i][x % s[i+1]]

→小さい方から見ればよかったかぁ…あと、xをすべてdpに乗せて調べ上げるという発想はなかった…
（無意識のうちに無理だろと思ってしまっていた）


"""
mod = 10**9+7
N,X = map(int,input().split())

S = list(map(int,input().split()))
S.sort()

dp = [[ i for i in range(X+1) ] for j in range(N+1)]

for i in range(N):

    for x in range(X+1):

        dp[i+1][x] = ((i) * dp[i][x] + dp[i][x % S[i]]) % mod


print (dp[N][X])

Solution (Python):
def main():
    import sys
    input = sys.stdin.readline

    N, x = map(int, input().split())
    S = list(map(int, input().split()))
    S.sort(reverse=True)
    mod = 1000000007

    ans = [0] * (10**5+1)
    ans[x] = 1
    for i, s in enumerate(S):
        for j in range(10**5+1):
            aj = ans[j]
            ans[j%s] = (ans[j%s] + aj)%mod
            ans[j] = (ans[j] + aj * (N-i-2))%mod
    ret = 0
    for k, a in enumerate(ans):
        ret = (ret + (a*k)%mod)%mod
    print(ret)



if __name__ == '__main__':
    main()

Solution (Python):
from collections import defaultdict

MOD = 10**9+7
n, x = map(int, input().split())
S = sorted(map(int, input().split()), reverse=True)
dp = [defaultdict(int) for _ in range(n+1)]
dp[0][x] = 1
for i in range(n):
  now = S[i]
  for j in dp[i]:
    dp[i+1][j%now] += dp[i][j]
    dp[i+1][j] += dp[i][j] * (n-1-i)
    dp[i+1][j%now] %= MOD
    dp[i+1][j] %= MOD
ans = 0
for k, v in dp[n].items():
  ans += k*v
  ans %= MOD
print(ans)

Solution (Python):
from math import factorial
from collections import defaultdict

MOD = 10**9+7
n, x = map(int, input().split())
S = sorted(map(int, input().split()), reverse=True)
dp = [defaultdict(int) for _ in range(n+1)]
dp[0][x] = 1
for i in range(n):
  now = S[i]
  p = pow(n-i, MOD-2, MOD)
  q = p*(n-1-i)%MOD
  for j in dp[i]:
    dp[i+1][j%now] += dp[i][j] * p
    dp[i+1][j] += dp[i][j] * q
    dp[i+1][j%now] %= MOD
    dp[i+1][j] %= MOD
ans = 0
for k, v in dp[n].items():
  ans += k*v
  ans %= MOD
ans *= (factorial(n) % MOD)
ans %= MOD
print(ans)

Solution (Python):
N, X = map(int, input().split())
S = [int(i) for i in input().split()]

mod = 10 ** 9 + 7

S.sort(reverse=True)

dp = [0] * (X + 1)
dp[X] = 1

for i in range(N) :
  ndp = [0] * (X + 1)
  for j in range(X + 1) :
    ndp[j%S[i]] += dp[j]
    ndp[j%S[i]] %= mod
    if i != N - 1 :
      ndp[j] += dp[j] * (N - 1 - i)
      ndp[j] %= mod
  dp = ndp

ret = 0
for i in range(X + 1) :
  ret += i * dp[i] % mod
  ret %= mod

print(ret)

Solution (Python):


n,x = map(int,input().split())
da = list(map(int,input().split()))
m = 100010
mod = 1000000007

dp = [0 for i in range(m)]
da.sort()
da.reverse()
dp[x] = 1
for i in range(n):
    for j in range(m):
        if j >= da[i]:
            dp[j%da[i]] += dp[j]
            dp[j] *= n-i-1
            if dp[j] >= mod: dp[j] %= mod
            if dp[j%da[i]] >= mod: dp[j] %= mod
        else:
            dp[j] *= n-i
            if dp[j] >= mod: dp[j]%=mod
ans = 0
for i in range(m):
    ans += dp[i] * i
    if ans >= mod: ans %= mod
print(ans)

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
def IR(n):
    return [I() for i in range(n)]
def LIR(n):
    return [LI() for i in range(n)]
def SR(n):
    return [S() for i in range(n)]
def LSR(n):
    return [LS() for i in range(n)]

sys.setrecursionlimit(1000000)
mod = 1000000007



def solve():
    n,x = LI()
    s = LI()
    s.sort()
    dp = [[0]*(x+1) for i in range(n+1)] 
    for j in range(x+1):
        dp[0][j] = j 
        for i in range(n):
            si = s[i]
            ni = i+1
            nd = dp[i][j%si] 
            nd += i*dp[i][j] 
            dp[ni][j] += nd
            if dp[ni][j] >= mod:
                dp[ni][j] %= mod
    print(dp[n][x])
    return


def solve_1():
    n,x = LI()
    s = LI()
    s.sort()
    
    dp = [j for j in range(x+1)]
    for i in range(n):
        si = s[i]
        ndp = [0]*(x+1)
        for j in range(x+1):
            nd = dp[j%si] 
            nd += i*dp[j] 
            ndp[j] += nd
            if ndp[j] >= mod:
                ndp[j] %= mod
        dp = [j for j in ndp]
    print(dp[x])
    return


if __name__ == "__main__":
    solve_1()

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
def IR(n):
    return [I() for i in range(n)]
def LIR(n):
    return [LI() for i in range(n)]
def SR(n):
    return [S() for i in range(n)]
def LSR(n):
    return [LS() for i in range(n)]

sys.setrecursionlimit(1000000)
mod = 1000000007



def solve():
    n,x = LI()
    s = LI()
    s.sort()
    dp = [[0]*(x+1) for i in range(n+1)] 
    for j in range(x+1):
        dp[0][j] = j 
        for i in range(n):
            si = s[i]
            ni = i+1
            nd = dp[i][j%si] 
            nd += i*dp[i][j] 
            dp[ni][j] += nd
            if dp[ni][j] >= mod:
                dp[ni][j] %= mod
    print(dp[n][x])
    return


if __name__ == "__main__":
    solve()

Solution (Python):
def solve():
  n,x = (int(i) for i in input().split())
  s = sorted(list(int(i) for i in input().split()))
  dp = [[0]*(x+1) for _ in range(n+2)]
  for i in range(x):
    dp[0][i] = i
  mod = 10**9+7
  for i in range(n):
    for j in range(x+1):
      dp[i+1][j] += dp[i][j%s[i]]
      dp[i+1][j] %= mod
      dp[i+1][j] += dp[i][j]*i
      dp[i+1][j] %= mod
  print(dp[n][x])
  
solve()

Solution (Python):
import sys
import math
input = sys.stdin.readline
MOD = 10**9 + 7

n, k = [int(item) for item in input().split()]
s = [int(item) for item in input().split()]
s.sort(reverse=True)
MAX_VAL = max(max(s), k)

dp = [[0] * (MAX_VAL + 1) for _ in range(n + 1)]
dp[0][k] = 1
for i, item in enumerate(s):
    for j in range(MAX_VAL + 1):
        if dp[i][j] == 0:
            continue
        
        dp[i+1][j % item] += dp[i][j]
        dp[i+1][j % item] %= MOD
        
        dp[i+1][j] += dp[i][j] * (n - 1 - i)
        dp[i+1][j] %= MOD

ans = 0
for i, item in enumerate(dp[n]):
    ans += i * item
    ans %= MOD
print(ans)

Solution (Python):


import sys
from functools import lru_cache
sys.setrecursionlimit(500000)

def make_modinv_list(n, mod=10**9+7):
    
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

def main():
    modinv = make_modinv_list(200)

    @lru_cache(maxsize=None)
    def solve(idx, X):
        
        if idx == N:
            return X
        
        res = 0
        
        res += modinv[N-idx] * solve(idx+1, X % S[idx]) % mod
        
        res += (N-idx-1) * modinv[N-idx] * solve(idx+1, X) % mod
        return res % mod

    mod = 10**9+7
    N, X = map(int, input().split())
    S = sorted(map(int, input().split()), reverse=True)
    ans = solve(0, X)
    for i in range(2, N+1):
        ans = ans * i % mod
    print(ans)

main()

Solution (Python):
mod = 10**9+7
N, X = map(int, input().split())

S = sorted(list(map(int, input().split())))
S_max = 10**5
dp = [[0]*(S_max+1) for _ in range(N+1)]

for x in range(S_max+1):
    dp[0][x] = x
for i in range(1, N+1):
    for x in range(S_max+1):
        dp[i][x] = (dp[i-1][x]*(i-1)+dp[i-1][x%S[i-1]])%mod
print(dp[N][X]%mod)