1、算法艺术与信息学竞赛标准课件,数论(二): 经典问题和算法刘汝佳,问题1: 大整数取余,给一个n位的大整数P和正整数m求P mod m的值输入int n, int mint digit：digit0为P第首位，digitn-1末位输出int d：P mod m的值,分析,公式一: (a*b) mod m = (a mod m) * (b mod m) mod m公式二: (a+b) mod m = (a mod m) + (b mod m) mod m注意: 在两个公式中, 都需要在最后对m取模,分析,由公式, 可以由前i-1位的余数di-1推出第i位的余数di = (di-1*10+digi
2、ti) mod m时间复杂度: O(n)每个di都只使用一次，空间复杂度O(1),d = 0;for(i = 0; i 1由于每个数只被用一次, 空间是O(1)的问题: dk-1*dk-1可能会overflow! (n=105时已经会) 解决方法: 使用更大的整数范围,分析,下面的数据类型LL取决于编译器. gcc使用long long, 而VC+使用_int64对于其他操作符, 把*d%p替换掉即可,LL ans = 1, d = a % p;do if(n ,问题3: 求1n的所有素数,求1n的所有素数输入int n输出bool isPrime：数i(1=i=n)是否为素数int prim
3、e：第i(1=i=(n+1)个素数int primeCount：素数总数,分析,假设要求1100的素数2是素数, 删除2*2, 2*3, 2*4, , 2*50第一个没被删除的是3, 删除3*3, 3*4, 3*5,3*33第一个没被删除的是5, 删除5*5, 5*6, 5*20得到素数p时, 需要删除p*p, p*(p+1), p*n/p, 运算量为n/p-p, 其中p不超过n1/2(想一想, 为什么),Eratosthenes的筛子,primeCount = 0;for(i = 2; i n; i+) isPrimei = true;for(i = 2; i n; i+) if(isPri
4、mei) primesprimeCount+ = i;if(i = n/i) for(j = i*i; j n. 更好的办法是递推求i2, 利用(i+1)2-i2=2i+1, 另设一个变量k=i2100, 1000, 10000109内的素数个数分别为：25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534用筛法一般最多用来计算107内的素数,优化,枚举过程也可以优化一下优化一：除了2以外偶数都不是素数，所以每次i可以增加2优化二：除了2、3以外，素数p除以6的余数只能是1或5，所以可以修改为每次交替增加2, 4时间优化并不明显, 但空间分别
5、缩小为原来的1/2和1/3,分析,时间复杂度显然为筛的时间复杂度+O(n)筛的过程不超过n/2+n/3+n/5+由公式得: 筛的过程是nlnlnn的, 总O(nloglogn)其中B1为Mertens常数0.2614972128,问题4: 素数判定,给定正整数p判断p是否为素数输入int p输出bool isPrime,算法一,朴素的枚举法, 枚举到n1/2为止任务: 统计1n的素数个数n=105时瞬间, n=106时需要数秒,for(i = 2; i * i = p; i+) if(p % i = 0) return false;return true;,优化,可以利用前面所说的模6优化,
6、速度为原来的若干倍. 优化后106次一秒之内出解, 107需要约10秒,if(p=2|p=3) return true;if(p%2=0|p%3=0) return false;for(i=5,k=4; i*i=p; i+=(k=6-k) if(p%i=0) return false;return true;,优化,也可以改成只用素数试除, 速度变快但速度并不是很明显(n=107时约两倍). 其中primes初始化为2和3, 随着循环测试的进行不断更新,for(i=0; primesi*primesi= p; i+) if(p % primesi = 0) return false;retur
7、n true;,Miller-Rabin测试,对于奇数n, 记n=2r*s+1, 其中s为奇数随机选a(1=a=n-1), n通过测试的条件是as1(mod n), 或者存在0=jb时a和b均为偶数, gcd(a,b)=2*gcd(a/2,b/2)a为偶数, b为奇数, gcd(a,b)=gcd(a/2,b)如果a和b均为奇数, gcd(a,b)=gcd(a-b,b)不需要除法, 只需减法和右移, 适合大整数,Euclid算法,测试: 求1=a, b=n的所有gcd之和n=1000时为4449880n=5000时为135637352需要的次数满足Lam:即steps= 1; c = 1; el
8、se c=0; if(!(b,问题6: 线性不定方程,给出整数a和b和c, 求出一组整数x, y, 满足ax + by = c输入int a, int b, int c输出int x, int y,分析,方程: ax+by=c设gcd(a,b)=d, 则如果c不是d的倍数, 一定无解, 因为方程左边为d的倍数否则可以令c=c*d, 当求出方程ax+by=d的任意解(x0, y0)后, (cx0, cy0)就是原方程的解, 因ax+by=acx0 + bcy0 =c(ax0+by0)=cd=c问题转化为求方程ax+by=d的一组解,分析,ax+by=d的一组解可以用以下算法求出:,gcd(int
9、 a, int b, int,证明,边界: b = 0时a*1+b*0 = a成立, 否则递归调用后, y和x满足by+a%b*x=d赋值后y = y-a/b*x, 计算ax+by即可注意a/b是整除运算, 它等于(a-a%b)/b, 故ax+by = ax+by- (a-a%b)/b*x)*b= ax+by- (a-a%b)*x = by+a%b*x = d,分析,如何通过这组解求出所有解?如果有两组解(x1, y1)和(x2, y2), 有ax1+by1=d, ax2+by2=d因此, a(x1-x2)=b(y2-y1), 消去公因子d得a(x1-x2)=b(y2-y1)因为a和b互素,
10、x1-x2必须是b的整数倍. 设x1-x2=kb, 则由等式得y2-y1=ka, 因此x1-x2必须取整数值，且可取任意整数值,问题7: 模线性方程,给出正整数a, b, n求解模线性方程axb(mod n)输入:Int a, int b, int n输出Int d: 解剩余系的个数int sol: 对0=id, xsoli(mod n)是解,分析,首先明确一点, 如果x是解(0=xn), 则对于任意整数k, x+kn也是解, 所以解应表示成一些剩余类xxiaxb(mod n)等价于存在整数y, 使得ax-ny=b这是一个线性同余方程, 首先判断d=(a,n)是不是b的约数, 如果是, 等价于
11、方程ax-ny=b, 相当于求解axb(mod n), (a,n)=1,分析,这个方程有两种解法直接解ax-ny=d, 再两边同时乘以b/d. 注意这只是一个解, 一共应该有d个, 间隔为n/d(因为对于模n/d来说解是唯一的)在模n的缩系中a是存在逆a-1的. 因此解为xa-1b(mod n), 同样扩展得d个解.在缩系中求逆也有两种方法求ax-ny=1的解, 和方法一等价利用欧拉定理a-1a(n)-1(mod n),分析,特别地, (a, n)=1时, 求a-1(mod n)程序如下:,int inverse(int a, int n) int d, x, y; gcd(a, n, d,
12、x, y); return x;,分析,方法一的程序如下,gcd(a, n, d, x, y);if(b%d) d = 0; / no answerelse sol0 = x * (b/d) % n; for(i = 1; i d; i+) soli = (soli-1 + n/d) % n;,分析,由刚才的分析, 方程axb(mod n)完全被转化为了xx0 (mod m0), 其中x0=sol0, m0=n/(a, n)下面考虑由形如xxi (mod mi)的方程所组成的方程组,问题8: 模线性方程组,有n个方程组成一个方程组xai (mod mi)其中所有模mi两两互素求任意一个解x输入
13、int n, int a, int m输出int x,分析,令M=m1m2mn中国剩余定理指出, 这样的解在0,M-1的范围内是唯一的, 设它为x, 则任意解都能写成x+kM的形式. 关键是要求出x记wi=M/mi，则(wi,mi)=1, 因此存在pi,qi使得wipi+miqi=1,分析,关键等式: wipi+miqi=1记ei=wipi，则j=i时ei1(mod mj) (等式两边取模mj)ji时ei0(mod mj) (ei中含有因子mj)因此e1a1+e2a2+enan就是一个解调整到范围0,M-1中程序实现: 一个一个方程考虑, 每次求出wi, pi, ei, 并累加eixi到结果中
14、,分析,注意调用gcd时不要把x覆盖掉了,M=1; for(i = 0; i n; i+) M*=mi;for(i = 0; i n; i+) w = M / mi; gcd(mi, w, dummy, dummy, y); x = (x + y*w*ai) % M;return (n+x%M)%M;,分析,顺便说一下, 对于任意的模m, 可以分解成素数幂乘积的形式, 则问题问题完全转化为本问题(需要先判断无解),问题9: 求欧拉函数,给整数n, 求(n)输入int n输出Int phi,分析,对于素数的幂pk, 和pk不互素的数一定有因子p, 即p, 2p, 3p, pk-1*p, 共pk-
15、1个. 因此(pk)=pk-pk-1=pk(1-1/p)因为是积性函数, 因此有公式需要得到素因数分解式. 注意每得到一个素因子p, 需要约去p的所有指数,分析,注意乘法总在除法之后, 因此不会overflow,phi = n;for(i=2,j=4;j1) phi = phi / n * (n-1);,问题10: 求1n的欧拉函数,给正整数n, 求1n每个数的欧拉函数值输入Int n输出Int phi,分析,n的每个素因子p对n的phi值都有(1-1/p)的贡献, 因此可以借助素数筛法实现, 但是必须从p开始而不是p*p开始, 因此更接近nlnlnn不需要额外的isPrime标志, 直接判断
16、目前的phi. 如果为0, 说明为素数先检查j的phi函数值. 如果为0, 说明从来没有被筛过, 赋初值j.筛的时候顺便把phij乘以(1-1/p), 注意先除再乘避免overflow,分析,可以顺便求出素数表, 只需要在if(!phii)时先进行primesprimeCount+ = i即可.,phi1 = 1;for(i = 2; i = n; i+) if(!phii) for(j=i; j=n; j+=i) if(!phij) phij = j; phij = phij / i * (i-1); ,问题11: 离散对数,给a, b, n, 求方程axb(mod n)的一个解(n为素数)
17、输入Int a, int b, int n输出Int x,分析,离散对数是一个困难的问题这里只考虑n为素数的情况. 由欧拉定理, 只需要检查x=0,1,2,n-1的情况, 时间复杂度为O(n)由于n为素数, 只要a不为0, 一定存在逆a-1首先检查前m项, 即a0, a1, a2, am-1模n的值是否为b, 并把ai mod n保存在ei里求出am的逆a-m,分析,下面检查am, am+1, a2m-1如果有解, 相当于存在i使得ei*amb(mod n)两边左乘v = a-m, 得到eiv*b(mod n)只需检查是否存在ei等于给定值事先给e排序, 每次可二分查找每处理一轮, 查找内容需
18、要乘以v,分析,需要O(mlogm)的时间计算e数组并排序, O(n/m*logm)的时间进行剩下的工作, 取m=n1/2, 则时间复杂度为O(n1/2logn)可以用sort+二分查找, 也可以用map下面的程序使用了一个map x, 其中xj代表满足ei=j的最小i(可能有多个i, 如a=1时),m = (int)ceil(sqrt(n);v = inverse(pow_mod(a, m, n), n);e = 1; x1 = 0;for(i = 1; i m; i+) e = e * (LL)a % n; if(!x.count(e) xe = i; for(i = 0; i m; i+) if(x.count(b) return i*m+xb; b=b * (LL)v % n; ,问题12: 模平方根,给整数a, n求满足x2a(mod n)的所有根输入Int a, int n输出Int d: 解的个数Int x: 所有解,分析,类似于实数方程, 可能有三种情况: 无解, 两等根和两异根下面介绍Shank算法,
