POJ 3696 The Luckiest number (整数的阶+欧拉定理)

题目：http://poj.org/problem?id=3696

题意：给你一个数L，求出只由数字8组成的数的最小长度x，使长度为x的只由8组成的数能整除L，无解输出0。

分析：只有8组成的数字可以表示为(10^x-1)/9*8。那么可以得到等式

    （10^x-1）/9*8=L*k,

   ->(10^x-1)*8=9*L*k

   ->(10^x-1)*8/gcd(9*L,8)=9*L/gcd(9*L)*k  而8/gcd(9*L,8)与9*L/gcd(9*L,8)互质

   ->9*L/gcd(9*L,8)|(10^x-1)

   ->(10^x-1)=9*L/gcd(9*L,8)*k'   令p=9*L/gcd(9*L,8)

   ->(10^x-1)=p*k'

   ->(10^x-1)=0 (mod p)

   ->10^x=1  (mod p)

   若同余式有解，则10^(x-1)*10=1 (mod p) 即gcd(10,p)==1 ，无解就是gcd(10,p)!=1。

   那么怎么求出最小的x？

   这里介绍整数的阶的定义及相关定理。

定义：

  设a和n是互质的正整数，使得a^x≡1(mod n)成立的最小的正整数x称为a模n的阶，记a模n的阶为ord n(a)。

定理：

  如果a和n是互质的整数且n>0，那么正整数x是同余式a^x≡1(mod n)的一个解,当且仅当(ord n(a))|x。

推论：

  如果a和n是互质的整数，且n>0，那么ord n(a)|φ(n)。

  那么题目就转变为求φ(p)的最小因子α，使得10^α=1 (mod p)。

  如何找出φ(p)的因子？在1~sqrt(φ(p))内枚举就行了。

代码：

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <vector>

using namespace std;
typedef long long LL;
const int maxn=50001;

LL nprime,prime[maxn];
bool isprime[maxn];

LL gcd(LL a,LL b)
{
	return b==0?a:gcd(b,a%b);
}

void doprime()
{
	nprime=0;
	LL i,j;
	for(i=1;i<maxn;i+=2)
		isprime[i]=true;
	isprime[1]=false;
	isprime[2]=true;
	for(i=2;i<maxn;i++)
		if(isprime[i])
		{
			prime[nprime++]=i;
			for(j=i*i;j<maxn;j+=i)
				isprime[j]=false; 
		} 
}

LL euler(LL L)
{
	LL i,j,t=L,ans=L;
	for(i=0;i<nprime && prime[i]*prime[i]<=L;i++)
	{
		if(t%prime[i]==0)
		{
			ans=ans-ans/prime[i];
			while(t%prime[i]==0)
				t/=prime[i];
		}
		if(t==1)
			break;
	}
	if(t>1)
		ans=ans-ans/t;
	return ans;
}

LL multi(LL a,LL b,LL mod)
{
	LL ans=0;
	while(b)
	{
		if(b&1)
		{
			ans=ans+a;
			if(ans>=mod)
				ans-=mod;
		}
		b>>=1;
		a=a<<1;
		if(a>=mod)
			a-=mod;
	}
	return ans;
}

LL my_pow(LL a,LL n,LL mod)
{
	LL ans=1;
	while(n)
	{
		if(n&1)
			ans=multi(ans,a,mod);
		n>>=1;
		a=multi(a,a,mod);
	}
	return ans;
}

int main()
{
	vector <LL > vt;
	doprime();
	LL L,oula,p,ncase=1;
	while(scanf("%lld",&L),L)
	{
		printf("Case %lld: ",ncase++);
		p=9ll*L/gcd(9ll*L,8);
		if(gcd(p,10)!=1)
		{
			printf("0\n");
			continue;
		 } 
		oula=euler(L);
		vt.clear();
		for(LL i=1;i*i<=oula;i++)
			if(oula%i==0)
			{
				vt.push_back(i);
				vt.push_back(oula/i);
			}
		sort(vt.begin(),vt.end());
		for(LL i=0;i<vt.size();i++)
			if(my_pow(10,vt[i],p)==1)
			{
				printf("%lld\n",vt[i]);
				break;
			 } 
	} 
	return 0;
}