hdu5525 Product 费马小定理

Product

 

 Accepts: 21

 

 Submissions: 171

 Time Limit: 6000/3000 MS (Java/Others)

 

 Memory Limit: 131072/131072 K (Java/Others)

问题描述

给n个数{A}_{1},{A}_{2}....{A}_{n}A​1​​,A​2​​....A​n​​，表示N=\prod_{i=1}^{n}{i}^{{A}_{i}}N=∏​i=1​n​​i​A​i​​​​。求N所有约数之积。

输入描述

输入有多组数据.
每组数据第一行包含一个整数n.(1\leq n\leq {10}^{5})(1≤n≤10​5​​)
第二行n个整数{A}_{1},{A}_{2}....{A}_{n}A​1​​,A​2​​....A​n​​,保证不全为0.(0\leq {A}_{i}\leq {10}^{5})(0≤A​i​​≤10​5​​).
数据保证 $\sum n\le 500000$.

输出描述

对于每组数据输出一行为答案对{10}^{9}+710​9​​+7取模的值.

输入样例

4
0 1 1 0
5
1 2 3 4 5

输出样例

36
473272463

把N化成 N=\prod_{i=1}^{k}{{p}_{i}}^{{a}_{i}}N=∏​i=1​k​​p​i​​​a​i​​​​ ，其中p为互不相等的质数，则含 {p}_{x}p​x​​ 个数为j的约数有 \frac{\prod_{i=1}^{k}({a}_{i}+1)}{{a}_{x}+1}​a​x​​+1​​∏​i=1​k​​(a​i​​+1)​​ 个,而j的取值范围是0到 {a}_{x}a​x​​ ,从而得到 {p}_{x}p​x​​ 在所有约数中出现了 \frac{(1+{a}_{x}){a}_{x}}{2}*\frac{\prod_{i=1}^{k}({a}_{i}+1)}{{a}_{x}+1}​2​​(1+a​x​​)a​x​​​​∗​a​x​​+1​​∏​i=1​k​​(a​i​​+1)​​ 次。根据费马小定理 {a}^{p-1}\equiv 1(mod p)a​p−1​​≡1(modp) ,统计 {a}_{x}a​x​​ 时可以对p-1取模，由于后面还要除2，所以先对2(p-1)取模。求 \frac{\prod_{i=1}^{k}({a}_{i}+1)}{{a}_{x}+1}​a​x​​+1​​∏​i=1​k​​(a​i​​+1)​​ 部分时不能取逆元，可以分别对前缀和后缀计算积。计算出N包含的所有质因子出现次数后用快速幂统计一遍，复杂度 $O(NlogN)$ .

#include<iostream>
#include<cstdio>
#include<bitset>
#include<vector>
#include<queue>
#include<set>
#include<map>
#include<cstring>
#include<cmath>
#include<algorithm>
#define ll long long
#define ull unsigned long long
#define debug puts("======");
#define inf (1<<30)
#define eps 1e-8
using namespace std;
const int maxn=100010;
ll mod=1000000007;
ll M=1000000006;
int prime[maxn];
bool vis[maxn];
int cnt;
void getPrime()
{
    int N=100000;
    int m=sqrt(N+0.5);
    memset(vis,0,sizeof(vis));
    for(int i=2;i<=m;i++) {
        if(vis[i]==0) {
            for(int j=i*i;j<=N;j+=i)
                vis[j]=1;
        }
    }
    cnt=1;
    for(int i=2;i<=N;i++) {
        if(!vis[i])
            prime[cnt++]=i;
    }
}
int n;
ll num[maxn];
int pos[maxn];
ll lef[maxn];
ll rig[maxn];
ll cal(ll a)
{
    if(a%2==0)
        return a/2%M*((a+1)%M)%M;
    else
        return (a+1)/2%M*(a%M)%M;
}
ll powMod(ll a,ll b)
{
    if(b==0) return 1;
    ll ans=powMod(a,b/2);
    ans=ans*ans%mod;
    if(b&1)
        ans=ans*a%mod;
    return ans;
}
int main()
{
    getPrime();
    for(int i=1;i<cnt;i++)
        pos[prime[i]]=i;
    int a;
    while(~scanf("%d",&n)) {
        memset(num,0,sizeof(num));
        for(int i=1;i<=n;i++) {
            scanf("%d",&a);
            int u=i;
            for(int j=1;j<cnt;j++) {
                if((ll)prime[j]*prime[j]>u)
                    break;
                while(u%prime[j]==0) {
                    u/=prime[j];
                    num[j]+=a;
                    num[j]%=(2*M);
                }
            }
            if(u>1) {
                num[pos[u]]+=a;
                num[pos[u]]%=(2*M);
            }
        }
        int m;
        for(int i=cnt-1;i>=1;i--) {
            if(num[i]>0) {
                m=i;
                break;
            }
        }
        lef[0]=rig[m+1]=1;
        for(int i=1;i<=m;i++)
            lef[i]=lef[i-1]*(num[i]+1)%M;
        for(int i=m;i>=0;i--)
            rig[i]=rig[i+1]*(num[i]+1)%M;
        ll ans=1;
        for(int i=1;i<=m;i++) {
            ll tot=cal(num[i]);
            tot=tot*lef[i-1]%M*rig[i+1]%M;
            ans=ans*powMod(prime[i],tot)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}