Problem:
 Score : 2718 points    Problem Statement There is a very long bench.
The bench is divided into M sections, where M is a very large integer.  Initially, the bench is vacant.
Then, M people come to the bench one by one, and perform the following action:   We call a section comfortable if the section is currently unoccupied and is not adjacent to any occupied sections.
If there is no comfortable section, the person leaves the bench.
Otherwise, the person chooses one of comfortable sections uniformly at random, and sits there.
(The choices are independent from each other).   After all M people perform actions, Snuke chooses an interval of N consecutive sections uniformly at random (from M-N+1 possible intervals), and takes a photo.
His photo can be described by a string of length N consisting of X and - : the i -th character of the string is X if the i -th section from the left in the interval is occupied, and - otherwise.
Note that the photo is directed.
For example, -X--X and X--X- are different photos.  What is the probability that the photo matches a given string s ?
This probability depends on M .
You need to compute the limit of this probability when M goes infinity.  Here, we can prove that the limit can be uniquely written in the following format using three rational numbers p, q, r and e = 2.718 \ldots (the base of natural logarithm):  p + \frac{q}{e} + \frac{r}{e^2}  Your task is to compute these three rational numbers, and print them modulo 10^9 + 7 , as described in Notes section.      Notes When you print a rational number, first write it as a fraction \frac{y}{x} , where x, y are integers and x is not divisible by 10^9 + 7 (under the constraints of the problem, such representation is always possible).
Then, you need to print the only integer z between 0 and 10^9 + 6 , inclusive, that satisfies xz \equiv y \pmod{10^9 + 7} .      Constraints  1 \leq N \leq 1000  |s| = N  s consists of X and - .         Input Input is given from Standard Input in the following format:  N  s       Output Print three rational numbers p, q, r , separated by spaces.        Sample Input 1 1
X      Sample Output 1 500000004 0 500000003  The probability that a randomly chosen section is occupied converge to \frac{1}{2} - \frac{1}{2e^2} .       Sample Input 2 3
---      Sample Output 2 0 0 0  After the actions, no three consecutive unoccupied sections can be left.       Sample Input 3 5
X--X-      Sample Output 3 0 0 1  The limit is \frac{1}{e^2} .       Sample Input 4 5
X-X-X      Sample Output 4 500000004 0 833333337  The limit is \frac{1}{2} - \frac{13}{6e^2} .       Sample Input 5 20
-X--X--X-X--X--X-X-X      Sample Output 5 0 0 183703705  The limit is \frac{7}{675e^2} .       Sample Input 6 100
X-X-X-X-X-X-X-X-X-X--X-X-X-X-X-X-X-X-X-X-X-X-X-X-X--X--X-X-X-X--X--X-X-X--X-X-X--X-X--X--X-X--X-X-X-      Sample Output 6 0 0 435664291

Solution (Python):
import sys
import numpy as np

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 1_000_000_007

lim = 1 << 10
fact = np.ones(lim, np.int64)
for n in range(1, lim):
    fact[n] = fact[n - 1] * n % MOD
fact_inv = np.ones(lim, np.int64)
fact_inv[-1] = pow(int(fact[-1]), MOD - 2, MOD)
for n in range(lim - 1, 0, -1):
    fact_inv[n - 1] = n * fact_inv[n] % MOD
inv = np.zeros(lim, np.int64)
inv[1:] = fact[:-1] * fact_inv[1:] % MOD

def integrate_poly_0x(F):
    F1 = np.zeros_like(F)
    F1[1:] = F[:-1] * inv[1:len(F)] % MOD
    return F1

def integrate_poly_x1(F):
    F1 = -integrate_poly_0x(F)
    F1[0] -= F1.sum() % MOD
    return F1 % MOD

def integrate_0x(a, F, G):
    """\int_0^x (ae^{-x} + F(x) + G(x)/e)"""
    F1 = integrate_poly_0x(F)
    G1 = integrate_poly_0x(G)
    F1[0] += a
    return -a, F1, G1

