Acdreamer博客数论学习

Day One

等比数列求和．

求解 Sn=(a+a2+...+an)modM

采用二分．n%2==0时 Sn=(1+a(n2))Sn2 ,n%2==1时， Sn=(1+an+12)Sn−12+an+12

#include<cstdio>
#include<iostream>
#define ll long long
using namespace std;
const int mod = 7;
ll pow(int a, int k){
ll ans = 1;
while(k){
if(k&1)ans=ans*a%mod;
a=a*a%mod;
k>>=1;
}
return ans;
}
ll S(int a, int n){
if(n==1)return a%mod;
if(n&1){
return ((1+pow(a,(n+1)/2))%mod*S(a,(n-1)/2)%mod + pow(a,(n+1)/2)%mod)%mod;
}else{
return (1+pow(a,n/2))%mod*S(a,n/2)%mod;
}
}
int main(){
int a,n;
while(~scanf("%d%d",&a,&n)){
printf("%lldn",S(a,n));
}
return 0;
}

题目：https://vjudge.net/problem/POJ-3233．

计算矩阵１－n次方和

int n,k,m;
int mod = 1e9+7;
#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
struct Matrix{
int mp[35][35];
Matrix(){
memset(mp,0,sizeof mp);
}
Matrix(int v){
memset(mp,0,sizeof mp);
for(int i = 1; i <= n; i++){
mp[i][i] = v;
}
}
};
Matrix multi(Matrix a, Matrix b){
Matrix c;
for(int i = 1; i <= n; i++){
for(int k = 1; k <= n; k++){if(a.mp[i][k])
for(int j = 1; j <= n; j++){if(b.mp[k][j])
c.mp[i][j] = c.mp[i][j] + a.mp[i][k]*b.mp[k][j];
c.mp[i][j] %= mod;
}
}
}
return c;
}
Matrix pls(Matrix a, Matrix b){
Matrix c;
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n; j++){
c.mp[i][j] = a.mp[i][j]+b.mp[i][j];
c.mp[i][j] %= mod;
}
}
return c;
}
Matrix powmul(Matrix a, int k){
Matrix c = Matrix(1);
for(int i = 1; i <= n; i++)
c.mp[i][i] = 1;
while(k){
if(k&1)c = multi(c,a);
a = multi(a,a);
k>>=1;
}
return c;
}
Matrix S(Matrix a, int k){
Matrix temp = Matrix(1);
if(k==1)return a;
if(k&1){
Matrix temp2 = powmul(a,(k+1)/2);
return pls(multi(pls(temp, temp2), S(a,(k-1)/2)), temp2);
}else{
Matrix temp2 = powmul(a,k/2);
return multi(pls(temp, temp2), S(a,k/2));
}
}
int main(){
while(~scanf("%d%d%d",&n,&k,&m)){
Matrix a;
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n; j++){
scanf("%d",&a.mp[i][j]);
}
}
mod = m;
a = S(a,k);
for(int i = 1 ; i <= n; i++){
for(int j = 1; j <= n; j++){
printf("%d",a.mp[i][j]);
if(j!=n)printf(" ");
}
puts("");
}
}
return 0;
}

素数定理

定理 gcd(an−1,am−1)=agcd(n,m)−1

在a>1,m,n>0时成立．

#include<cstdio>
#include<iostream>
#define ll long long
using namespace std;
int mod = 1e9+7;
ll gcd(ll a, ll b){
return b?gcd(b,a%b):a;
}
ll pow(ll a, ll k){
ll ans = 1;
while(k){
if(k&1)ans=ans*a%mod;
a = a*a%mod;
k>>=1;
}
return ans;
}
int main(){
int a,m,n,k1,t;scanf("%d",&t);
while(t--){
scanf("%d%d%d%d",&a,&m,&n,&k1);
mod = k1;
ll ans = pow(a,gcd(n,m))-1+mod;
printf("%lldn",ans%mod);
}
return 0;
}

定理扩展 a>b,gcd(a,b)=1,gcd(am−bm,an−bn)=agcd(m,n)−bgcd(m,n)

定理 G=gcd(C1n,C2n,C3n,...,Cnn−1)

