Problem:
 Score : 1000 points    Problem Statement Inspired by the tv series Stranger Things , bear Limak is going for a walk between two mirror worlds.  There are two perfect binary trees of height H ,
    each with the standard numeration of vertices from 1 to 2^H-1 .
The root is 1 and the children of x are 2 \cdot x and 2 \cdot x + 1 .  Let L denote the number of leaves in a single tree, L = 2^{H-1} .  You are given a permutation P_1, P_2, \ldots, P_L of numbers 1 through L .
It describes L  special edges that connect leaves of the two trees.
There is a special edge between vertex L+i-1 in the first tree and vertex L+P_i-1 in the second tree.   drawing for the first sample test, permutation P = (2, 3, 1, 4) , special edges in green  Let's define product of a cycle as the product of numbers in its vertices.
Compute the sum of products of all simple cycles that have exactly two special edges , modulo (10^9+7) .  A simple cycle is a cycle of length at least 3, without repeated vertices or edges.      Constraints  2 \leq H \leq 18  1 \leq P_i \leq L where L = 2^{H-1}  P_i \neq P_j (so this is a permutation)         Input Input is given from Standard Input in the following format (where L = 2^{H-1} ).  H  P_1  P_2  \cdots  P_L       Output Compute the sum of products of simple cycles that have exactly two special edges.
Print the answer modulo (10^9+7) .        Sample Input 1 3
2 3 1 4      Sample Output 1 121788  The following drawings show two valid cycles (but there are more!) with products 2 \cdot 5 \cdot 6 \cdot 3 \cdot 1 \cdot 2 \cdot 5 \cdot 4 = 7200 and 1 \cdot 3 \cdot 7 \cdot 7 \cdot 3 \cdot 1 \cdot 2 \cdot 5 \cdot 4 \cdot 2 = 35280 .
The third cycle is invalid because it doesn't have exactly two special edges.        Sample Input 2 2
1 2      Sample Output 2 36  The only simple cycle in the graph indeed has two special edges, and the product of vertices is 1 \cdot 3 \cdot 3 \cdot 1 \cdot 2 \cdot 2 = 36 .       Sample Input 3 5
6 14 15 7 12 16 5 4 11 9 3 10 8 2 13 1      Sample Output 3 10199246

Solution (Python):


H=int(input())


P=list(map(lambda x:int(x)-1,input().split()))
mod=10**9+7
inv=pow(2,mod-2,mod)

L=2**(H-1)

base_t=[1 for i in range(1+1<<H)]
base_t[0]=0
for i in range(2,1+1<<H):
    base_t[i]=i*base_t[i>>1]
    base_t[i]%=mod

base_s=[(i*pow(base_t[i]**2,mod-2,mod))%mod for i in range(1+1<<H)]

ans=0
def solve(i):
    global ans
    if i>=L:
        i-=L
        val=(base_t[i+L]*base_t[P[i]+L])%mod
        val2=pow(val,2,mod)
        id=P[i]+L
        res=[[(P[i]+L)>>(H-1-j),val,val2] for j in range(H)]
        return res
    L1=solve(2*i)
    L2=solve(2*i+1)
    Li,ttt=merge(L1,L2)
    ttt=(ttt*(base_s[i]-base_s[i>>1]))%mod
    ans=(ans+ttt)%mod
    return Li

def merge(L1,L2):
    res=[]
    L1.append([10**18,0])
    L2.append([10**18,0])
    L1=L1[::-1]
    L2=L2[::-1]
    ttt=0
    if L1[-1][0]<=L2[-1][0]:
        res.append(L1.pop())
    else:
        res.append(L2.pop())
    while L1 and L2:
        if L1[-1][0]<=L2[-1][0]:
            if res and res[-1][0]==L1[-1][0]:
                res[-1][1]+=L1[-1][1]
                res[-1][1]%=mod
                res[-1][2]+=L1[-1][2]
                res[-1][2]%=mod
                L1.pop()
            else:
                j,val,val2=res[-1]
                tmp=(base_s[j]-base_s[j>>1]) %mod
                c=val**2 %mod
                c=(c-val2)%mod
                tmp =(tmp*c) %mod
                ttt=(ttt+tmp)%mod
                res.append(L1.pop())
        else:
            if res and res[-1][0]==L2[-1][0]:
                res[-1][1]+=L2[-1][1]
                res[-1][1]%=mod
                res[-1][2]+=L2[-1][2]
                res[-1][2]%=mod
                L2.pop()
            else:
                j,val,val2=res[-1]
                tmp=(base_s[j]-base_s[j>>1]) %mod
                c=val**2 %mod
                c=(c-val2)%mod
                tmp =(tmp*c) %mod
                ttt=(ttt+tmp)%mod
                res.append(L2.pop())
    res.pop()
    return res,ttt

