Problem:
 Score : 1200 points    Problem Statement One day Mr. Takahashi picked up a dictionary containing all of the N ! permutations of integers 1 through N .
The dictionary has N ! pages, and page i ( 1 ≤ i ≤ N! ) contains the i -th permutation in the lexicographical order.  Mr. Takahashi wanted to look up a certain permutation of length N in this dictionary, but he forgot some part of it.  His memory of the permutation is described by a sequence P_1, P_2, ..., P_N .
If P_i = 0 , it means that he forgot the i -th element of the permutation; otherwise, it means that he remembered the i -th element of the permutation and it is P_i .  He decided to look up all the possible permutations in the dictionary.
Compute the sum of the page numbers of the pages he has to check, modulo 10^9 + 7 .      Constraints  1 ≤ N ≤ 500000  0 ≤ P_i ≤ N  P_i ≠ P_j if i ≠ j ( 1 ≤ i, j ≤ N ), P_i ≠ 0 and P_j ≠ 0 .       Partial Score  In test cases worth 500 points, 1 ≤ N ≤ 3000 .         Input The input is given from Standard Input in the following format:  N  P_1  P_2  ...  P_N       Output Print the sum of the page numbers of the pages he has to check, as modulo 10^9 + 7 .        Sample Input 1 4
0 2 3 0      Sample Output 1 23  The possible permutations are [ 1 , 2 , 3 , 4 ] and [ 4 , 2 , 3 , 1 ].
Since the former is on page 1 and the latter is on page 22 , the answer is 23 .       Sample Input 2 3
0 0 0      Sample Output 2 21  Since all permutations of length 3 are possible, the answer is 1 + 2 + 3 + 4 + 5 + 6 = 21 .       Sample Input 3 5
1 2 3 5 4      Sample Output 3 2  Mr. Takahashi remembered all the elements of the permutation.       Sample Input 4 1
0      Sample Output 4 1  The dictionary consists of one page.       Sample Input 5 10
0 3 0 0 1 0 4 0 0 0      Sample Output 5 953330050

Solution (Python):

"""

https://atcoder.jp/contests/code-festival-2016-qualc/tasks/codefestival_2016_qualC_e

ページの合計に与える寄与を考えればよい
ある位置iにある数字の与える寄与は
(自分より右にある数字の数)! * 自分より右にあるより小さい数字の数

元からおかれている数字同士に関しての寄与はBITで終わり
元からおかれている数字 & 自由な数字の寄与は
順列の総数 * 自由な数字の中で自分より小さいものの数 * 右にある空欄の数 / 全ての空欄の数 * fac[N-1-i]
順列の総数 * 自由な数字の中で自分より大きいものの数 * 左にある空欄の数 / 全ての空欄の数 * fac[N-1-i]
で求まる

自由同士は、半分になるはずなので
fac[N-1-i] * 右にある0の数 // 2
それより左にあるfacの平均を出す？→avrage
順列の総数 * 自由な数字の中で自分より大きいものの数 * 左にある空欄の数 / 全ての空欄の数 * average

"""

mod = 10**9+7

def inverse(a,mod): 
    return pow(a,mod-2,mod)





def modfac(n, MOD):
 
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv
    return factorials, invs

def modnCr(n,r,mod,fac,inv): 
    return fac[n] * inv[n-r] * inv[r] % mod

def bitadd(a,w,bit): 
 
    x = a
    while x <= (len(bit)-1):
        bit[x] += w
        x += x & (-1 * x)
 
def bitsum(a,bit): 
 
    ret = 0
    x = a
    while x > 0:
        ret += bit[x]
        x -= x & (-1 * x)
    return ret

import sys
from sys import stdin
N = int(stdin.readline())
a = list(map(int,stdin.readline().split()))

BIT = [0] * (N+1)

fac,inv = modfac(N+10,mod)
zeronum = 0
 
ans1 = 0
RZ = [0] * N
app = [False] * (N+1)
app[0] = True

for i in range(N-1,-1,-1):

    if a[i] != 0:
        ans1 += fac[N-1-i] * ( bitsum(a[i],BIT) )
        ans1 %= mod
        bitadd(a[i],1,BIT)
        app[a[i]] = True
    else:
        zeronum += 1
    RZ[i] = zeronum

