AtCoder Beginner Contest 239 Ex

https://atcoder.jp/contests/abc239/tasks/abc239_h

例えば、N=4, M=5で考えてみましょう。
1回サイコロを振って、目が1ならそのまま、2なら次以降が2より大きければ終わり、3なら次以降が1より大きければ終わり、4と5も同様です。 E_mをmになる回数の期待値とすると、

 \displaystyle E_5 = \frac{1}{5}(1 + E_5) + \frac{1}{5}(1 + E_2) + \frac{1}{5}(1 + E_1) + \frac{1}{5}(1 + E_1)

これを整理すると、

 \displaystyle E_5 = \frac{1}{4}(5 + E_2 + E_1 + E_1)

一般には、

 \displaystyle E_m = \frac{1}{N-1}(N+1 + \sum_{n=2}^{N}{E_{\lfloor \frac{m}{n} \rfloor}})

となります。計算量は O(N\sqrt{M})ありそう間に合いません。
ただ、最後の2項はこの場合、同じですね。例えば、N=100, M=100なら、n=51からn=100まで同じです。これは、よくあるパターンで、

https://atcoder.jp/contests/abc230/tasks/abc230_e

と同じです。m/nの最小値に注意するとよいです。

# coding: utf-8
# Dice Product 2

from itertools import *


#################### library ####################

def read_tuple():
    return tuple(map(int, raw_input().split()))

# ax = by + c
def linear_diophantine(a, b, c):
    def f(a, b, c):
        if a == 1:
            return (b + c, 1)
        elif a == -1:
            return (-b - c, 1)
        
        d = b / a
        r = b % a
        t = f(r, a, -c)
        return (t[1] + d * t[0], t[0])
    
    return f(a, b, c)

def inverse(n, p):
    x, y = linear_diophantine(n, p, 1)
    return x


#################### naive ####################

def F_naive(N, M):
    memo = { }
    def E(m):
        if m == 0:
            return 0
        elif m in memo:
            return memo[m]
        
        s = N
        for n in range(2, N+1):
            s += E(m/n)
        s = s * inv % D
        memo[m] = s
        return s
    
    inv = inverse(N-1, D)
    return E(M)


#################### process ####################

def read_input():
    N, M = read_tuple()
    return (N, M)

def F(N, M):
    memo = { }
    def E(m):
        if m == 0:
            return 0
        elif m in memo:
            return memo[m]
        
        s = 0
        k0 = max(1, m/N)
        for k in count(k0):
            if m/k <= k:
                break
            a, b = m/(k+1), m/k
            b = min(N, b)
            s += (b - a) * E(k)
        
        for n in range(2, m/k+1):
            s += E(m/n)
        
        s = (N + s) * inv % D
        memo[m] = s
        return s
    
    inv = inverse(N - 1, D)
    return E(M)


#################### main ####################

D = 10**9+7
N, M = read_input()
print F(N, M)