hdu3307 欧拉函数

 题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3307

Description has only two Sentences

Time Limit: 3000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 1733    Accepted Submission(s): 525

Problem Description

a_n = X*a_n-1 + Y and Y mod (X-1) = 0.
Your task is to calculate the smallest positive integer k that a_k mod a₀ = 0.

 

Input

Each line will contain only three integers X, Y, a₀ ( 1 < X < 2³¹, 0 <= Y < 2⁶³, 0 < a₀ < 2³¹).

 

Output

For each case, output the answer in one line, if there is no such k, output "Impossible!".

 

Sample Input

2 0 9

 

Sample Output

1

 

 a_n = X*a_n-1 + Y ，给你X,Y,a₀，叫你求出最小的整数k(k>0)使得 a_k% a₀=0。根据公式我们可以求出a_n= a₀*X^n-1+(x⁰+x¹+x²+……+x^n-1)Y。由等比数列前项和公式可得

 a_n=a₀*X^n-1+(xⁿ-1)*Y/(x-1)。

所以只需令(x^k-1)*Y/(x-1)%a₀=0，求出k的值。

令Y'=Y/(x-1),d=gcd(Y/(x-1),a₀)

Y'/=d,a₀/=d;

此时Y'与a₀互质

(x^k-1)*Y/(x-1)%a₀=0

等价于

(x^k-1)%a₀=0

x^k%a₀=1

而这就符合欧拉定理a^φ(n)Ξ 1(mod n)

x^φ(a₀)Ξ 1(mod a₀)

如果gcd(x,a₀)!=1则k无解

否则k为phi(a₀)除以满足条件的phi(a₀)的因子

#include<iostream>
using namespace std;
#define ll long long
ll zys[1000][2],cnt;
ll gcd(ll a,ll b)
{
    return b?gcd(b,a%b):a;
}
ll phi(ll x)
{
    ll ans=x;
    for(ll i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            ans=ans/i*(i-1);
            while(x%i==0)x/=i;
        }
    }
    if(x>1)
    ans=ans/x*(x-1);
    return ans;
}
void fj(ll ans)
{
    for(ll i=2;i*i<=ans;i++)
    {
        if(ans%i==0)
        {
            zys[cnt][0]=i;
            zys[cnt][1]=0;
            while(ans%i==0)
            {
                ans/=i;
                zys[cnt][1]++;
            }
            cnt++;
        }
    }
    if(ans>1)
    {
        zys[cnt][0]=ans;
        zys[cnt++][1]=1;
    }
}
ll poww(ll a,ll b,ll mod)
{
    ll ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
int main()
{
    ll x,y,a,c,d;
    while(cin>>x>>y>>a)
    {
        if(x==1)
        {
            if(gcd(y,a)==1)cout<<a<<endl;
            else cout<<a/gcd(y,a);
            continue;
        }
        c=y/(x-1);
        if(c%a==0)
        {
            cout<<1<<endl;
            continue;
        }
        d=gcd(c,a);
        if(gcd(x,a/d)!=1)cout<<"Impossible!"<<endl;
        else
        {
            a/=d;
            ll ans=phi(a);
            cnt=0;
            fj(ans);
            for(int i=0;i<cnt;i++)
            {
                for(int j=1;j<=zys[i][1];j++)
                {
                    if(poww(x,ans/zys[i][0],a)==1)ans/=zys[i][0];
                }
            }
            cout<<ans<<endl;
        }
    }
    return 0;
}