bzoj3209 花神的数论题

http://www.lydsy.com/JudgeOnline/problem.php?id=3209

题意：求sum(1)*sum(2)*...*sum(n)，其中sum(x)表示x的二进制表达中1的数量。答案模10^7+7。n<=10^15。

考虑枚举sum(x)的值为t，然后求有多少个不大于n的正整数的二进制表达恰好有t个1，设这个答案为f(t)。

这一部分可以用数位dp来实现。例如要求解1...12，可以将其分成1...7,8...11,12...12这三层。每一层中，满足范围内的所有数的前缀若干位固定，而后缀若干位（不妨设位len位）可以取遍所有2^len种状态，即后缀这些位的每一位均可取0或1。设固定的前缀已经包含了k个1，则这一层对f(i)(k<=i<=len+k)的贡献是C(len,i-k)，其中C()是组合数。

对每层分别计算贡献后，答案即为1^f(1)*2^f(2)*...*m^f(m)，其中m是n的二进制表达的位数。这一步用快速幂实现。

需要注意的是，由于f()的值是作为乘方运算中的指数来对答案做贡献，因此不能对10^7+7取模，而应该对phi(10^7+7)=9988440取模。预处理组合数时也要注意应该对9988440取模。

时间复杂度O((logn)^2)。

代码：

#include<cstdio>
#include<cstring>
#define rpt(i,l,r) for(i=l;i<=r;i++)
#define rpd(i,r,l) for(i=r;i>=l;i--)
#define ll long long
#define p 10000007
#define ph 9988440
int C[50][50];
int i,j,k;
int f[50];
ll n,ans;
ll mi(ll x,ll y){
    if(y==0) return 1;
    ll temp=mi(x,y>>1);
    (temp*=temp)%=p;
    if(y&1) (temp*=x)%=p;
    return temp;
}
int main(){
    rpt(i,0,49){
        C[i][0]=C[i][i]=1;
        rpt(j,1,i-1) C[i][j]=(C[i-1][j-1]+C[i-1][j])%ph;
        rpt(j,i+1,49) C[i][j]=0;
    }
     
    scanf("%lld",&n);
    n++;k=0;
    memset(f,0,sizeof(f));
    rpd(i,49,0){
        if(n&(1LL<<i)){
            rpt(j,k,49) (f[j]+=C[i][j-k])%=ph;
            k++;
        }
    }
    ans=1;
    rpt(i,1,49) (ans*=mi(i,f[i]))%=p;
    printf("%d",(int)ans);
}