数学小专题

 

正文

简单复习一下数学的几个模板。
总感觉数学还是要多复习,不然很容易忘记。

简单数论

O(n)判定质数

不用多讲。

Miller-Rabin算法

bool mr(ll x,ll b)
{
    if(qp(b,x-1,x)!=1)
        return false;
    ll k=x-1;   
    while((k&1)==0)
    {
        k>>=1;
        ll d=qp(b,k,x);//快速幂,计算b^k mod x
        if(d!=1 && d!=x-1)
            return false;
        if(d==x-1)
            return true;
    }
    return true;
}
bool mr(ll x)
{
    if(x==46856248255981)
        return false;
    if(x==2||x==3||x==7||x==61||x==24251)
    return true;
    return mr(x,2)&&mr(x,3)&&mr(x,7)&&mr(x,61)&&mr(x,24251);
}

对于一个数x,我们选一个质数b作为基底,依次计算bx1,bx12,bx122,。注意:

  • 在第一次计算时不考虑bx1x1的情况。
  • 遇到bx12kx1或者x12k为奇数的情况就直接判定正确。

注意5个素因子积:2,3,5,7,24251。那个很强的伪素数46856248255981我倒觉得真心没必要背。

约数O(N)筛法

又叫试除法。每次围绕N对称地找到两个因子p,q,把这两个因子加入约数列表中。注意当p2=N时要特判一下,不然会加多了。

倍数法

和线性筛的代码比较类似,适合一次筛多个质数的问题。

gcd

gcd(a,b)=gcd(b,amodb)

埃式筛

比较好理解。

线性筛

for(rg int i = 1; i <= N; ++ i)
{
    if(v[i] < 0)
    {
        v[i] = i;
        prime[++ pn] = i;
    }
    for(rg int j = 1; j <= pn; ++ j)
    {
        if(prime[j] * i > N || prime[j] > v[i])
            break;
        v[i * prime[j]] = prime[j];
    }
}

暂时还没有试过这个代码。从理论上来说应该是没错的。

进阶数论

冷知识

设一个数的分解式:

N=i=1kpiαi

我们可以把指数提取出来:

N=(α1,α2,)

它可以帮助我们从全新的角度观察数论。其实说白了就是把数论中的运算变成对数形式。

举个例子,求两个数的gcd:

gcd((αi),(βi))=(min{αi,βi})

这样我们就可以把方程gcd(x,a)=b变成若干个整数不等式了。

再举个例子,整除关系:

(αi)|(βi)αiβi

这样就也把整除关系变成若干个不等式了。从本质上来说,<Z,|>本身就是一个偏序集,因此才会有这样的特性。

可以结合这道题看一下:CSP-S 2009 Hankson的趣味题

重要函数的计算与预处理

欧拉函数ϕ(n)

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

这个是直接根据定义式计算的。
它的预处理版本:

void init()
{
    for(rg int i = 2; i <= n; ++ i)
        phi[i] = i;
    for(rg int i = 2; i <= n; ++ i)
        if(phi[i] == i)
            for(int j = i; j <= n; j += i)
                phi[j] = phi[j] / i * (i - 1);
}

当然,还有基于线性筛的版本:

void init()
{
    for(rg int i = 2; i <= N; ++ i)
    {
        if(v[i] == 0)
        {v[i] = i, prime[++ pn] = i, phi[i] = i - 1;}
        for(rg int j = 1; j <= pn; ++ j)
        {
            if(prime[j] * i > N || prime[j] > v[i])
                break;
            v[i * prime[j]] = prime[j];
            phi[i * prime[j]] = phi[i] * (i % prime[j] == 0 ? prime[j] : prime[j] - 1);
        }
    }
}

这里用到了ϕ(i)的两个性质,这里也不多赘述。

莫比乌斯函数μ(n)

for(rg int i = 1; i <= N; ++ i)
    mu[i] = 1, v[i] = 0;
for(rg int i = 2; i <= N; ++ i)
{
    if(v[i])
        continue;
    mu[i] = -1;
    for(rg int j = (i << 1); j <= N; j += i)
    {
        v[j] = 1;
        if((j / i) % i == 0)
            mu[j] = 0;
        else
            mu[j] *= -1;
    }
}

顺带一提,学长zsy讲过一个非常有趣的理解方式。这里不妨粘贴之。


我们都知道莫比乌斯反演:(简便起见,这里直接写成卷积形式)

f=Igg=μf

从之前的素因子式N=(αi)的角度理解,(αi)其实是一个高维向量,而g(αi)就是这个向量的点权,f(αi)就是g的前缀和。联想到二维前缀和ai=SiSi1,莫比乌斯反演事实上就是对于高维前缀和f(αi)的差分变换。