solve(1)

print(ans*inv %mod)

Solution (Python):
import os
import sys

import numpy as np


def solve(inp):
    def mod_pow(x, a, MOD):
        ret = 1
        cur = x
        while a:
            if a & 1:
                ret = ret * cur % MOD
            cur = cur * cur % MOD
            a >>= 1
        return ret

    h = inp[0]
    ppp = inp[1:]
    n = 1 << (h - 1)
    ppp += n - 1
    MOD = 10 ** 9 + 7

    cumprods = np.ones(2 * n, dtype=np.int64)
    cumprods_rev = np.ones(n, dtype=np.int64)
    cumprods_through = np.zeros(n, dtype=np.int64)
    for i in range(2, 2 * n):
        cumprods[i] = cumprods[i >> 1] * i % MOD
    for i in range(n):
        cumprods_rev[i] = mod_pow(cumprods[i], MOD - 2, MOD)
        cumprods_through[i] = cumprods[i + n] * cumprods[ppp[i]] % MOD

    cumprods_from_tree1 = np.zeros(2 * n, dtype=np.int64)
    ans = 0

    for lca in range(1, n):
        d = lca
        digit = h
        while d:
            d >>= 1
            digit -= 1

        leftmost_leaf_in_left_subtree = lca << digit
        leftmost_leaf_in_right_subtree = ((lca << 1) + 1) << (digit - 1)
        rightmost_leaf_in_right_subtree = (lca + 1) << digit

        rev = cumprods_rev[lca >> 1]
        for leaf in range(leftmost_leaf_in_left_subtree, leftmost_leaf_in_right_subtree):
            v = ppp[leaf - n]
            cp = cumprods_through[leaf - n] * rev % MOD
            while v > 1:
                cumprods_from_tree1[v] += cp * cumprods_rev[v >> 1] % MOD
                v >>= 1

        rev = cumprods_rev[lca]
        for leaf in range(leftmost_leaf_in_right_subtree, rightmost_leaf_in_right_subtree):
            v = ppp[leaf - n]
            cp = cumprods_through[leaf - n] * rev % MOD
            while v > 1:
                ans += cumprods_from_tree1[v ^ 1] % MOD * cp % MOD * cumprods_rev[v >> 2] % MOD
                v >>= 1
            ans %= MOD

        for leaf in range(leftmost_leaf_in_left_subtree, leftmost_leaf_in_right_subtree):
            v = ppp[leaf - n]
            while cumprods_from_tree1[v] != 0:
                cumprods_from_tree1[v] = 0
                v >>= 1

    return ans


if sys.argv[-1] == 'ONLINE_JUDGE':
    from numba.pycc import CC

    cc = CC('my_module')
    cc.export('solve', '(i8[:],)')(solve)
    cc.compile()
    exit()

if os.name == 'posix':
    
    from my_module import solve
else:
    from numba import njit

    solve = njit('(i8[:],)', cache=True)(solve)
    print('compiled', file=sys.stderr)

inp = np.fromstring(sys.stdin.read(), dtype=np.int64, sep=' ')
ans = solve(inp)
print(ans)

Solution (Python):
import sys

h, *ppp = map(int, sys.stdin.buffer.read().split())
n = 1 << (h - 1)
ppp = [p + n - 1 for p in ppp]
MOD = 10 ** 9 + 7

cumprods = [1] * (2 * n)  
cumprods_rev = [1] * n  
cumprods_through = [0] * n  
for i in range(2, 2 * n):
    cumprods[i] = cumprods[i >> 1] * i % MOD
for i in range(n):
    cumprods_rev[i] = pow(cumprods[i], MOD - 2, MOD)
    cumprods_through[i] = cumprods[i + n] * cumprods[ppp[i]] % MOD



cumprods_from_tree1 = [0] * (2 * n)

ans = 0

