第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· ·