天天看点

szu 寒训复习day #4数论入门详解[修改]用Mardown再改了一下

整数的取余运算

定义:带余除法 设 a , b a,b a,b是整数,且 b > 0 b>0 b>0,则存在非负整数 q , r q,r q,r ,使得

a = b q + r a=bq+r a=bq+r

且 0 ≤ r < b 0≤r<b 0≤r<b ,称 q 为商, r 为余数

显然带余除法中的商和余数都是唯一的,在下文中将商记为 a / b a/b a/b ,将余数记为 a % b a\%b a%b,“/”与“%”的运算优先级与乘除法相同,当然,在C语言中二者分别对应 a / b a/b a/b 与 a % b a\%b a%b

1.1.1 定理 整数的加、减、乘对于取余的右分配律

( a + b ) % c = ( a % c + b % c ) % c (a+b)\%c=(a\%c+b\%c)\%c (a+b)%c=(a%c+b%c)%c

( a − b ) % c = ( a % c − b % c + c ) % c (a-b)\%c=(a\%c-b\%c+c)\%c (a−b)%c=(a%c−b%c+c)%c

( a b ) % c = ( a % c ) ( b % c ) % c (ab)\%c=(a\%c)(b\%c)\%c (ab)%c=(a%c)(b%c)%c

注意,在计算减法时,通常需要加 c ,防止变成负数,因为在取模运算中是没有负数的

在计算乘法时,如果 c 较大(但不超过64位整数范围),可以使用快速乘法进行计算,原理与快速幂运算类似,复杂度为 O(logb)

下面瞅代码

ll fastMul(ll a,ll b,ll p){
    a%=p;
    ll ans=0;
    while(b>0){
        if(b&1)ans=(ans+a)%p;
        b>>=1;
        a=(a+a)%p;
    }
    return ans;
}
           

如果需要更快的快速乘法,可以使用 long double 数据类型进行计算,复杂度为 O(1)

(背反正我也不知道是啥qwq)

ll modMul(ll a,ll b,ll p){
    if (p<=1000000000) return a*b%p;
    else if (p<=1000000000000ll) return (((a*(b>>20)%p)<<20)+(a*(b&((1<<20)-1))))%p;
    else {
        ll d=(ll)floor(a*(long double)b/p+0.5);
        ll ret=(a*b-d*p)%p;
        if (ret<0) ret+=p;
        return ret;
    }
}
           

乘法逆元

虽然取余运算对于“+”、“-”、“×”不难,但通常情况下

a / b % c ≠ ( a % c ) / ( b % c ) % c a / b \% c ≠ (a\%c) / (b\%c) \%c a/b%c​=(a%c)/(b%c)%c

如何计算 a / b % c a/b\%c a/b%c ?我们有时可以找一个乘法逆元 b − 1 b^{-1} b−1 ,使得 b ∗ b − 1 ≡ 1 ( m o d c ) b*b^{-1}≡1(mod c) b∗b−1≡1(modc) ,那么就有

a / b % c = a b − 1 % c a/b\%c=ab^{-1}\%c a/b%c=ab−1%c

如果 c 是素数,有下面的定理

定理 费马小定理 设 b 是一个整数,c 是一个素数,二者互质,那么

b c − 1 ≡ 1 ( m o d c ) b^{c-1}≡1(mod c) bc−1≡1(modc)

将上式改写一下,得到

b ∗ b ( c − 2 ) ≡ 1 ( m o d c ) b*b^{(c-2)}≡1(mod c) b∗b(c−2)≡1(modc)

因此取 b − 1 = b c − 2 b^{-1}=b^{c-2} b−1=bc−2 即可,一般需要用快速幂计算,但要注意,与除数不能为0类似,要保证 b % c ≠ 0 b\%c≠0 b%c​=0

ll inv(ll a,ll p){
    return fpow(a,p-2);
}
ll fpow(ll a, ll b)
{
    ll ans = 1;
    while(b)
    {
       if(b & 1) ans = ans * a % p;
       a = a * a % p;
       b >>= 1;
    }
    return ans % p;
}
           

2. 最大公因数与最小公倍数

2.1 最大公因数

顾名思义,最大公因数就是公因数中最大的那个,我们记 a,b 的最大公因数为 gcd(a,b) ,有如下性质

2.1.1 性质 gcd(a,b)=gcd(b,a)

2.1.2 性质 gcd(a,b)=gcd(a-b,b)(a≥b)

2.1.3 性质 gcd(a,b)=gcd(a%b,b)

2.1.4 性质 gcd(a,b,c)=gcd(gcd(a,b),c)

