最大公约数 [动态规划]

最大公约数


color{red}{正解部分}

最优情况中, gcdgcd 每次变化都是 质因数分解中的一个质因子次数减少 11, 所以 fmax(n)f_{max}(n)[1,n][1, n] 中 质因数分解质数次幂的和 的 最大值 加上 11.

首先 fmax(n)f_{max}(n) 分解质因数后一定是 2x3y2^x3^y 的形式, 其中 xN,y{0,1}x in N, y in {0, 1}
假设存在 t>3t > 3, 使得 fmax(n)f_{max}(n) 分解质因数后为 t2x3yt2^x3^y, 则 tt 可以替换为 2?2^?,答案不会更差, 因此 fmax(n)=log2n+1f_{max}(n) = lfloor log_2n floor + 1 .

接下来计算 f(p)=fmax(n)f(p) = f_{max}(n)pp 的个数, 考虑 从左往右 填数, 前面的数字一定是 2x3y2^x3^y 的形式,
F[i,x,y]F[i, x, y] 表示 填了前 ii 个数字, 前 ii 个数字的 gcdgcd2x3y2^x3^y 的方案数,

初值

  • F[1,log2N,0]=1F[1, lfloor log_2 N floor, 0] = 1
  • F[1,log2N1,1]=1              (2log2N1×3N)F[1, lfloor log_2 N floor - 1, 1] = 1 (2^{lfloor log_2 N floor - 1} imes3 le N)

转移 时使得 xxyy 至多其中一个减一,

  • x,yx, y 不变: F[i,x,y] +=F[i1,x,y]×(N2x3y(i1))F[i, x, y] += F[i-1, x, y] imes left( lfloor frac{N}{2^x3^y} floor -(i-1) frac{}{} ight)

  • xx11:     F[i,x,y] +=F[i1,x+1,y]×(N2x3yN2x+13y) F[i, x, y] + = F[i-1, x+1,y] imes left( lfloor frac{N}{2^x3^y} floor - lfloor frac{N}{2^{x+1}3^y} floor ight)

  • yy11:     F[i,x,y] +=F[i1,x,y+1]×(N2x3yN2x3y+1) F[i, x, y] + = F[i-1, x, y+1] imesleft( lfloor frac{N}{2^x3^y} floor - lfloor frac{N}{2^x3^{y+1}} floor ight)

最后答案即为 F[N,0,0]F[N, 0, 0] .


color{red}{实现部分}

#include<bits/stdc++.h>
#define reg register

const int mod = 1e9 + 7;
const int maxn = 1e6 + 5;

int N;
int F[maxn][22][2];

int Ksm(int a, int b){ int s=1; while(b){ if(b&1)s=1ll*s*a%mod; a=1ll*a*a%mod;b>>=1; } return s; }

int main(){
        scanf("%d", &N); int t = log2(N);
        F[1][t][0] = 1; if(Ksm(2, t-1)*3 <= N) F[1][t-1][1] = 1; 
        for(reg int i = 2; i <= N; i ++)
                for(reg int x = 0; x <= t; x ++)
                        for(reg int y = 0; y <= 1; y ++){
                                int &cur = F[i][x][y], t1 = Ksm(2,x)*Ksm(3,y), t2 = Ksm(2,x+1)*Ksm(3,y), t3 = Ksm(2,x)*Ksm(3,y+1);
                                if(N/t1 > i-1) cur = (cur + 1ll*F[i-1][x][y]*(N/t1-i+1)%mod) % mod;
                                if(x != t) cur = (cur + 1ll*F[i-1][x+1][y]*(N/t1 - N/t2)%mod) % mod;
                                if(y != 1) cur = (cur + 1ll*F[i-1][x][y+1]*(N/t1 - N/t3)%mod) % mod;
                        }
        printf("%d
", F[N][0][0]);
        return 0;
}
原文地址:https://www.cnblogs.com/zbr162/p/11822426.html