ACM必做50题的解题-数论

poj 1061 青蛙的约会

这题用的是解线性同余方程的方法，之前没接触到过，搜索资料后看到一个人的博客里面讲的很好就拷贝过来了。内容如下：
对于方程：ax≡b(mod m)，a，b，m都是整数，求解x 的值，这个就是线性同余方程。
符号说明：
mod表示：取模运算
ax≡b(mod m)表示：(ax - b) mod m = 0，即同余
gcd(a,b)表示：a和b的最大公约数
求解ax≡b(mod n)的原理：对于方程ax≡b(mod n)，存在ax + by = gcd(a,b)，x，y是整数。而ax≡b(mod n)的解可以由x，y来堆砌。具体做法如下：

第一个问题：求解gcd(a,b)

定理一：gcd(a,b) = gcd(b,a mod b)

欧几里德算法

int Euclid(int a,int b)
{
if(b == 0)
return a;
else
return Euclid(b,mod(a,b));
}
附：取模运算

int mod(int a,int b)
{
if(a >= 0)
return a % b;
else
return a % b + b;
}
第二个问题：求解ax + by = gcd(a,b)
定理二：ax + by = gcd(a,b)= gcd(b,a mod b) = b * x' + (a mod b) * y'
= b * x' + (a - a / b * b) * y'
= a * y' + b * (x' - a / b * y')
= a * x + b * y
则：x = y'

y = x' - a / b * y'

以上内容转自/redcastle/blog/item/934b232dbc40d336349bf7e4.html


由这个可以得出扩展的欧几里德算法：

int exGcd(int a, int b, int &x, int &y)
{
if(b == 0)
{
x = 1;
y = 0;
return a;
}
int r = exGcd(b, a % b, x, y);
int t = x;
x = y;
y = t - a / b * y;
return r;
}

代码：

#include
#include
#include
#include
using namespace std;

__int64 mm,nn,xx,yy,l;
__int64 c,d,x,y;

__int64 modd(__int64 a, __int64 b)
{
if(a>=0)
return a%b;
else
return a%b+b;
}

__int64 exGcd(__int64 a, __int64 b)
{
if(b==0)
{
x=1;
y=0;
return a;
}
__int64 r=exGcd(b, a%b);
__int64 t=x;
x=y;
y=t-a/b*y;
return r;
}

int main()
{
scanf("%I64d %I64d %I64d %I64d %I64d",&xx,&yy,&mm,&nn,&l);
if(mm>nn) //分情况
{
d=exGcd(mm-nn,l);
c=yy-xx;
}
else
{
d=exGcd(nn-mm,l);
c=xx-yy;
}
if(c%d != 0)
{
printf("Impossible\n");
return 0;
}
l=l/d;
x=modd(x*c/d,l); ///取模函数要注意
printf("%I64d\n",x);

system("pause");
return 0;
}



POJ 1142 SmithNumber

题意：寻找最接近而且大于给定的数字的SmithNumber

什么是SmithNumber？

用sum(int)表示一个int的各位的和，那一个数i如果是SmithNumber,则sum（i） = sigma( sum(Pj )),Pj表示i的第j个质因数。例如4937775= 3*5*5*65837，4+9+3+7+7+7+5 = 42，3+5+5+（6+5+8+3+7） = 42，所以4937775是SmithNumber。

思路：要寻找大于给定数字且最接近给定数字的SmithNumber，只要

将给定数字不断的加1，判断它是否是SmithNumber就行了，如果是SmithNumber就立即输出。

但是如何判断是否是SmithNumber呢？首先就是要对数字进行质因数分解。质因数分解要保证因子都是质数。这里先介绍一个如何判断一个数int是否是质数呢，如果对于这个数，i = 2.....sqrt(int)都不是它的约数，那int就是一个质数。所以我们可以将i初始化为2，然后不断递增i，看i是否是int的一个约数，如果是的话那i就是int的一个质因数（因为这个数是从2开始第一个可以整除int的数，它不可能是一个可以分解的合数，否则，它的约数在它之前就整除int），然后将int除以该质因数，重置i为2，重新对int判断它是否是质数。这样最后剩下的int一定是一个质数，从而对int进行了质因数分解

