小 C 的数学题

题目描述

已知某函数 f(n)f(n)的定义如下：

f(n)={1,n≤2f(n−1)+f(n−2)+⌊log2n⌋,n≥3{\displaystyle f(n)={\begin{cases}1,&{n\le2}\\ f(n-1)+f(n-2)+\lfloor log_2n\rfloor,&{n\ge 3}\end{cases}}}

求函数 f(n)f(n)的值，为了防止数据过大，你只需要输出答案对 109+710^9+7 取模之后的结果。

输入

最多十行，每行一整数 nn，1≤n<1091\le n < 10^9

输出

f(n)f(n)

样例输入

666
666666
666666666

样例输出

581824003
310868947
609519880

分析

乍一看是矩阵快速幂，但是随着 nn 的增大，⌊log2n⌋\lfloor log_2n\rfloor也会随之改变。

\begin{pmatrix}f_n\\f_{n-1}\\1\end{pmatrix}={\begin{pmatrix}1&1&\lfloor log_2n\rfloor \\1&0&0\\0&0&1\end{pmatrix}\begin{pmatrix}f_{n-1}\\f_{n-2}\\1\end{pmatrix}

变化中蕴含着不变性，在区间 n∈[2p,2p+1)n\in [2^{p},2^{p+1}) 内，⌊log2n⌋\lfloor log_2n\rfloor 的值不变。因此我们只需要将 nn 分成若干个这样的区间，进行分段的矩阵快速幂即可。

代码

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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include<iostream>  
#include<cmath>
#include<cstring>  
using namespace std;
typedef long long ll;
const int N=4;
const ll mod=1e9+7;
ll n;
struct matrix//定义矩阵
{
	ll m[N][N];
	matrix()//构造函数 
	{
		memset(m,0,sizeof(m));
	}
};
//定义矩阵乘法
matrix operator *(const matrix &x,const matrix &y)
{
	matrix res;
	int i,j,k;
	for(i=1;i<=3;i++)
	{
		for(j=1;j<=3;j++)
		{
			for(k=1;k<=3;k++)
			{
				res.m[i][j]+=x.m[i][k]*y.m[k][j];
				res.m[i][j]%=mod;
			}
		}
	}
	return res;
}
//矩阵快速幂
matrix quick_matrix_pow_mod(matrix a,ll b)
{
	matrix res;
	int i;
	for(i=1;i<=3;i++) res.m[i][i]=1;//单位矩阵 
	while(b)
	{
		if(b&1) res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
matrix t,f;
void solve()
{
	//t为系数矩阵
	t.m[1][1]=1;
	t.m[1][2]=1;
	t.m[2][1]=1;
	t.m[3][3]=1;
	f.m[1][1]=f.m[2][1]=f.m[3][1]=1;
	ll i;
	for(i=3;i<=n;)
	{
		ll l,r,now=(ll)(log((double)i)/log((double)2));
		l=i;
		r=min(n,(1LL<<(now+1))-1);//分段
		t.m[1][3]=now;//更改系数矩阵
		matrix temp=quick_matrix_pow_mod(t,r-l+1);
		f=temp*f;//得到新的f
		i=r+1;//进入下一段
	}
	printf("%lld\n",f.m[1][1]);//输出
}
int main()
{ 
	while(~scanf("%lld",&n)) solve();
	return 0;
}