for lca in range(1, n):
    digit = h - lca.bit_length()
    leftmost_leaf_in_left_subtree = lca << digit
    leftmost_leaf_in_right_subtree = ((lca << 1) + 1) << (digit - 1)
    rightmost_leaf_in_right_subtree = (lca + 1) << digit

    rev = cumprods_rev[lca >> 1]
    for leaf in range(leftmost_leaf_in_left_subtree, leftmost_leaf_in_right_subtree):
        v = ppp[leaf - n]
        cp = cumprods_through[leaf - n] * rev % MOD
        while v > 1:
            cumprods_from_tree1[v] += cp * cumprods_rev[v >> 1] % MOD
            v >>= 1

    rev = cumprods_rev[lca]
    for leaf in range(leftmost_leaf_in_right_subtree, rightmost_leaf_in_right_subtree):
        v = ppp[leaf - n]
        cp = cumprods_through[leaf - n] * rev % MOD
        while v > 1:
            ans = (ans + cumprods_from_tree1[v ^ 1] * cp * cumprods_rev[v >> 2]) % MOD
            v >>= 1

    for leaf in range(leftmost_leaf_in_left_subtree, leftmost_leaf_in_right_subtree):
        v = ppp[leaf - n]
        while cumprods_from_tree1[v] != 0:
            cumprods_from_tree1[v] = 0
            v >>= 1

print(ans)

Solution (Python):
import sys
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

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

MOD = 1_000_000_007

@njit((i8, i8), cache=True)
def mpow(a, n):
    p = 1
    while n:
        if n & 1:
            p = p * a % MOD
        a = a * a % MOD
        n >>= 1
    return p


@njit((i8, ), cache=True)
def fact_table(N):
    """Compute 0!, 1!, ..., (N-1)! and their inverses modulo MOD. 

    Parameters
    ----------
    N : int
        length of the output array. 

    Returns
    -------
    fact : np.ndarray
        0!, 1!, ..., (N-1)!

    fact_inv : np.ndarray
        inverse of 0!, 1!, ..., (N-1)! modulo MOD. 
    """
    fact = np.empty(N, np.int64)
    fact[0] = 1
    for n in range(1, N):
        fact[n] = n * fact[n - 1] % MOD
    fact_inv = np.empty(N, np.int64)
    fact_inv[N - 1] = mpow(fact[N - 1], MOD - 2)
    for n in range(N - 1, 0, -1):
        fact_inv[n - 1] = fact_inv[n] * n % MOD
    return fact, fact_inv

@njit((i8, i8[:]), cache=True)
def main(H, P):
    fact, fact_inv = fact_table(1 << 20)
    inv = np.zeros(1 << 20, np.int64)
    inv[1:] = fact_inv[1:] * fact[:-1] % MOD
    
    f = np.ones_like(P)
    f_inv = np.ones_like(P)
    g = np.ones_like(P)
    g[0] = 0
    for v in range(2, 1 << H):
        parent = v >> 1
        f[v] = f[parent] * parent % MOD
        f_inv[v] = f_inv[parent] * inv[parent] % MOD
        g[v] = f_inv[v] * f_inv[v] % MOD * inv[v] % MOD
    
    dp1 = np.zeros(1 << H, np.int64)  
    dp2 = np.zeros(1 << H, np.int64)  
    ws = np.zeros(1 << H, np.int64)
    ptr = 0
    ans = 0
    for x1 in range(1, 1 << H):
        l, r = x1, x1 + 1
        while 2 * l < 1 << H:
            l, r = 2 * l, 2 * r
        for v in range(l, r):
            pv = f[v] * v % MOD * f_inv[x1] % MOD  
            w = P[v]
            pw = pv
            while w:
                pw = pw * w % MOD
                if dp1[w] == 0:
                    ws[ptr] = w
                    ptr += 1
                dp1[w] += pw
                dp2[w] += pw * pw % MOD
                w >>= 1
        for x2 in ws[:ptr]:
            dp1[x2] %= MOD
            dp2[x2] %= MOD
            x = (dp1[x2]**2 - dp2[x2]) % MOD
            x = x * inv[2] % MOD
            coef = (f[x1] * f[x2] % MOD)**2 % MOD
            x = x * coef % MOD
            p1, p2 = x1 >> 1, x2 >> 1
            coef = g[x1] * g[x2] % MOD
            coef -= g[p1] * g[x2] % MOD
            coef -= g[x1] * g[p2] % MOD
            coef += g[p1] * g[p2] % MOD
            coef %= MOD
            ans += coef * x % MOD
        for x2 in ws[:ptr]:
            dp1[x2] = 0
            dp2[x2] = 0
        ptr = 0
    return ans % MOD

H = int(readline())
P = np.zeros(1 << H, np.int64)
P[1 << H - 1:] = (np.array(read().split(), np.int64) - 1) | 1 << H - 1

print(main(H, P))