def integrate_x1(a, F, G):
    """\int_x^1 (ae^{-x} + F(x) + G(x)/e)"""
    F1 = integrate_poly_x1(F)
    G1 = integrate_poly_x1(G)
    G1[0] -= a
    return a, F1, G1

def integrate_01(a, F, G):
    a, F, G = integrate_x1(a, F, G)
    F[1:] = 0
    G[1:] = 0
    F[0] += a
    return 0, F, G

def one_to_zero(dp0, dp1, dp2):
    dp0 = integrate_0x(*dp2)
    dp1 = integrate_x1(*dp2)
    dp2 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
    return dp0, dp1, dp2

def zero_to_one(dp0, dp1, dp2):
    a0, F0, G0 = integrate_01(*dp0)
    a1, F1, G1 = integrate_x1(*dp1)
    dp0 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
    dp1 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
    dp2 = [a0 + a1, F0 + F1, G0 + G1]
    return dp0, dp1, dp2

def zero_to_zero(dp0, dp1, dp2):
    a0, F0, G0 = integrate_01(*dp0)
    dp0 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
    dp1 = [a0, F0, G0]
    dp2 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
    return dp0, dp1, dp2

def integrate_poly_01_with_exp(F):
    
    a, b = 0, 0
    F = F.copy()
    while len(F) > 0:
        a += F[0]
        b -= F.sum() % MOD
        F = F[1:] * np.arange(1, len(F)) % MOD
    return a, b

def integrate_01_with_exp(a, F, G):
    
    ans = np.zeros(3, np.int64)
    if a & 1:
        a += MOD
    ans[0] += a // 2
    ans[2] -= a // 2
    a, b = integrate_poly_01_with_exp(F)
    ans[0] += a
    ans[1] += b
    a, b = integrate_poly_01_with_exp(G)
    ans[1] += a
    ans[2] += b
    return ans % MOD

def get_ans(dp0, dp1, dp2):
    ans = np.zeros(3, np.int64)

    a, F, G = integrate_01(*dp0)
    ans[0] += a + F[0]
    ans[1] += G[0]
    a, F, G = integrate_01(*dp1)
    ans[0] += a + F[0]
    ans[1] += G[0]
    ans -= integrate_01_with_exp(*dp1)
    ans += integrate_01_with_exp(*dp2)
    ans %= MOD
    return ans

N = int(readline())
S = [1 if x == 'X' else 0 for x in readline().rstrip().decode()]
N += 10

dp0 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
dp1 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]
dp2 = [0, np.zeros(N, np.int64), np.zeros(N, np.int64)]

if S[0] == 1:
    dp2[0] = 1
else:
    dp0[0] = -1
    dp0[1][0] = 1
    dp1[0] = 1
for s, t in zip(S, S[1:]):
    if (s, t) == (1, 1):
        print(0, 0, 0)
        exit()
    elif (s, t) == (1, 0):
        dp0, dp1, dp2 = one_to_zero(dp0, dp1, dp2)
    elif (s, t) == (0, 1):
        dp0, dp1, dp2 = zero_to_one(dp0, dp1, dp2)
    elif (s, t) == (0, 0):
        dp0, dp1, dp2 = zero_to_zero(dp0, dp1, dp2)
print(*get_ans(dp0, dp1, dp2))

Solution (Python):
from copy import *
N=int(input())
S=input()
DP=[[[[0]*2 for i in range(2)] for j in range(2)] for k in range(N)]
mod=10**9+7
for i in range(N):
  for j in range(2):
    for k in range(2):
      for l in range(2):
        DP[i][j][k][l]=[[],[],0]
if S[0]=='X':
  DP[0][0][1][0]=[[],[],1]
else:
  DP[0][0][0][1]=[[],[],1]
  DP[0][1][1][1]=[[1],[],mod-1]