2.1.5 性质 gcd(ka,kb)=k gcd(a,b)

2.2 辗转相除法

根据2.1.3 性质,得到辗转相除法的参考代码模板

typedef long long ll;
ll gcd(ll a, ll b){
    return b? gcd(b,a%b)  : a;
}
           

注意当 b≠0 时,返回值为 gcd(b,a%b) 而不是 gcd(a%b,b) ,否则会不断递归导致栈溢出

2.3 最小公倍数

最小公倍数就是公倍数中最小的那个,我们记 a,b 的最小公倍数为 lcm(a,b) ,有如下性质

2.3.1 性质

lcm(a,b)=ab/(gcd(a,b))

下面是最小公倍数的参考代码模板

ll lcm(ll a, ll b){
    return a/gcd(a,b)*b;
}
           

注意: 是先除后乘,避免在中间过程中数据超出64位整数的范围

2.4 扩展欧几里得算法

考虑二元一次不定方程

a x + b y = c ax+by=c ax+by=c

其中 a,b,c 是已知的正整数,如何求出方程的解呢?

**2.4.1 定理 上述方程有解的充要条件是 g c d ( a , b ) ∣ c gcd(a,b)|c gcd(a,b)∣c (c 是 gcd(a,b) 的倍数)

可以理解为, g c d ( a , b ) gcd(a,b) gcd(a,b) 是 a x + b y ax+by ax+by 可以表示出的最小正整数

裴蜀定理 : 对任意一对正整数 a , b a ,b a,b一定存在非零 x , y x ,y x,y使得 a ∗ x + b ∗ y = g c d ( a , b ) a * x + b * y = gcd(a,b) a∗x+b∗y=gcd(a,b); a ∗ x + b ∗ y = d a * x + b * y = d a∗x+b∗y=d 则 d d d就是 a a a和 b b b的最大公约数的倍数 ;

1.设 a x 1 + b y 1 = g c d ( a , b ) , b x 2 + ( a % b ) y 2 = g c d ( b , a % b ) ; ax_1+by_1=gcd(a,b), bx_2+(a\%b)y_2=gcd(b,a\%b); ax1​+by1​=gcd(a,b),bx2​+(a%b)y2​=gcd(b,a%b);

2.由 g c d ( a , b ) = g c d ( b , a % b ) gcd(a,b)=gcd(b,a\%b) gcd(a,b)=gcd(b,a%b),可得: a x 1 + b y 1 = b x 2 + ( a % b ) y 2 ax_1+by_1=bx_2+(a\%b)y_2 ax1​+by1​=bx2​+(a%b)y2​;

3.即: a x 1 + b y 1 = b x 2 + ( a − ⌊ a b ⌋ ) ∗ b ) y 2 = a y 2 + b x 2 − ⌊ a b ⌋ ∗ b y 2 ; ax_1+by_1=bx_2+(a-\lfloor {a\over b} \rfloor) * b)y_2 =ay_2+bx_2-\lfloor {a\over b} \rfloor * by_2; ax1​+by1​=bx2​+(a−⌊ba​⌋)∗b)y2​=ay2​+bx2​−⌊ba​⌋∗by2​;

即: a x 1 + b y 1 = a y 2 + b ( x 2 − ⌊ a b ⌋ ∗ y 2 ) ax_1+by_1=ay_2 + b(x_2-\lfloor {a\over b} \rfloor* y_2) ax1​+by1​=ay2​+b(x2​−⌊ba​⌋∗y2​)

求解递归上传下的规律 根据恒等定理,对应项相等,得: x 1 = y 2 ; y 1 = x 2 − ⌊ a b ⌋ ∗ y 2 x_1=y_2; y_1=x_2-\lfloor {a\over b} \rfloor * y_2 x1​=y2​;y1​=x2​−⌊ba​⌋∗y2​;

每次x y都会倒过来一次 这样我们就得到了:x1,y1的值基于x2,y2,所以我们可以通过递归求解**

下面是参考代码模板:

ll ext_gcd(ll a,ll b,ll& x,ll& y){
    if (!b)
    {
        x = 1;y = 0;
        return a;
    }
    ll d = ext_gcd(b, a%b, y, x);
     y -= a / b * x;
    return d;
}
           

扩展欧几里得算法板子题

在求逆元时,要找到 b − 1 b^{-1} b−1 使得 b ∗ b − 1 ≡ 1 ( m o d c ) b*b^{-1}≡1(mod c) b∗b−1≡1(modc) ,实质上是求解方程 b x + c y = 1 bx+cy=1 bx+cy=1 中的 x ,因此可以用扩展欧几里得算法来求逆元,当然只有 g c d ( b , c ) = 1 gcd(b,c)=1 gcd(b,c)=1 时才有解,否则逆元不存在