mEX = [0] * (N+1) 
for i in range(1,N+1):
    if app[i] == False:
        mEX[i] = mEX[i-1] + 1
    else:
        mEX[i] = mEX[i-1]

ans1 *= fac[zeronum]

ans2 = 0
ans3 = 0
tmpinv = inverse(zeronum,mod)
tmpsum = 0
for i in range(N):
    if a[i] != 0:
        ans2 += fac[zeronum] * mEX[a[i]] * RZ[i] * tmpinv * fac[N-1-i]
        ans2 += fac[zeronum] * (zeronum-mEX[a[i]]) * tmpinv * tmpsum
        ans2 %= mod
    else:
        ans3 += fac[N-1-i] * fac[zeronum] * (RZ[i]-1) * inverse(2,mod)
        tmpsum += fac[N-1-i]

print (ans1 , ans2 , ans3 , fac[zeronum],file=sys.stderr)
print ((ans1 + ans2 + ans3 + fac[zeronum]) % mod)

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


M = 10**9+7
n = int(input())
p = list(map(int, input().split()))
g1 = [None]*(n+10)
v = 1
g1[0] = 1
for i in range(1,len(g1)):
    v *= i
    v %= M
    g1[i] = v
    
def init(bit, values):
    for i,v in enumerate(values):
        add(bit,i+1,v)

def query(bit,i):
    res = 0
    while i > 0:
        res += bit[i]
        i -= i&(-i)
    return res


def add(bit,i,x):
    if i==0:
        raise RuntimeError
    while i <= len(bit)-1:
        bit[i] += x
        i += i&(-i)
    return


index = [i for i,num in enumerate(p) if num==0]
used = [False]*n
for i in range(n):
    if p[i]>0:
        used[p[i]-1] = True
nl = []
for i in range(n):
    if not used[i]:
        nl.append(i+1)
nl.sort()
scores = [0]*(n+1) 
v = len(nl)
cur = 0
for i in range(1,n+1):
    if nl and cur<len(nl) and nl[cur]<i:

        cur += 1
        v -= 1
    scores[i] = v
ans = 0
bit = [0]*(n+1)
s = set(index)
val = 0
pp = len(index) 
ss = sum(nl)
num = 0
inv2 = pow(2, M-2, M)
for i in range(n):
    if i in s:
        ans += (g1[pp-1] * (ss - pp - val - pp*(num)*inv2)) * g1[n-i-1]
        num += 1
    else:
        ans += (g1[pp] * (p[i]-query(bit, p[i])-1) - g1[pp-1]*num*(len(nl)-scores[p[i]])) * g1[n-i-1]
        add(bit, p[i], 1)
        val += scores[p[i]]
    ans %= M

print((ans + g1[pp])%M)

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


M = 10**9+7
n = int(input())
p = list(map(int, input().split()))
g1 = [None]*(n+10)
v = 1
g1[0] = 1
for i in range(1,len(g1)):
    v *= i
    v %= M
    g1[i] = v
    
def init(bit, values):
    for i,v in enumerate(values):
        add(bit,i+1,v)

def query(bit,i):
    res = 0
    while i > 0:
        res += bit[i]
        i -= i&(-i)
    return res


def add(bit,i,x):
    if i==0:
        raise RuntimeError
    while i <= len(bit)-1:
        bit[i] += x
        i += i&(-i)
    return


index = [i for i,num in enumerate(p) if num==0]
used = [False]*n
for i in range(n):
    if p[i]>0:
        used[p[i]-1] = True
nl = []
for i in range(n):
    if not used[i]:
        nl.append(i+1)
nl.sort()
scores = [0]*(n+1) 
v = len(nl)
cur = 0
for i in range(1,n+1):
    if nl and cur<len(nl) and nl[cur]<i:

        cur += 1
        v -= 1
    scores[i] = v
ans = 0
bit = [0]*(n+1)
s = set(index)
val = 0
pp = len(index) 
ss = sum(nl)
num = 0
inv2 = pow(2, M-2, M)
for i in range(n):
    if p[i]==0:
        ans += (g1[pp-1] * (ss - pp - val - pp*(num)*inv2)) * g1[n-i-1]
        num += 1
    else:
        ans += (g1[pp] * (p[i]-query(bit, p[i])-1) - g1[pp-1]*num*(len(nl)-scores[p[i]])) * g1[n-i-1]
        add(bit, p[i], 1)
        val += scores[p[i]]
    ans %= M