然后就很简单的将数字各质因数的各位加起来，看和是否等于该数字的各位和，如果相等那它可能就是SmithNumber，为什么说只是可能呢，因为这个数可能是质数，但是质数不是SmithNumber。



#include

#include
int Sum( int number )
{
int sum = 0;
while( number != 0 )
{
sum += number % 10;
number /= 10;
}
return sum;
}

bool SmithNumber( int number )
{
int i = 2;
int temp = number;
int sumOfNumber = Sum( number );
int sum = 0;
while( i <= (int)sqrt( (double)number ) )
{
if ( number % i ==0 )
{
sum += Sum( i );
number /= i;
i = 2;
}
else
{
++i;
}

//以上的代码做了无谓的计算，可用下面的代码，更新于20090904

//while( number % i == 0 )

//{

// sum += sum(i);

// number /= i;

//}

// ++i;
}
sum += Sum( number );
if ( sum == sumOfNumber && number != temp )
{
return true;
}
else
{
return false;
}
}


int main()
{
while( true )
{
int num;
scanf("%d",&num );
if ( num == 0 )
{
break;
}
while( true )
{
if ( SmithNumber(++num))
{
printf("%d\n", num);
break;
}
}
}
return 0;
}

ACM——POJ 2262（Goldbach's Conjecture）

题目地址：/JudgeOnline/problem?id=2262
题目思路：对于任何一个偶数 n，从 x = 1 和 y = n - 1 开始，看 x、y 是否是质数，不是则 x += 2, y += 2
这题需要开很大的内存，我 RE n 次，居然是因为数组开太小了，我这题耗时不是很理想，但我的耗内存
在我看到的中是最小的，所以看来 OJ 上的题只要能开内存的就尽量开。估计我这题用栈耗时了。

#include
#include
#include
#include
#include
#include

using namespace std;

// 判断是否为质数的函数
int IsPrime ( int x )
{
int i;
if( x < 2 )

return 0;

for( i = 2; i <= (int) ( sqrt( (double)x + 0.5 ) ); i++ )
if( x % i == 0)
return 0;

return 1;
}

int main()
{
// 输入数，输入数 / 2 向上延伸，输入数 / 2 向下延伸，输入数 / 2
int m_Input, m_Num_Max, m_Num_Min, m_InputToTwo;

// 总体输出
char m_Output[ 1000000 ];
memset( m_Output, 0, 1000000 );

// 标识 m_Output 的 Pos
int m_Output_Pos = 0;

// 是否找到标识
bool b_Find;

// 栈
stack m_Stack;

// 临时数
int m_Value_Top;

// 循环输入
while ( scanf( "%d", &m_Input ) && ( m_Input != 0 ) )
{
b_Find = true;

// m_Input 肯定是一个偶数
m_InputToTwo = m_Input / 2;

// 置值
m_Num_Max = m_Input - 1;
m_Num_Min = 1;

// 寻找，直至都为素数 或者 找不到 为止
while ( ( !IsPrime( m_Num_Max ) ) || ( !IsPrime( m_Num_Min ) ) )
{
// 否则，前进 & 后退 2 格
m_Num_Max -= 2;
m_Num_Min += 2;

// 如果发生如下情况，则输出 "Goldbach's conjecture is wrong."
if ( ( m_Num_Max < m_InputToTwo ) || ( m_Num_Min > m_InputToTwo ) )
{
char* m_TempChar = "Goldbach's conjecture is wrong.\n";
strcat( m_Output, m_TempChar );
b_Find = false;
m_Output_Pos += strlen( m_TempChar );
break;
}
}

// 如果找到了
if ( b_Find )
{
// 将 m_Input 转换为字符串存入 m_Output
while ( m_Input != 0 )
{
m_Value_Top = m_Input % 10;
m_Stack.push( m_Value_Top );
m_Input /= 10;
}

while ( !m_Stack.empty() )
{
m_Value_Top = m_Stack.top();
m_Output[ m_Output_Pos++ ] = m_Value_Top + 48;
m_Stack.pop();
}

// 加入 =
m_Output[ m_Output_Pos++ ] = ' ';
m_Output[ m_Output_Pos++ ] = '=';
m_Output[ m_Output_Pos++ ] = ' ';

// 将 m_Num_Min 转换为字符串存入 m_Output
while ( m_Num_Min != 0 )
{
m_Value_Top = m_Num_Min % 10;
m_Stack.push( m_Value_Top );
m_Num_Min /= 10;
}

while ( !m_Stack.empty() )

{
m_Value_Top = m_Stack.top();
m_Output[ m_Output_Pos++ ] = m_Value_Top + 48;
m_Stack.pop();
}

// 加入 +
m_Output[ m_Output_Pos++ ] = ' ';
m_Output[ m_Output_Pos++ ] = '+';
m_Output[ m_Output_Pos++ ] = ' ';

// 将 m_Num_Max 转换为字符串存入 m_Output
while ( m_Num_Max != 0 )
{
m_Value_Top = m_Num_Max % 10;
m_Stack.push( m_Value_Top );
m_Num_Max /= 10;
}

while ( !m_Stack.empty() )
{
m_Value_Top = m_Stack.top();
m_Output[ m_Output_Pos++ ] = m_Value_Top + 48;
m_Stack.pop();
}

// 加入 \n
m_Output[ m_Output_Pos++ ] = '\n';
}
}

// 输出
printf( "%s", m_Output );

system( "pause" );
return 0;
}


