数论

最新推荐文章于 2023-10-06 01:41:12 发布

北弑鸣狮

最新推荐文章于 2023-10-06 01:41:12 发布

阅读量258

点赞数 1

分类专栏：数论文章标签：数论

本文链接：https://blog.csdn.net/qq_38890926/article/details/81513315

版权

数论专栏收录该内容

5 篇文章 0 订阅

订阅专栏

GCD问题和LCM问题

两个数的gcd

g=gcd(a,b) 则g<=a && g<=b

ll gcd(ll a,ll b)
{
    return b==0?a:gcd(b,a%b);
}

由该算法我们可以知道，给定a,b由a,b进行加减得出来的数都是gcd(a,b)的倍数，并且能得到任意倍数

两个数的lcm

l=gcd(a,b) 则l>=a && l>=b

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

gcd循环(辗转相除法的原理)

如果gcd(n,m)=g,则g=gcd(n,m)=gcd(n,m+i*n)=gcd(n+i*m,m)=gcd(1,g);其逆运算很关键，即我们可以用加法性质来简化运算较大结果的部分。

质因数大于都m的数

UVa11440

一个数的所有质因数都大于m,这种条件一定要转化到互质gcd=1才好用结论处理，因此，我们认为该数x和小于m的所有质数p互质，也即是 $gcd(x,p_{1}p_{2}..p_{n}=1)$ ,由于我们为了方便寻找这样的数，我们能想到互质定理(欧拉定理)，由于小于等于m的数都是由以上质数组成，故gcd(x,m!)=1也是成立的。也即是我们可以利用欧拉定理求得小于等于m!的x的个数。再利用gcd循环求得n!中x的个数即可。

扩展欧几里德

ax+by=gcd(x,y),已知a和b的情况下，求取gcd(a,b)，x，y的值

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

欧拉函数

定义： $m=p_{1}^{e_{1}}p_{2}^{e_{2}}....p_{n}^{e_{n}}$ 其中p1,p2...pn都是质数， $\varphi (m)=m*\prod (p_{i}-1)/p_{i}$ ，欧拉函数表示不超过m,且和m互质的数的个数。

//打出n的表 o(nlog(n))
int eular[maxn];
void get_eular()   //phi
{
    for(int i=1;i<maxn;i++)
        eular[i]=i;

    for(int i=2;i<maxn;i++)
        if(eular[i]==i)
            for(int j=i;j<maxn;j=j+i)
                eular[j]=eular[j]/i*(i-1);
}


//欧拉筛打表o(n)
int tot=0;
int phi[maxn],prime[maxn];
bool isPrime[maxn];
void getphi()
{
    for(int i=1;i<maxn;i++)isPrime[i]=true;
    phi[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(isPrime[i]==true)prime[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot;j++)
        {
            if(i*prime[j]>=maxn)break;
            isPrime[i*prime[j]]=false;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }else phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
}

//直接获得某个值的欧拉函数值 o(sqrt(n))
int eular(int temp)
{
    int result=temp;
    int sqrtt=sqrt(temp);
    for(int i=2;i<=sqrtt;i++)
    {
        if(temp%i==0)
        {
            result=result/i*(i-1);
            while(temp%i==0)      //剔除相同因数
                temp=temp/i;
        }
    }
    if(temp!=1)      //剩下一个质数
        result=result/temp*(temp-1);

    return result;
}

结论：

结论1：m是质数，则φ(m)=m-1,由文自定义容易获得。

结论2：对于和m互质的数x，x^φ(m) mod m =1;费马小定理的一般形式 x^(p-1) mod p=1 ( x<p,p为质数)

故：当 b很大时 $x^{b}=x^{b\bmod \varphi (m) } (\bmod m)$ ,可以将b缩小范围

结论3：欧拉函数是积性函数，则对于任意互质的整数a和b有性质f(ab)=f(a)f(b)

结论4：当n为奇数时，。从定义容易看出，因为n中没有2的因子。

原根

x是奇素数p的原根，当且仅当{x^i %p | 1<=i<=p-1} 等于集合{ 1，...，p-1}。

结论：由于p是质数，他有Phi(p-1)个原根

求取a关于范围b的互质个数

UVA - 10214

由于欧拉函数只能求取a关于范围1-a中和a互质的结果，因此，为了求取范围1-b中的结果，我们如下考虑gcd的循环：

将b分解成1-a a+1-2*a 2a+1-3a... na+1-b的部分，其中完整的部分有b/a个，最后一部分有a-b/a*a个结果，单独运算即可

//a在范围1-b上的互质的个数
int eular(int a,int b)
{
    int ans=b/a*eular[a];

    for(int i=1;i<=b%a;i++)
        if(gcd(a,i)==1)
            ans++;

    return ans;
}

Farey级数

F(n)(n>=2)表示有一系列不能约分的分数a/b（0<a<b<=n,且gcd(a,b)=1）组成，求Fn。

由于gcd(a,b)=1，因此我们可以知道就是对于每个b，其中多少个a与之互质，也即是b的欧拉函数。

故Farey[i]=Farey[i-1]+eular[i]; Farey[2]=eular[2]=1;

图像形式，也即是：看下半部分的，从原点射出一条线段，碰到的第一个点。加上(1,1)(0,1)(1,0)三个点，因为gcd(x,y)=1

利用莫比乌斯反演可更快得出结果:大致如下 $g(1)=\sum_{i=1}^{n}mu[i]*(n/i)*(n/i)+2$

质数问题

定义：除了1和本身以外不含有其他因数的数。也即是对于x来说i取(1,x-1)，gcd(i,x)=1。

质数判定

二次探测：

如果x是一个质数，0<m<x,则方程m^2%x=1，只有m=1或m=(x-1)的解

Millar_rabin+二次探测做法：

记x=(2^j)*m+1,则a^(x-1)%x ?=1,即变成了检测a^(2^j*m)%x是否等于1，也即是(a^2)^j * a^m %x ?= 1

//0(sqrt(n))算法
bool isprime(ll x)
{
    ll sqrtx=sqrt(x);
    for(int i=2;i<=sqrtx;i++)
        if(x%i==0)
        return false;
    return true;
}


// Miller-Rabin测试
ll multi(ll a,ll b,ll n)   //大数乘法 a*b mod n  目标乘法变加法，极端消耗时间
{
    a=a%n;b=b%n;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res=res+a;
            if(res>=n)res=res-n;
        }
        a=a<<1;
        if(a>=n)a=a-n;
        b=b>>1;
    }
    return res;
}
ll qpow(ll a,ll b,ll mod)   //超大数快速幂
{
    ll temp=1;
    while(b)
    {
        if(b&1)
            temp=multi(temp,a,mod);
        a=multi(a,a,mod);
        b=b>>1;
    }
    return temp;
}
inline ll random(ll n)
{
    return (ll)(rand()%(n-1)+1);
}
bool witness(ll a,ll n,ll x,int t)
{
    ll r=qpow(a,x,n);
    ll last=r;
    for(int i=1;i<=t;i++)
    {
        r=multi(r,r,n);
        if(r==1 && last!=1 && last!=n-1)return true;
        last=r;
    }
    if(r!=1)return true;
    return false;
}
bool Miller_rabin(ll n)
{
    if(n==2 || n==3 || n==5) return true;
    if((n<2) || (n&1==0) || (n%3==0) || (n%5==0))
        return false;

    ll x=n-1,t=0;
    while((x&1)==0)x=x>>1,t++;
    for(int k=0;k<10;k++)
        if(witness(random(n),n,x,t))   //如果rand的值过大
            return false;
    return true;
}

