数论(二)

Miller-Rabin测试的代码
• 以下代码适合于大于7, 小于109的数(如果要 扩展到2.5*1013, 则需要单独判断 3215031751, 并避免LL乘法溢出)
if(!miller_rabin(n, 2)) return false; if(!miller_rabin(n, 3)) return false; if(!miller_rabin(n, 5)) return false; if(!miller_rabin(n, 7)) return false; return true;
问题5: 最大公约数
• 给a, b, 求a和b的最大公约数gcd(a, b) • 输入
– int a – int b
• 输出
– int d
分析
• 方法一: 先分解素因数, 然后求最大公约数 • 方法二: (Euclid算法)利用公式gcd(a, b)=gcd(b, a mod b), 时间复杂度为O(logb) • 方法三: (二进制算法) gcd(a,a)=a, a>b时
分析
• ax+by=d的一组解可以用以下算法求出:
gcd(int a, int b, int& d, int& x, int& y){ if(!b){ d = a; x = 1; y = 0; } else{ gcd(b, a%b, d, y, x); y -= x*(a/b); }
证明
• • • • 边界: 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*y’ = d
for(i = 2; i * i <= p; i++) if(p % i == 0) return false; return true;
优化
• 可以利用前面所说的模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;
• 得到素数p时, 需要删除p*p, p*(p+1), … p*[n/p], 运算量为[n/p]-p, 其中p不超过 n1/2(想一想, 为什么)
Eratosthenes的筛子
primeCount = 0; for(i = 2; i < n; i++) isPrime[i] = true; for(i = 2; i < n; i++) if(isPrime[i]){ primes[primeCount++] = i; if(i <= n/i) for(j = i*i; j <= n; j += i) isPrime[j] = false; }
– 预先判断是否能被2, 3, 5, 7整除(而不是在每个 测试中做), n=107减小到14秒, n=108约2.5分钟
• 任务二: 测试109~109+106的数
– 答案: 共48155个素数. 比较试除法和MillerRabin测试的速度
用a测试n的代码
int r = 0, s = n - 1, j; if(!(n%a)) return false; while(!(s&1)){ s >>= 1; r++; } LL k = pow_mod(a, s, n); if(k == 1) return true; for(j = 0; j < r; j++, k = k * k % n) if(k == n - 1) return true; return false;
n1的情况时间复杂度由于n为素数只要a不为0一定存在逆a1值是否为b并把amodn保存在ei里如果有解相当于存在i使得事先给e排序每次可二分查找每处理一轮查找内容需要乘以v分析需要om的时间计算e数组onmlogm的时间进行剩下的工作则时间复杂度为on12logn可以用sort二分查找也可以用map其中xj代表满足eifxcountbreturn问题12
Miller-Rabin测试
• 对于奇数n, 记n=2r*s+1, 其中s为奇数 • 随机选a(1<=a<=n-1), n通过测试的条件是
– as≡1(mod n), 或者 – 存在0<=j<=r-1使得a(2^j)*s≡-1(mod n)
• 素数对于所有a通过测试, 合数通过测试的 概率不超过1/4 • 注意: 先要判断此数本身是否为a的倍数
问题3: 求1~n的所有素数
• 求1~n的所有素数 • 输入
– int n
• Байду номын сангаас出
– bool isPrime[]：数i(1<=i<=n)是否为素数 – int prime[]：第i(1<=i<=π(n+1))个素数 – int primeCount：素数总数
分析
• 假设要求1~100的素数
– 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
优化
• 枚举过程也可以优化一下
– 优化一：除了2以外偶数都不是素数，所以每 次i可以增加2 – 优化二：除了2、3以外，素数p除以6的余数只 能是1或5，所以可以修改为每次交替增加2, 4
• 时间优化并不明显, 但空间分别缩小为原来 的1/2和1/3
分析
• 时间复杂度显然为筛的时间复杂度+O(n) • 筛的过程不超过n/2+n/3+n/5+… • 由公式
《算法艺术与信息学竞赛》 标准课件
数论(二): 经典问题和算法 刘汝佳
问题1: 大整数取余
• 给一个n位的大整数P和正整数m • 求P mod m的值 • 输入
– int n – int digit[]：digit[0]为P第首位，digit[n-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位的余数d[i-1]推出第i 位的余数d[i] = (d[i-1]*10+digit[i]) mod m • 时间复杂度: O(n) • 每个d[i]都只使用一次，空间复杂度O(1)
• 得: 筛的过程是nlnlnn的, 总O(nloglogn) • 其中B1为Mertens常数0.2614972128…
问题4: 素数判定
• 给定正整数p • 判断p是否为素数 • 输入
– int p
• 输出
– bool isPrime
算法一
• 朴素的枚举法, 枚举到n1/2为止 • 任务: 统计1~n的素数个数 • n=105时瞬间, n=106时需要数秒
– 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之和
二进制算法
int t = 1, c, d; while(a != b){ if(a < b) swap(a, b); if(!(a&1)){ a >>= 1; c = 1; }else c=0; if(!(b&1)){ b >>= 1; d = 1; }else d=0; if(c && d) t <<= 1; else if(!c && !d) a -= b; } return t * a;
分析
• 设si为数(nknk-1…ni), 则nk = sk & 1 • 从右往左递推
– dk = (dk-1 * dk-1) mod p – sk = sk-1 >> 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 & 1) ans = ans * d % p; d = (d * d) % p; }while(n >>= 1);
– n=1000时为4449880 – n=5000时为135637352
• 需要的次数满足[Lamé]:
• 即steps<=4.785lgN + 1.6723 • Lamé定理: 算法的最坏情况为计算gcd(Fn, Fn-1) • Heilbronn定理: 平均次数12ln2/π2lgn=0.843lgn
优化
• 也可以改成只用素数试除, 速度变快但速度 并不是很明显(n=107时约两倍). • 其中primes初始化为2和3, 随着循环测试的 进行不断更新
for(i=0; primes[i]*primes[i]<= p; i++) if(p % primes[i] == 0) return false; return true;
Miller-Rabin测试
• • • • 能通过a=2的最小合数是2047=23*89 能通过a=2, 3的最小合数是1373653 能通过a=2, 3, 5的最小合数是 能通过a=2, 3, 5, 7的, 2.5*1013以内唯一的 合数是3215031751
实现与实验
• 计算出k = as mod n是首先判断k是否为1或 者-1, 然后j每次加1, 相当于每次t平方, 而不 需要每次都调用幂取模的函数重算一次 • 任务一: 统计1~n的素数个数
Euclid算法
• 递归形式和迭代形式(效率基本相当) int gcd(int a, int b){ return (b? gcd(b, a%b) : a); } int gcd(int a, int b){ int t; while(b){ t = a % b; a = b; b = t; } return a; }
d = 0; for(i = 0; i < n; ++i) d = (d * 10 + digit[i]) % m
问题2: 模取幂
• 给出a, n, p, 求an mod p的值 • 输入
– int a – int n – int p
• 输出
– int d
分析
• 算法一: 利用公式 an mod p = (an-1 mod p) * a mod p • 时间复杂度: O(n) • 算法二: 把n写成二进制n = (nknk-1…n0)2 • 设di为a的2i次方模p的结果, 则 an mod p = (dk^nk) * (dk-1^nk-1) * … • 即把nk为1的那些dk乘起来 • 时间复杂度: O(logn)
问题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)后, (c’x0, c’y0)就是原方程的解, 因ax+by=ac’x0 + bc’y0 =c’(ax0+by0)=c’d=c • 问题转化为求方程ax+by=d的一组解
分析
• 问题：i*j可能超界！ • 解决办法： 改用除法判断，并预先判断是 否有i2>n. 更好的办法是递推求i2, 利用 (i+1)2-i2=2i+1, 另设一个变量k=i2 • 100, 1000, 10000…109内的素数个数分别 为：25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534 • 用筛法一般最多用来计算107内的素数