1. n为素数，答案为n；
2. n有多个素因子，答案为１；
3. n只有一个素因子,答案就是该素因子(其实包括了第一种情况

题目：https://vjudge.net/problem/HDU-2582

#include<cstdio>
#include<iostream>
#include<cstring>
#define ll long long
#define maxn 1000100
using namespace std;
bool judge[maxn];
int prime[maxn/10];
ll biao[maxn];
int tot = 0;
void init(){
memset(prime,0,sizeof prime);
memset(judge,0,sizeof judge);
for(int i = 2; i < maxn; i++){
if(!judge[i]){
for(int j = i+i; j < maxn; j+=i){
judge[j]=1;
}
prime[tot++]=i;
}
}
fill(biao,biao+maxn,1);
for(int i = 0; i < tot; i++){
for(ll j = prime[i]; j < maxn; j*=prime[i]){
biao[j]=prime[i];
}
}
for(int i = 4; i < maxn; i++){
biao[i]+=biao[i-1];
}
}
int main(){
int n;
init();
while(~scanf("%d",&n)){
printf("%lldn",biao[n]);
}
return 0;
}

定理 Fn为斐波那契额数列，gcd(Fm,Fn)=Fgcd(m,n)

http://acm.nyist.net/JudgeOnline/problem.php?pid=468

此处需要学会一种新的判断素数的方式叫做Miller_Rabin素数测试．在DayTwo中有．

#include<stdio.h>
#include<algorithm>
#include<iostream>
using namespace std;
long long multi(long long a, long long b, long long mod){
long long temp = a, sum = 0;
while(b){
if(b & 1) sum = (sum + temp) % mod;
temp = (temp + temp) % mod;
b = b >> 1;
}
return sum;
}
long long Modular_exponent(long long a, long long x, long long mod){
long long t = a % mod, r = 1;
while(x){
if(x & 1) r = multi(r, t, mod);
t = multi(t, t, mod); x = x >> 1;
}
return r;
}
bool Miller_Rabin(long long n){
long long time = 20;
if(n < 2) return false;
if(n == 2) return true;
bool flag = false;
for(long long k = 0; k <= time; ++k){
flag = false;
long long d = n - 1, r = 0, t, a = rand() % (n - 2) + 2;
while((d & 1) == 0){
d = d>>1;
r++;
}
t = Modular_exponent(a, d, n);
if(t == 1 || t == n-1) {flag = true; continue;}
for(long long i = 1; i < r; i++){
t = multi(t, t, n);
if(t == 1) {flag = false; return flag;}
if(t == n-1) {flag = true; break;}
}
if(!flag) break;
}
return flag;
}
int main()
{
long long n;
while(scanf("%lld",&n)!=EOF)
puts((n == 4 || (n != 2 && Miller_Rabin(n))) ? "Yes" : "No");
return 0;
} 

定理：给定两个互素的正整数Ａ和Ｂ，那么他们最大不能组合的数为AB-A-B,不能组合的个数为num= (A−1)(B−1)2

ax+by=c的存在非负解的条件就是c>AB-A-B.

HDU - 1792

#include<cstdio>
#include<iostream>
using namespace std;
int main(){
int a,b;
while(~scanf("%d%d",&a,&b)){
printf("%d %dn",a*b-a-b,(a-1)*(b-1)/2);
}
return 0;
}

定理 ∑gcd(i,N)=∑dϕ(Nd)

数论专题里也有一题如是说．

https://vjudge.net/problem/POJ-2480

#include<cstdio>
#include<cstring>
#include<iostream>
#define ll long long
#define maxn 50000
using namespace std;
bool judge[maxn];
int prime[maxn];
int tot;
void init(){
tot = 0;
memset(judge,0,sizeof judge);
memset(prime, 0, sizeof prime);
for(int i = 2; i < maxn; i++){
if(!judge[i]){
for(int j = i+i;
j < maxn; j+=i){
judge[j] = 1;
}
prime[tot++]=i;
}
}
}
int main(){
init();
ll n;
while(~scanf("%lld",&n)){
ll ans = n;
for(int i = 0; i < tot && prime[i]*prime[i] <= n; i++){
if(n%prime[i]==0){
int cot = 0;
while(n%prime[i]==0){cot++;n/=prime[i];}
ans += ans * cot * (prime[i]-1)/prime[i];
}
}
if(n!=1)ans += ans*(n-1)/n;
printf("%lldn",ans);
}
}

剩下没有题目的定理

结论：

Day Two

卢卡斯－莱默检验法 检验梅森素数

#include<cstdio>
#include<iostream>
#define ll long long
using namespace std;
ll T,p;
ll multi(ll a, ll b, ll m){
ll ans = 0;
while(b){
if(b&1)ans = (ans+a)%m;
a = (a+a)%m;
b>>=1;
}
return ans;
}
int main(){
scanf("%lld",&T);
while(T--){
scanf("%lld",&p);
ll n = ((ll) 1<<p)-1;
ll r = 4;
if(p==2)puts("yes");
else{
while((--p)!=1)r=(multi(r,r,n)-2+n)%n;
if(r%n==0)puts("yes");
else puts("no");
}
}
}

Miller-Rabin素数测试

​ 费马小定理：对于素数p和任意整数a，有

　　伪素数：如果n是一个正整数，如果存在和n互素的正整数a满足 an−1≡1(modn) ，我们说n是基于a的伪素数。如果一个数是伪素数，那么它几乎肯定是素数。

​ 序列:
​ a^m%n; a^(2m)%n; a^(4m)%n; …… ;a^(m*2^q)%n把上述测试序列叫做Miller测试,关于Miller测试,有下面的定理:

定理:若n是素数,a是小于n的正整数,则n对以a为基的Miller测试,结果为真.Miller测试进行k次,将合数当成素数处理的错误概率最多不会超过4^(-k).

　　Miller-Rabin测试：不断选取不超过n-1的基b(s次)，计算是否每次都有

​ 二次探测定理 ：p为素数，０

#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<iostream>
#define ll long long
using namespace std;
const int Times = 10;
ll multi(ll a, ll b ,ll m){
ll ans = 0;
a%=m;
while(b){
if(b&1)ans = (ans+a)%m;
a = (a+a)%m;
b >>= 1;
}
return ans;
}
ll quick_mod(ll a, ll b, ll m){
ll ans = 1;
a %= m;
while(b){
if(b&1)ans = multi(ans, a, m);
a = multi(a,a,m);
b>>=1;
}
return ans;
}
bool Miller_Rabin(ll n){
if(n==2)return true;
if((n<2) || !(n&1))return false;
ll m = n-1;
int k = 0;
while((m&1)==0){
k++;
m>>=1;
}
for(int i = 0; i < Times; i++){
ll a = rand()%(n-1)+1;
ll x = quick_mod(a,m,n);
ll y = 0;
for(int j = 0; j < k; j++){
y = multi(x,x,n);
if(y==1 && x!=1 && x!= n-1)return false;
x = y;
}
if(y!=1)return false;
}
return true;
}
int main(){
int T;
scanf("%d",&T);
while(T--){
ll n;
scanf("%lld",&n);
if(Miller_Rabin(n))puts("Yes");
else puts("No");
}
return 0;
}

例题：[http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1186

java写大数：

import java.math.BigInteger;
import java.util.Random;
import java.util.Scanner;
/**
* Created by huchi on 17-7-18.
*/
public class test {
public static final int Times = 10;
public static BigInteger quick_mod(BigInteger a, BigInteger b, BigInteger m){
BigInteger ans = BigInteger.ONE;
a = a.mod(m);
while(!b.equals(BigInteger.ZERO)){
if(b.mod(BigInteger.valueOf(2)).equals(BigInteger.ONE)){
ans = (ans.multiply(a)).mod(m);
b = b.subtract(BigInteger.ONE);
}
b = b.divide(BigInteger.valueOf(2));
a = (a.multiply(a)).mod(m);
}
return ans;
}
public
static boolean Miller_Rabin(BigInteger n){
if(n.equals(BigInteger.valueOf(2)))return true;
if(n.equals(BigInteger.ONE))return false;
if((n.mod(BigInteger.valueOf(2))).equals(BigInteger.ZERO))return false;
BigInteger m = n.subtract(BigInteger.ONE);
BigInteger y = BigInteger.ZERO;
int k = 0;
while((m.mod(BigInteger.valueOf(2))).equals(BigInteger.ZERO)){
k++;
m = m.divide(BigInteger.valueOf(2));
}
Random d = new Random();
for(int i = 0; i < k; i++){
int t = 0;
//随机数
if(n.compareTo(BigInteger.valueOf(10000))==1){
t = 10000;
}else{
t = n.intValue()-1;
}
int a = d.nextInt(t) + 1;
BigInteger x = quick_mod(BigInteger.valueOf(a), m, n);
//测试
for(int j = 0; j < k; j++){
y = (x.multiply(x)).mod(n);
if(y.equals(BigInteger.ONE)
&& !(x.equals(BigInteger.ONE))
&& !(x.equals(n.subtract(BigInteger.ONE)))){
return false;
}
x = y;
}
if(!(y.equals(BigInteger.ONE)))return false;
}
return true;
}
public static void main(String[] args){
Scanner in = new Scanner(System.in);
while(in.hasNextBigInteger()){
BigInteger n = in.nextBigInteger();
if(Miller_Rabin(n))System.out.println("Yes");
else System.out.println("No");
}
}
}

扩展欧几里得

二元一次方程

a>b时，b==0时，x=1,y=0,　否则设 bx2+(a mod b)y2=gcd(b,a mod b)=gcd(a,b)

又有 ax1+by1=bx2+(a−[ab]b)y2

则有 x1=y2,y1=x2−[ab]y2 ,　用递推实现．

void e_gcd(int a, int b, int &d, int &x, int &y){
if(!b){
d=a,x=1,y=0;
return;
}else{
e_gcd(b,a%b,d,y,x);
y -= (a/b)*x;
}
}

代表的是ax+by=gcd(a,b)的解x,y,和d为gcd(a,b);，得到的x,y中a>b则x为正．在gcd(a,b)==1的时候，可以领y=(y%a+a)％a把y变成正的，再根据等式把x变成负的就行．

例题：https://vjudge.net/problem/POJ-2142

#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;
int gcd(int a, int b){
return b?gcd(b,a%b):a;
}
void e_gcd(int a, int b, int &d, int &x, int &y){
if(!b){
d=a,x=1,y=0;
return;
}else{
e_gcd(b,a%b,d,y,x);
y -= (a/b)*x;
}
}
int main(){
int a,b,c;
while(~scanf("%d%d%d",&a,&b,&c) && (a+b+c)){
int x, y, d, tx, ty;
int gd = gcd(a,b);
a /= gd; b/= gd; c/= gd;
e_gcd(a,b,d,x,y);
//获取x非负的时候的情况
tx = x*c;
tx = (tx%b + b)%b;
ty = (c-tx*a)/b;
if(ty < 0)ty = -ty;
//获取y非负的情况
y = y*c;
y = (y%a + a)%a;
x = (c-b*y)/a;
if(x<0)x =-x;
if(x+y>tx+ty){
x=tx;
y=ty;
}
printf("%d %dn",x,y);
}
return 0;
}

https://vjudge.net/problem/HIT-2815

题目没读清，特判就写错了，扩展欧几里得和上面一道题一样的用．

#include<cstdio>
#include<iostream>
#define ll long long
using namespace std;
ll gcd(ll a, ll b){
return b?gcd(b, a%b):a;
}
void e_gcd(ll a, ll b, ll &d, ll &x, ll &y){
if(!b){
d = a; x = 1; y = 0;
return ;
}else{
e_gcd(b,a%b,d,y,x);
y -= (a/b)*x;
}
}
int main(){
ll t,a,b,d,x,y;scanf("%lld",&t);
while(t--){
scanf("%lld%lld",&a,&b);
if(a*b==0){
if(a==1 || b==1){
puts("1");
}else{
puts("-1");
}
}else if(a==1 || b==1){
if(max(a,b)==2)
puts("1");
else puts("2");
}
else if(gcd(a,b)!=1){
puts("-1");
}else{
e_gcd(a,b,d,x,y);
ll tx,ty;
tx = (x%b+b)%b;
ty = (1-a*x)/b;
if(ty<0)ty = -ty;
y = (y%a+a)%a;
x = (1-b*y)/a;
if(x<0)x=-x;
if(x+y>tx+ty){
x = tx;
y = ty;
}
printf("%lldn",x+y-1);
}
}
return 0;
}

毕达哥拉斯三元组的解

x2+y2=z2 ，只考虑本原解．不妨假设x为偶数,y,z为奇数.这样就必有 x=2mn,y=m2−n2,z=m2+n2 ．

例题：https://vjudge.net/problem/POJ-1305

输出两个数据，一个是z不大于n的本原解的个数，和本原解中没有用到的数的个数．

#include<cstdio>
#include<cstring>
#include<iostream>
#define maxn 1000111
using namespace std;
int biao[maxn];
bool flag[maxn];
int gcd(int a, int b){
return b?gcd(b,a%b):a;
}
int main(){
int n;
memset(flag,0,sizeof flag);
for(int i = 2; i*i < maxn; i++){
for(int j = 1; j < i; j++){
if(gcd(i,j)!=1)continue;
if((i&1)&&(j&1))continue;
if(i*i+j*j<maxn)biao[i*i+j*j]++;
}
}
for(int i = 1; i < maxn; i++){
biao[i]+=biao[i-1];
}
while(~scanf("%d",&n)){
memset(flag,0,sizeof flag);
int ans = 0;
for(int i = 2; i*i < maxn; i++){
for(int j = 1; j < i; j++){
if(gcd(i,j)!=1)continue;
if((i&1)&&(j&1))continue;
for(int k = 1; k*(i*i+j*j)<=n; k++){
flag[2*i*j*k]=flag[k*(i*i-j*j)]=flag[k*(i*i+j*j)]=1;
}
}
}
for(int i = 1; i <= n; i++)
if(flag[i])ans++;
printf("%d %dn",biao[n],n-ans);
}
return 0;
}

例题二:

利用毕达哥拉斯三元组＋欧拉函数＋容斥原理

https://vjudge.net/problem/HDU-3939???

寻找满足的i,j满足gcd(i,j)=1,dfs()容斥原理来找有多少互斥的数．

这里需要判断i<＝j和i>j的情况．其中,i<=j时，i为偶数,则直接加phi[i]，i为奇数时,由于j必须为偶数,不能用phi[i]，所以用容斥原理，一共i>>1个数(晒去了２的倍数)．大于时同理．

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#define ll long long
#define maxn 1001000
using namespace std;
bool judge[maxn];
int prime[maxn/10];
int phi[maxn];
int check[40];
int tot,num;
ll ans;
void init(){
tot = 0;
memset(judge, 0, sizeof judge);
for(int i = 2; i < maxn; i++){
if(!judge[i]){
for(int j = i+i; j < maxn; j+=i){
judge[j]=1;
}
prime[tot++]=i;
}
}
for(int i = 1; i < maxn; i++)phi[i]=i;
for(int i = 2; i < maxn; i+=2)phi[i]>>=1;
for(int i = 3; i < maxn; i+=2)if(phi[i]==i){
for(int j = i; j < maxn; j+=i){
phi[j] = phi[j] - phi[j]/i;
}
}
}
void prime_check(int n){
num = 0;
if(!judge[n]){
check[num++]=n;
return ;
}
for(int i = 0; i < tot && n > 1; i++){
if(n%prime[i]==0){
check[num++]=prime[i];
while(n%prime[i]==0) n/=prime[i];
if(n>1 && !judge[n]){
check[num++]=n;
return ;
}
}
}
}
//容斥原理计算无关的数fun(j)=j-(j/c1+j/c2+j/c3)+(j/(c1*c2)+j/(c1*c3)+j/(c2*c3))-(j/(c1*c2*c3))。
void dfs(int k, int r, int s, int n){
if(k==num){
if(r&1) ans -= n/s;
else ans+= n/s;
return;
}
dfs(k+1, r, s, n);
dfs(k+1, r+1, s*check[k], n);
}
ll L;
int main(){
init();
int t;scanf("%d",&t);
while(t--){
ans = 0;
scanf("%lld",&L);
int m = (int)sqrt(1.0*L);
for(int i = m; i > 0; i--){
int p = (int)sqrt(L - (ll)i*i);
//需要寻找gcd(i,j)=1的个数
if(i<=p){
if(i&1){
prime_check(i);
dfs(0,0,1,i>>1);
}else ans+=phi[i];
}else{
// i>p时，j所取的范围在[1,p]
prime_check(i);
if(i&1)dfs(0,0,1,p>>1);
else dfs(0,0,1,p);
}
}
printf("%lldn",ans);
}
}

欧拉函数与欧拉定理

定理一：

定理二: p是素数，有 ϕ(pa)=pa−pa−1 ,容斥原理．

定理三: n=pα11pα22...pαkk时，ϕ(n)=n(1−1p1)(1−1p2)...(1−1pk) 　定理１＋３

定理四: ∑d|nϕ(d)=n 莫比乌斯反演．

定理五: (a,m)=1,则 x≡baϕ(m)−1是同于方程ax≡b(mod m)

定理六: n>2时， ϕ(n) 为偶数

定理三求解欧拉函数：

int phi(int n){
int ans = n;
for(int i = 2; i*i <= n; i++){
if(n%i==0){
ans = ans - ans/i;
while(n%i==0) n/=i;
}
}
if(n>1)ans = ans-ans/n;
return ans;
}

容斥原理，筛法去减不满足的．

#define maxn 100000
int p[maxn];
void phi_table(){
memset(p,0,sizeof p);
for(int i = 1; i < maxn; i++)p[i]=i;
for(int i = 2; i < maxn; i+=2)p[i]>>=1;
for(int i = 3; i < maxn; i+=2){
if(p[i]==i){
for(int j = i; j < maxn; j+=i)
p[j]=p[j]-p[j]/i;
}
}
}

二次筛法

这道题有类似的在kuangbin的数论基础专题里有出现过．

例题：https://vjudge.net/problem/POJ-2689

注意1　2这个样例．．．还有晒的时候不能把质数也晒了．．．．

#include<cstdio>
#include<iostream>
#include<cstring>
#define maxn 50000
int prime[maxn];
bool judge2[21*maxn];
bool judge1[maxn];
using namespace std;
int tot=0;
void init(){
memset(judge1,0,sizeof judge1);
for(int i = 2; i<maxn; i++){
if(!judge1[i]){
prime[tot++] = i;
for(int j = i+i; j < maxn; j+=i){
judge1[j]=true;
}
}
}
}
int main(){
init();
unsigned int u,v;
while(~scanf("%d%d",&u,&v)){
if(u==1)u=2;
memset(judge2,0,sizeof judge2);
for(int i = 0; i < tot; i++){
if(prime[i]>v)break;
unsigned int j=u/prime[i]*prime[i];
if(u%prime[i]==0 && u!=prime[i])judge2[0]=1;
while(j<=u)j+=prime[i];
if(j==prime[i])j+=j;
for(;j<=v;j+=prime[i]){
judge2[j-u]=1;
}
}
int st,stlen=maxn,ed,edlen=0;
int cot=0,len=0;
for(int i = 0; i <= v-u; i++){
if(!judge2[i]){
if(len &&len<stlen){
st = i-len;
stlen = len;
}
if(len && len>edlen){
ed = i-len;
edlen = len;
}
cot++;
len = 0;
}
len ++ ;
}
if(cot!=2 && st==ed){
ed = st+stlen;
}
if(cot>1)printf("%d,%d are closest, %d,%d are most distant.n",u+st,u+st+stlen,u+ed,u+ed+edlen);
else puts("There are no adjacent primes.");
}
return 0;
}

逆元

正整数a,m，如果有 ax≡1(mod m) ,这个同于方程中x的最小正整数解叫做a模m的逆元．

逆元一般用扩展欧几里得算法来求．m为素数时，由费马小定理得到逆元为 am−2modm

逆元常见问题

已知b|a，求解 ans=abmod m

常见解法:扩展欧几里得．m是素数时，直接用费马小定理．但他们都需要a与m互素．

除法取模求逆元的通用方法： ans=a mod(mb)b

求 AB 的所有因子和对9901取余的结果．

考虑

所有的因子和：

可以用二分求等比数列和，也可以用等比数列公式(需要用逆元)直接求.分别写在了Method_one和Method_two中．

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#define ll long long
#define maxn 10000
using namespace std;
ll mod =9901;
bool judge[maxn];
ll check[maxn][2];
int prime[maxn/2];
int tot,num;
ll a,b;
void init(){
tot=0;
memset(judge, 0, sizeof judge);
for(int i = 2; i < maxn; i++){
if(!judge[i]){
for(int j = i+i; j < maxn; j+=i){
judge[j]=1;
}
prime[tot++]=i;
}
}
}
ll multi(ll a, ll b, ll m){
ll ans = 0;
a %= m;
while(b){
if(b&1)ans = (ans+a)%m;
a = (a+a)%m;
b >>= 1;
}
return ans;
}
ll largequickmod(ll p, ll k, ll m){
ll ans = 1;
while(k){
if(k&1)ans = multi(ans,p,m);
p = multi(p,p,m);
k>>=1;
}
return ans;
}
ll quickmod(ll p, ll k){
ll ans=1;
p%=mod;
while(k){
if(k&1)ans = ans*p%mod;
p = (p*p)%mod;
k>>=1;
}
return ans;
}
ll S(ll p, ll k){
if(p==0)return 0;
if(k==0)return 1;
ll temp = quickmod(p,(k>>1))%mod;
if(k&1){
return S(p,k>>1)%mod*(1+temp*p%mod)%mod;
}else{
return (S(p,(k>>1)-1)%mod*(1+temp*p%mod)%mod+temp)%mod;
}
}
void Method_one(){
ll ans = 1;
for(int i = 0; i < num; i++){
ans = (ans * S(check[i][0], b*check[i][1]))%mod;
}
if(a==0)ans=0;
printf("%lldn",ans);
}
void Method_two(){
ll ans = 1;
for(int i = 0; i < num; i++){
ll M = mod*(check[i][0]-1);
ans *= (largequickmod(check[i][0], check[i][1]*b+1, M)+M-1)/(check[i][0]-1);
ans %= mod;
}
printf("%lldn",ans);
}
int main(){
init();
while(~scanf("%lld%lld",&a,&b)){
ll temp = a;
num=0;
memset(check, 0, sizeof check);
for(int i = 0; i < tot; i++){
if(temp<=1)break;
if(temp%prime[i]==0){check[num++][0]=prime[i];}
while(temp%prime[i]==0){
temp/=prime[i];
check[num-1][1]++;
}
}
if(temp>1){
check[num][0]=temp;
check[num++][1]=1;
}
//Method_one();
Method_two();
}
return 0;
}

1→p模p的所有逆元，这里p为奇质数．用快速幂的复杂度Ｏ(plog(p)),用递推式可以优化到Ｏ(p)：

http://www.lydsy.com/JudgeOnline/problem.php?id=2186

例题：1-N!中与M!互质的数的个数，M<=N;M!|N!, ．当n是m的倍数时，1-n中与m互素的个数为 nmϕ(m)

题目中　 N!M!ϕ(M!)=N!M!M!(1−1p1)(1−1p2)...(1−1pk) .　　

时间卡的特别紧．bool开的就超时了．inv[i]算逆元．multi[maxn]算阶乘取模，要求p与m互质，否则需要用除法取模．

#include<cstdio>
#include<bitset>
#define ll long long
#define maxn 10000005
using namespace std;
ll mod;
bitset<maxn> judge;
void init(){
judge.set();
for(int i = 4; i < maxn; i+=2)judge[i]=0;
for(int i = 3; i < maxn; i+=2){
if(judge[i]){
for(int j = i+i; j < maxn; j+=i){
judge[j]=0;
}
}
}
}
ll multi[maxn];
int inv[maxn];
ll ans[maxn];
int main(){
init();
int T,R;
scanf("%d%d",&T,&R);
mod = R;
multi[0]=1;
for(int i = 1; i < maxn; i++){
multi[i]=multi[i-1]*i%mod;
}
inv[1]=1;
for(int i = 2; i < mod && i < maxn; i++){
inv[i]=(mod-mod/i)*inv[mod%i]%mod;
}
ans[1]=1;
for(int i = 2; i < maxn; i++){
if(judge[i]){
ans[i] = ans[i-1]*(i-1)%mod;
ans[i] = ans[i]*inv[i%mod]%mod;
}else{
ans[i] = ans[i-1];
}
}
int n,m;
while(T--){
scanf("%d%d",&n,&m);
ll anss = 1LL*multi[n]*ans[m]%mod;
printf("%lldn",anss);
}
return 0;
}

欧拉筛

void getphi()
{
int i,j;
phi[1]=1;
for(i=2;i<=N;i++)//相当于分解质因式的逆过程

{
if(!mark[i])
{
prime[++tot]=i;//筛素数的时候首先会判断i是否是素数。

phi[i]=i-1;//当 i 是素数时 phi[i]=i-1

}
for(j=1;j<=tot;j++)
{
if(i*prime[j]>N)
break;
mark[i*prime[j]]=1;//确定i*prime[j]不是素数

if(i%prime[j]==0)//接着我们会看prime[j]是否是i的约数

{
phi[i*prime[j]]=phi[i]*prime[j];break;
}
else
phi[i*prime[j]]=phi[i]*(prime[j]-1);//其实这里prime[j]-1就是phi[prime[j]]，利用了欧拉函数的积性

}
}
}

最后

以上就是幸福哑铃为你收集整理的Acdreamer博客数论学习(1)Acdreamer博客数论学习Day OneDay Two的全部内容，希望文章能够帮你解决Acdreamer博客数论学习(1)Acdreamer博客数论学习Day OneDay Two所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。