Project Euler 759

Project Euler 759

题目

A squared recurrence relation

The function f is defined for all positive integers as follows:

f(1)=1f(2n)=2f(n)f(2n+1)=2n+1+2f(n)+1nf(n)

It can be proven that f(n) is integer for all values of n.

The function S(n) is defined as S(n)=i=1nf(i)2.

For example, S(10)=1530 and S(102)=4798445.

Find S(1016). Give your answer modulo 1000000007.

解决方案

暴力枚举出 f 的前几项,在 OEIS 中查询得到结果为 A245788

b(n) 表示 n 的二进制表示中有多少个 1,那么按照题意,f(n)=nb(n)

S 进行改写:

S(n)=i=1log2(n+1)(i2mn,b(m)=im2)

也就是说,问题转化成枚举 i,求 n 以内的所有数中,满足 b(m)=i 的所有 m2 之和。

n 写成一个长度为 l 的二进制数字符串 d1d2dl,也就是有 l=log2(n+1)

中间状态 c1(i,j)(1il,0j9i) 分别表示如下含义:有多少个 i有前导 0 的数,其二进制有 j 1,并等于由字符串 d1d2di 表示的数。令状态 s1(i,j) 表示 c1(i,j) 中的所有数之和。令状态 t1(i,j) 表示 c1(i,j) 中的所有数的平方和

不难发现,如果 j=k=1idk,那么 c1(i,j)=1s1(i,j) 的值则为 d1d2di 这一个字符串的二进制数,并且 t1(i,j)=s1(i,j)2;否则 c1(i,j)=s1(i,j)=t1(i,j)=0

状态 c1,s1,t1 的意义在于规定了答案严格在界限上的情况。而接下来定义的状态 c0,s0,t0 则是严格不在界限上的情况。

令状态 c0(i,j) 表示有多少个 i有前导 0 的数,其数位和为 j,并严格小于由字符串 d1d2di 表示的数。令状态 s0(i,j) 表示 c0(i,j) 中的所有数之和,t0(i,j) 表示 c0(i,j)。中的所有数的平方和。

那么可以先写出 c0 的状态转移方程:

c0(i,j)={1ifj=00else ifi=1c0(i1,j1)+c0(i1,j)+[di=1]c1(i1,j)else

其中,[] 表示示性函数,如果 [] 内的表达式为真,那么值为 1,否则为 0。方程的前两项,如果之前的状态已经是严格小于 d1d2di1 的数,那么接下来无论在后面拼接一个 0 还是一个 1,都是严格小于的;对于第三项而言,如果之前取的值都是严格等于的,那么接下来如果 di=1,那么可以拼接一个 0,从而达到严格小于的情况。

根据上面的思想,不难写出 s0(i,j) t0(i,j)

s0(i,j)={0ifi=1j=02s0(i1,j1)+c0(i1,j1)+2s0(i1,j)+[di=1]2s1(i1,j)else

t0(i,j)={0ifi=1j=04t0(i1,j1)+4s0(i1,j1)+c0(i1,j1)+4t0(i1,j)+[di=1]2s1(i1,j)else

需要注意的是,对于 t0 的状态转移方程最后一行的前三项,系数是 4,4,1,来自于平方式 (2n+1)2=4n2+4n+1

计算完成后,得到最终答案为:

S(n)=i=1li2(t0(l,i)+t1(l,i))

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
N = 10 ** 16
mod = 1000000007
d = list(map(int, list(bin(N)[2:])))
l = len(d)
c = [[[0, 0] for _ in range(l + 1)] for _ in range(l)]
s = [[[0, 0] for _ in range(l + 1)] for _ in range(l)]
t = [[[0, 0] for _ in range(l + 1)] for _ in range(l)]
c[0][1][1] = s[0][1][1] = t[0][1][1] = 1
c[0][0][0] = 1
cnt = 1
for i in range(1, l):
cnt += d[i]
c[i][cnt][1] = 1
s[i][cnt][1] = s[i - 1][cnt - d[i]][1] * 2 + d[i]
# 因为d[i]=0,1,d[i] ** 2 = d[i]
t[i][cnt][1] = t[i - 1][cnt - d[i]][1] * 4 + s[i - 1][cnt - d[i]][1] * d[i] * 4 + d[i]
c[i][0][0] = 1
s[i][0][0] = t[i][0][0] = 0
for j in range(1, l + 1):
c[i][j][0] = c[i - 1][j - 1][0] + c[i - 1][j][0]
s[i][j][0] = s[i - 1][j - 1][0] * 2 + c[i - 1][j - 1][0] + s[i - 1][j][0] * 2
t[i][j][0] = t[i - 1][j - 1][0] * 4 + s[i - 1][j - 1][0] * 4 + c[i - 1][j - 1][0] + t[i - 1][j][0] * 4
if d[i] == 1:
c[i][j][0] += c[i - 1][j][1]
s[i][j][0] += s[i - 1][j][1] * 2
t[i][j][0] += t[i - 1][j][1] * 4

ans = 0
for i in range(l + 1):
ans += sum(t[-1][i]) * i * i
ans %= mod
print(ans)