您的位置:首页 > 其它

数论整理DAY1

2016-11-28 22:59 176 查看

整除性

1.如果m>0且比值n/m是一个整数,那么就是m整除n(或者说n被m整除),记作m|n.

2.两个数的最大公因子(greatest common divisor)是能整除它们的最大整数.记作gcd(a,b).

3.最小公倍数(least common multiple),即最小的能同时被两数整除的数,记作lcm(a,b).


最大公约数

求法:辗转相除法(自行百度).容易证出求最小公倍数的复杂度是log级别的.


int gcd(int a,int b){
return b?gcd(b,a%b):a;
}


最小公倍数

lcm(a,b)=a∗bgcd(a,b)

扩展欧几里得

解方程ax+by=1,gcd(a,b)=1.即ax+by=1=gcd(a,b)

而gcd(a,b)=gcd(b,a%b).那么bx+(a%b)y=gcd(b,a%b)

最终有:gcd(1,0)=1∗x+0∗y=1,此组解x=1y=0

再返回迭代回初始的解。

ax+by=gcd(a,b)=gcd(b,a%b)=bx′+(a%b)y′=bx′+(a−a÷b∗b)y′

bx′+(a−a÷b∗b)y′=ay′+b(x′−a÷b∗y′)

x=y′,y=x′−a÷b∗y′,迭代回去就好了.

同余方程

求a∗xmodb=1的最小正整数解,题目中有gcd(a,b)=1

var a,b,x,y:longint;
function gcd(a,b:longint):longint;
var t,tp:longint;
begin
if b=0 then begin
x:=1;
y:=0;
exit(a);
end;
tp:=gcd(b,a mod b);
t:=x;
x:=y;
y:=t-a div b*y;
exit(tp);
end;
begin
read(a,b);
gcd(a,b);
if x<0 then inc(x,b);
writeln(x);
end.


对于gcd(a,b)>1的来说,直接同时除以gcd(a,b),然后处理a’x+b’y=1即可。

事实上,这里所求的x即a模b意义下的逆元.这一点在处理带除法的的模运算的时候会用到.

即abmodp=a∗b−1,b−1即b模p下的逆,可用于加速运算等等.

hdu5976

这里用到逆元加速运算。具体方法这里不再赘述,请看题解

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <bitset>
using namespace std;
typedef long long LL;
const LL mod=1e9+7;
const LL maxn=1e5+5;
LL A[maxn],s[maxn];
LL Inv[maxn];
void getINV(){
Inv[1]=1;
for(LL i=2; i<maxn; i++)
Inv[i] = (mod-mod/i)*Inv[mod%i]%mod;
}
void init(){
s[1]=0;
for(LL i=2;i<maxn;i++)
s[i]=s[i-1]+i;
A[0]=1;
for(LL i=1;i<maxn;i++)
A[i]=A[i-1]*i%mod;
getINV();
}

int main(){
init();
int T;
cin>>T;
while(T--){
int pos;
LL x,sum;
scanf("%lld",&x);
LL n=(LL)(sqrt(2*x+2));
for(int i=n;i>=1;i--)
if(s[i]<=x){
pos=i;
break;
}
if(x==1) { puts("1"); continue; }
if(s[pos]==x) { printf("%lld\n",A[pos]); continue; }
if(s[pos]+pos-1>=x){
LL tot=s[pos]+pos-1-x;
sum=(A[pos+1]*Inv[tot+2])%mod;
}
else{
LL tot=s[pos+1]-2+pos-1-x;
sum=((A[pos+2]*Inv[2])%mod)*Inv[tot+3]%mod;
}
printf("%lld\n",sum);
}
return 0;
}


费马小定理

gcd(a,p)=1则ap−1modp=1

若a,p互质且a

int inv(int a,int m){
if(a==1) return 1;//是1不是l
return inv(m%a,m)*(m-m/a)%m;
}


素数

素数筛(判断素数)

素数筛选,判断小于MAXN的数是不是素数.notprime是一张表,为false表示是素数,true表示不是素数.


const int MAXN=1000001;
bool notprime[MAXN];
void init(){
memset(notprime,false,sizeof(notprime));
notprime[0]=notprime[1]=true;
for(int i=2;i<MAXN;i++) if(!notprime[i]){
if(i>MAXN/i)continue;
for(int j=i*i;j<MAXN;j+=i) notprime[j]=true;
}
}


素数筛(找出素数)

素数筛选小于等于MAXN的素数