POJ 2407 Relatives

这题从题意可以看出就是求比从1~n - 1从有几个数和n没有公共因子， 通常的算法很简单就能够想到， 我开始也是按通常的做法写了一个， 觉得

可能会TLE, 果不其然， 后来上网查了一下， 知道了欧拉函数， 这个就是求比n小的数中与n互质（也就是没有公共因子）的算法， 看来还是那些经典的算法效率比较高， 比纯用暴力试探高多了...

欧拉函数描述如下：

利用欧拉函数和它本身不同质因数的关系，用筛法计算出某个范围内所有数的欧拉函数值。

欧拉函数和它本身不同质因数的关系：欧拉函数ψ（Ｎ）＝Ｎ｛∏ｐ｜Ｎ｝（１－１/ｐ）。（Ｐ是数Ｎ的质因数）

如：

ψ（１０）＝１０×（１－１/２）×（１－１/５）＝４；

ψ（３０）＝３０×（１－１/２）×（１－１/３）×（１－１/５）＝８；

ψ（４９）＝４９×（１－１/７）＝４２。

注意的是P是N的质因子， 这里求质因子还是不能够用常规的判断这个数是不是质数， 这样的话可能还会TLE, 网上学到他们用的一个while() 循环，感觉还挺巧的， 学习了...

#include
#include
int enlerFun(int n)
{
int count = n;
int i = 2;
for(; i<=n; i++)
if(n % i == 0)
{
count -= count / i;
while(n % i == 0)
n /= i;
}
return count;
}
int main()
{
int inputVal = 0;
int count = 0;

while(scanf("%d", &inputVal) && inputVal != 0)
{
count = enlerFun(inputVal);

printf("%d\n", cou

nt);
}

return 0;
}

POJ 1811 Prime Test

Miller Robin素性测试 + Pollard rho寻找素因子
Miller Robin 和 Pollard rho的理论想非常强，细节这里就不说了，可以参考
算法导论第31章

#include
#include
#include
#define MAX_L 64 //最长位数
#define TIMES 8 //miller robin素性测试的测试次数
#define MAX_VAL (pow(2.0, 60)) //定义最大值
#define CVAL 200
using namespace std;

//最小的素因子
__int64 minFactor;

