第16章 数论问题 if(!(n&1))return 0; u64 k = 0, i, j, m, a; m = n - 1;
while(m%2==0)m=(m>>1),k++; for(i=0;i
a = square_multiply(random()%(n-1)+1, m, n);//平方乘 if(a==1)continue; for(j=0;j
if(j
return 1; }
//Pollard p,只找出一个因子。 u64 gcd(u64 a,u64 b){ if(b==0) return a; else return gcd(b,a%b); }
//用公式f(x) = x^2 + 1检验碰撞。 u64 f(u64 x, u64 n ){
return (by(x,x,n)+1)%n; }
//分解不到，return 0 u64 Pollard(u64 n){
if(n<=2)return 0;
if(!(n&1))return 2; //必不可少 u64 i, p, x,xx;
for(i=1;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; }
243·
·
ACM程序设计培训教程
u64 prime(u64 a){
if(Miller_Rabin(a))return 0;
u64 t = Pollard(a), p,d;d=a/t;if(d
int main(){
u64 l, a, t; limit = 1;
limit = limit<<63; //动态化分段使用 scanf(\
while(t--){scanf(\
if(Miller_Rabin(a))printf(\ else printf(\ }
return 0; }
16.4 其他—Pell方程
16.4.1 〖案例6〗Smith问题
计算机将人们从繁重的计算中间，在计算机发明之前，人们必须靠手算，而且容易错。但现在，借助计算机，大量的计算也是瞬间而就。
但今天Smith教授遇到一个问题，他需要知道下列形式方程的整数解：
x2?Ny2?1
其中N是一个正整数。但是对于很多N，教授的破计算机似乎找不到解。 “这就奇怪了”Smith教授说“一定是我的程序哪里出错了”。你，教授的助手，需要解决这个问题。
有多个测试序列。每个测试序列一行一个整数N (1 < =N < =10^8)。
对于每个测试序列，输出方程的在10^1000内的最小整数解X和Y，如果没有这样的解，输出一行\。
52
244· ·
第16章 数论问题 256 991
649 90
No solution!
379516400906811930638014896080 12055735790331359447442538767
形如x2?Ny2?1的方程是所谓的Pell方程的特例，而起最小整数解可以用连分数法得到。也就是将N化成如下的形式：
N?a0?a1?x?a0?y1a1?1??11?2?111?16，而
11 ??an1如果有ak?2a0,且
1ak?1，则（x,y）就是最小整数解。
14?3?例如
15?3?411?112?1，所以（15，4）就是方程
x2?14y2?1的最小整数解。
当然N是完全平方数时方程没有整数解，因为整数平方差不可能为1。 为了求得N的近似值，使用以下的递推方程：
a0??N?p0?a0p1?a1p0?1pn?anpn?1?pn?2 q0?1q1?a1qn?anqn?1?qn?2这里的?x?表示下取整。计算另外两个数：
·245·
ACM程序设计培训教程
P0?0P1?a0Pn?an?1Qn?1?Pn?2Q0?1n?Pn2Qn?Qn?1?a?Pn?an??0?Qn??可以证明数列满足下列两个重要性质：
pnqn?1?pn?1qn???1?22pn?Nqn???1?n?1n?1Qn?1
所以方程的最小正整数解为
??p,q??x,y???rr??p2r?1,q2r?1?r奇数r偶数
#include
#define Base 1000000000 #define Cap 200
typedef __int64 hugeint;
struct bignum {
int Len; int Data[Cap]; bignum() : Len(0) {}
bignum(const bignum& V) : Len(V.Len){ memcpy(Data, V.Data, Len * sizeof *Data);}
bignum(int V) : Len(0) { for (; V > 0; V /= Base) Data[Len++] = V % Base; } bignum& operator=(const bignum& V) { Len = V.Len;
memcpy(Data, V.Data, Len * sizeof *Data); return *this; } int& operator[](int Index) { return Data[Index]; } int operator[](int Index) const { return Data[Index]; } };
246· ·