2.4.4 推论 逆元的存在性充要条件是 gcd(b,c)=1

扩展欧几里得算法解决线性同余方程

szu 寒训复习day #4数论入门详解[修改]用Mardown再改了一下

线性同余方程组

szu 寒训复习day #4数论入门详解[修改]用Mardown再改了一下

中国剩余定理不互质情况

x ≡ ( a 1 m o d m 1 ) 和 x ≡ ( a 2 m o d m 2 ) x \equiv (a_1 mod m_1)和x \equiv (a_2 mod m_2) x≡(a1​modm1​)和x≡(a2​modm2​);

变形一下

x = k 1 a 1 + m 1 ; x = k 2 a 2 + m 2 x = k_1a_1 + m_1;x = k_2a_2 + m_2 x=k1​a1​+m1​;x=k2​a2​+m2​

联立方程式

k 1 ∗ a 1 + m 1 = k 2 ∗ a 2 + m 2 k_1 * a_1 + m_1 = k_2 * a_2 + m_2 k1​∗a1​+m1​=k2​∗a2​+m2​

相同项移到一边

k 1 ∗ a 1 − k 2 ∗ a 2 = m 1 − m 2 k_1 * a_1 - k_2 * a_2 = m_1 - m_2 k1​∗a1​−k2​∗a2​=m1​−m2​

给定一个 a 1 a 2 a_1 a_2 a1​a2​(相当于 a ,b) m 1 − m 2 m_1 - m_2 m1​−m2​ (相当于)d

转化为扩展欧几里得算法 ( m 1 − m 2 ) ∣ g c d ( a 1 , a 2 ) (m_1 - m_2)| gcd(a_1,a_2) (m1​−m2​)∣gcd(a1​,a2​)

假设现在我们求出来了 k 1 和 k 2 k_1和k_2 k1​和k2​

a ∗ x + b ∗ y = d a * x + b * y = d a∗x+b∗y=d其中 d d d是最大公约数