//(1)计算a * b mod n, 思路: 利用b的二进制表示进行拆分计算
//(2)例如: b = 1011101那么a * b mod n = (a * 1000000 mod n + a * 10000 mod n + a * 1000 mod n + a * 100 mod n + a * 1 mod n) mod n
//(3)思路就是上面描述的那样, 那么可以用从低位往高位遍历b, 并用a来记录当前位为1的值，每次遇到b当前位为
//1就将结果值加上a并 mod n，然后a 要乘以2
__int64 multAndMod(__int64 a, __int64 b, __int64 n)
{
a = a % n;
__int64 res = 0;

while(b)
{
//当前位为1
if(b & 1)
{
//加上当前权位值
res += a;
//相当于mod n
if(res >= n) res -= n;
}
//乘以2，提高一位
a = a<<1;
//mod n
if(a >= n) a -= n;
b = b >> 1;
}

return res;
}

//(1)计算a ^ b mod n, 思路: 和上面类似，也是利用b的二进制表示进行拆分计算
//(2)例如: b = 1011101那么a ^ b mod n = [(a ^ 1000000 mod n) * (a ^ 10000 mod n) * (a ^ 1000 mod n) * (a ^ 100 mod n) * (a ^ 1 mod n)] mod n
//(3)思路就是上面描述的那样, 那么可以用从低位往高位遍历b, 并用a来记录当前位为1的值，每次遇到b当前位为
//1就将结果乘上a并 mod n，然后a 要乘以a以提升一位
__int64 modAndExp(__int64 a, __int64 b, __int64 n)
{
a = a % n;
__int64 res = 1;
while(b >= 1)
{
//遇到当前位为1，则让res * 当前a并mod n
if(b & 1)
res = multAndMod(res, a, n);
//a * a以提升一位
a = multAndMod(a, a, n);
b = b >> 1;
}
return res;
}

//MillerRobin素性测试，true:素数，flase:合数
bool millerRobin(__int64 a, __int64 n)
{
__int64 u = 0, cur = n - 1;
int t = 0;
bool find1 = false;
while(cur != 0)
{
if(!find1)
{
int pb = cur % 2;
if(pb == 0) t++;
else find1 = true;
}
if(find1)
break;
cur = cur / 2;
}
u = cur;

cur = modAndExp(a, u, n);
__int64 now;
for(int p = 1; p <=

t; p++)
{
now = modAndExp(cur, 2, n);
if(cur != 1 && now == 1 && cur != n - 1)
{
//printf("%d %d\n", cur, now);
return false;
}
cur = now;
}
if(cur != 1)
{
//printf("a:%I64d u:%I64d n:%I64d val:%I64d\n", a, u, n, start);
return false;
}
//printf("a:%I64d u:%I64d n:%I64d val:%I64d\n", a, u, n, start);
return true;
}

//利用Miller Robin对n进行n次素性测试
bool testPrime(int times, __int64 n)
{
if(n == 2) return true;
if(n % 2 == 0) return false;

__int64 a; int t;
srand(time(NULL));
for(t = 1; t <= times; t++)
{
a = rand() % (n - 1) + 1;
if(!millerRobin(a, n)) return false;
}
return true;
}

__int64 gcd(__int64 a, __int64 b)
{
if(b == 0) return (a);
return gcd(b, a % b);
}

__int64 PollardRho(__int64 n, int c)
{
int i = 1;
srand(time(NULL));
__int64 x = rand() % n;
__int64 y = x;
int k = 2;
while(true)
{
i = i + 1;
x = (modAndExp(x, 2, n) + c) % n;
__int64 d = gcd(y - x, n);
if(1 < d && d < n) return d;
if(y == x) return n; //重复了, 说明当前x下无解，需要重新启动PollardRho
if(i == k)
{
y = x;
k = k * 2;
}
}
}

