Luogu | P1762 偶数 | 数论 / 动态规划 | Lucas 定理 + 数位

Description

给定一个正整数 $N$ ，请输出杨辉三角形前 $N$ 行的偶数个数$\mod 1000003$ 的结果。

Input Format

一行一个正整数 $N$ 。

Output Format

一行一个正整数，表示结果。

Input Sample

Sample 01

6

Output Sample

Sample 01

6

Data Range

$ 1 \le N \le 10^{15}$

Solution

根据 Lucas 定理，当 p 为质数时有：$\binom{n}{m} \equiv \binom{n \div p}{m\div p} \times \binom{n \mod p}{m \mod p} (\mod p)$。

这里的 p = 2 ，那么可以发现最后 $\binom{n}{m} \mod 2$ 的值其实为 n 和 m 二进制按位分解的组合数的乘积 $\mod 2$ 。

又因为有

$$\binom{0}{0} = 1$$

$$\binom{0}{1} =0$$

$$\binom{1}{0} = 1$$

$$\binom{1}{1} = 1$$

那么可以发现，在上述四式中，m & n == m 等价于 $\binom{n}{m} = 1$。

这意味着，在原式中有 m & n ==m 等价于 $\binom{n}{m} \mod 2 =1$。

则杨辉三角形的第 $x$ 行的奇数个数为$\sum_{i=1}^{x} ((i \ \& \ x) == i) $。对于 $x$ ，对 $x$ 进行二进制分解，假设可以得到 $k$ 个 1 ，那么这 $k$ 个位置的 1 对应的组合数中的 $m$ 的该位上的数字可以为 1 或者 0 两种选法，可得 $x$ 对答案的贡献为 $(\sum_{i=1}^{x} ((i \ \& \ x) == i)) = 2^k$ 。

考虑计算出 ≤ N 的所有数字对于答案的贡献。这个过程可以用 数位 dp 计算出 $\sum$ **≤ N 的数字中所有数字的二进制分解的 1 的个数为 k 的数字的个数 × 2k ** 完成。

注意中途的 long long。

Code

#include<cstdio>
#include<cstdlib>
#include<cstring>

#define min(x,y) ((x)<(y)?(x):(y))
#define max(x,y) ((x)>(y)?(x):(y))
#define swap(x,y) {int t=x; x=y,y=t;}
#define wipe(x,y) memset(x,y,sizeof(x))
#define dbgIn(x) freopen(x".in","r+",stdin)
#define rep(x,y,z) for(int x=y,I=z;x<=I;++x)
#define dbgOut(x) freopen(x".out","w+",stdout)
#define file(x) freopen(x".in","r+",stdin); freopen(x".out","w+",stdout)

#define mo 1000003

inline void Read(int &x){
    x=0; char ch=0; while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-48,ch=getchar(); return;
}

int len;

long long N;

int tmp[65];
int res[65];
int sum[65];

long long F[65][65];

void Decompose(long long X){
    while(X){
        tmp[++len]=X%2;
        X/=2;
    }
    rep(i,1,len)
        res[i]=tmp[len-i+1];
    return;
}

void Prec(){
    F[0][0]=1;
    rep(i,1,63){
        rep(j,0,i){
            F[i][j]=F[i-1][j];
            if(j>=1)
                F[i][j]=(F[i][j]+F[i-1][j-1])%mo;
        }
    }
}

long long Calc(long long X){
    long long ret=0;
    Decompose(X);
    rep(i,1,len)
        sum[i]=sum[i-1]+res[i];
    rep(i,1,len){
        rep(j,0,res[i]-1){
            rep(k,0,len-i)
                ret=(ret+(F[len-i][k]%mo)*((long long)1<<(k+sum[i-1]))%mo)%mo;
        }
    }
    return ret;
}

int main(){
    Prec();
    scanf("%lld",&N);
    printf("%lld\n",((((((N+1)%mo)*(N%mo)/2)-Calc(N)))%mo+mo)%mo);
    return 0;
}