质数筛

// o(nloglog(n))
bool isPrime[maxn];
int prime[maxn];
int primeNum;

void getPrime()
{
    primeNum=0;
    for(int i=1;i<maxn;i++)
        isPrime[i]=true;

    for(int i=2;i<maxn;i++)
    {
        if(isPrime[i]==true)
        {
            prime[++primeNum]=i;
            for(ll j=(ll)i*(ll)i;j<maxn;j=j+i)
                isPrime[j]=false;
        }
    }
}

// 欧拉筛 O(n)    基于欧拉筛可以简化很多筛法
int prime[maxn],tot=0;
bool isPrime[maxn];
void getPrime()
{
    memset(isPrime,true,sizeof(isPrime));
    for(int i=2;i<maxn;i++)
    {
        if(isPrime[i])prime[++tot]=i;
        for(int j=1;j<=tot;j++)
        {
            if(i*prime[j]>=maxn)break;
            isPrime[i*prime[j]]=false;
            if(i%prime[j]==0) break;
        }
    }
}

质数筛-函数法

质数分解

正常分解：

        int sqrtx=sqrt(x);
        for(int i=2;i<=sqrtx;i++)
        {
            while(x%i==0)
                x=x/i;
        }
        if(x!=1)
        {
            /**/
        }

先质数筛再分解，节约时间.按照选出来的数进行分解，对于一个数x，只需要筛选出sqrt(x)就好：

bool isPrime[maxn];
int prime[maxn];
int primeNum;

void getPrime()
{
    primeNum=0;
    for(int i=1;i<maxn;i++)isPrime[i]=true;

    for(int i=2;i<maxn;i++)
    {
        if(isPrime[i]==true)
        {
            prime[++primeNum]=i;
            for(ll j=(ll)i*(ll)i;j<maxn;j=j+i)
                isPrime[j]=false;
        }
    }
}

int main()
{
    int sqrtx=sqrt(x);
    for(int i=1;prime[i]<=sqrtx;i++)
    {
        int p=prime[i];
        while(x%p==0)x=x/p;
    }
    if(x!=1)
    {
        /**/
    }
}

Pollard_Rho分解

对于较大一点的数x来说，我们需要求他的sqrt(x)来进行因式分解也是一件比较麻烦的事情，因此我们考虑基于类似Miller_Rabbin的概率算法来处理相应问题。

处理思路是，先用Miller_Rabbin算法判定x是否是一个质数，如果是则进行Pollard分解

//如果太慢 记得把qpow 换成直接快速幂