void getSmallest(__int64 n, int c)
{
if(n == 1) return;
//判断当前因子是否为素数
if(testPrime(TIMES, n))
{
if(n < minFactor) minFactor = n;
return;
}
__int64 val = n;
//循环，知道找到一个因子
while(val == n)
val = PollardRho(n, c--);
//二分
getSmallest(val, c);
getSmallest(n / val, c);
}
int main()
{
int caseN;
__int64 n;
scanf("%d", &caseN);
while(caseN--)
{
scanf("%I64d", &n);
minFactor = MAX_VAL;
if(testPrime(TIMES, n)) printf("Prime\n");
else
{
getSmallest(n, CVAL);
printf("%I64d\n", minFactor);
}
}
return 0;
}



poj 2447 RSA

密码学:RSA公钥密码

#include
#include
#include
#include
using namespace std;
typedef unsigned __int64 u64;
#define MAX 30
#define MAXN 5
u64 len, dig, limit;
u64 factor[MAXN];
u64 mod(u64 a, u64 b, u64 n){
if(!a)return 0;
else return ( ((a&dig)*b)%n + (mod(a>>len,b,n)<}
u64 by(u64 a, u64 b, u64 n){
u64 p;
p = 8, len = 61;
while(pp<<=4;
len -=4;
}

dig = ((limit/p)<<1) - 1;
return mod(a,b,n);
}
u64 random(){//产生随机数
u64 a;
a = rand();
a *= rand();
a *= rand();
a *= rand();
return a;
}
u64 square_multiply(u64 x, u64 c, u64 n){//(x^c)%n快速算法
u64 z=1;
while(c){
if(c%2==1)z = by(z,x,n);
x = by(x,x,n);
c=(c>>1);
}
return z;
}
bool Miller_Rabin(u64 n){//Miller_Rabin素数测试
if(n<2)return false;
if(n==2)return true;
if(!(n&1))return false;
u64 k = 0, i, j, m, a;
m = n - 1;
while(m%2==0)m=(m>>1),k++;
for(i=0;ia = square_multiply(random()%(n-1)+1, m, n);//平方乘
if(a==1)continue;
for(j=0;jif(a==n-1)break;
a = by(a,a,n);
}
if(jreturn false ;
}
return true;
}
u64 gcd(u64 a,u64 b){
if(a==0) return b;
else return gcd(b%a,a);
}
u64 f(u64 x, u64 n ){
return (by(x,x,n)+1)%n;
}
u64 Pollard(u64 n){
if(n<=2)return 0;
if(!(n&1))return 2;
u64 i, p, x,xx;
for(i=0;ix = random()%n; //或者直接用 x = i
xx = f(x,n);
p = gcd((xx+n-x)%n , n);
while(p==1){
x = f(x,n);
xx = f(f(xx,n),n);
p = gcd((xx+n-x)%n,n)%n;
}
if(p)return p;
}
return 0;
}
u64 prime(u64 a){
if(Miller_Rabin(a))return a;
u64 t = Pollard(a);
if(t!=0)
{return prime(t);}

}
u64 Euclid(u64 a,u64 b,__int64 &x,__int64 &y)
{
if(b==0)
{
x=1,y=0;
return a;
}
u64 t,d;
if(b!=0)
d=Euclid(b,a%b,x,y);
t=x;
x=y;
if(b!=0)
y=t-a/b*y;
return d;
}
int main(){
u64 c,e,n,i,j,m,t,n0,T,ans,l;
__int64 T0,x,y,d;
limit = 1;
limit = limit<<63;
while(scanf("%I64u%I64u%I64u",&c,&e,&n)!=EOF){
m=0;n0=n;
while(n%2==0){n/=2;factor[m++]=2;}
while(n>1){
if(Miller_Rabin(n))break;
t = prime(n);
factor[m++] = t;
if(t!=0)
n/=t;
}
if(n>1)factor[m++]=n;
//for(l=0;lT0=(__int64)(factor[0]-1)*(factor[1]-1);
T=(u64)T0;
Euclid(e,T,x,y);
d=x;
//printf("%I64d\n",d);
//while(d

<0)d+=T0;
d=(d%T0+T0)%T0;
//d=(__int64)d;
// printf("%I64d %I64d\n",d,T0);
ans=square_multiply(c,(u64)d,n0);
printf("%I64u\n",ans);

}
}