HDU 5407 CRB and Candies 组合数+逆元

【写于2016-08-05，现在从hexo搬运过来】

对于正整数n，计算

n是1e6级别的

做法比较多的，不过都是一个道理。
有的是用kummer定理，那个看得挺奥义的。。
其实这个基于一个公式就可以直接做了

打表功力足够的话可能可以打表看出这个结论，不过。。好像不容易
这里有一篇论文证明这个结论，分析性比较强：这里
这样我们就转化为求1e6规模下的多个数lcm的问题了

#####我们先考虑任意n个数字的lcm，其最大值为N
每个数都有若干质因子，现假定这n个数一共涉及到了若干有限个质数因子，且这些质因子都不大于N，进而最小公倍数等于全体出现一次的质因子之积全体出现两次的质因子之积… 全体出现了n次的质因子之积
有点绕。。换句话说，对于一个质因子p，考虑最大的正整数k，s.t. p^k<=n+1
答案乘上p^k，对所有N以内的质数因子，都这么做一遍，最后的答案即为lcm

#####对于1~n+1的lcm
n+1以内的所有质数都会以质因子形式出现，按上述做法做一遍即可了

当然最后还得除以n+1,取模的意义下，乘其逆元即可，这里采用的ex_gcd求逆元
以后总结一下逆元的3种写法。。目前好像就知道三种

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<algorithm>
#include<functional>
#include<iostream>
#include<cmath>
#include<string>
#include<cctype>
#include<stack>
#include<queue>
#include<set>
#include<sstream>
#include<map>
#include<ctime>
using namespace std;
#define LL long long
#define ULL unsigned long long
//#define mod 100
#define ms(a,x) memset(a,x,sizeof(a));
#define lg2 0.30102999566
const double eps=1e-2;
const double epss=1e-5;
const double pi=acos(-1);
const int maxm = 100000+10;
const int maxn = 1000000+10;
const int inf = 0x7f7f7f7f;
const int mod = 1e9+7;
//HDU 5407 solution 1
//ex_gcd_inverse & a theory.
LL prime[maxn];
bool check[maxn]={false};
int index;
void init()//get prime in 1~1e6
{
    index=0;
    for(int i=2;i<=maxn;i++)
    {
        if(!check[i])prime[index++]=i;
        for(int j=i*2;j<=maxn;j+=i)check[j]=true;
    }
}
LL ex_gcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0){
        x=1,y=0;
        return a;
    }
    LL d=ex_gcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
LL inv(LL a,LL n)//a 's inverse mod n;
{
    LL x,y;
    LL gcd=ex_gcd(a,n,x,y);
    return (x%n+n)%n;
}
LL quick_mod(LL a,LL b)
{
    LL ans=1;
    while(b)
    {
        if(b&1) ans=(ans*a)%mod;
        a=(a*a)%mod;
        b>>=1;
    }
    return ans;
}
int T;
LL n;
int main()
{
    init();
    //for(int i=0;i<15;i++)printf("%d ",prime[i]);
    scanf("%d",&T);
    while(T--)
    {
        scanf("%I64d",&n);
        int pos=upper_bound(prime,prime+index,n)-prime+1;
        LL ans=1;
        //cout<<"***"<<pos<<endl;
        for(int i=0;i<pos;i++)
        {
            LL p=(LL)(log(n+1)/log(prime[i]));
            //cout<<"**"<<p<<endl;
            ans=(ans*quick_mod(prime[i],p)%mod)%mod;
        }
        ans=ans*inv(n+1,mod) %mod;
        printf("%I64d\n",ans);
    }
}