ll multi(ll a,ll b,ll n)   //大数乘法 a*b mod n  目标乘法变加法，极端消耗时间
{
    a=a%n;b=b%n;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res=res+a;
            if(res>=n)res=res-n;
        }
        a=a<<1;
        if(a>=n)a=a-n;
        b=b>>1;
    }
    return res;
}
ll qpow(ll a,ll b,ll mod)   //超大数快速幂
{
    ll temp=1;
    while(b)
    {
        if(b&1)
            temp=multi(temp,a,mod);
        a=multi(a,a,mod);
        b=b>>1;
    }
    return temp;
}
inline ll random(ll n)
{
    return (ll)(rand()%(n-1)+1);
}
bool witness(ll a,ll n,ll x,int t)
{
    ll r=qpow(a,x,n);
    ll last=r;
    for(int i=1;i<=t;i++)
    {
        r=multi(r,r,n);
        if(r==1 && last!=1 && last!=n-1)return true;
        last=r;
    }
    if(r!=1)return true;
    return false;
}
bool Miller_rabin(ll n)
{
    if(n==2 || n==3 || n==5) return true;
    if((n<2) || (n&1==0) || (n%3==0) || (n%5==0))
        return false;

    ll x=n-1,t=0;
    while((x&1)==0)x=x>>1,t++;
    for(int k=0;k<10;k++)
        if(witness(random(n),n,x,t))   //如果rand的值过大
            return false;
    return true;
}
/*      
ll gcd(ll a,ll b)   //分解出只有正数的快速方法
{
    if(a==0)return 1;
    if(a<0)return gcd(-a,b);
    while(b)
    {
        ll t=a%b;a=b;b=t;
    }
    return a;
}
*/
ll gcd(ll a, ll b)
{
    return b==0?a:gcd(b,a%b);
}
ll Pollard_rho(ll n,ll c)
{
    ll x=random(n);
    ll y=x,i=1,k=2;
    int num=0;
    while(true)
    {
        i++;
        x=(multi(x,x,n)+c)%n;
        ll d=gcd(y-x+n,n);  //由于用的是递归gcd，没判断正负，需要传入正值+n会影响较多速度，gcd中处理较好
        if(d!=1 && d!=n)return d;
        if(y==x)return n;
        if(i==k)
        {
            y=x;
            k=k<<1;
        }
    }
}
void find_factor(ll n)
{           //递归进行质因数分解N
    if(Miller_rabin(n))  //寻找到质因数
    {
        minnum=min(minnum,n);       //寻找到一个质因数n，对其进行合适作操作
        return;
    }
    ll p=n;
    while (p>=n) p=Pollard_rho(p,random(n));
    find_factor(p);
    find_factor(n/p);
}

常见结论

1.在10^20的范围类，两个质数之间的差值不会超过1000；

MOD问题

mod求和

UVa1363

给定n,k，求 $\sum_{i=1}^{n} k\bmod i$

我们将k mod i 进行分组，按照k/i=d的结果进行分组，我们即可看到规律。

当d=0时，k mod i=k.

当d=1时，i取(k/2+1,k)

当d=2时，i取(k/3+1,k/2)

当为d时，i取(k/(d+1)+1,k/d)

当d>=sqrt(k)时，i取(k/(d+1)+1,k/d),即该范围不多于1个值，可能无法取值。

故分三段取值d=0 d取(1,sqrt(k)) d>sqrt(k) 第一段直接得结论，第二段用等差数列得结果，第三段枚举

逆元问题

我们知道乘法和触发本质是一样的，因此在mod k 的前提下，我们需要求取（a * b / c ）mod k 的值，其实就是在求取((a mod k ) *(b mod k ) *( 1/c mod k ) )mod k 而 (1/c mod k)其实就是c在mod k下的逆元叫做inv(c)

//扩展欧几里德求取a在m下的逆元
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y=y-(a/b)*x;
    return d;
}
ll inv(ll a,ll m)
{
    ll x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}

如果k为质数，我们可以使用费马定理，利用快速幂处理

ll qpow(ll a,ll n,ll mod)
{
    ll temp=1;
    while(n)
    {
        if(n&1)temp=(temp*a)%mod;
        a=(a*a)%mod;
        n=n>>1;
    }
    return temp;
}
ll inv(ll a,ll m)
{
    return qpow(a,m-2,m);
}

逆元打表

ll inv[maxn];
void initInv()
{
    inv[0]=inv[1]=1;
    for(int i=2;i<mod+10;i++)inv[i]=(((-mod/i)*inv[mod%i])%mod+mod)%mod;
}

积性函数

定义：

积性函数：对于任意互质的整数a和b有性质f(ab)=f(a)f(b)的数论函数。

完全积性函数：对于任意整数a和b有性质f(ab)=f(a)f(b)的数论函数。

φ(n) －欧拉函数

μ(n) －莫比乌斯反演，关于非平方数的质因子数目

gcd(n,k) －最大公因子，当k固定的情况

d(n) －n的正因子数目

σ(n) －n的所有正因子之和

莫比乌斯反演

//单独求解
int mu(int x)
{
    int cnt=1;
    int sq=sqrt(x);
    for(int i=2;i<=x&&i<=sq;i++)
    {
        if(x%i!=0)continue;
        x=x/i;
        if(x%i==0)return 0;
        cnt=-cnt;
    }
    if(x!=1)cnt=-cnt;
    return cnt;
}


//O(n)
int mu[maxn];
bool vis[maxn];
int prime[maxn],tot;
void getMu()
{
    memset(vis,0,sizeof(vis));
    mu[1]=1;
    tot=0;
    for(int i=2;i<maxn;i++)
    {
        if(vis[i]==false)
        {
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot;j++)
        {
            if(i*prime[j]>=maxn)break;
            int tmp=i*prime[j];
            vis[tmp]=true;
            if(i%prime[j]==0)
            {
                mu[tmp]=0;
                break;
            }
            else mu[tmp]=-mu[i];
        }
    }
}

多项式问题

差分

给定多项式 $f(n)=a_{n}n^{b{n}}+a_{n-1}n^{b_{n-1}}+...+a_{0}$ , 当n可以取1.....n等值，使得计算出来的结果满足一定的条件，由于在多项式中的系数an可能比较大，使得结果不方便计算，或者由于n的取值比较大使得我们的for循环模拟比较大，我们不妨考虑差分的思想，想一想f(n)-f(n-1)满足什么样的结论。

结论

结论1： $f(n)=c_{n-1}f(n-1)+c_{n-2}f(n-2)+....c_{n-k-1}f(n-k-1) + (k=bn)$ 可从方程观查看未知数个数判断出结论.

例：uva 1069

给出形如 (P(n))/D 的式子，求出这个式子关于正n取值整数是否全是整数.

由于n的取值很大，我们可以考虑从差分入手，每个值都能被D整除，则他们的差分也能被D整除，则差分的差分也能被D整除，故我们能够得出一个结论：当n取0-bn时，对于每一个f(n)/D结果都是整数，我们就可以得到一个结论对于n去任意正整数f(n)/D为整数。

组合数问题

组合数公式： c(p,q)=p! / ((p-q)! * q!)

组合数除法

计算C（p， q） / C(s, r) uva 10375

利用组合数公式进行运算算则尽量少的运算次数，即q>p/2时，p=p-q

运算时乘大数除大数乘小数除小数防止精度越界。

int p,q,r,s;
int main()
{
    while(scanf("%d %d %d %d",&p,&q,&r,&s)!=EOF)
    {
        if(q>p/2)q=p-q;
        if(s>r/2)s=r-s;
        double ans=1.00;
        if(q>s)
        {
            for(int i=1;i<=s;i++)
                ans=ans*(p+1-i)/(r+1-i);
            for(int i=s+1;i<=q;i++)
                ans=ans*(p+1-i)/i;
        }
        else
        {
            for(int i=1;i<=q;i++)
                ans=ans*(p+1-i)/(r+1-i);
            for(int i=q+1;i<=s;i++)
                ans=ans/(r+1-i)*i;
        }
        printf("%.5f\n",ans);
    }
    return 0;
}

组合数取模

组合数取模也即是C(m,n)%p,其中p为质数，怎么求，有公式可知： $C(m,n)\%p=C(m/p,n/p)*C(m\%p,n\%p)\%p$ .

该定理又叫lucas定理。

这里写图片描述

故代码如下：(p较小，且p为质数)

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y=y-(a/b)*x;
    return d;
}
inline ll inv(ll a,ll m)
{
    ll x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}
ll C(ll n,ll m,ll p)
{
    if(m>n)return 0;
    ll up=1,down=1;
    for(int i=n-m+1;i<=n;i++)up=up*i%p;
    for(int i=1;i<=m;i++)down=down*i%p;
    return up*inv(down,p)%p;
}
ll Lucas(ll n,ll m,ll p)
{
    if(m==0)return 1;
    return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}

当p不为质数时，我们如何处理呢？

我们知道对p进行分解可以得到p=p1^a1 * p2 ^a2 *.....*pn^an, 且每一项之间两两互质；故如果求得X%(p1^a1)=b1,X%(p2^a2)=b2,.......,X%(pn^an)=bn,我们是不是可以用CRT(因为各项互质)将上式合并得到关于X%p=b的结果b；

假如ai都为1，即p为square-free-number，则求解方式直接采用质数的卢卡斯定理即可.

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y=y-(a/b)*x;
    return d;
}
inline ll inv(ll a,ll m)
{
    ll x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}
ll C(ll n,ll m,ll p)
{
    if(m>n)return 0;
    ll up=1,down=1;
    for(int i=n-m+1;i<=n;i++)up=up*i%p;
    for(int i=1;i<=m;i++)down=down*i%p;
    return up*inv(down,p)%p;
}
ll Lucas(ll n,ll m,ll p)
{
    if(m==0)return 1;
    return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}

ll multi(ll a,ll b,ll n)
{
    a=(a%n+n)%n;b=(b%n+n)%n;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res=res+a;
            if(res>=n)res=res-n;
        }
        a=a<<1;b=b>>1;
        if(a>=n)a=a-n;
    }
    return res;
}
ll crt(ll *a,ll *m,int num)
{
    ll lcm=1;
    for(int i=0;i<num;i++)lcm=lcm*a[i];
    ll ans=0;
    for(int i=0;i<num;i++)
       ans=ans+multi(multi(inv(lcm/a[i],a[i]),lcm/a[i],lcm),m[i],lcm);
    return (ans+lcm)%lcm;
}

int num;
ll a[20],b[20];
ll exLucas(ll n,ll m,ll p)
{
    for(int i=0;i<p;i++)
    {
        scanf("%lld",&a[i]);
        b[i]=Lucas(n,m,a[i]);
    }
    ll ans=crt(a,b,p);
    return ans;
}

如果p不是square-free-number，CRT部分只需各数互质，不会出现问题，但是lucas定理部分由于pi^ai不是质数,需要进行改动：

计算C(n,m)%p^t。我们知道，C(n,m)=n!/m!/(n-m)!，若我们可以计算出n!%p^t，我们就能计算出m!%p^t以及(n-m)!%p^t。我们不妨设x=n!%p^t,y=m!%p^t,z=(n-m)!%p^t,那么答案就是x*reverse(y,p^t)*reverse(z,p^t)(reverse(a,b)计算a对b的乘法逆元)。那么下面问题就转化成如何计算n!%p^t。比如p=3,t=2,n=19

n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19
=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6)=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6)
根据这个例子发现，求解n!可以分为3部分：第一部分是pipi的幂的部分，也就是3636即pi⌊npi⌋pi⌊npi⌋，可以直接求解；第二部分是一个新的阶乘，也就是6!即⌊npi⌋!⌊npi⌋!，可以递归下去求解；第三部分是除前两部分之外剩下的数
考虑第三部分如何求解
发现第三部分在模pikipiki意义下是以pikipiki为周期的，即(1∗2∗4∗5∗7∗8)≡(10∗11∗13∗14∗16∗17)(modpiki)(1∗2∗4∗5∗7∗8)≡(10∗11∗13∗14∗16∗17)(modpiki)，所以只求pikipiki长度的即可；但是还剩下一个孤立的19，可以发现剩下孤立的数长度不会超过pikipiki，只需要暴力求解即可
e.最后一个问题是对于求出的m!%pikim!%piki和(n−m)!%piki(n−m)!%piki有可能与pikipiki不互质，无法求逆元
所以要将m!%pikim!%piki和(n−m)!%piki(n−m)!%piki中质因子pipi先全部除去，求出逆元后再全部乘回去
计算n!中质因子p的个数x的公式为x=⌊np⌋+⌊np2⌋+⌊np3⌋+...x=⌊np⌋+⌊np2⌋+⌊np3⌋+...
递推式也可以写为f(n)=f(⌊np⌋)+⌊np⌋

ll qpow(ll a,ll b,ll m)
{
    ll tmp=1;
    while(b)
    {
        if(b&1)tmp=tmp*a%m;
        a=a*a%m;
        b=b>>1;
    }
    return tmp;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y=y-(a/b)*x;
    return d;
}
ll inv(ll a,ll m)
{
    ll x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}
ll fac(ll n, ll p, ll pk)
{
    if(!n)return 1;
    ll res=1;
    if(n/pk)
    {
        for(ll i=2;i<=pk;i++)if(i%p)res=res*i%pk;
        res=qpow(res,n/pk,pk);
    }
    for(ll i=2;i<=n%pk;i++)if(i%p)res=res*i%pk;
    return res*fac(n/p,p,pk)%pk;
}
ll Lucas(ll n,ll m,ll pi,ll pk)
{
    if(m>n)return 0;
    ll a=fac(n,pi,pk),b=fac(m,pi,pk),c=fac(n-m,pi,pk),k=0;
    for(ll i=n;i!=0;i=i/pi)k=k+i/pi;
    for(ll i=m;i!=0;i=i/pi)k=k-i/pi;
    for(ll i=n-m;i!=0;i=i/pi)k=k-i/pi;
    return a*inv(b,pk)%pk*inv(c,pk)%pk*qpow(pi,k,pk)%pk;
}
ll multi(ll a,ll b,ll n)
{
    a=(a%n+n)%n;b=(b%n+n)%n;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res=res+a;
            if(res>=n)res=res-n;
        }
        a=a<<1,b=b>>1;
        if(a>=n)a=a-n;
    }
    return res;
}
ll crt(ll *a,ll *m,int num)
{
    ll lcm=1;
    for(int i=0;i<num;i++)lcm=lcm*a[i];
    ll ans=0;
    for(int i=0;i<num;i++)
       ans=ans+multi(multi(inv(lcm/a[i],a[i]),lcm/a[i],lcm),m[i],lcm);
    return (ans+lcm)%lcm;
}

int tot;
ll a[20],b[20],px[maxn];
void divide(ll x)
{
    tot=0;
    ll sq=(ll)sqrt(x*1.0+2.0);
    for(int i=2;i<=sq&&x!=1;i++)
    {
        if(x%i==0)
        {
            px[tot]=i,a[tot]=1;
            while(x%i==0)a[tot]=a[tot]*i,x=x/i;
            tot++;
        }
    }
    if(x!=1)px[tot]=a[tot++]=x;
}
ll exLucas(ll n,ll m,ll p)
{
    divide(p);
    for(int i=0;i<tot;i++)b[i]=Lucas(n,m,px[i],a[i]);
    ll ans=crt(a,b,tot);
    return ans;
}

直角三角形

1.毕达哥拉斯定理

所有形如 $a^{2}+b^{2}=c^{2}$ 的数，有且仅有一下原型数：

令 $x=m^{2}-n^{2},y=2mn,z=m^{2}+n^{2}$ ,其中 m>n,gcd(m,n)=1,m (mod 2)!=n(mod2)

暴力搜索得：

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

int maxs=sqrt(l);
for(int i=1;i<=maxs;i++)
{
    int nn=i*i;
    for(int j=i+1;j<=maxs;j++)
    {
        int mm=j*j;
        if()
        {
            //超过，退出条件
        }
        if(i%2!=j%2 && gcd(i,j)==1)
        {
            int x=mm-nn;
            int y=2*i*j;
            int z=mm+nn;
            /*****/
            num=num++;
        }
    }
}

佩尔方程

佩尔方程：即形如 $x^{2}-dy^{2}=1$ 形式的方程；

佩尔方程的解法：即先计算出最小解： $x=sqrt(1+dy^{2}) y=1.2.3..$ ，其他所有的解可以用递推形式写出来

$\binom{x_{n}}{y_{n}}=\binom{x_{n-1}\quad dy_{n-1}}{y_{n-1}\quad x_{n-1}} \binom{x_{1}}{y_{1}}$ ,写成幂的形式，也即是

$\binom{x_{n}}{y_{n}}=\binom{x_{1}\quad dy_{1}}{y_{1}\quad x_{1}}^{n-1} \binom{x_{1}}{y_{1}}$

int n,k;
int main()
{
    while(cin>>n>>k)
    {
        int x=0,y;           //寻找特解
        for(y=1;y<5000;y++)
        {
            x=sqrt(1+n*y*y);
            if(x*x-n*y*y==1)
                break;
        }
        if(y==5000)
        {
            //no answer
        }
                //快速幂运算  计算第k个解
        Mat e;
        e.m[1][1]=x;e.m[1][2]=n*y;
        e.m[2][1]=y;e.m[2][2]=x;
        e.row=2,e.col=2;
        Mat res=mat_pow(e,k-1);
        int ans=(res.m[1][1]*x+res.m[1][2]*y)%mod;
        printf("%d\n",ans);
    }
    return 0;
}

中国剩余定理

给出几个线性同余方程：求最小的总的个数

一元线性同余方程组

其中m1,m2,mn相互互质；

容易知道，当求得一个结果xtemp的时候，令lcm=m1*m2*..*mn,则结果为x=xtemp+k*lcm。

其中最小的那个值为 x=(xtemp%lcm+lcm)%lcm;

引论：对于一个式子 x%a=m,我们很容易写出 x=(a*k+1)*m，本质上即找到mod a为1的数，使其乘余数m即可。

对于多个式子，本质上是这样多个式子的和，即xtemp=(a1*k1+1)*m1+...+(an*kn+1)*mn

不妨记为：xtemp=x1*m1+...+xn*mn;

其中xi表示mod ai为1的数，且为ki*π(aj)(j!=i),，即xi=ki*(lcm/ai).

故可以得到方程： ki*(lcm/ai)=1(mod ai)

故可以的得到式子：(lcm/ai)*ki+ai*y=1

即可以用exgcd求出ki的值从而得到xi的值，最终得到xtemp的结果。

代码过程：按照上述结果，先求出lcm，再根据每一个a[i]，求出(lcm/ai)*ki+ai*y=1中ki的值,最后求取ki*(lcm/ai)*mi的和

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y=y-(a/b)*x;
    return d;
}
inline ll inv(ll a,ll m)
{
    ll x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}
ll crt(ll *a,ll *m,int num)
{
    ll lcm=1;
    for(int i=0;i<num;i++)lcm=lcm*a[i];
    ll ans=0;
    for(int i=0;i<num;i++)ans=(ans+inv(lcm/a[i],a[i])*(lcm/a[i])*m[i])%lcm; //溢出取模
    return (ans+lcm)%lcm;
}

扩展中国剩余定理

扩展中国剩余定理即a1,a2,a3之间不一定互质的情况出现。

这时候需要对所给出的式子进行两个两个的组合，

考虑一下式子：

x ≡ m1 mod a1

x ≡m2 mod a2

所以有x= a1*k1+m1 = a2*k2+m2

即 a1*k1 = a2*k2+m2-m1

由于a1与a2不一定互质，故令d=gcd(a1,a2)。要求(m2-m1)%d=0;如果不为0，说明无法得到结果。

上式可以化成(a1/d)*k1=(m2-m1)/d+a2/d*k2;

即 (a1/d)*k1=c(mod (a2/d)); 因为c不一定是1，故求出a1/d关于a2/d的逆元，也即是 k1=c*inv(a1/d) (mod a2/d)

使用一次exgcd可以将k1temp的结果计算出来。

因为k1每变化a2/d的值，k2就会变化a1/d的值，则k1=k1temp+(a2/d)*y.

所以x=a1*(k1temp+(a2/d)*y )+m1; 展开可得 x=(a1*k1temp+m1) +(a1*a2/d)*y;

也即是x=(a1*k1temp+m1) (mod (a1*a2/d));

也就是说，找到了这么一个数，他满足12等式计算出来的结果，并且每(a1*a2/d)的差都满足该等式。

合并到最后得到x=m (mod a);也就是x=k*a+m；

代码过程：

提前定义一个大A和大M，表示x=M(mod A);

对于每一组a1，a2，求出其gcd值d,再代入exgcd(a1/d,a2/d,x,y),求得的x即为a1/d在a2/d下的逆元，乘以c也即是k1的值。

就可以得到x=(a1*k1temp+m1) (mod (a1*a2/d))；

最后x=(m%a+a)%a;

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y=y-(a/b)*x;
    return d;
}
ll excrt(ll *a,ll *m,int num)
{
    ll A=a[0],M=m[0];
    for(int i=1;i<num;i++)
    {
        ll x,y;
        ll d=exgcd(A,a[i],x,y);    //求两个式子a的gcd值，同时求出inv(a1/d)(mod a2/d)的结果

        ll c=m[i]-M;    //判断两个式子组成的结果是否正确
        if(c%d)return -1;
        c=c/d,A=A/d,a[i]=a[i]/d;
        ll k=((x*c)%a[i]+a[i])%a[i];//计算c*inv(a1/d)(mod a2/d)，求出k值并防止x过大或者过小
        M=M+A*d*k;      //获得新式子的M
        A=A*a[i]*d;     //获得新式子的A
        M=M%A;          //防止M过大过小
    }
    ll ans=M>0?M:M+A;    //当所求数大于0采用。
    //ll ans=(M%A+A)%A;  //M=A的情况哦，可以为0次时使用！！！！ 不是输出0
    return ans;
}

BSGS

BSGS算法即处理a^x≡b(mod p) (p为质数)的算法。

已知a,b,p的情况，求对应满足成立式子的最小x。

由费马定理可以知道，a^1 %p...a^(p-1)%p会遍历正好一遍0-p-1 因此x必定在0-p-1之间。

由此创建了bsgs算法，令x=m*i-t 其中m是个定值m=sqrt(p+0.5);即利用平方分割的思想进行处理，在平方分割的每个桶中，遍历寻找结果。

a^(m*i-t)≡b(mod p) 等价于 ((a^m)^i)≡b*a^t (mod p), 令c=a^m ，得c^i≡b*a^t (mod p)由此我们可以将每一个t对应的[1,m]对应的值存入map中，需要的t值越大越好，因为用的减法，故每次计算出的值，直接map就好，不用判断。

    ll m=ceil(sqrt(p+0.5)); //获取n的平方根
    hashmap.clear();   //map 放里面不用每次清空 浪费空间，节约时间
    ll a_t=b;  //利用连乘求值
    for(int t=1;t<=m;t++)
    {
        a_t=(a_t*a)%p;
        hashmap.add(a_t,t); //计算b*a^t(mod p) 反向映射，用结果映射t t为1的时候计算结果为b*a mmap[b*a]=t;
    }

让i从[0,m]遍历，去寻找对应的值是否存在，第一个找到的值对应的i，t就是我们需要的值。

typedef long long ll;

ll qpow(ll a,ll b,ll n) //快速幂，用于求a^m 以及求逆元
{
    ll temp=1;
    while(b)
    {
        if(b&1)temp=(temp*a)%n;
        a=(a*a)%n;
        b=b>>1;
    }
    return temp;
}

ll bsgs(ll a,ll b,ll p)
{
    ll m=ceil(sqrt(p+0.5)); //获取n的平方根
    hashmap.clear();   //map 放里面不用每次清空 浪费空间，节约时间
    ll a_t=b;  //利用连乘求值
    for(int t=1;t<=m;t++)
    {
        a_t=(a_t*a)%p;
        hashmap.add(a_t,t); //计算b*a^t(mod p) 反向映射
    }
    ll c=qpow(a,m,p);     //获取c
    ll c_i=1;
    for(int i=1;i<=m;i++)   //计算
    {
        c_i=(c_i*c)%p;
        ll num=hashmap.get(c_i);
        if(num!=-1)return m*i-num;
    }
    return -1;
}

int main()
{
    ll p,b,n;
    while(scanf("%lld %lld %lld",&p,&b,&n)!=EOF)
    {
        ll result=bsgs(b,n,p);

        if(result==-1)
            printf("no solution\n");
        else
            printf("%lld\n",result);
    }
    return 0;
}

扩展BSGS

扩展BSGS即在取mod的时候，p不一定是质数。

优先判断，当b==1的时候，即a^0==1 返回0；

即a^x≡b(mod p) 可以写成 a^x=kp+b, 由于a和p不互质，故存在d=gcd(a,p);使得a/d * a^(x-1)=k*(p/d)+b/d;成立，因为三个部分都是整数，故b/d一定是整数。如果不是，说明不存在这样的式子可以成立。

记tmp=a/d,p=p/d,b=b/d 故上式变成了 a/d * a^(x-1)=k*p+b; 即tmp*a^(x-1)=b(mod p),同理继续上述操作，直到d=gcd(a,p)=0；

有tmp=a/d1 * a/d2 *......*a/dk p=p/d1/d2../dk b=b/d1/d2.../dk a^(x-k)

故上式变成 tmp*a^(x-k)=b(mod p);

在这个计算过程中如果发现了某一个tmp==b 说明了a^(x-k)==1,故直接返回k即可。

    ll tmp=1,k=0;   //tmp被分割出来的部分的值 k进行了多少次gcd分割
    while(true)
    {
        ll d=gcd(a,p);
        if(d==1)break;
        if(b%d) return -1;
        k++;
        b=b/d;p=p/d;
        tmp=(a/d*tmp)%p;
        if(tmp==b)return k;
    }

处理完后，即可发现上式变成了a^(x-k)=b/tmp(mod p);

这里需要处理b/tmp(mod p),由于先需要处理b/tmp的最小公倍数，再将化简后的结果用exgce求取tmp关于p的逆元，因为p不一定是质数，故不能用费马定理求逆元。用exgcd来求逆元，

    b=(b*inv(tmp,p))%p;

用一次常规bsgs即可

    ll result=bsgs(a,b,p);
    if(result==-1)return -1;
    else return result+k;

总的代码

struct HashNode
{
    ll key;
    ll val;
    int next;
};
const int hashsize=31643,maxnode=100000;
struct Hash
{
    HashNode node[maxnode];
    int head[maxnode];
    int tot;
    void clear()
    {
        memset(head,-1,sizeof(head));
        tot=0;
    }
    void add(ll key,ll val)
    {
        ll x=key%hashsize;
        node[tot].key=key;
        node[tot].val=val;
        node[tot].next=head[x];
        head[x]=tot++;
    }
    ll get(ll key)
    {
        ll x=key%hashsize;
        for(int i=head[x];i!=-1;i=node[i].next)
            if(node[i].key==key) return node[i].val;
        return -1;
    }
};
Hash hashmap;

ll qpow(ll a,ll b,ll n) //快速幂，用于求a^m 以及求逆元
{
    ll temp=1;
    while(b)
    {
        if(b&1)temp=(temp*a)%n;
        a=(a*a)%n;
        b=b>>1;
    }
    return temp;
}
//扩展欧几里德求取a在m下的逆元
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    ll d=a;
    if(b==0)
    {
        x=1;y=0;
    }
    else
    {
        d=exgcd(b,a%b,y,x);
        y=y-(a/b)*x;
    }
    return d;
}
ll inv(ll a,ll m)
{
    ll x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}
ll bsgs(ll a,ll b,ll p)
{
    ll m=ceil(sqrt(p+0.5)); //获取n的平方根
    hashmap.clear();   //map 放里面不用每次清空 浪费空间，节约时间
    ll a_t=b;  //利用连乘求值
    for(int t=1;t<=m;t++)
    {
        a_t=(a_t*a)%p;
        hashmap.add(a_t,t); //计算b*a^t(mod p) 反向映射，用结果映射t t为1的时候计算结果为b*a mmap[b*a]=t;
    }
    ll c=qpow(a,m,p);     //获取c
    ll c_i=1;
    for(int i=1;i<=m;i++)   //计算
    {
        c_i=(c_i*c)%p;
        ll num=hashmap.get(c_i);
        if(num!=-1)return m*i-num;
    }
    return -1;
}
ll gcd(ll a ,ll b)
{
    return (b==0?a:gcd(b,a%b));
}
ll exbsgs(ll a,ll b,ll p)
{
    if(b==1)return 0; //特殊情况先考虑
    ll tmp=1,k=0;   //tmp被分割出来的部分的值 k进行了多少次gcd分割
    while(true)
    {
        ll d=gcd(a,p);
        if(d==1)break;
        if(b%d) return -1;
        k++;
        b=b/d;p=p/d;
        tmp=(a/d*tmp)%p;
        if(tmp==b)return k;
    }
    b=(b*inv(tmp,p))%p;
    ll result=bsgs(a,b,p);
    if(result==-1)return -1;
    else return result+k;
}
//  538 23454 14704
int main()
{
    ll k,p,n;
    while(scanf("%lld %lld %lld",&k,&p,&n)!=EOF)
    {
        if(n>p)
        {
            puts("Orz,I can’t find D!");
            continue;
        }
        ll result=exbsgs(k,n%p,p);
        if(result==-1)printf("Orz,I can’t find D!\n");
        else printf("%lld\n",result);
    }
    return 0;
}

伯努利数

给定k，计算i的k次方求和转换为多项式求和的系数

北弑鸣狮

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
数论

GCD问题和LCM问题两个数的gcdg=gcd(a,b) 则g<=a && g<=bll gcd(ll a,ll b){ return b==0?a:gcd(b,a%b);}由该算法我们可以知道，给定a,b由a,b进行加减得出来的数都是gcd(a,b)的倍数，并且能得到任意倍数两个数的lcml=gcd(a,b) 则l>=a &...
复制链接

扫一扫