本次学习的主要内容:
- 3.1素数运算的实验范例
- 3.2求解不定方程和同余方程的实验范例
- 3.3 特殊的同余式
- 3.4 积性函数的实验范例
- 3.5 高斯素数的实验范例
素数运算部分
主要是一些使用素数筛法和蛮力方法的一些题目
E - How many prime numbers
题意:给n个数,判断里面素数的个数。
直接使用试除法判断是否是素数,然后直接对素数个数计数
代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
ll exgcd(ll a,ll b,ll& x,ll& y) {
ll d=a;
if(b) {
d=exgcd(b,a%b,y,x);
y-=a/b*x;
} else x=1,y=0;
return d;
}
ll cal(ll a,ll b,ll c) {
ll x,y;
ll gcd=exgcd(a,b,x,y);
if(c%gcd!=0)return -1;
x=x*(c/gcd);
y=y*(c/gcd);
ll ans=1e15+7,t;
t=(y-x)*gcd/(a+b);
t=t-3;
for(int i=0; i<7; i++,t++) {
ll x1=x+b*t/gcd;
ll y1=y-a*t/gcd;
if(x1*y1>0)
ans=min(ans,max(abs(x1),abs(y1)));
else ans=min(ans,abs(x1-y1));
}
return ans;
}
int main() {
int cas;
scanf("%d",&cas);
while(cas--) {
ll A,B,a,b;
scanf("%lld%lld%lld%lld",&A,&B,&a,&b);
ll s=cal(a,b,abs(A-B));
printf("%lld\n",s);
}
return 0;
}
L - Goldbach’s Conjecture
题意:哥德巴赫猜想 任意大于 2 的偶数可以被分成两个素数,本题是大于4.
使用筛法筛出所有素数后,直接枚举每个素数 i 看 i - n是否是素数。如果是,直接可以输出。
代码:
#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <vector>
#define ll long long
using namespace std;
int n,m;
int main() {
while(scanf("%d",&n) != EOF) {
int cnt = 0;
for(int i = 0 ; i < n; i++){
int flag = 0;
scanf("%d",&m);
for(int i = 2; i <= sqrt(m); i++){
if( m % i == 0){
flag = 1;
break;
}
}
if(flag == 0)
cnt++;
}
printf("%d\n",cnt);
}
return 0;
}
M - Summation of Four Primes
题意:给一个数 n,将它分解成四个素数的和。
根据L题提到的哥德巴赫猜想,大于2的偶数能够分成两个素数的和 。根据输入的奇偶性,可以直接确定定死前两个素数 若输入是奇数,则定为2 3 ? ? 将原数减去5 就是偶数, 若是偶数,则定为2 2 ? ? ,将原数减去4还是偶数。
代码:
#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <vector>
#define ll long long
using namespace std;
const int Max = 1e7;
int n,m,f[Max];
vector<int> prime;
void init() {
for(int i = 2; i < Max; i++) f[i] = 1;
for(int i = 2; i < Max; i++) {
for(ll j = (ll)i * i; j < Max; j += i) {
f[j] = 0;
}
}
for(int i = 2; i < Max; i++)
if(f[i]) prime.push_back(i);
}
int main() {
init();
int size = prime.size();
while(scanf("%d",&n) != EOF) {
if(n < 8) {
printf("Impossible.\n");
} else if( n == 8) {
printf("2 2 2 2\n");
} else if(n == 9) {
printf("2 2 2 3\n");
} else if(n == 10) {
printf("2 2 3 3\n");
} else if(n == 11) {
printf("2 3 3 3");
} else if(n == 12) {
printf("3 3 3 3");
} else {
if(n % 2 == 0) {
printf("2 2");
n -= 4;
for(int i = 0 ; i < size; i++) {
if(f[n - prime[i]]) {
printf(" %d %d\n",prime[i],n-prime[i]);
break;
}
}
} else {
printf("2 3");
n -= 5;
for(int i = 0 ; i < size; i++) {
if(f[n - prime[i]]) {
printf(" %d %d\n",prime[i],n-prime[i]);
break;
}
}
}
}
}
return 0;
}
N - Digit Primes
题目:给一个区间[L,R],找出区间中满足 n 是素数,同时 n 的每一个数位的和也是素数的条件的数n。
此题要注意直接暴力求区间个数一定会超时,采用类似前缀和的思想。从1开始记录到当前位有多少满足条件的数的个数dp[i], dp[r] - dp[l-1] 就是答案。
代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
const int Max = 1e6 + 10;
int n,t1,t2;
int f[Max],dp[Max]={0};
void init() {
for(int i = 2; i < Max; i++) f[i] = 1;
for(int i = 2; i < Max; i++) {
for(ll j = (ll)i * i; j < Max; j += i) {
f[j] = 0;
}
}
dp[1] = 0;
for(int i = 2; i < Max; i++){
if(f[i]){
int sum = 0,tmp = i;
while(tmp){
int r = tmp % 10;
sum += r;
tmp /= 10;
}
if(f[sum]){
dp[i] = dp[i-1] + 1;
}else
dp[i] = dp[i-1];
}else
dp[i] = dp[i-1];
}
}
int main(){
init();
// for(int i = 2; i < 100; i++)
// printf(" %d: %d ",i,dp[i]);
scanf("%d",&n);
for(int i = 0 ; i < n; i++){
scanf("%d %d",&t1,&t2);
printf("%d\n",dp[t2]-dp[t1-1]);
}
return 0;
}
O - Prime Gap
题意:给一个数n,如果n是素数直接输出0,不然输出n左右各一个离得最近的素数的差。
直接依据题意进行暴力就过了,枚举素数。
代码:
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#define ll long long
using namespace std;
const int MAX=1299999;
int a[MAX],k;
bool vis[MAX];
void solve() {
memset(vis,1,sizeof(vis));
for(int i=2; i<=1299709; i++) {
if(vis[i]==1) {
for(int j=2; j*i<=1299709; j++) {
vis[i*j]=0;
}
}
}
for(int i=1; i<=1299709; i++) {
if(vis[i]==1) {
a[k]=i;
k++;
}
}
}
int main() {
int n;
solve();
while(scanf("%d",&n) & n) {
if(vis[n]==1) {
printf("0\n");
} else {
for(int i=1; i<k; i++) {
if(a[i]>n&&a[i-1]<n) {
printf("%d\n",a[i]-a[i-1]);
break;
}
}
}
}
return 0;
}
P - Happy 2006
题意:给两个数 k 和 m,找出第k大的和m互素的数。
因为k给的范围很到,如果直接找的话,根本无法打表。所以这里使用了一些小技巧。
用到了一个推论:
欧几里德算法GCD(a,b) = GCD(b,a mod b)
可以推出GCD(a * t + b,a ) = GCD(a, (a * t + b) mod a) = GDC(a,b)。
同时如果a与b互素,则b * t+a与b也一定互素;如果a与b不互素,则b * t+a与b也一定不互素。
与m互素的数 对m取模具有周期性。如下图
代码:
#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <vector>
#define ll long long
using namespace std;
const int Max = 1e6 + 10;
int k,m;
ll a[Max],cnt = 0;
int gcd(int a,int b) {
return b == 0 ? a : gcd(b,a % b);
}
int main() {
while(scanf("%d %d",&m,&k) != EOF) {
cnt = 0;
memset(a,0,sizeof(a));
for(int i = 1; i <= m; i++) {
if(gcd(i,m) == 1)
a[cnt++] = i;
}
printf("%lld\n",a[(k - 1) % cnt] + (k-1) / cnt * m);
}
return 0;
}
求解不定方程部分
本部分主要是使用扩展欧几里得算法来进行不定方程的求解
我认为这部分最重要的就是要搞清楚扩展欧几里得到底求出来的是什么 和 扩展欧几里得原理。
1.扩展欧几里得可以 求出不定方程ax + by = c的一组特解(x0,y0)。
2.原理
Ⅰ.吴永辉老师PPT
设d= GCD(a, b),a’=a DIV d,b’=b DIV d,并且c’=c DIVd。则不定方程ax+by=c可以被等价地写为a’x+b’y=c’,GCD(a’,b’)==1。如果c不是d的倍数,则不存在(x,y)满足不定方程。
采用扩展的欧几里得算法求解a’x+b’y=1,(x’, y’)是整数根。设x0=x’ * c’,y0=y’ * c’,则(x0,y0)是ax+by=c的一个特解,也就是说,a * x0 + b * y0=c。
所以,a(x0+b)+b(y0-a)=c,a(x0+2b)+b(y0-2a)=c,…,a(x0+kb)+b(y0-ka) = c,k是整数。所以,不定方程ax+by=c的通解是x= x0+kb,y= y0-ka,k是整数。
Ⅱ.董永建等《信息学奥赛一本通·提高篇》
当gcd(a,b)整除c时,设g = gcd(a,b),a’ = a/g ,b’ = b/g,c’ = c/g。
我们可以用exgcd(a’,b’)求出不定方程a’x’ + b’y’ = 1的整数解。 那么a’c’x’ + b’c’y’ = c’
a’gc’x’ + b’gc’y’ = c’g 即ac’x’ + bc’y’ = c,所以,x0 = c’x’,y0 =c’y’是方程的一组解。
因为不定方程a’x + b’y = c’ 等价于 同余方程 a‘x ≡ c’(mod b’),所以 x
是模b’同余的一个剩余类, 所以该不定方程的通解为 x = x0 + b’ * t ,y = y0 - a’ * t(t∈ Z)
c不是gcd(a,b)的倍数时没有上述步骤,方程无解。
A - The Balance
题意:有一个天平和多个质量为 a 和 b的砝码,求称出某一固定重量药物最小要多少个a砝码和b砝码。
思路很简单,设需要x个a砝码,y个b砝码,ax + by = 药物,求出 |x| + |y|的最小值。
最小值是个困难
用扩展欧几里德去求根,求出x,y然后分成两组,求出x可以根据方程推出y,求出y可以根据方程推出x.然后取两组相加的最小值即可。
其意义解释如下:
假设将物品放在天平的右边,对两种情况求解:
1)求x是作为解的最小正整数,即a毫克的砝码放在天平左边的最优解(如果y<0,则b毫克的砝码放天平的右边);
2)求y是作为解的最小正整数,即b毫克的砝码放在天平左边的最优解(如果x<0,则a毫克的砝码放天平的右边);
最后,两者中|x|+|y|小的就是结果。
代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
ll exgcd(ll a,ll b,ll &x,ll &y){
if(b == 0){
x = 1;
y = 0;
return a;
}
ll t = exgcd(b,a%b,x,y);
ll x0 = x,y0 = y;//x0 和 y0 是下一次 递归中 x 和 y 的解
x = y0;
y = x0 - a / b * y0;
return t;
}
int main(){
ll a,b,d,x,y;
while(scanf("%lld %lld %lld",&a,&b,&d) && (a + b + d != 0)){
ll gcd = exgcd(a,b,x,y);
a /= gcd;
b /= gcd;
d /= gcd;
exgcd(a,b,x,y);
ll tx = x * d;//求出x的解,然后根据方程推y
tx = (tx % b + b) % b;
ll ty = (d - a*tx ) / b;
if(ty < 0) ty = -ty;
ll yt = y * d;//求出y的解,然后根据方程推x
yt = (yt % a + a) % a;
ll xt = (d - b*yt ) / a;
if(xt < 0) xt = -xt;
if(xt + yt < tx + ty){//比较
tx = xt;
ty = yt;
}
printf("%lld %lld\n",tx,ty);
}
return 0;
}
B - One Person Game
题意:一个人从数轴上A点出发要到B点,他每次可以选择向左右任意走a步,b步,c步 共六种动作,且a+b = c 问A到B最少要多少步。
首先想到搜索,由于范围给的是2的-31次方到2的31次方,放弃。
设ax+by=|B-A|。如果|B-A|不是GCD(a, b)的倍数,则ax+by=|B-A|无解,否则用扩展的欧几里得算法求解。设 (x0, y0)是ax+by=B-A的一个解,则不定方程ax+by=c的通解是x=x0+kb,y=y0-ka,k是整数。
因为c可以表示为a+b,如果x==y,则同向行走a或b的步数相同,合并为c,x就是步数;
如果x不等于y,且xy>0,则同向行走a或b,步数是max(x, y),其中min(x, y)步行走c;
否则,如果x不等于y,且xy<0,则逆向行走a或b,步数是|x|+|y|。
最后,求最少步数。
同时也可以通过画一个图来理解这个题:
(来自litble博客)
通过扩欧可以得到一个特解x0和y0,然后通解可以表示为x=x0+at,y=y0−bt
可以在坐标系上画出两条直线l1,l2(注意上式中的x不同于坐标系的x轴,x轴准确的说应该是t轴!),并且斜率是一正一负(a,b同正,题目已说)
所以呢,取一个t,如果x和y异号,答案就是你画一条平行于y轴的直线,它和l1,l2的交点间的距离。如果同号就是在这两条直线上横坐标为t的离x轴远的那个点离x轴的距离,那么意会可得,l1,l2交点处取到最优解,由于交点处的t不一定是整数,所以我们要做一些处理,考虑交点附近的整点。
代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
ll exgcd(ll a,ll b,ll& x,ll& y) {
ll d=a;
if(b) {
d=exgcd(b,a%b,y,x);
y-=a/b*x;
} else x=1,y=0;
return d;
}
ll solve(ll a,ll b,ll c) {
ll x,y;
ll gcd=exgcd(a,b,x,y);
if(c%gcd!=0)return -1;
x=x*(c/gcd);
y=y*(c/gcd);
ll ans=1e15+7,t;
t=(y-x)*gcd/(a+b);
t=t-3;
for(int i=0; i<7; i++,t++) {
ll x1=x+b*t/gcd;
ll y1=y-a*t/gcd;
if(x1*y1>0)
ans=min(ans,max(abs(x1),abs(y1)));
else ans=min(ans,abs(x1-y1));
}
return ans;
}
int main() {
int cas;
scanf("%d",&cas);
while(cas--) {
ll A,B,a,b;
scanf("%lld%lld%lld%lld",&A,&B,&a,&b);
ll s=solve(a,b,abs(A-B));
printf("%lld\n",s);
}
return 0;
}
求解同余方程
C - Raising Modulo Numbers
题意:
使用同余理论的运算规则求解:
- (a + b) % p = (a % p + b % p) % p (1)
- (a ^ b) % p = ((a % p)^b) % p (4)
代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <string>
#define ll long long
using namespace std;
const int Max = 45010;
int cas,m,h,a[Max],b[Max];
int qpow(int a,int b,int mod){
int ans = 1;
while(b){
if(b & 1){
ans = ans * a;
ans = ans % mod;
}
a *= a;
a %= mod;
b >>= 1;
}
return ans;
}
int main(){
scanf("%d",&cas);
while(cas--){
int res = 0;
scanf("%d %d",&m,&h);
for(int i = 1; i <= h; i++){
scanf("%d %d",a+i,b+i);
res += (qpow(a[i]%m,b[i],m)%m);
}
res %= m;
printf("%d\n",res);
}
return 0;
}