g1=[1,1]
g2=[1,1]
inv=[0,1]
c=[1,2]
for i in range(2,N+5):
  g1.append((g1[-1]*i)%mod)
  inv.append((-inv[mod%i]*(mod//i))%mod)
  g2.append((g2[-1]*inv[-1])%mod)
  c.append((c[-1]+g2[i])%mod)

def add(x,y):
  for j in range(2):
    for i in range(len(y[j])-len(x[j])):
      x[j].append(0)
    for i in range(len(y[j])):
      x[j][i]=(x[j][i]+y[j][i])%mod
  x[2]=(x[2]+y[2])%mod

def integrate(p):
  p[2]=-p[2]%mod
  for j in range(2):
    p[j].append(0)
    for i in range(len(p[j])-2,-1,-1):
      p[j][i+1]=p[j][i]*inv[i+1]%mod
      p[j][i]=0

def fx0(p):
  integrate(p)
  p[0][0]=(p[0][0]-p[2])%mod

def f1x(p):
  integrate(p)
  for j in range(2):
    for i in range(len(p[j])):
      p[j][i]=-p[j][i]%mod
  p[2]=-p[2]%mod
  for i in range(2):
    p[i][0]=-sum(p[i])%mod
  p[1][0]=(p[1][0]-p[2])%mod

for i in range(N-1):
  if S[i+1]=='X':
    add(DP[i][1][1][0],DP[i][1][1][1])
    fx0(DP[i][1][1][0])
    add(DP[i][1][1][1],DP[i][1][0][1])
    add(DP[i][1][1][1],DP[i][0][1][1])
    add(DP[i][1][1][1],DP[i][0][0][1])
    f1x(DP[i][1][1][1])
    for j in range(2):
      add(DP[i+1][0][1][0],DP[i][1][1][j])
  else:
    t=[[],[],0]
    add(DP[i][0][1][0],DP[i][0][1][1])
    add(DP[i][1][1][0],DP[i][1][1][1])
    add(DP[i+1][1][1][1],DP[i][0][1][0])
    fx0(DP[i+1][1][1][1])
    add(t,DP[i][0][1][0])
    add(t,DP[i][1][1][0])
    f1x(t)
    fx0(DP[i][1][1][0])
    add(DP[i+1][0][0][1],DP[i][1][1][0])
    add(DP[i+1][0][0][1],t)

def f10(p,o,e):
  r=[0]*3
  q=deepcopy(p)
  integrate(q)
  r[0]=(sum(q[0])-q[2])*o%mod
  r[1]=(sum(q[1])+q[2])*o%mod
  for j in range(2):
    for i in range(len(p[j])):
      r[j]=(r[j]+p[j][i]*g1[i]*e)%mod
      r[j+1]=(r[j+1]-p[j][i]*g1[i]*c[i]*e)%mod
  r[0]=(r[0]+p[2]*inv[2]*e)%mod
  r[2]=(r[2]-p[2]*inv[2]*e)%mod
  return r

add(DP[-1][0][1][1],DP[-1][1][1][1])
add(DP[-1][0][1][0],DP[-1][1][1][0])
add(DP[-1][0][0][1],DP[-1][1][0][1])
X=[f10(DP[-1][0][1][1],1,0),f10(DP[-1][0][1][0],0,1),f10(DP[-1][0][0][1],1,mod-1)]
A=[sum(map(lambda x:x[i],X))%mod for i in range(3)]
print(*A)

Solution (Python):



mod=10**9+7
def value(obj):
    if type(obj)==modint:
        tmp=obj.number
    else:
        tmp=obj
    return tmp%mod
class modint():
    number=0
    def __init__(self,num):
        self.number=num%mod
    def __add__(self,other):
        tmp=value(other)
        return modint(self.number+tmp)
    def __iadd__(self,other):
        tmp=value(other)
        self.number+=tmp
        self.number%=mod
        return self
    def __sub__(self,other):
        tmp=value(other)
        return modint(self.number-tmp)
    def __isub__(self,other):
        tmp=value(other)
        self.number-=tmp
        self.number%=mod
        return self
    def __mul__(self,other):
        tmp=value(other)
        return modint(self.number*tmp)
    def __imul__(self,other):
        tmp=value(other)
        self.number*=tmp
        self.number%=mod
        return self
    def __truediv__(self,other):
        tmp=value(other)
        return self*pow(tmp,mod-2,mod)
    def __itruediv(self,other):
        tmp=value(other)
        self.number=self/tmp
        self.number%=mod
        return self
    def __repr__(self):
        return str(self.number)
    def __eq__(self,other):
        return value(self)==value(other)
    def __ne__(self,other):
        return value(self)!=value(other)
    def __str__(self):
        return str(self.number)
    def __int__(self):
        return self.number
    def __hash__(self):
        return self.number

class polynomial():
    sequence=[0]
    def __init__(self,seq):
        self.sequence=[]
        for i in seq:
            self.sequence.append(value(i))
    def get(self,x):
        res=modint(0)
        tmp=modint(1)
        for a in self.sequence:
            res+=tmp*a
            tmp*=x
            if tmp==0:
                break
        return res
    def __add__(self,other):
        res=[]
        for i,a in enumerate(self.sequence):
            try:
                res[i]=a
            except IndexError:
                res.append(modint(0))
                res[i]=a
        for i,a in enumerate(other.sequence):
            try:
                res[i]+=a
            except IndexError:
                res.append(modint(0))
                res[i]+=a
        
        return polynomial(res)
    def __iadd__(self,other):
        return self+other
    def __sub__(self,other):
        res=[]
        for i,a in enumerate(self.sequence):
            try:
                res[i]=a
            except IndexError:
                res.append(modint(0))
                res[i]=a
        for i,a in enumerate(other.sequence):
            try:
                res[i]-=a
            except IndexError:
                res.append(modint(0))
                res[i]-=a
        return polynomial(res)
    def __isub__(self,other):
        return self-other
    def indef_integral(self):
        res=[modint(0)]
        for i,a in enumerate(self.sequence):
            res.append(modint(a)/(i+1))
        return polynomial(res)
    def integrate(self,flag):
        
        
        if flag==0:
            return self.indef_integral()
        else:
            tmp=self.indef_integral()
            return polynomial([tmp.get(1)])-tmp
    def __repr__(self):
        return str(self.sequence)
    def __str__(self):
        return str(self.sequence)

class func():
    P=polynomial([modint(0)])
    Q=polynomial([modint(0)])
    C=modint(0)
    def __init__(self,poly1,poly2,cnst):
        self.P=poly1
        self.Q=poly2
        self.C=cnst
    def __add__(self,other):
        PP=self.P+other.P
        QQ=self.Q+other.Q
        CC=self.C+other.C
        return func(PP,QQ,CC)
    def __iadd__(self,other):
        return self+other
    def __sub__(self,other):
        PP=self.P-other.P
        QQ=self.Q-other.Q
        CC=self.C-other.C
        return func(PP,QQ,CC)
    def __isub__(self,other):
        return self-other
    def dintegrate(self,flag):
        if flag==0:
            PP=self.P.integrate(0)+polynomial([self.C])
            QQ=self.Q.integrate(0)
            CC=modint(-1)*self.C
        else:
            PP=self.P.integrate(1)
            QQ=self.Q.integrate(1)+polynomial([modint(-1)*self.C])
            CC=self.C
        return func(PP,QQ,CC)
MAX_N=1009
N=int(input())
s=input()
fact=[modint(0) for i in range(MAX_N+1)]
inv=[modint(0) for i in range(MAX_N+1)]
fact[0]=modint(1)
for i in range(MAX_N):
    fact[i+1]=fact[i]*(i+1)
inv[-1]=modint(1)/fact[-1]
for i in range(MAX_N-1,-1,-1):
    inv[i]=inv[i+1]*(i+1)
dp=[[[[func(polynomial([modint(0)]),polynomial([modint(0)]),modint(0)) for l in range(2)]for k in range(2)]for j in range(2)] for i in range(MAX_N+1)]

if s[0]=="X":
    dp[0][0][1][0]=func(polynomial([modint(0)]),polynomial([modint(0)]),modint(1))
else:
    dp[0][0][0][1]=func(polynomial([modint(0)]),polynomial([modint(0)]),modint(1))
    dp[0][1][1][1]=func(polynomial([modint(1)]),polynomial([modint(0)]),modint(-1))
for i in range(1,N):
    if s[i]=="X":
        dp[i][0][1][0]+=dp[i-1][1][1][0].dintegrate(0)
        dp[i][0][1][0]+=dp[i-1][1][1][1].dintegrate(0)
        dp[i][0][1][0]+=dp[i-1][0][0][1].dintegrate(1)
        dp[i][0][1][0]+=dp[i-1][0][1][1].dintegrate(1)
        dp[i][0][1][0]+=dp[i-1][1][0][1].dintegrate(1)
        dp[i][0][1][0]+=dp[i-1][1][1][1].dintegrate(1)
    else:
        dp[i][0][0][1]+=dp[i-1][1][1][0].dintegrate(0)
        dp[i][0][0][1]+=dp[i-1][1][1][1].dintegrate(0)
        dp[i][1][1][1]+=dp[i-1][0][1][0].dintegrate(0)
        dp[i][1][1][1]+=dp[i-1][0][1][1].dintegrate(0)
        dp[i][0][0][1]+=dp[i-1][0][1][0].dintegrate(1)
        dp[i][0][0][1]+=dp[i-1][0][1][1].dintegrate(1)
        dp[i][0][0][1]+=dp[i-1][1][1][0].dintegrate(1)
        dp[i][0][0][1]+=dp[i-1][1][1][1].dintegrate(1)
a=modint(0)
b=modint(0)
c=modint(0)

p = dp[N-1][0][1][1].P.integrate(0)
a+= p.get(1)-p.get(0)
p = dp[N-1][0][1][1].Q.integrate(0)
b+= p.get(1)-p.get(0)
a+= dp[N-1][0][1][1].C
b-= dp[N-1][0][1][1].C
p = dp[N-1][1][1][1].P.integrate(0)
a+= p.get(1)-p.get(0)
p = dp[N-1][1][1][1].Q.integrate(0)
b+= p.get(1)-p.get(0)
a+= dp[N-1][1][1][1].C
b-= dp[N-1][1][1][1].C
for x in range(2):
    for i,v in enumerate(dp[N-1][x][1][0].P.sequence):
        a+=modint(1)*v*fact[i]
        for j in range(i+1):
            b-=modint(1)*v*fact[i]*inv[j]
    for i,v in enumerate(dp[N-1][x][1][0].Q.sequence):
        b+=modint(1)*v*fact[i]
        for j in range(i+1):
            c-=modint(1)*v*fact[i]*inv[j]
    a+=dp[N-1][x][1][0].C/2
    c-=dp[N-1][x][1][0].C/2

p = dp[N-1][0][0][1].P.integrate(0)
a+= p.get(1)-p.get(0)
p = dp[N-1][0][0][1].Q.integrate(0)
b+= p.get(1)-p.get(0)
a+= dp[N-1][0][0][1].C
b-= dp[N-1][0][0][1].C
p = dp[N-1][1][0][1].P.integrate(0)
a+= p.get(1)-p.get(0)
p = dp[N-1][1][0][1].Q.integrate(0)
b+= p.get(1)-p.get(0)
a+= dp[N-1][1][0][1].C
b-= dp[N-1][1][0][1].C
for x in range(2):
    for i,v in enumerate(dp[N-1][x][0][1].P.sequence):
        a-=modint(1)*v*fact[i]
        for j in range(i+1):
            b+=modint(1)*v*fact[i]*inv[j]
    for i,v in enumerate(dp[N-1][x][0][1].Q.sequence):
        b-=modint(1)*v*fact[i]
        for j in range(i+1):
            c+=modint(1)*v*fact[i]*inv[j]
    a-=dp[N-1][x][0][1].C/2
    c+=dp[N-1][x][0][1].C/2

print(a,b,c)

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return input().rstrip()
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
fc = [1] * 1010
ifc = [1] * 1010
Bn = [0] * 1010

def ad(x, y=0):
    return (x + y) % mod

def mlt(x, y=1):
    return (x * y) % mod
    
def pwr(x, y):
    return pow(x, y, mod)

def inv(x):
    return pwr(x, mod - 2)

def bn(x):
    if x > 0 and Bn[x] == 0: Bn[x] = ad(1, mlt(x, bn(x - 1)))
    return Bn[x]


def solve():
    for i in range(1, 1010): fc[i] = fc[i - 1] * i
    for i in range(1010): ifc[i] = inv(fc[i])
    num = 1
    dp = [[[0] * 1010 for _ in range(2)] for k in range(1010)]
    dpc = [0] * 1010
    f = [False] * 1010
    n = II()
    s = S()
    if n == 1 and s == "-":
        print("{} 0 {}".format((mod + 1) // 2, (mod + 1) // 2))
        return
    q = 0
    r = n
    n = -1
    if s[q] == "-":
        q += 1
    if s[q] == "-":
        q += 1
        n += 1
        num += 1
    if s[q] == "-":
        print("0 0 0")
        return
    for i in range(q, r):
        if s[i] == "X":
            if num == 0:
                print("0 0 0")
                return
            if n >= 0:
                f[n] = num - 1
            n += 1
            num = 0
        else:
            if num > 1:
                print("0 0 0")
                return
            num += 1
    if num == 2:
        f[n] = 1
        n += 1
    dpc[0] = 1
    for i in range(n):
        if f[i]:
            for j in range(2):
                for k in range(2 * n + 2):
                    kari = inv(k * k + 3 * k + 2)
                    dp[i + 1][j][0] = ad(dp[i + 1][j][0], mlt(dp[i][j][k], kari))
                    dp[i + 1][j][1] = ad(dp[i + 1][j][1], mlt(dp[i][j][k], mod - kari))
            dp[i + 1][1][0] = ad(dp[i + 1][1][0], dpc[i])
            dp[i + 1][1][1] = ad(dp[i + 1][1][1], mod - dpc[i])
        else:
            for j in range(2):
                for k in range(2 * n + 2):
                    kari = inv(k + 1)
                    dp[i + 1][j][0] = ad(dp[i + 1][j][0], mlt(dp[i][j][k], kari))
                    dp[i + 1][j][1] = ad(dp[i + 1][j][1], mlt(dp[i][j][k], mod - kari))
                    dp[i + 1][j][k + 2] = ad(dp[i + 1][j][k + 2], mlt(dp[i][j][k], inv(k * k + 3 * k + 2)))
            dp[i + 1][1][1] = ad(dp[i + 1][1][1], dpc[i])
            dp[i + 1][1][0] = ad(dp[i + 1][1][0], mod - dpc[i])
            dpc[i + 1] = ad(dpc[i + 1], dpc[i])
    ans = [0] * 3
    for j in range(2):
        for k in range(2 * n + 2):
            ans[j] = ad(ans[j], mlt(dp[n][j][k], fc[k]))
            ans[j + 1] = ad(ans[j + 1], mlt(dp[n][j][k], mod - ad(fc[k], bn(k))))
    kari = (mod + 1) // 2
    ans[0] = ad(ans[0], mlt(dpc[n], kari))
    ans[2] = ad(ans[2], mlt(dpc[n], mod - kari))
    print("{} {} {}".format(ans[0], ans[1], ans[2]))
    return



if __name__ == '__main__':
    solve()