print((ans + g1[pp])%M)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
sys.setrecursionlimit(10 ** 7)

import itertools

MOD = 10**9 + 7



N,*A = map(int,read().split())

A_rev = A[::-1]

class BIT():
    def __init__(self, max_n):
        self.size = max_n + 1
        self.tree = [0] * self.size
        
    def get_sum(self,i):
        s = 0
        while i:
            s += self.tree[i]
            i -= i & -i
        return s
 
    def add(self, i, x):
        while i < self.size:
            self.tree[i] += x
            i += i & -i

right_zero = [0] * (N+1)
for i,x in enumerate(A_rev):
    right_zero[i+1] = right_zero[i] 
    if x == 0:
        right_zero[i+1] += 1

zero_cnt = right_zero[-1]
inv = pow(zero_cnt,MOD-2,MOD)

rest_smaller = [1] * (N+1)
for x in A:
    rest_smaller[x] = 0
rest_smaller = list(itertools.accumulate(rest_smaller))

right_smaller_prob = [(MOD+1)//2] * N
for i,x in enumerate(A_rev):
    if x != 0:
        right_smaller_prob[i] = rest_smaller[x] * inv

right_smaller_filled = [0] * N
bit = BIT(max_n = N)
p = 0
for i,x in enumerate(A_rev):
    if x != 0:
        right_smaller_filled[i] = bit.get_sum(x)
        bit.add(x,1)
        p += zero_cnt-rest_smaller[x]
    else:
        right_smaller_filled[i] = p * inv

fact = [1] * (N+1)
for n in range(1,N+1):
    fact[n] = fact[n-1] * n % MOD

answer = (1+sum(y*(x*p+z) for x,p,y,z in zip(right_zero,right_smaller_prob,fact,right_smaller_filled))) * fact[zero_cnt] % MOD
print(answer)

Solution (Python):
class Bit:
    def __init__(self, n):
        sz = 1
        while n >= sz:
            sz *= 2
        self.size = sz
        self.data = [0]*sz

    def sum(self, i):
        assert i > 0
        s = 0
        while i > 0:
            s += self.data[i]
            i -= i & -i
        return s

    def add(self, i, x):
        assert i > 0
        while i < self.size:
            self.data[i] += x
            i += i & -i

import bisect
N = int(raw_input())
a = [int(x)-1 for x in raw_input().split()]
weight = [1] * N
M = 10**9+7
def factorial(n):
    acc = 1
    for i in xrange(1, 1 + n):
        acc = (acc * i) % M
    return acc
def inv(n):
    return pow(n, M-2, M)
for i in xrange(N-1):
    weight[N-2-i] = (weight[N-1-i] * (i+1)) % M
missing = list(set(xrange(N)) - set(a))
missing.sort()
n_missing = len(missing)
sum_missing = sum(missing)
avg_missing = sum_missing * inv(n_missing)
n_present = N - n_missing
n_rows = factorial(n_missing)
acc = 0
i_missing = 0
i_present = 0
missing_lt_tot = 0

present_values_bit = Bit(N+5)
inv_n_missing = inv(n_missing)
inv_2 = inv(2)
for i, (w, x) in enumerate(zip(weight, a)):
    if x == -1:
        ai = avg_missing
        bi = (i_missing * inv_2 + i_present - missing_lt_tot * inv_n_missing) % M
        i_missing += 1
    else:
        ai = x
        
        bi = present_values_bit.sum(x + 1)
        
        present_values_bit.add(x + 2, 1)
        '''
        missing_lt = 0
        for j in missing:
            if j < a[i]:
                missing_lt += 1
        '''
        missing_lt = bisect.bisect(missing, a[i])
        missing_lt_tot += missing_lt
        
        bi += missing_lt * inv_n_missing * i_missing
        bi %= M
        i_present += 1
    
    acc = (acc + (((ai - bi) * n_rows) % M) * w) % M
print (acc + n_rows) % M