const int MAXN=10001;
int prime[MAXN];
void init(){
memset(prime,0,sizeof(prime));
for(int i=2;i<MAXN;i++){
if(!prime[i])prime[++prime[0]]=i;
for(int j=1;j<=prime[0]&&prime[j]<=MAXN/i;j++){                       prime[prime[j]*i]=1;
if(i%prime[j]==0) break;
}
}
}


分解质因子

void solve(){
for(int i=1;i<=prime[0];i++){//prime[]为上段中的prime
top++;
if(tmp%prime[i]!=0) continue;
while(tmp%prime[i]==0){
s[top]++;
tmp/=prime[i];
if(tmp==1) break;
}
}
if(tmp>1){
top++;
s[top]++;//s存质数,有必要的话再开个数组存因子.
}
}


MR算法

假如p是质数,且gcd(a,p)=1,那么ap−1≡1(mod p)。即假如p是质数,且a,p互质,那么ap−1modp=1。(费马小定理)

该定理的逆命题是不一定成立的,但大多数情况是成立的。

网上套来的代码:

typedef long long llt;
//利用二进制计算a*b%mod
llt multiMod(llt a,llt b,llt mod){
llt ret = 0LL;
a %= mod;
while( b ){
if ( b & 1LL ) ret = ( ret + a ) % mod, --b;
b >>= 1LL;
a = ( a + a ) % mod;
}
return ret;
}

//计算a^b%mod
llt powerMod(llt a,llt b,llt mod){
llt ret = 1LL;
a %= mod;
while( b ){
if ( b & 1LL ) ret = multiMod(ret,a,mod),--b;
b >>= 1LL;
a = multiMod(a,a,mod);
}
return ret;
}

//Miller-Rabin测试,测试n是否为素数
bool Miller_Rabin(llt n,int repeat){
if ( 2LL == n || 3LL == n ) return true;
if ( !( n & 1LL ) ) return false;

//将n分解为2^s*d
llt d = n - 1LL;
int s = 0;
while( !( d & 1LL ) ) ++s, d>>=1LL;

srand((unsigned)time(0));
for(int i=0;i<repeat;++i){//重复repeat次
llt a = rand() % ( n - 3 ) + 2;//取一个随机数,[2,n-1)
llt x = powerMod(a,d,n);
llt y = 0LL;
for(int j=0;j<s;++j){
y = multiMod(x,x,n);
if ( 1LL == y && 1LL != x && n-1LL != x ) return false;
x = y;
}
if ( 1LL != y ) return false;
}
return true;
}


Meisell-Lehmer

求N以内素数个数,只知道模板,原理不懂.不献丑了.

hdu 5901

//Meisell-Lehmer
//G++ 218ms 43252k
#include<cstdio>
#include<cmath>
using namespace std;
#define LL long long
const int N = 5e6 + 2;
bool np
;
int prime
, pi
;
int getprime()
{
int cnt = 0;
np[0] = np[1] = true;
pi[0] = pi[1] = 0;
for(int i = 2; i < N; ++i)
{
if(!np[i]) prime[++cnt] = i;
pi[i] = cnt;
for(int j = 1; j <= cnt && i * prime[j] < N; ++j)
{
np[i * prime[j]] = true;
if(i % prime[j] == 0)   break;
}
}
return cnt;
}
const int M = 7;
const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
int phi[PM + 1][M + 1], sz[M + 1];
void init()
{
getprime();
sz[0] = 1;
for(int i = 0; i <= PM; ++i)  phi[i][0] = i;
for(int i = 1; i <= M; ++i)
{
sz[i] = prime[i] * sz[i - 1];
for(int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
}
}
int sqrt2(LL x)
{
LL r = (LL)sqrt(x - 0.1);
while(r * r <= x)   ++r;
return int(r - 1);
}
int sqrt3(LL x)
{
LL r = (LL)cbrt(x - 0.1);
while(r * r * r <= x)   ++r;
return int(r - 1);
}
LL getphi(LL x, int s)
{
if(s == 0)  return x;
if(s <= M)  return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
if(x <= prime[s]*prime[s])   return pi[x] - s + 1;
if(x <= prime[s]*prime[s]*prime[s] && x < N)
{
int s2x = pi[sqrt2(x)];
LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
for(int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
return ans;
}
return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
}
LL getpi(LL x)
{
if(x < N)   return pi[x];
LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
return ans;
}
LL lehmer_pi(LL x)
{
if(x < N)   return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(sqrt3(x));
LL sum = getphi(x, a) +(LL)(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++)
{
LL w = x / prime[i];
sum -= lehmer_pi(w);
if (i > c) continue;
LL lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
}
return sum;
}
int main()
{
init();
LL n;
while(~scanf("%lld",&n))
{
printf("%lld\n",lehmer_pi(n));
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: