欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

数论相关

程序员文章站 2022-07-09 10:19:14
...

数论相关

线性筛相关


线性筛相比传统的筛法要快的多,因为我们保证了每个数一定都是由它最小的一个质因子筛出来的,因此平摊的时间复杂度大概在O(n)左右,可以看作是线性的复杂度。所以叫做线性筛。
线性筛的应用可以说是十分的广泛了,积性函数可以通过它迅速求得。例如欧拉函数,莫比乌斯函数,约数和等。

莫比乌斯函数及莫比乌斯反演

首先莫比乌斯函数的定义为

μ(1)=1
一个数p,如果p能拆分为k个互不相同的质数的乘积,μ(p)=k。其他情况μ(p)=0。

莫比乌斯反演的问题主要适用于当形如

F(1)=f(1)
F(2)=f(1)+f(2)
F(3)=f(1)+f(3)
F(4)=f(1)+f(2)+f(4)

其中的f(x)是我们问题的答案,但是我们没有办法直接去进行求解。所以我们通过以下性质求F(x),从而求得f(x)。

f(1)=F(1)
f(2)=F(2)-F(1)
f(3)=F(3)-F(1)
f(4)=F(4)-F(2)

另外反演可以分为两种形式

F(n)=∑ f(d) {d/n} -> f(n)=∑ μ(d)F(n/d) {d/n}
F(n)=∑ f(d) {n/d} -> f(n)=∑ μ(n/d)F(d) {n/d}

第一种形式是f(n)的约数和,第二种则是倍数和。
也就是说当我们遇到一个问题的解为f(x),我们可以通过求x的约数和或者倍数和,然后通过反演求得。

此外我们求莫比乌斯函数值的方法,是基于莫比乌斯函数是积性函数的性质通过线性筛求得。

积性函数:对于两个互质的数x,y f(xy)=f(x)f(y)

线性筛的时间复杂度平摊下来是O(n)的,因为我们保证了每一个不是素数的数都是被它的最小质因数所筛去的。

//线性筛
for(int i=2;i<=MAXN;i++)
{
    if(!vis[i])
        prime.push_back(i);
    for(int j=0;j<prime.size();j++)
    {
        int t=prime[j]*i;
        if(t>MAXN)
            break;
        vis[t]=1;
        if(i%prime[j]==0)
            break;
    }
}

//基于积性函数的性质求莫比乌斯函数值
miu[1]=1;
for(int i=2;i<=MAXN;i++)
{
    if(!vis[i])
    {
        miu[i]=-1;
        prime.push_back(i);
    }
    for(int j=0;j<prime.size();j++)
    {
        int t=prime[j]*i;
        if(t>MAXN)
            break;
        vis[t]=1;
        if(i%prime[j]==0)//这里可以看作i和质数p两个没有互质,因而并没有满足积性函数的性质,所以特判。
        {
            miu[t]=0;
            break;
        }
        miu[t]=miu[i]*miu[prime[j]];
    }
}

欧拉函数

欧拉函数的值是比这个数小的且与它互质的数的个数。它是个积性函数,暂不证明…
对于任意质数p有两个性质
1.φ(p)=p-1
2.φ(p^k)=(p-1)*p^(k-1)
所以我们通过线性筛算Phi时就比较方便。
对于质数当然Phi值就是它本身减一,如果是个合数我们考虑两种情况。第一种,i%p!=0这说明i和p互质,因此φ(ip)=φ(i)(p1)。第二种情况就是两个不互质了,这个适合就不能直接算了。我们可以通过它的一些性质,转换一下。
φ(i*p)=φ(i/p^k)φ(p^(k+1))=φ(i/p^k)*p^k(p-1)=φ(i/p^k)*φ(p^k)*p=φ(i)*p

小结

我平时线性筛的应用不是很多,只是应用了筛素数。此次重新复习将一些积性函数的筛法复习,但应该仍然有遗漏,待第二次补充。

高斯消元


高斯消元可以理解为将方程组转换为矩阵,通过简单变化进行消元,达到解方程组的目的。
遇到这类题目,就我目前做过的题来看,大都可以看出如何列出方程组,一些题的方程组或许会特别恶心。但是只要认真仔细地处理一定可以做出来。坚信这一点!

高斯约当消元

高斯约到消元将系数矩阵消到只有对角线上有值的情况,这样我们进行求解的时候就会轻松很多。

//此为浮点数版本,可以不输出小数点后面的部分达到输出整数的目的
void gauss()
{
    int i,j,r,c;
    int maxr;
    for(r=0,c=0;r<n&&c<m;r++,c++)
    {
        maxr=r;
        for(i=r+1;i<n;i++)
            if(fabs(A[i][c])>fabs(A[maxr][c]))
                maxr=i;
        if(fabs(A[maxr][c])<EPS)
        {
            r--;
            continue;
        }
        if(maxr!=r)
            for(j=c;j<=m;j++)
                swap(A[r][j],A[maxr][j]);
        for(i=0;i<n;i++)
        {
            if(i!=r&&fabs(A[i][c])>EPS)
            {
                for(j=m;j>=c;j--)
                    A[i][j]-=A[r][j]/A[r][c]*A[i][c];
            }
        }
    }
    ran=r;//ran存储为*变元的数量
}

//判断是否有解,以及标记*变元部分以及求得答案矩阵
bool check()
{
    int i, j;
    int cnt, pos;
    for(i = ran; i<n;i++)
        if(fabs(A[i][m])> EPS)
            return false;
    for(i = 0; i < m; i++)
        Free[i] = true;
    for(i=ran-1; i >= 0; i--)
    {
        cnt = 0;
        for(j = 0; j < m; j++)
            if(fabs(A[i][j]) > EPS && Free[j])
                cnt++, pos = j;
        if(cnt == 1)
        {
            Free[pos] = false;
            X[pos] = A[i][m] / A[i][pos];
        }
    }
    return true;
}

行列式

行列式可以说是高斯消元的一种应用。我们求行列式的值就需要应用高斯消元,利用矩阵的初等变换,将矩阵消成只有上三角有值的情况,然后将对角线上的数相乘,便得到行列式的值。应为在矩阵的变换中,我们没有取模运算,我们利用一下两条性质可以得到“取模”这个运算。
1.任意两行交换,行列式的值变号
2.任意一行加减另一行的倍数,行列式的值不变。
这样我们便可以利用整除求得余数了。

long long guass(int n,int m)//返回的是行列式的值
{
    long long ans=1LL;
    int r=1,c=1;
    for(;c<=n&&r<=m;r++,c++)
    {
        int k=r;
        while(k<=m&&!a[k][c])
            k++;
        if(k>m)
            return 0;
        if(r!=k)//两行交换,行列式的值不变
            swa(r,k,n),ans=-ans;
        for(int i=k+1;i<=m;i++)
        {
            while(a[i][c])//此处利用了辗转相除法,迅速地将a[c][r]这个位置化为0.
            {
                long long t=a[i][c]/a[r][c];
                for(int j=c;j<=n;j++)
                    a[i][j]=(a[i][j]-a[r][j]*t)%MOD;
                if(!a[i][c])
                    break;
                swa(i,r,n);
                ans=-ans;
            }
        }
    }
    if(r!=m+1)
        return 0;
    for(int i=1;i<=n;i++)//最终结果就是对角线相乘。
    {
        ans*=a[i][i];
        ans%=MOD;
    }
    return (ans+MOD)%MOD;//注意在此处ans可能为负的,所以取模的时候将其转换为正数。
}

异或方程组

异或方程组主要应用在一些状态只有两个,且每次只改变一个或多个位置的状态。我们根据题目的不同进行不同的处理。例如开关灯问题,很明显,可以列出异或方程,然后通过高斯约当消元,可以解出答案。我们只需要在约当消元的时候,将模运算改为异或运算就可以了。

小结

高斯消元类的问题大概都比较明显了,根据题目列出方程然后求解。相对来说比较直接(暴力求解…)。主要是在矩阵的变换中可能会出现一些纰漏,注意细节。

辗转相除法(欧几里得算法)


辗转相除法

在一些题中,有很大一部分是有关最大公约数的,所以我们需要快速地求到最大公约数。这里就要用到辗转相除法。

gcd(a,b)表示a和b的最大公约数
则有gcd(a,b)=gcd(b,a%b)

证明我不知道,所以只会用结论。但似乎好像在数论课的课件上讲过,有时间可以去翻翻。

//极其简单的代码
int gcd(int a,int b)
{
    if(b==0)
        return a;
    return gcd(a,b);
}

拓展欧几里得算法

普通的欧几里得算法我们用的并不多,更多是用拓展欧几里得算法解不定方程或模线性方程。
对于不定方程ax+by=gcd(a,b)=1转移成用拓展欧几里得求解的过程(即a,b**互质**的情况):
数论相关
所以我们只需要考虑一种特殊情况即x=1,y=0通过这个解向上回代,可以求得一组解(x0,y0)。

x=x0+(b/d)k y=y0-(a/d)k k是整数

通过以上的证明我们就可以写出拓展欧几里得的代码

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

那么问题来了,如果a和b**不互质**的情况那么大致相同。
我们将原来的方程转化为两个系数是互质的情况就可以处理了。
大概就是这样的
数论相关

最后一个问题,模线性方程组怎么解。
方程
数论相关
数论相关

小结

拓展欧几里得的应用比较广泛,应该熟练掌握,尤其是模线性方程,应该特别注意推导过程以及在代码中的实现。

大素数的判定与大整数的质因子分解


米勒拉宾素数测试

定理:对一个大于2的素数p,方程x^2≡1(mod p),的解只能为x=1或x=-1
通过这一个定理,我们可以知道如果方程x^2≡1(mod p)的解不为1或-1,那么p就不是素数了。
因此可以由此推得算法。若一个数p是质数,则p-1是个偶数,可以分解为p-1=2^k*d,都费马小定理可以得到a^(p-1)≡1(mod p),将之前的式子带入,得到(a^(2^(s-1)*d))^2≡1(mod p)。
上面的式子继续推导将会得到,a^d≡1(mod p),a^(2^r*d)≡-1(mod p)上面的两个式子至少有一个是成立的话,那么p就有可能是个素数。但是你可能会注意到这个a的选取对于这个式子的成立似乎是有一定影响的。也就是说有可能p是个合数是,这个式子依旧成立,则称p骗过了基a的测试。
因此为了避免这种情况发生,我们要多次选择a,进行测试,这样p骗过测试的几率会小得多。在这里你可以选择使用随机数进行,也可以选择用几个比较小的素数进行。根据实战,大概用二十以内的素数进行测试可以满足竞赛所需的数据范围了。

bool MillerTest(long long d,long long n,int k)
{
    long long a=prime[k];
    long long x=PowMod(a,d,n);
    if(x==1||x==n-1)
        return 1;
    while(d!=n-1)
    {
        x=(x*x)%n;
        d*=2;
        if(x==n-1)
            return 1;
    }
    return 0;
}

bool isP(long long n)
{
    if(n<=1||n==4)
        return 0;
    if(n<=3)
        return 1;
    long long d=n-1;
    while(d%2==0)
        d/=2;
    for(int k=0;k<8&&prime[k]<n;k++)
        if(!MillerTest(d,n,k))
            return 0;
    return 1;
}

Pollard‘s rho素因数分解算法

对于一些较大的数,我们需要对其进行质因数分解时,如果用朴素的算法枚举1~sqrt(n)的话是远远不能达到我们对时间的需求的,所以我们需要一个快速地分解办法。
此算法注意是通过生日悖论来提高效率的(过程略,详参课件),随机选取两个数,它们两个的差是否是n的因数。

long long PR(long long n)
{
    if(n==1||IsP(n))//判断是否是质数
        return n;
    if(n%2LL==0)
        return 2LL;
    long long x=(rand()-1)%n+1LL;
    long long y=x;
    long long c=(rand()-2)%n+2LL;
    long long d=1LL;
    while(d==1LL)
    {
        x=(AddMod(x,x,n)+c+n)%n;
        y=(AddMod(y,y,n)+c+n)%n;
        y=(AddMod(y,y,n)+c+n)%n;
        d=gcd(myabs(x-y),n);//如果d为零,说明x和y相等,说明已经产生了环,或正好求得了一个因数。
    }
    long long r=n/d;
    if(r==1)
        return PR(d);
    return min(PR(d),PR(r));//求最小的质因数
}

小结

这两个算法在一些比较特殊的题目中可能会遇到,平时遇到的可能性比较小,但仍旧应该掌握!

相关标签: c++ 数论