hdu 5407 CRB and Candies（素数筛选法，除法取模（乘法逆元））

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=5407

解题思路：

官方题解：

The problem is just to calculate $g(N)=LCM(C(N,0),C(N,1),...,C(N,N))$.

Introducing function $f(n)=LCM(1,2,...,n)$, the fact $g(n)=f(n+1)/(n+1)$ holds.

We calculate $f(n)$ in the following way.

$f(1)=1$.

If n\ =p^{k}n =p​k​​ then $f(n)=f(n-1)\times p$, else $f(n)=f(n-1)$.

Time complexity:$O(N\cdot logN)$

如果不会做，可以取个巧，已知前n项，可以在这个网站找找规律再做。。。http://oeis.org/?language=english

• 求
  LCM((n0),(n1)…(nn))
• 令
  g(n)=LCM((n0),(n1)…(nn))
  以及
  f(n)=LCM(1,2,…n)
  则有
  g(n)=f(n+1)n+1

• 显然 f(1)=1 ，考虑 f(n) 的质因数分解形式，知道

  f(n)=⎧⎩⎨f(n−1) * p,f(n−1),if n=pkif n!=pk
  其中 p 表示一个质数，即判断 n 是否能表达成某个质数的整数 k 次方.

• 最后的除以 n+1 部分需要用逆元处理一下

证明过程如下：

http://www.zhihu.com/question/34859879/answer/60168919

http://arxiv.org/pdf/0906.2295v2.pdf

首先用素数筛选法求出1～1000000之间的素数，然后再计算f[n]，但是计算f(n)的时候不能直接除以(n+1)，应乘以其逆元. 

下面有两种求逆元的模板：

1.利用扩展欧几里德算法求逆元

ll extend_gcd(ll a,ll b,ll &x,ll &y){
    if(a == 0 && b == 0)
        return -1;//无最大公约数
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    ll d = extend_gcd(b,a%b,y,x);
    y -= a/b*x;
    return d;
}

//*********求逆元素*******************
//ax = 1(mod n)
ll mod_reverse(ll a,ll n)
{
    ll x,y;
    ll d = extend_gcd(a,n,x,y);
    if(d == 1)
        return (x%n+n)%n;
    else
        return -1;
}

2.利用快速幂求逆元

ll pow_mod(ll x, int n) {
    ll ret = 1;
    while (n) {
        if (n&1)
            ret = ret * x % MOD;
        x = x * x % MOD;
        n >>= 1;
    }
    return ret;
}

ll mod_reverse(ll x) {
    return pow_mod(x, MOD-2);
}

AC代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;

typedef long long ll;
const int MOD = 1000000007;
const int N = 1000005;
ll f[N];
int nprime;
int vis[N];
int prime[80000];

void getprime(){
    nprime = 0;
    memset(vis,0,sizeof(vis));
    memset(prime,0,sizeof(prime));
    for(int i = 2; i <= N-5; i++){
        int t = (N-5)/i;
        for(int j = 2; j <= t; j++){
            vis[i*j] = 1;
        }
    }
    for(int i = 2; i <= N-5; i++){
        if(!vis[i])
            prime[nprime++] = i;
    }
    memset(vis,0,sizeof(vis));
    for(int i = 0; i < nprime; i++){
        ll a = prime[i];
        ll b = a;
        for(; a < N; a*=b)
            vis[a] = b;
    }
}

void init(){
    getprime();
    f[1] = 1;
    for(int i = 2; i <= N-4; i++){
        if(vis[i])
            f[i] = f[i-1]*vis[i]%MOD;
        else
            f[i] = f[i-1];
        f[i] %= MOD;
    }
}
/*
ll extend_gcd(ll a,ll b,ll &x,ll &y){
    if(a == 0 && b == 0)
        return -1;//无最大公约数
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    ll d = extend_gcd(b,a%b,y,x);
    y -= a/b*x;
    return d;
}

//*********求逆元素*******************
//ax = 1(mod n)
ll mod_reverse(ll a,ll n)
{
    ll x,y;
    ll d = extend_gcd(a,n,x,y);
    if(d == 1)
        return (x%n+n)%n;
    else
        return -1;
}
*/

ll pow_mod(ll x, int n) {
    ll ret = 1;
    while (n) {
        if (n&1)
            ret = ret * x % MOD;
        x = x * x % MOD;
        n >>= 1;
    }
    return ret;
}

ll mod_reverse(ll x) {
    return pow_mod(x, MOD-2);
}

int main(){
    init();
    int T;
    scanf("%d",&T);
    while(T--){
        int n;
        scanf("%d",&n);
        //ll ni = mod_reverse(n+1,MOD);
        ll ni = mod_reverse(n+1);
        printf("%lld\n",f[n+1]*ni%MOD);
    }
    return 0;
}