結果

問題 No.665 Bernoulli Bernoulli
ユーザー maspymaspy
提出日時 2020-01-03 02:42:40
言語 Python3
(3.12.2 + numpy 1.26.4 + scipy 1.12.0)
結果
AC  
実行時間 618 ms / 2,000 ms
コード長 3,503 bytes
コンパイル時間 163 ms
コンパイル使用メモリ 13,184 KB
実行使用メモリ 46,740 KB
最終ジャッジ日時 2024-05-02 07:03:41
合計ジャッジ時間 12,994 ms
ジャッジサーバーID
（参考情報） 		judge2 / judge1
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 565 ms
46,224 KB
testcase_01 AC 569 ms
46,204 KB
testcase_02 AC 597 ms
46,204 KB
testcase_03 AC 603 ms
46,596 KB
testcase_04 AC 615 ms
46,212 KB
testcase_05 AC 601 ms
46,720 KB
testcase_06 AC 601 ms
46,464 KB
testcase_07 AC 595 ms
46,200 KB
testcase_08 AC 593 ms
46,488 KB
testcase_09 AC 601 ms
46,720 KB
testcase_10 AC 603 ms
46,352 KB
testcase_11 AC 618 ms
46,224 KB
testcase_12 AC 601 ms
46,464 KB
testcase_13 AC 598 ms
46,356 KB
testcase_14 AC 596 ms
46,736 KB
testcase_15 AC 593 ms
46,344 KB
testcase_16 AC 611 ms
46,740 KB
testcase_17 AC 598 ms
46,720 KB
testcase_18 AC 593 ms
46,464 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

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

import numpy as np

N,K = map(int,read().split())

MOD = 10 ** 9 + 7

def make_power(a, L, MOD=MOD):
    B = L.bit_length()
    x = np.empty((1<<B), np.int64)
    x[0] = 1
    for n in range(B):
        x[1<<n:1<<(n+1)] = x[:1<<n] * a % MOD
        a *= a; a %= MOD
    return x[:L]
def cumprod(A, MOD = MOD):
    L = len(A); Lsq = int(L**.5+1)
    A = np.resize(A, Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        A[:,n] *= A[:,n-1]; A[:,n] %= MOD
    for n in range(1,Lsq):
        A[n] *= A[n-1,-1]; A[n] %= MOD
    return A.ravel()[:L]
def make_fact(U, MOD = MOD):
    x = np.arange(U, dtype = np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact,fact_inv

def fft_convolve(f, g, MOD = MOD):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    fft = np.fft.rfft; ifft = np.fft.irfft
    Lf = len(f); Lg = len(g); L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fl = f & (1 << 15) - 1; fh = f >> 15
    gl = g & (1 << 15) - 1; gh = g >> 15
    conv = lambda f,g: ifft(fft(f,fft_len) * fft(g,fft_len))[:L]
    x = conv(fl, gl) % MOD
    y = conv(fl+fh, gl+gh) % MOD
    z = conv(fh, gh) % MOD
    a, b, c = map(lambda x: (x + .5).astype(np.int64), [x,y,z])
    return (a + ((b - a - c) << 15) + (c << 30)) % MOD

def fps_inv(F,MOD=MOD):
    L = len(F)
    i = (L-1).bit_length()
    N = 1 << i
    F = np.resize(F,N)
    n = 1
    x = pow(int(F[0]),MOD-2,MOD)
    R = np.array([x],np.int64)
    while n < N:
        RF = fft_convolve(R,F[:n+n])[:n+n]
        RRF = fft_convolve(R,-RF)[:n+n]
        RRF[:n] += R+R
        R = RRF; R %= MOD
        n += n
    return R[:L]

def fps_exp(h,MOD=MOD):
    # 面倒なので2べきに
    L = len(h)
    i = (L-1).bit_length()
    N = (1 << i)
    h = np.resize(h,N)
    dh = np.empty_like(h)
    dh[0:N-1] = h[1:N] * np.arange(1,N,dtype=np.int64) % MOD
    f = np.zeros_like(h); g = np.zeros_like(h)
    f[:2] = np.array([1,h[1]],np.int64); g[0] = 1; m = 2
    while m <= N//2:
        fg = fft_convolve(f[:m],g[:m],MOD)[:m]
        fgg = fft_convolve(fg,g[:m],MOD)[:m]
        g[:m] *= 2; g[:m] -= fgg; g %= MOD
        q = dh[:m-1]
        s = np.zeros(m+m,np.int64); s[1:m] = f[1:m] * np.arange(1,m,dtype=np.int64) % MOD
        r = fft_convolve(f[:m],q,MOD)
        s[1:1+len(r)] -= r
        s[0:m] += s[m:m+m]
        s = s[:m] % MOD
        t = fft_convolve(g[:m],s,MOD)[:m]
        t *= inv[m:m+m]; t %= MOD
        u = h[m:m+m] - t; u %= MOD
        v = fft_convolve(f,u,MOD)[:m]
        f[m:m+m] += v
        m *= 2
    return f[:L]

def Bernoulli(N,MOD = MOD):
    den = fact_inv[1:N+1]
    B = fps_inv(den)
    B *= fact[:len(B)]; B %= MOD
    return B

def kth_power_sum_polynomial(k,MOD=MOD):
    B = Bernoulli(k+1,MOD)
    S = np.zeros(k+2,np.int64)
    C = fact[k] * fact_inv[1:k+2] % MOD * fact_inv[:k+1][::-1] % MOD
    S[1:] = B[::-1] * C
    S[(K&1)::2] *= (-1)
    return S%MOD

fact, fact_inv = make_fact(10**5)

S = kth_power_sum_polynomial(K)

N %= MOD # 今回の状況だと分母にMODが入らないので正当
power = make_power(N,len(S))

answer = (S * power % MOD).sum() % MOD
print(answer)
0