考虑平面上的所有整点，若$x,y,x+y$都是完全平方数，那么$(x,y)$不可通过，否则可以通过。例如$(9, 16)$不可通过，而$(0, 4), (3,1), (9,4)$都可以通过。

现考察从点$(0, 0)$到点$(n, n)$的合法路径，合法路径是指每次只往上或往右移动一格，且只经过可通过的点的路径。

记$P(n)$为合法路径路径总数，可以验证$P(5)=252, P(16)=596994440, P(1000) \bmod 10^9+7=341920854$。

求$P(10^7) \bmod 10^9+7$。

本题难度:

解答

若所有的点都可通过，那么路径总数是$\binom{2n}{n}$（从2n步中选出n步向上）。

存在不可通过的点时，反复使用动态规划和容斥原理来计算出从原点到各不可通过的点的合法路径数，最后从总路径中减去这些数量即可。

不可通过的点用勾股数的欧几里得公式很容易生成。

最终结果是$299742733$。

import fractions

target=10**7
mod=1000000007
f=[1]
for k in range(1,2*target+1):
    f.append((f[-1]*k)%mod)

g=[1]*(2*target+1)
g[-1]=pow(f[2*target],mod-2,mod)
for k in range(2*target,1,-1):
    g[k-1]=(g[k]*k)%mod

p=[]
for m in range(2,2000):
    for n in range(1,m):
        u,v=m*m-n*n,2*m*n
        if fractions.gcd(u,v)==1:
            a,b=u,v
            while a*a<=target and b*b<=target:
                p.append((a*a,b*b))
                p.append((b*b,a*a))
                a+=u
                b+=v
p.sort()

s=f[2*target]*g[target]*g[target]
d=[0]*(len(p)+1)

for k,(x0,y0) in enumerate(p):
    d[k]=(f[x0+y0]*g[x0]*g[y0])%mod
    for j in range(k):
        x1,y1=p[j][0],p[j][1]
        if y1<=y0:
            d[k]=(d[k]-d[j]*f[x0-x1+y0-y1]*g[x0-x1]*g[y0-y1])%mod
    s-=d[k]*f[2*target-x0-y0]*g[target-x0]*g[target-y0]

print s%mod