{ x 0 = k 1 + k ( a 2 / g c d ( a 1 , a 2 ) ) y 0 = k 2 − k ( a 1 / g c d ( a 1 , a 2 ) ) \begin{cases} x_0 = k1 + k (a2 / gcd(a1,a2))\\ y_0 = k2 - k (a1 / gcd(a1,a2)) \end{cases} {x0​=k1+k(a2/gcd(a1,a2))y0​=k2−k(a1/gcd(a1,a2))​

上面是 k 1 ∗ a 1 − k 2 ∗ a 2 = m 1 − m 2 k1 * a1 - k2 * a2 = m1 - m2 k1∗a1−k2∗a2=m1−m2通解公式

x = k 1 a 1 + m 1 x = k1 a1 + m1 x=k1a1+m1;

带入公式

x = ( k 1 + k ∗ ( a / b ) ) ∗ a 1 + m 1 ; x = (k1 + k * (a / b)) * a1 + m1; x=(k1+k∗(a/b))∗a1+m1;

x = a 1 ∗ k 1 + m 1 + k ∗ [ a 1 , a 2 ] x = a1 * k1 + m1 + k * [a1,a2] x=a1∗k1+m1+k∗[a1,a2];//把两个不定方程合并

x = x 0 + k ∗ a x = x0 + k * a x=x0+k∗a则对比一下一开始的 x = k 1 a 1 + m 1 x = k1a1 + m1 x=k1a1+m1 由对应项相同可得 m 1 = x 0 , a = l c a [ a 1 , a 2 ] m1 = x0,a = lca[a1,a2] m1=x0,a=lca[a1,a2] 一直递归迭代下去就可以求解咯

中国剩余模板题传送门

看一下核心部分

1.先通过扩展欧几里得算法求出最大公约数d

2.判断一下a[i] - a[0] 是否为d的倍数如果不是就直接break掉了该方程组无解

3.如果有解对于k1 * a1 - k2 * a2 = m1 - m2 我们通过exgcd求解的是最小公倍数所以对于k1 我们要扩展到相同倍数

4.通过上面我们知道k1+a2/d为一组解则(k1%t+t)%t是为了求最小正整数解

5.最后一步更新x=k1a1+m1

m1=a1*k1+m1;a1=lcm(a1,a2)

下面a和m反过来了

for(int i = 1; i < n; ++ i)
     	{
     		LL d = exgcd(m[0],m[i],k1,k2);
     		if((a[i] - a[0]) % d) 
     		{
     		   flag = true;
				break;	
		    }
		   //现在我们要求的是 a2 - a1 等式两边乘以一个倍数
          k1 *= (a[i] - a[0]) / d;//
		  LL t = m[i] / d;
		  k1 = (k1 % t + t) % t;
		  a[0] = m[0] * k1 + a[0];
		  m[0] = abs(m[i] / d * m[0]);
		}
           

还有一点这道题有个坑就是它的余数有全为0的情况就说明它们是等比数列

我们就要特判一下假如a[0] == 0的话我们就要从求遍开始没考虑wa了好多发后来自己看来一下边界情况qwq

#include <iostream>
#include <cmath>
using namespace std;
typedef long long LL;
LL m[110];
LL a[110];

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

int main()
{
  int T;
  cin >> T;
  int C = 1;
  while(T --)
  {
     	int n;
     	cin >> n;
     	for(int i = 0; i < n; ++ i)
     	  cin >> m[i];
     	for(int j = 0; j < n; ++ j)
     	  cin >> a[j];
     	  bool flag = false;
     	  LL k1, k2;//表示 x ,y 取什么值时为gcd 
     	  int t = a[0];
     	for(int i = 1; i < n; ++ i)
     	{
     		LL d = exgcd(m[0],m[i],k1,k2);
     		if((a[i] - a[0]) % d) 
     		{
     		   flag = true;
				break;	
		    }
		   //现在我们要求的是 a2 - a1 等式两边乘以一个倍数
          k1 *= (a[i] - a[0]) / d;
		  LL t = m[i] / d;
		  k1 = (k1 % t + t) % t;
		  a[0] = m[0] * k1 + a[0];
		  m[0] = abs(m[i] / d * m[0]);
		}
		
		if(a[0] == 0)
		{
			a[0] = 1;
		  for(int i = 0; i < n; ++ i)
		    a[0] = a[0] * m[i] / exgcd(a[0],m[i],k1,k2);
 		} 
     	if(flag) printf("Case %d: -1\n",C ++);
     	else printf("Case %d: %lld\n",C ++ ,a[0]);
   }	
} 
           

4. 素数

素数是只有1和本身两个因数的数,1不是素数

4.1 素数的判断

bool isPrime(ll n){
    if(n==1) return false;
    for(ll i = 2; i * i <= n; i ++)
        if(n % i == 0)return false;
    return true;
}
           

这是最简单的素数的判断的参考代码模板,复杂度为 O(√n)

原理其实很简单,对于一个大于1的整数,如果 x 是它的一个大于 √n 的因子,那么 n/x 是它的小于 √n 的因子

在大多数情况下,这种判断方式的复杂度已经足够小了,如果要追求更高的效率,可以考虑 kn+i 法

一个大于1的整数如果不是素数,那么一定有素因子,因此在枚举因子时只需要考虑可能为素数的因子即可。 kn+i 法即枚举形如 kn+i 的数,例如取 k=6那么 6n+2,6n+3,6n+4,6n+6 都不可能为素数,只需要枚举形如 6n+1,6n+5 的数即可,这样复杂度降低了 2/3 。下面的模板是 kn+i 法 k=30 的版本

bool isPrime(ll n){
    if(n == 2 || n == 3 || n == 5)return 1;
    if(n % 2 == 0 || n % 3 == 0|| n % 5 == 0||n == 1)return 0;
    ll c=7,a[8]={4,2,4,2,4,6,2,6};
    while(c * c <= n)
     for(auto i : a)
    {
      if(n % c == 0)return 0;
      c += i;
    }
    return 1;
}
           

如果 n 极大,可以使用素数测试算法,素数测试算法可以通过控制迭代次数来间接控制正确率,常用的有下面的Miller-Rabin方法

想要深入解看这篇博客

传送门

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
LL num[3]={2,7,61};
LL mul(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;
}
LL powm(LL a,LL b,LL m){
    LL ans=1;
    while(b){
        if(b&1)ans=mul(ans,a,m);
        a=mul(a,a,m);
        b>>=1;
    }
    return ans;
}
bool test(LL m,LL a,LL d){
    if(!(m&1))return false;
    while(!(d&1))d>>=1;
    LL t=powm(a,d,m);
    while(d!=m-1&&t!=m-1&&t!=1){
        t=mul(t,t,m);
        d<<=1;
    }
    return t==m-1||(d&1);
}
bool isprime(LL m){
    for(int i=0;i<3;i++){
        if(m==num[i])return true;
        if(!test(m,num[i],m-1))return false;
    }
    return true;
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    LL N;
    cin>>N;
    if(N==1) cout<<"No\n";
    else{
        if(isprime(N))cout<<"Yes\n";
        else cout<<"No\n";
    }
return 0;}
           

4.2 素数筛

如果要求出不超过 n 的所有素数,素数筛是最好的选择,下面是一种朴素的筛法

void getPrime(bool p[],int n){
    for(int i=1;i<=n;i++)p[i]=true;
    p[1]=false;
    for(int i=2;i<=n;i++){
        if(p[i]){
            for(int j=i+i;j<=n;j+=i)p[j]=false;
        }
    }
}
           

这种方法的原理是从小到大将素数的倍数筛掉,复杂度为 O(nlogn) ,注意到每个合数如果有多个素因子,那么就会被重复筛掉,造成复杂度的浪费,因此,用下面的方法可以保证每个合数只被它最小的素因子筛掉一遍,以 O(n) 的复杂度解决上述问题

欧拉筛

ll getPrime(ll n,bool vis[],ll prime[]){
    ll tot=0;
    for(ll i=1;i<=n;i++)vis[i]=0;
    for(ll i=2;i<=n;i++){
        if(!vis[i])prime[tot++]=i;
        for(ll j=0;j<tot;j++){
            if(prime[j]*i>n)break;
            vis[prime[j]*i]=1;
            if(i%prime[j]==0)break;
        }
    }
    return tot;
}

           

高斯消元

上学期刚学线性代数现在还热乎qwq

szu 寒训复习day #4数论入门详解[修改]用Mardown再改了一下

想起上学期算线性方程就想吐太容易算错了qwq

现在福利来了写一个场程序自动计算线性方程组就是高斯消元的办法

下面我们来看一下代码的思路

前置知识:初等行(列)变换

1.把某一行乘一个非0的数 (方程的两边同时乘上一个非0数不改变方程的解)交换某两行 (交换两个方程的位置)

2.把某行的若干倍加到另一行上去 (把一个方程的若干倍加到另一个方程上去)

接下来,运用初等行变换,把增广矩阵,变为阶梯型矩阵

将多元组的系数抽出来

形成系数矩阵

将系数矩阵进行初等行变换

线性代数

等价变化

1.完美阶梯型 唯一解

2.0 = 非0 无解;

3.0 = 0 无穷多组解

高斯消元

枚举每一列 c

1.找到当前这一列绝对值最大的一行

2.将该行换到最上面

3.将该行第一个数变成1

4.将下面所有行的第c列消成0

5.第一个方程处理完就固定

继续迭代

#include <iostream>
#include <algorithm>
#include <cmath>

using namespace std;
const int N = 110;
const double eps = 1e-10;
int n;
double a[N][N];

int gauss()
{
    int c, r;//c 是枚举列 r 是枚举行
    for( c = 0, r = 0; c < n; ++ c)//1先枚举列
    {
        int t = r;
        for(int i = r; i < n; ++ i)//2枚举行找到该列中行的绝对值的最大值
          if(fabs(a[i][c]) - fabs(a[t][c]) > eps)
             t = i;
             
          if(fabs(a[t][c]) < eps) continue;//如果最大值为0就说明这一行的数值全部为0 就continue
          
          for(int i = c; i <= n; ++ i)//3将该行和最顶端的交换(除固定的外)
          swap(a[t][i],a[r][i]);
         
          for(int i = n; i >= c; -- i) a[r][i] /= a[r][c];
          //4将该行首项元素变为 1
          for(int i = r + 1; i < n; ++ i)//然后再枚举每一行将下面的非0行全部消成0
            if(fabs(a[i][c]) > eps)//如果该行的首不为0
              for(int j = n; j >= c; -- j)
                a[i][j]  -= a[r][j] * a[i][c];
        r ++;//统计非0行
    }
    if(r < n)//如果非0行的行数小于n
    {
        for(int i = r; i < n; ++ i)
         if(fabs(a[i][n]) > eps)//如果出现0!= 0的情况就是无解
            return 2;
            
            return 1;//否者就是有无穷多组解
    }
    for(int i = n; i >= 0; -- i)//如果是有唯一解
      for(int j = i + 1; j < n; ++ j)//按照我们手算的顺序从下到上去消元
        a[i][n] -= a[i][j] * a[j][n];
    return 0;
}

int main()
{
     cin >> n;
     for(int i = 0; i < n;  ++ i)
      for(int j = 0; j < n + 1; ++ j)
        cin >> a[i][j];
        
        int t = gauss();
        if(t == 0)
         for(int i = 0; i < n; ++ i) printf("%.2lf\n",a[i][n]);
        else if(t == 1) puts("Infinite group solutions");
        else puts("No solution");
        
        return 0;
}
           

继续阅读