实际上,我之前补充过一个冷知识:广义的莫比乌斯反演。对于偏序关系,定义它的莫比乌斯函数:

μ(x,y)={1x=yxuyμ(x,u)xy0else

则有:

f(x)=0uxg(u)g(x)=0uxμ(u,x)f(u)

这实际上说明了莫比乌斯反演实质上就是一个局部有限偏序集上的高维前缀和差分变换。这也进一步说明了莫比乌斯函数在容斥原理和数论上的链接作用。


模算术

正如你所见,amodb这个函数在设计时出现了一些问题,使得两个字母之间间隔太大。因此,下文中均用C++写法a%b代表amodb。同理,为了方便书写,我们统一用a/b代表ab

欧拉定理

若正整数a,p互质,则:

aϕ(p)1(modp)

对于任意正整数b,还有:

abab%ϕ(p)(modp)

如果ap不一定互质,则应该使用下面这个公式:

ab{ab%ϕ(p)+ϕ(p)b>ϕ(p)abbϕ(p)(modp)

证明均略去。

裴蜀定理和exgcd

对于整数方程ax+by=c,存在整数解(x,y)的充要条件是gcd(a,b)|c。它的另一个版本是ax+by=gcd(a,b)一定有整数解。
它的证明是基于gcd(a,b)=gcd(b,a%b)的。在递归终点b=0时,显然有x=1,y=0。否则,
根据数学归纳法,假设下面这个方程成立:

bx+(a%b)y=gcd(b,a%b)

根据模数的另一种定义:

a%b=ab(a/b)

我们可以得到:

bx+(ab(a/b))y=gcd(b,a%b)=gcd(a,b)

整理得到:

ay+b(x(a/b)y)=gcd(a,b)

因此,我们令x=y,y=(x(b/a)y),就可以得到 ax+by=gcd(a,b)了。
在程序中可以这样简写:

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

一次同余方程FCE

一般形式为axb(modp)。可以转换成ax+py=b的形式,然后用exgcd解出(x,y)
假设gcd(a,p)b,那么方程无解。
注意一下得到特解x0后,方程的通解xx0(modpgcd(a,p))

逆元

a在模p意义下的逆元为a1,则:aa11(modp)。这是一个FCE,直接求解就可以了。

值得注意的是逆元的几种非常特殊的求法:

  1. 如果p是质数,那么a的逆元就是ap2(modp)

  2. 线性求逆元:i1(pp/i)×(p%i)1(modp)

  3. 线性求阶乘逆元:1(i+1)!×(i+1)1i!(modp)

稍微注意一下,在用阶乘求逆元时可能会出现p|i!的情况,此时阶乘算出来都是0

中国剩余定理

设方程组:

{xa1(modp1)xa2(modp2)xak(modpk)

p1,p2,,pk两两互质时,设P=i=1kpi,则方程一定有整数解:

xi=1kai(Ppi)(Ppi)1(modP)

其中(Ppi)1(Ppi)1(modpi)

假设p1,p2,,pk不满足两两互质,那么应该用扩展版本。

xk1是通过EXCRT计算出的前k1个方程的特解,而Pk1=i=1k1pi。为了让它满足第k个方程,有:

xk1+λPk1ak(modpk)

已知xk1,Pk1,ak,pk,这个方程就是一个关于λ的FCE,可以直接计算。当然,如果它无解,那整个方程组就都无解了。新的特解xk=xk1+λPk1通过加上或减去Pk就是前k个方程的特解。
另外有个小技巧:求乘积Pk可以改为求前k个模数的lcm。粘贴核心代码:

int main()
{
	n=qr(1);
	RP(i,1,n)
		m[i]=qr(1ll),a[i]=qr(1ll);//x = a[i] (mod m[i])
	
	ans=0,M=1;//M = prod
	RP(i,1,n)
	{
		ll p=FCE(M,a[i]-ans,m[i]);
		ans=ans+p*M;
		M*=m[i]/gcd(M,m[i]);
		ans=((ans%M)+M)%M;
	}
	printf("%lld",((ans%M)+M)%M);
	
	return 0;
}

高次同余方程(指数型)

求解形如axb(modp)类型的,关于x的方程。其中有a,p互质,x非负。
这里没有一般的解析方法。但是由于a,p互质,我们可以自由计算逆元,进行和a有关的乘除法。
x=ϵtη,其中t=p0η<t。事先声明,这样的记法其实有点“混用”的意味,但事实上相比普通的拉丁字母,它更有助于记忆理解。
方程变成了:

aϵtηb(modp)

即:

aϵtbaη(modp)

对于所有的baη,0η<t,我们将它插入一个Hash表中。枚举ϵ的所有取值0,1,,t,依次判断在Hash表中是否有对应的baη,然后拼凑得到答案。
这种方法叫做Baby step,Giant step。是一种类分块的优化枚举。

简单的线性代数

矩阵可以看成一个二维数表,可以用两个紧挨的下标ij来标识位置(i,j)处的元素。
矩阵加法就是对应的元素相加,即:Cij=Aij+Bij
矩阵乘法相对复杂。只有形如A(n×k)B(k×m)才能相乘得到一个形如C(n×m)的矩阵。标准的公式:

Cij=λ=1kAiλBλj

可以理解为:Cij=A的第i行行向量B的第j列列向量。正所谓“行列式”,它说明行在前,列在后。

矩阵乘法的应用

对于一个k阶的齐次线性递推关系:

xn=a1xn1+a2xn2++akxnk

由于是线性递推关系,我们可以把所有的状态xnxnk记录下来,封装到一个向量中去。即:

[a1a2ak1ak100001000010][xnxn1xn2xnk]=[xn+1xnxn1xnk+1]

此时这个k阶方矩阵:

T=[a1a2ak1ak100001000010]

又叫作转移矩阵。
如果知道这个递推关系,想要从x0,x1,,xk开始递推c次,这个过程就等价于对这个列向量左乘c次,即Tc。由于矩阵乘法满足结合律,这个过程可以用矩阵快速幂优化到O(k3logc)
如果上面不是齐次递推关系,而是带了一个常数p(如an=2an+3),该怎么办呢?可以直接把常数接在状态向量的最下面,转移矩阵赋值Tk+1,k+1=1,原封不动地转移到新状态去就可以了。
当然,直接代入矩阵乘法的公式,对着公式依次赋值,才是最稳妥的做法。

当然,矩阵乘法的应用不止于此。考虑这样一个问题:

给定一个带边权无向连通图,求过k个点(不算起点),从1走到n的最短路?

设邻接矩阵Eij表示点ij初始时的距离。如果两点无边相连,则赋值。由于+×都具有交换律和结合律,而min+同样具有交换律和结合律,我们可以把矩阵乘法+,×重定义成min,+。也就是说,每一次“乘法”:

Eij=minλ=1k{Eiλ+Eλj}

因此最终的答案就是E1,nk。可以用矩阵快速幂优化到O(n3logk)

高斯消元

为节省空间,这里直接接一个链接过去。

线性基

狭义的线性基仅指异或空间下的线性基。其将大小为N的二进制集合,转换为log值域大小级别的集合,大大降低了时间和空间复杂度。但是线性基往往会破坏顺序结构:你可以得到集合中若干数异或后的最大值,最小值,但却无法具体得知你到底异或了哪些数。

线性基可以看成数表{p0,p1,,pk},其中pi表示为1的最高位为i的数。它可能是由原来若干个数异或得到的。

当你想新插入一个数a到线性基中,你需要从高到低枚举每一个已有的基底pi。如果这个基底还没有数插入,你当然可以选择将这个位置加入a。但如果这一位上有数了,为了提高存储效率,你需要令a=axorpi从而消去这一位的值,然后继续检查低位的基底。

举个例子。这个线性基:

100010010010000110000001

就是一个不错的线性基。而:

100010010010011100000110000101

则非常糟糕。最高位为2,4的基底发生了重复,这两对基底应该分别异或起来,再插入这个线性基中。
粘贴一下核心代码:

void bits(ll x)
{
	DRP(i,62,0)
	{
		if(!(x >> (ll)i))
			continue;
		if(!p[i])
		{p[i] = x; break;}
		x ^= p[i];
	}
}
int main()
{
	scanf("%d", &n);
	RP(i,1,n)
		scanf("%lld", &a[i]), bits(a[i]);
	DRP(i,62,0)
		if((ans ^ p[i]) > ans)
			ans ^= p[i];
	cout<<ans;
	return 0;
}
posted @   LinearODE  阅读(354)  评论(0编辑  收藏  举报
编辑推荐:
· .NET 原生驾驭 AI 新基建实战系列:向量数据库的应用与畅想
· 从问题排查到源码分析:ActiveMQ消费端频繁日志刷屏的秘密
· 一次Java后端服务间歇性响应慢的问题排查记录
· dotnet 源代码生成器分析器入门
· ASP.NET Core 模型验证消息的本地化新姿势
阅读排行:
· ThreeJs-16智慧城市项目(重磅以及未来发展ai)
· .NET 原生驾驭 AI 新基建实战系列(一):向量数据库的应用与畅想
· Ai满嘴顺口溜,想考研?浪费我几个小时
· Browser-use 详细介绍&使用文档
· 软件产品开发中常见的10个问题及处理方法
点击右上角即可分享
微信分享提示