2020寒假专题训练：数论

Day1 20.02.11

一、素数判定：Miller-Rabbin

①费马小定理：$a^{p-1} equiv 1 pmod p$（$p$为素数）.

因此可以选取若干个 $a$，对 $p$ 进行判定.  满足条件的 $p$ 为素数的概率在 $frac{3}{4}$ 左右.

②二次探测定理：$x^2 equiv 1 pmod p$ 当 $p$ 为大于 $2$ 的素数时仅有两个解 $x_1=1,x_2=p-1$.

因此考虑选取 $9$ 至 $12$ 个素数，先使用费马小定理判定，再对 $a^{frac{p-1}{2^k}}$ 使用二次探测判定.

代码如下：

inline long long mul ( long long x,long long y,long long mod ) { return (long long)((__int128)x*y%mod); }
inline long long power ( long long x,long long y,long long mod )
{
    long long z=1;
    for ( ;y;y>>=1,x=mul(x,x,mod) ) if ( y&1 ) z=mul(z,x,mod);
    return z;
}
const long long Prime[13]={1,2,3,5,7,11,13,17,19,61,2333,4567,24251};
long long Pow[100],ans;
inline bool Miller_Rabbin ( long long p )
{
    if ( p==1 ) return false;
    for ( int i=1;i<=12;i++ )
    {
        if ( p==Prime[i] ) return true;
        if ( !(p%Prime[i]) or power(Prime[i]%p,p-1,p)!=1 ) return false;
        long long x=p-1;
        while ( !(x&1) )
        {
            x>>=1;
            long long res=power(Prime[i]%p,x,p);
            if ( res!=1 and res!=p-1 ) return false;
            if ( res==p-1 ) break;
        }
    }
    return true;
}

二、素数筛法：线性筛（欧拉筛）

①埃氏筛：枚举每个素数 $p$，将 $p$ 的倍数标记为非素数. 可以证明效率为 $O(n log log n)$.

②欧拉筛：考虑对埃氏筛进行优化. 不难发现，埃氏筛的问题在于对多因子的合数进行了多次标记.

考虑只对每个数标记一次. 令每个数只被其最小的素因子筛到. 考虑转化为枚举除去最小素因子的剩下的因子.

对于每个数枚举小于其的素数. 将这两个数的积标记为非素数. 当枚举的素数是当前数的因子是退出循环即可.

效率 $O(n)$. 代码如下：

bool flag[100000000];
std::vector<int> pr;
inline void solve ( int n )
{
    for ( int i=2;i<=n;i++ )
    {
        if ( !flag[i] ) pr.push_back(i);
        for ( int j:pr )
        {
            if ( i*j>n ) break;
            flag[i*j]=true;
            if ( !(i%j) ) break;
        }
    }
}

线性筛还可以用来求积性函数. 后文 Day4 部分有提及.

三、质因数分解：Pollard-Rho

前置知识：素数判定 Miller-Rabbin.

引理：生日悖论：$23$ 人中，存在两个人在同一天生日的概率超过 $frac{1}{2}$.

考虑找到一个数 $a$ 使得 $1 < gcd(a,n) < n$. 则有：$gcd(a,n) mid n, frac{n}{gcd(a,n)} mid n$.

考虑用随机的方式找到 $a$. 构造函数 $f(x)=(x^2+c) mod n$. 其中 $c$ 为一个固定的随机数.

先任意取 $x_1$，然后令 $x_i=f(x_{i-1})(i geq 2)$ 即可.

观察生成的随机数列，发现分布较为均匀.

但注意到一个问题：生成的数列在一定长度后会陷入循环，形成 $ ho$ 形数列.

因此每次判断 $d=gcd(lvert x_i-x_1 vert,n)$ 是否符合要求. 当 $x_i=x_1$ 时退出循环. 由引理知这样一个环的期望长度为 $O(sqrt n)$.

考虑对求解的过程进行优化. 由于计算 $gcd$ 的时间较慢，考虑每隔一段时间统一计算. 采用倍增的思想，统计 $res=(prod lvert x_i - x_1  vert) mod n$.

当 $res=0$ 时出现环. 每隔 $2^k$ 次计算一次 $gcd$，重新初始化即可. 代码如下：

inline long long work ( long long x )
{
    long long s=0,t=0,c=1LL*rand()*rand()*rand()%(x-1)+1;
    for ( int step=1;;step<<=1,s=t )
    {
        long long res=1;
        for ( int i=1;i<=step;i++ )
        {
            t=(mul(t,t,x)+c)%x;
            res=mul(res,std::abs(t-s),x);
            if ( !(i%127) )
            {
                long long d=std::__gcd(res,x);
                if ( d!=1 ) return d;
            }
        }
        long long d=std::__gcd(res,x);
        if ( d!=1 ) return d;
    }
    return -1;
}
inline void Pollard_Rho ( long long x )
{
    if ( x==1 or x<=ans ) return;
    if ( Miller_Rabbin(x) ) { ans=std::max(ans,x);return; }
    long long y=x;
    while ( y>=x ) y=work(x);
    while ( !(x%y) ) x/=y;
    Pollard_Rho(x);Pollard_Rho(y);
}

四、BSGS/exBSGS

1. BSGS

求解方程：$y^x equiv z pmod p$（$p$ 为素数）.

令 $x=am-b(0 leq b leq m-1)$，则有：$y^{am} equiv z imes y^b pmod p$.

等式右边有 $m$ 种取值，用 hash_table 或 unordered_map 记录.

等式左边枚举 $a$，查询右边是否出现过该值即可.

当 $m=sqrt p$ 时效率最优，效率为 $O(sqrt p)$.

实际写代码时由于数据原因，常将 $m$ 与一个常量（例如 $1000$，视数据而定）取 $min$.

代码如下：

std::unordered_map<int,int> map;
inline int BSGS ( int y,int z,int p )
{
    int m=(int)sqrt(p)+1;
    if ( z%p==1 ) return 0;
    if ( !(y%p) ) return -1;
    map.clear();
    for ( int i=0,t=z;i<m;i++,t=1LL*t*y%p ) map[t]=i;
    for ( int i=1,v=power(y,m,p),t=v;i<=m+1;i++,t=1LL*t*v%p ) if ( map.count(t) ) return i*m-map[t];
    return -1;
}

2. exBSGS

求解方程：$y^x equiv z pmod p$（$p$ 为任意正整数）.

将等式展开有：$y imes y^{x-1} + k imes p = z(k in mathbb{Z})$.

由裴属定理（Day2 内容）知，当且仅当 $d=gcd(y,k) mid z$ 时有解.

若有解则有：$frac{y}{d}  imes y^{x-1} + k imes frac{p}{d} = z$.

不断递归求解，最终得到 $gcd(y,p)=1$ 时即可用 BSGS 求解.

设递归 $c$ 次，$g=prodlimits_{i=1}^{c} d_i$.

则有：$y^{x-c} imes frac{y^c}{g} equiv frac{z}{g} pmod {frac{p}{g}}$.

代码如下：

std::unordered_map<int,int> map;
inline int exBSGS ( int y,int p,int z )
{
    if ( z==1 ) return 0;
    int c=0,g=1;
    for ( int d=std::__gcd(y,p);d!=1;d=std::__gcd(y,p) )
    {
        if ( z%d ) return -1;
        c++;z/=d;p/=d;g=1LL*g*y/d%p;
        if ( z==g ) return c;
    }
    int m=sqrt(p)+1;map.clear();
    for ( int i=0,t=z;i<m;i++,t=1LL*t*y%p ) map[t]=i;
    for ( int i=1,v=power(y,m,p),t=1LL*g*v%p;i<=m+1;i++,t=1LL*t*v%p ) if ( map.count(t) ) return i*m-map[t]+c;
    return -1;
}

---

---

Day2 20.02.12

一、裴蜀定理

裴蜀定理：方程 $sumlimits_{i=1}^{n} a_i imes x_i = b$ 有整数解的充要条件为 $gcd(a_1,a_2,dots,a_n mid b$.

证明：

①必要性：

令 $d=gcd(a_1,a_2,dots,a_n)$，则 $forall i in [1,n], d mid a_i$.

$ herefore d mid sumlimits_{i=1}^{n} a_i imes x_i$，即 $gcd(a_1,a_2,dots,a_n) mid b$.

②充分性：

考虑数学归纳法.

Ⅰ. $n=2$ 时成立（详见 exgcd 部分）.

Ⅱ. 假设 $n=k$ 时成立，考虑 $n=k+1$ 时，

设 $S=sumlimits_{i=1}{k} a_i,d=gcd(a_1,a_2,dots,a_k)$.

则方程 $S imes X + a_{k+1} imes y = gcd(d,a_{k+1})$ 有整数解.

即存在 $(X,Y)$ 使得 $sumlimits_{i=1}^{k+1} a_i imes x_i = gcd(a_1,a_2,dots,a_{k+1})$ 一组解为 $x_i=X(i in [1,k]),x_{k+1}=Y$.

$ herefore n=k+1$ 时成立.

综上，证毕.

二、欧拉定理

欧拉定理：$a^{varphi(m)} equiv 1 pmod m(a perp m)$.

证明：考虑所有与 $m$ 互素的数 $x_i(i in [1,varphi(m)]$，令 $p_i=ax_i$，不难得到 ${p_1 mod m,p_2 mod m,dots,p_{varphi(m)} mod m}$ 与 ${x_1,x_2,dots,x_{varphi(m)}$ 一一对应.

累乘得：$prodlimits_{i=1}^{varphi(m)} p_i =prodlimits_{i=1}^{varphi(m)} x_i $，则 $a^{varphi(m)} equiv 1 pmod m$.

扩展欧拉定理：当 $b geq varphi(m)$ 时，$a^b equiv a^{b mod varphi(m) + varphi(m) pmod m$.

证明考虑 $m$ 的每个素因子即可，此处不予给出.

三、exgcd

求解方程：$ax+by=c$ 或 $ax equiv c pmod b$.

当 $gcd(a,b) mid c$ 时原方程有解. 考虑先构造出方程 $ax+by=gcd(a,b)$ 的解，再令 $x'=frac{c}{gcd(a,b)} imes x$ 即可，$y$ 同理.

设：$egin{cases}ax_1+by_1=gcd(a,b)\bx_2+(amod b)y_2=gcd(b,amod b)end{cases}$.

又 $gcd(a,b)=gcd(b,a mod b)$ 得，$ax_1+by_1=bx_2+(a mod b)y_2$.

将 $amod b$ 展开并整理得：$ax_1+by_1=ay_2+b(x_2-lfloor frac{a}{b} floor y_2)$.

不断递归处理即可. 边界条件为当 $b=0$ 时，$x=1,y=0$. 代码如下：

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

四、CRT/exCRT

求解方程组：$egin{cases}xequiv a_1 pmod {n_1}\xequiv a_2 pmod {n_2}\ dots \xequiv a_k pmod {n_k}end{cases}$. 其中 $forall i eq j, n_i perp n_j$.

CRT（中国剩余定理）：

① 令 $n=prodlimits_{i=1}^{k} n_i$.

② 对于第 $i$ 个方程：

Ⅰ. 计算 $m_i=frac{n}{n_i}$.

Ⅱ. 计算 $m_i$ 在 $mod n_i$ 意义下的逆元 $m_i^{-1}$. 即方程 $am_i+bn_i=1$ 的最小正整数解 $a$，可使用 exgcd 求解.

Ⅲ. 令 $c_i=m_im_i^{-1}$（不对 $n_i$ 取模）.

③ 则唯一解为 $x equiv sumlimits_{i=1}^{k} a_ic_i pmod n$.

证明直接代入即可. 代码如下：

inline long long CRT ( int n,long long *a,long long *m )
{
    long long M=1,res=0;
    for ( int i=1;i<=n;i++ ) M*=m[i];
    for ( int i=1;i<=n;i++ )
    {
        long long w=M/m[i],x,y;exgcd(w,m[i],x,y);x=(x%m[i]+m[i])%m[i];
        res=(res+mul(mul(a[i],w,M),x,M)+M)%M;
    }
    return res;
}

求解方程组：$egin{cases}xequiv a_1 pmod {n_1}\xequiv a_2 pmod {n_2}\ dots \xequiv a_k pmod {n_k}end{cases}$. 其中 $n_i$ 为任意正整数.

exCRT：

考虑 $egin{cases} xequiv a_1 pmod {n_1}\xequiv a_2 pmod {n_2} end{cases}$. 有：$egin{cases} x=pn_1+a_1 \ x=qn_2+a_2 end{cases}$.

则 $n_1p-n_2q=a_2-a_1$. 用 exgcd 求出一组 $p,q$，则有：$x equiv n_1p+a_1 pmod {operatorname{lcm} (n_1,n_2)}$. 两两合并即可. 

代码如下：

inline long long exCRT ( void )
{
    long long m=a[1],ans=r[1],x=0,y=0;
    for ( int i=1;i<=n;i++ )
    {
        long long B=((r[i]-ans)%a[i]+a[i])%a[i],gcd=exgcd(m,a[i],x,y);
        x=mul(x,B/gcd,a[i]);ans+=m*x;m*=a[i]/gcd;ans=(ans+m)%m;
    }
    return (ans%m+m)%m;
}

五、Lucas/exLucas

Lucas 定理：对于素数 $p$ 有：$C_n^m mod p = C_{lfloor n/p floor}^{lfloor m/p floor} imes C_{n mod p}^{m mod p} mod p$.

上式的含义即为 $p$ 进制下 $n,m$ 的拆位. 证明可用二项式定理求解，此处略去. 效率 $O(p log_p n)$. 一般 $p$ 为 $10^7$ 及以下级别的素数，如 $10000019$.

 代码如下：

inline long long Lucas ( int x,int y )
{
    if ( x<y ) return 0;
    else if ( x<p ) return fct[x]*inv[y]*inv[x-y]%p;
    else return Lucas(x/p,y/p)*Lucas(x%p,y%p)%p;
}

exLucas：求 $C_n^m mod p$，其中 $p$ 为任意正整数.

考虑对 $p$ 进行分解. 转化为每个在模素因数的幂意义下的组合数取余. 分别解出后用 CRT 合并.

现在考虑 $C_n^m mod p^k$. 将组合数展开为阶乘形式，分别计算.

对于 $n!$，考虑提取所有 $p$ 的倍数的因子，以 $n=19,p=3$ 为例有：

$19!=1 imes 2 imes dots imes 19=1 imes 2 imes 4 imes 5 imes dots imes 16 imes 17 imes 19 imes (3 imes 6 imes dots imes 18)$.

前面一部分整段的直接算，最后一段每个单独计算即可. 后面部分提取出一个 $p$ 后转化为更小的阶乘，递归求解.

代码如下：

inline long long FCT ( long long n,long long x,long long p )
{
    if ( !n ) return 1;
    long long s=1;
    for ( long long i=1;i<=p;i++ ) if ( i%x ) s=mul(s,i,p);
    s=power(s,n/p,p);
    for ( long long i=n/p*p+1;i<=n;i++ ) if ( i%x ) s=mul(s,i,p);
    return mul(s,FCT(n/x,x,p),p);
}
inline long long INV ( long long a,long long b ) { long long x,y;exgcd(a,b,x,y);return (x%b+b)%b; }
inline long long multiLucas ( long long n,long long m,long long x,long long p )
{
    int cnt=0;
    for ( long long i=n;i;i/=x ) cnt+=i/x;
    for ( long long i=m;i;i/=x ) cnt-=i/x;
    for ( long long i=n-m;i;i/=x ) cnt-=i/x;
    return mul(mul(power(x,cnt,p),FCT(n,x,p),p),mul(INV(FCT(m,x,p),p),INV(FCT(n-m,x,p),p),p),p);
}
inline long long exLucas ( long long n,long long m,long long p )
{
    int cnt=0;long long mod[20]={0},a[20]={0};
    for ( long long i=2;i*i<=p;i++ ) if ( !(p%i) )
    {
        mod[++cnt]=1;
        while ( !(p%i) ) mod[cnt]*=i,p/=i;
        a[cnt]=multiLucas(n,m,i,mod[cnt]);
    }
    if ( p>1 ) mod[++cnt]=p,a[cnt]=multiLucas(n,m,p,p);
    return CRT(cnt,a,mod);
}

---

---

Day3 20.02.13

一、常系数一阶线性递推

递推式：$a_i=ka_{i-1}+b(k eq 0)$.

通项公式：$a_n=egin{cases}a_0+nb&k=1\a_0k^n+frac{b}{k-1}k^n-frac{b}{k-1}&k eq 1end{cases}$.

证明：构造 $a_i+x_0=k(a_{i-1}+x_0)$ 后运用等比数列相关知识即可.

【例题】DTOJ 1378 随机数生成器

【题目链接】http://59.61.75.5:8018/problem/1378

【题解】

直接分类讨论即可. 解方程时用 BSGS 求解.

【代码】

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
inline int power ( int x,int y,int mod )
{
    int z=1;
    for ( ;y;y>>=1,x=1LL*x*x%mod ) if ( y&1 ) z=1LL*z*x%mod;
    return z;
}
std::unordered_map<int,int> map;
inline int BSGS ( int y,int z,int p )
{
    int m=(int)sqrt(p)+1;map.clear();
    if ( z%p==1 ) return 1;
    for ( int i=0,t=z;i<m;i++,t=1LL*t*y%p ) map[t]=i;
    for ( int i=1,v=power(y,m,p),t=v;i<=m+1;i++,t=1LL*t*v%p ) if ( map.count(t) ) return (i*m-map[t]+1)%p;
    return -1;
}
signed main()
{
    for ( int T=read();T--; )
    {
        int p=read(),a=read(),b=read(),X1=read(),t=read(),ans=-1;
        if ( X1==t ) ans=1;
        else if ( a==1 )
        {
            if ( b ) ans=(1LL*(t-X1+p)*power(b,p-2,p)+1)%p;
        }
        else if ( a==0 )
        {
            if ( b==t ) ans=2;
        }
        else
        {
            int res=1LL*b*power((a-1+p)%p,p-2,p)%p;
            ans=BSGS(a,1LL*(t+res)%p*power((X1+res)%p,p-2,p)%p,p);
        }
        if ( !ans ) ans+=p;
        printf("%d
",ans);
    }
    return 0;
}

二、常系数二阶齐次线性递推

递推式：$a_i=pa_{i-1}+qa_{i-2}(p,q eq 0)$.

通项公式：设 $alpha,eta$ 为方程 $x^2-px-q=0$ 两根，则 $a_n=Aalpha^n+Beta^n$. 其中 $A,B$ 为常数，可用待定系数法求解.

证明：构造 $a_n - alpha a_{n-1}= eta(a_{n-1}-alpha a_{n-2})$ 即可.

注意到 $alpha,eta$ 常为无理数，可用二次剩余或矩阵维护.

【例题】DTOJ 4558「JLOI2015」有意义的字符串

【题目链接】http://59.61.75.5:8018/problem/4558

【题解】

考虑共轭. 构造方程 $x^2+frac{d-b^2}{4}x+b=0$. 设两根为 $p,q$. 则题目要求为 $p^n$.

考虑 $f_n=p^n+q^n=(p+q)(p^{n-1}+q^{n-1})-pq(p^{n-2}+q^{n-2})=(p+q)f_{n-1}-pqf_{n-2}$.

直接计算即可. 又 $q^n$ 较小，对其讨论即可.

【代码】

#include<bits/stdc++.h>
const unsigned long long mod=7528443412579576937llu;
unsigned long long b,d,n,ans;
inline unsigned long long mul ( unsigned long long a,unsigned long long b )
{
    unsigned long long res=0;
    while ( b )
    {
        if ( b&1 ) res=(res+a)%mod;
        a=(a+a)%mod;b>>=1;
    }
    return res;
}
struct matrix
{
    unsigned long long a[2][2];
    inline friend matrix operator * ( const matrix &A,const matrix &B )
    {
        matrix C;memset(C.a,0,sizeof(C.a));
        for ( int i=0;i<2;i++ ) for ( int k=0;k<2;k++ ) for ( int j=0;j<2;j++ ) C.a[i][j]=(C.a[i][j]+mul(A.a[i][k],B.a[k][j]))%mod;
        return C;
    }
    inline friend matrix operator ^ ( matrix A,unsigned long long n )
    {
        matrix B;B.a[0][0]=B.a[1][1]=1;B.a[0][1]=B.a[1][0]=0;
        while ( n )
        {
            if ( n&1 ) B=B*A;
            A=A*A;n>>=1;
        }
        return B;
    }
}A;
signed main()
{
    scanf("%llu%llu%llu",&b,&d,&n);
    if ( n==0 ) ans=2;
    else if ( n==1 ) ans=b;
    else A.a[0][0]=b,A.a[0][1]=mul((d-b*b%mod+mod)%mod,mul((mod+1)>>1,(mod+1)>>1)),A.a[1][0]=1,A=A^(n-1),ans=(mul(A.a[0][0],b)+mul(A.a[0][1],2))%mod;
    if ( b*b!=d and !(n&1) ) ans=(ans-1+mod)%mod;
    return !printf("%llu
",ans);
}

三、Fibonacci 循环节

显然 Fibonacci 数列在 $mod p$ 意义下有循环节，且循环节长度不大于 $p^2$.

由于证明较为复杂，此处仅给出结论：

设 $G(p)$ 为 $mod p$ 意义下循环节. 对 $p$ 分解有： $p=prodlimits_{i=1}^{k} p_i^{a_i}$.

则 $G(p)=prodlimits_{i=1}^{k} G(p_i^{a_i})=prodlimits_{i=1}^{k} G(p_i) imes p_i^{a_i-1}$.

对于 $G(p_i)$ 又有：若 $p_i^2 equiv 1 pmod 5$，则 $G(p_i) mid p-1$. 否则 $G(p_i) mid 2(p+1)$.

又实际上无需确切算出最小循环节，用 $p-1$ 及 $2(p+1)$ 代替即可.

【例题】DTOJ 4688 迫害 DJ

【题目链接】http://59.61.75.5:8018/problem/4688

【题解】

迭代 $k$ 次，每次求出循环节即可.

【代码】

#include<bits/stdc++.h>
inline long long read ( void )
{
    long long x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
long long a,b,pr[1000000],cnt[1000000],tot;
std::unordered_map<long long,long long> map;
inline long long mul ( long long x,long long y,long long p ) { return x*y-(long long)((long double)x*y/p)*p; }
std::pair<long long,long long> calc ( long long n,long long p )
{
    if ( n==-1 ) return std::make_pair(1,0);
    if ( n==0 ) return std::make_pair(0,1);
    if ( n%2==1 )
    {
        std::pair<long long,long long> res=calc(n-1,p);
        return std::make_pair(res.second,(res.first+res.second)%p);
    }
    else
    {
        std::pair<long long,long long> res=calc(n/2,p);
        return std::make_pair(mul((2*res.second-res.first+p)%p,res.first,p),(mul(res.first,res.first,p)+mul(res.second,res.second,p))%p);
    }
}
inline long long get ( long long n,long long p )
{
    std::pair<long long,long long> res=calc(2*n-1,p);
    return (mul(a,res.first,p)+mul((b-a),res.second,p))%p;
}
inline long long G ( long long p )
{
    if ( p==2 ) return 3;
    if ( p==3 ) return 8;
    if ( p==5 ) return 20;
    return ( p%5==4 or p%5==1 ) ? p-1 : 2*(p+1) ;
}
inline long long find_loop ( long long n )
{
    tot=0;
    for ( long long i=2;i*i<=n;i++ ) if ( !(n%i) )
    {
        pr[++tot]=i;cnt[tot]=0;
        while ( !(n%i) ) cnt[tot]++,n/=i;
    }
    if ( n>1 ) pr[++tot]=n,cnt[tot]=1;
    long long ans=1;
    for ( int i=1;i<=tot;i++ )
    {
        long long res=G(pr[i]);
        for ( int j=2;j<=cnt[i];j++ ) res*=pr[i];
        ans=ans/std::__gcd(ans,res)*res;
    }
    return ans;
}
long long n,k,Mod[150];
signed main()
{
    for ( int T=read();T--; )
    {
        a=read();b=read();n=read();k=read();Mod[1]=read();
        for ( int i=2;i<=k;i++ ) Mod[i]=find_loop(Mod[i-1]);
        n%=find_loop(Mod[k]);
        for ( int i=k;i;i-- ) n=get(n,Mod[i]);
        printf("%lld
",n);
    }
    return 0;
}

【例题】DTOJ 3098 斐波那契数

【题目链接】http://59.61.75.5:8018/problem/3098

【题解】

预处理两个循环节，每次迭代即可.

【代码】

#include<bits/stdc++.h>
inline long long read ( void )
{
    long long x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
inline long long read ( long long mod )
{
    long long x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=((x<<1)+(x<<3)+(ch^48))%mod;
    return x;
}
long long a,b,pr[1000000],cnt[1000000],tot;
std::unordered_map<long long,long long> map;
inline long long mul ( long long x,long long y,long long p )
{
    long long z=0;
    for ( ;y;y>>=1,x=(x+x)%p ) if ( y&1 ) z=(z+x)%p;
    return z;
}
std::pair<long long,long long> calc ( long long n,long long p )
{
    if ( n==-1 ) return std::make_pair(1,0);
    if ( n==0 ) return std::make_pair(0,1);
    if ( n%2==1 )
    {
        std::pair<long long,long long> res=calc(n-1,p);
        return std::make_pair(res.second,(res.first+res.second)%p);
    }
    else
    {
        std::pair<long long,long long> res=calc(n/2,p);
        return std::make_pair(mul((2*res.second-res.first+p)%p,res.first,p),(mul(res.first,res.first,p)+mul(res.second,res.second,p))%p);
    }
}
inline long long get ( long long n,long long p )
{
    std::pair<long long,long long> res=calc(n,p);
    return res.first;
}
inline long long G ( long long p )
{
    if ( p==2 ) return 3;
    if ( p==3 ) return 8;
    if ( p==5 ) return 20;
    return ( p%5==4 or p%5==1 ) ? p-1 : 2*(p+1) ;
}
inline long long find_loop ( long long n )
{
    tot=0;
    for ( long long i=2;i*i<=n;i++ ) if ( !(n%i) )
    {
        pr[++tot]=i;cnt[tot]=0;
        while ( !(n%i) ) cnt[tot]++,n/=i;
    }
    if ( n>1 ) pr[++tot]=n,cnt[tot]=1;
    long long ans=1;
    for ( int i=1;i<=tot;i++ )
    {
        long long res=G(pr[i]);
        for ( int j=2;j<=cnt[i];j++ ) res*=pr[i];
        ans=ans/std::__gcd(ans,res)*res;
    }
    return ans;
}
long long n,Mod[150];
signed main()
{
    Mod[1]=1000000007;Mod[2]=find_loop(Mod[1]);Mod[3]=find_loop(Mod[2]);
    for ( int T=read();T--; ) printf("%lld
",get(get(read(Mod[3]),Mod[2]),Mod[1]));
    return 0;
}

---

---

Day4 20.02.14

一、狄利克雷卷积

定义：数论函数 $f,g$ 的狄利克雷卷积为：$f*g(n)=sumlimits_{d mid n} f(d) imes g(frac{n}{d})$.

显然狄利克雷卷积具有非常多的性质，如交换律、结合律、分配律等.

定义：元函数 $epsilon(n)=[n=1]$. 则有：$epsilon*f=f$.

定义如下函数：$mathbf{id}(n)=n,mathbf{1}(n)=1$.

定义逆元：若 $f*g=epsilon$，则称 $g$ 为 $f$ 的逆元.

定义积性函数：若数论函数 $f$ 满足 $f(nm)=f(n)f(m)(nperp m)$，则 $f$ 为积性函数. 若 $f(nm)=f(n)f(m)$，则 $f$ 为完全积性函数.

显然根据定义可以得到，积性函数的狄利克雷卷积也是积性函数.

根据积性函数的特性，可以与线性筛结合预处理. 在筛的时候顺路处理函数值即可.

显然欧拉函数 $varphi$ 及上述提到的数论函数均为积性函数.

结论：$mathbf{id}=varphi * mathbf{1}$. 证明根据定义.

二、莫比乌斯反演

定义 $mathbf{1}$ 的逆元为 $mu$. 则有：$epsilon=mu * mathbf{1}$.

则有：若 $g=f*mathbf{1}$，则 $g*mu = f*mathbf{1}*mu$，即 $f=g*mu$.

莫比乌斯反演定理：若 $g(n)=sumlimits_{dmid n} f(d)$，则 $f(n)=sumlimits_{dmid n} g(d)mu(frac{n}{d})$.

关于 $mu$ 有：$mu(n)=egin{cases}1&n=1\(-1)^k&n=prodlimits_{i=1}^k p_i\0&otherwiseend{cases}$.

另一方向，即枚举的 $d$ 改为 $n$ 的倍数，也有类似结论的反演定理.

具体为：若 $g(n)=sumlimits_{nmid d} f(d)$，则 $f(n)=sumlimits_{nmid d} g(d)mu(frac{d}{n})$.

这个式子常常作为容斥的思路出现，往往答案需要求的为 $f(1)$，而 $f$ 本身不容易求出，$g$ 容易求出，则有 $f(1)=sumlimits_{i=1}^{n}g(i)mu(i)$.

【例题】DTOJ 4698 亚特兰大

【题目链接】http://59.61.75.5:8018/problem/4690

【题解】

显然 $gcd$ 为 $x$ 的路径条数不好求，但 $gcd$ 为 $x$ 的倍数的可以求出.

考虑反演，用第二个公式可得：$f(1)=sumlimits_{i=1}^{n}g(i)mu(i)$.

现在考虑如何求 $g$. 显然对于 $g(d)$，将所有边权为 $d$ 的边提出处理即可. 具体处理方式可以用可持久化并查集实现.

【代码】

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
const int maxn=100000+10;
const int maxm=1000000+10;
int n,m=1000000,Q,mu[maxm];bool flag[maxm],vis[200];
struct edge { int u,v,w,t; } E[maxn];
std::vector<edge> G[maxm];
std::pair<int,int> q[200];
std::vector<int> pr,B;
int f[maxn],siz[maxn],W[200][200],map[maxn],tot;
long long Ans[200],res;
inline void R ( int x ) { f[x]=x;siz[x]=1; }
inline int find ( int x ) { return f[x]==x ? x : find(f[x]); }
std::stack<std::pair<int,int>> st;
std::set<int> A;
inline void merge ( int u,int v )
{
    u=find(u);v=find(v);
    res-=1LL*siz[u]*(siz[u]-1)/2;
    res-=1LL*siz[v]*(siz[v]-1)/2;
    siz[u]+=siz[v];f[v]=u;
    res+=1LL*siz[u]*(siz[u]-1)/2;
    st.push(std::make_pair(u,v));
}
inline void split ( void )
{
    auto p=st.top();st.pop();
    int u=p.first,v=p.second;
    res-=1LL*siz[u]*(siz[u]-1)/2;
    siz[u]-=siz[v];f[v]=v;
    res+=1LL*siz[v]*(siz[v]-1)/2;
    res+=1LL*siz[u]*(siz[u]-1)/2;
}
signed main()
{
    mu[1]=1;
    for ( int i=2;i<=m;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),mu[i]=-1;
        for ( int j:pr )
        {
            if ( i*j>m ) break;
            flag[i*j]=true;mu[i*j]=-mu[i];
            if ( !(i%j) ) { mu[i*j]=0;break; }
        }
    }
    n=read();
    for ( int i=1;i<n;i++ ) E[i].u=read(),E[i].v=read(),E[i].w=read(),E[i].t=-1;
    Q=read();
    for ( int i=1,k,x;i<=Q;i++ ) k=read(),x=read(),E[k].t=0,q[i]=std::make_pair(k,x),A.insert(k);
    for ( int k:A ) B.push_back(k),map[k]=tot++;
    for ( int k=0;k<tot;k++ ) for ( int i=0;i<=Q;i++ ) W[k][i]=E[B[k]].w;
    for ( int i=1;i<=Q;i++ ) for ( int j=i;j<=Q;j++ ) W[map[q[i].first]][j]=q[i].second;
    for ( int i=1;i<n;i++ ) if ( E[i].t==-1 ) for ( int j=1;j*j<=E[i].w;j++ ) if ( !(E[i].w%j) )
    {
        if ( mu[j] ) G[j].push_back(E[i]);
        if ( j*j!=E[i].w and mu[E[i].w/j] ) G[E[i].w/j].push_back(E[i]);
    }
    for ( int j=0;j<=Q;j++ ) for ( int i=0;i<tot;i++ )
    {
        int u=E[B[i]].u,v=E[B[i]].v;
        for ( int k=1;k*k<=W[i][j];k++ ) if ( !(W[i][j]%k) )
        {
            if ( mu[k] ) G[k].push_back((edge){u,v,W[i][j],j});
            if ( k*k!=W[i][j] and mu[W[i][j]/k] ) G[W[i][j]/k].push_back((edge){u,v,W[i][j],j});
        }
    }
    for ( int x=1;x<=m;x++ ) if ( mu[x] )
    {
        res=0;
        for ( edge e:G[x] ) R(e.u),R(e.v);
        for ( int i=0,j;i<(int)G[x].size();i=j )
        {
            for ( j=i;j<(int)G[x].size() and G[x][j].t==G[x][i].t;j++ ) merge(G[x][j].u,G[x][j].v);
            if ( ~G[x][i].t )
            {
                Ans[G[x][i].t]+=mu[x]*res;
                for ( int k=j-1;k>=i;k-- ) split();
            }
            else
            {
                for ( int k=0;k<=Q;k++ ) vis[k]=false;
                for ( int k=j;k<(int)G[x].size();k++ ) vis[G[x][k].t]=true;
                for ( int k=0;k<=Q;k++ ) if ( !vis[k] ) Ans[k]+=mu[x]*res;
            }
        }
    }
    for ( int i=0;i<=Q;i++ ) printf("%lld
",Ans[i]);
    return 0;
}

【例题】求 $sumlimits_{i=1}^n sumlimits_{j=1}^m gcd(i,j)$.

【题解一】

设 $N=min(m,n)$，考虑每个数的贡献：

$$egin{aligned} sumlimits_{i=1}^n sumlimits_{j=1}^m gcd(i,j) &=sumlimits_{i=1}^N sumlimits_{j=1}^N gcd(i,j)
\ &= sumlimits_{d=1}^N sumlimits_{i=1}^{lfloor frac{n}{d} floor} sumlimits_{j=1}^{lfloor frac{m}{d} floor} [gcd(i,j)==1]
\ &=sumlimits_{d=1}^N d sumlimits_{i=1}^{lfloor frac{N}{d} floor} mu(i) lfloor frac{n}{id} floor lfloor frac{m}{id} floor
\&=  sumlimits_{d=1}^N  d sumlimits_{d mid T} mu(frac{T}{d}) lfloor frac{n}{T} floor lfloor frac{m}{T} floor
\&=sumlimits_{T=1}^{N} sumlimits_{d mid T} dmu(frac{T}{d}) lfloor frac{n}{T} floor lfloor frac{m}{T} floor
\&=sumlimits_{T=1}^{N} varphi(T)  lfloor frac{n}{T} floor lfloor frac{m}{T} floor end{aligned}$$

这里应用到了一个重要结论：$varphi=mathbf{id}*mu$，证明由狄利克雷卷积部分最后给出的一个式子得到.

【题解二】

$$egin{aligned} sumlimits_{i=1}^n sumlimits_{j=1}^m gcd(i,j) &= sumlimits_{i=1}^n sumlimits_{j=1}^m sumlimits_{dmid i, dmid j} varphi(d)
\&=sumlimits_{d=1}^{min(n,m)} varphi(d) lfloor frac{n}{d} floor lfloor frac{m}{d} floor end{aligned}$$

---

---

Day5 20.02.15

一、整除分块（数论分块）

在数论的题目中，常常出现 $sumlimits_{i=1}^{n} f(i) lfloor frac{n}{i} floor$ 形式的式子（如 Day4 莫比乌斯反演例题）。

考虑到 $lfloor frac{n}{i} floor$ 的取值只有 $sqrt n$ 种（易证），若 $f$ 的前缀和能够求出，则可以 $O(sqrt n)$ 计算上式.

代码如下：

inline long long solve ( int n )
{
    long long res=0;
    for ( int i=1,j;i<=n;i=j+1 ) j=n/(n/i),res+=(calc(j)-calc(i-1))*(n/i);
    return res;
}

二、杜教筛

令 $S(n)=sumlimits_{i=1}^{n} f(i)$.

则有：$sumlimits_{i=1}^{n} f*g(i)=sumlimits_{i=1}^{n} sumlimits_{xy=i}f(x)g(y)=sumlimits_{y=1}^{n}g(y)sumlimits_{xy leq n} f(x)=sumlimits_{y=1}^{n}g(y)S(lfloor frac{n}{y} floor)$.

因此可得 $g(1)S(n)=sumlimits_{i=1}^{n} f*g(i)-sumlimits_{y=2}^{n}g(y)S(lfloor frac{n}{y} floor)$.

因此若 $f*g$ 及 $g$ 的前缀和容易求出，则可以预处理 $n^{frac{2}{3}}$ 的部分，剩余的 $O(sqrt n)$ 计算（数论分块+记忆化搜索），效率 $O(n^{frac{2}{3}})$（为什么是这个效率我也不知道）.

【例题】luogu P4213 【模板】杜教筛（Sum）

【题解】$mathbf{id}=varphi * mathbf{1},epsilon=mu * mathbf{1}$. 杜教筛.

【代码】

#include<bits/stdc++.h>
#define print(n) std::cerr<<#n<<'='<<n<<'
'
inline int read ( void )
{
    int n=0;char ch;bool flag=true;
    while ( !isdigit(ch=getchar()) ) if ( ch=='-' ) flag=false;
    for ( n=ch^48;isdigit(ch=getchar()); ) n=(n<<1)+(n<<3)+(ch^48);
    return flag ? n : -n ;
}
const int m=1700000;
int pr[m+10],tot;
long long mu[m+10],phi[m+10];
bool flag[m+10];
std::unordered_map<int,long long> summu,sumphi;
inline int sum_mu ( int n )
{
    if ( n<=m ) return mu[n];
    if ( summu.count(n) ) return summu[n];
    long long ans=1;
    for ( register int i=2,j;i<=n;i=j+1 ) j=n/(n/i),ans-=(j-i+1)*sum_mu(n/i);
    return summu[n]=ans;
}
inline long long sum_phi ( int n )
{
    if ( n<=m ) return phi[n];
    if ( sumphi.count(n) ) return sumphi[n];
    long long ans=1LL*n*(n+1)/2;
    for ( register int i=2,j;i<=n;i=j+1 ) j=n/(n/i),ans-=(j-i+1)*sum_phi(n/i);
    return sumphi[n]=ans;
}
int main()
{
    mu[1]=phi[1]=1;
    for ( register int i=2;i<=m;i++ )
    {
        if ( !flag[i] ) pr[++tot]=i,phi[i]=i-1,mu[i]=-1;
        for ( register int j=1;j<=tot && i*pr[j]<=m;j++ )
        {
            flag[i*pr[j]]=true;
            if ( !(i%pr[j]) ) { mu[i*pr[j]]=0;phi[i*pr[j]]=phi[i]*pr[j];break; }
            mu[i*pr[j]]=-mu[i];phi[i*pr[j]]=phi[i]*(pr[j]-1);
        }
    }
    for ( register int i=2;i<=m;i++ ) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
    for ( int T=read(),n;T;T-- ) n=read(),printf("%lld %d
",sum_phi(n),sum_mu(n));
    return 0;
}

 

---

---

Day6 20.02.16

【例题】http://59.61.75.5:8018/problem/4704

【题解】https://www.cnblogs.com/RenSheYu/p/12266604.html

【例题】https://www.luogu.com.cn/problem/P3768

【题解】

$$egin{aligned}sumlimits_{i=1}^n sumlimits_{j=1}^n ijgcd(i,j) &=  sumlimits_{i=1}^n sumlimits_{j=1}^n ij sumlimits_{d mid i,d mid j} varphi(d)
\&=sumlimits_{d=1}^n d^2varphi(d) (sumlimits_{i=1}^{lfloor frac{n}{d} floor} i) (sumlimits_{j=1}^{lfloor frac{n}{d} floor} j)
\&=sumlimits_{d=1}^n d^2varphi(d) (S(lfloor frac{n}{d} floor))^2 end{aligned}$$

其中 $S(n)=sumlimits_{i=1}^{n} i=frac{n(n+1)}{2}$. 设 $f(n)=n^2varphi(n),g(n)=n^2,h(n)=sumlimits_{d=1}^n d^2varphi(d)$. 则有：

$$h(n)=sumlimits_{i=1}^{n} f*g(i)-sumlimits_{y=2}^{n}g(y)S(lfloor frac{n}{y} floor)$$

$$f*g(i)=sumlimits_{dmid i} f(d)g(lfloor frac{n}{d} floor)=i^2 sumlimits_{dmid i} varphi(d)=i^3$$

$$h(n)=sumlimits_{i=1}^{n} i^3-sumlimits_{y=2}^{n}y^2S(lfloor frac{n}{y} floor)$$

杜教筛.

【代码】

#include<bits/stdc++.h>
const int maxm=10000000+10;
int m=10000000,p,val[maxm],inv2,inv6;
long long n,phi[maxm];
bool flag[maxm];
std::vector<int> pr;
std::map<long long,int> map;
inline int power ( int x,int y )
{
    int z=1;
    for ( ;y;y>>=1,x=1LL*x*x%p ) if ( y&1 ) z=1LL*z*x%p;
    return z;
}
inline int S1 ( long long n ) { return n%p*((n+1)%p)%p*inv2%p; }
inline int S2 ( long long n ) { return n%p*((n+1)%p)%p*((2*n+1)%p)%p*inv6%p; }
inline int S3 ( long long n ) { return 1LL*S1(n)*S1(n)%p; }
inline int h ( long long n )
{
    if ( n<=m ) return val[n];
    if ( map.count(n) ) return map[n];
    int res=S3(n);
    for ( long long i=2,j;i<=n;i=j+1 ) j=n/(n/i),res=(res-1LL*h(n/i)*(S2(j)-S2(i-1)+p)%p+p)%p;
    return map[n]=res;
}
signed main()
{
    scanf("%d%lld",&p,&n);int ans=0;
    inv2=power(2,p-2);inv6=power(6,p-2);
    phi[1]=1;
    for ( int i=2;i<=m;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),phi[i]=i-1;
        for ( int j:pr )
        {
            if ( i*j>m ) break;
            flag[i*j]=true;phi[i*j]=phi[i]*(j-1);
            if ( !(i%j) ) { phi[i*j]=phi[i]*j;break; }
        }
    }
    for ( int i=1;i<=m;i++ ) val[i]=(val[i-1]+1LL*i*i%p*phi[i])%p;
    for ( long long i=1,j;i<=n;i=j+1 ) j=n/(n/i),ans=(ans+1LL*(h(j)-h(i-1)+p)*S1(n/i)%p*S1(n/i))%p;
    return !printf("%d
",ans);
}

【例题】http://59.61.75.5:8018/problem/1551

【题解】

显然 $O(abc)$ 枚举. 考虑线性求 $mathbf{d}(n)$.

考虑约数个数定理，有约数个数为素因数幂的指数加一之积.

记 $num(i)$ 为 $i$ 包含的最小素因数的幂的指数. 线性筛处理即可. 细节参见代码.

【代码】

#include<bits/stdc++.h>
const unsigned int mod=1073741824;
const int maxn=8000000+10;
int a,b,c,n,num[maxn];
unsigned d[maxn],ans;
bool flag[maxn];
std::vector<int> pr;
signed main()
{
    scanf("%d%d%d",&a,&b,&c);n=a*b*c;
    d[1]=1;
    for ( int i=2;i<=n;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),num[i]=1,d[i]=2;
        for ( int j:pr )
        {
            if ( i*j>n ) break;
            flag[i*j]=true;num[i*j]=1;d[i*j]=d[i]*d[j];
            if ( !(i%j) ) { num[i*j]=num[i]+1;d[i*j]=d[i]/(num[i]+1)*(num[i*j]+1);break; }
        }
    }
    for ( int i=1;i<=a;i++ ) for ( int j=1;j<=b;j++ ) for ( int k=1;k<=c;k++ ) ans=(ans+d[i*j*k])%mod;
    return !printf("%u
",ans);
}

【例题】http://59.61.75.5:8018/problem/2304

【题解】

设 $f(n)=sumlimits_{i=1}^n lfloor frac{n}{i} floor=sumlimits_{i=1}^n mathbf{d}(i)$.

$$egin{aligned}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} mathbf{d}(ij)
&= sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} sumlimits_{amid i} sumlimits_{b mid j} [gcd(a,b)==1]
\&= sumlimits_{a=1}^n sumlimits_{b=1}^m [gcd(a,b)==1] lfloor frac{n}{a}  floor lfloor frac{m}{b}  floor
\&= sumlimits_{d=1}^{min(n,m)} mu(d) (sumlimits_{a=1}^{lfloor frac{n}{d} floor} lfloor frac{n}{ad} floor)(sumlimits_{b=1}^{lfloor frac{m}{d} floor} lfloor frac{m}{bd} floor)
\&= sumlimits_{d=1}^{min(n,m)} mu(d)f(lfloor frac{n}{d} floor)f(lfloor frac{m}{d} floor)end{aligned}$$

约数个数可以线性筛处理（同上题），数论分块即可.

【代码】

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
const int maxn=50000+10;
int a,b,c,mu[maxn],num[maxn];
long long d[maxn];bool flag[maxn];
std::vector<int> pr;
signed main()
{
    mu[1]=d[1]=1;
    for ( int i=2;i<=50000;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),num[i]=1,d[i]=2,mu[i]=-1;
        for ( int j:pr )
        {
            if ( i*j>50000 ) break;
            flag[i*j]=true;num[i*j]=1;d[i*j]=d[i]*d[j];mu[i*j]=-mu[i];
            if ( !(i%j) ) { num[i*j]=num[i]+1;d[i*j]=d[i]/(num[i]+1)*(num[i*j]+1);mu[i*j]=0;break; }
        }
    }
    for ( int i=2;i<=50000;i++ ) mu[i]+=mu[i-1],d[i]+=d[i-1];
    for ( int T=read();T--; )
    {
        int n=read(),m=read();long long ans=0;
        for ( int i=1,j;i<=std::min(n,m);i=j+1 )
        {
            j=std::min(n/(n/i),m/(m/i));
            ans+=1LL*d[n/i]*d[m/i]*(mu[j]-mu[i-1]);
        }
        printf("%lld
",ans);
    }
    return 0;
}

---

---

Day7 20.02.17

【例题】https://www.luogu.com.cn/problem/P1829

【题解】

设 $n < m, S(n) = sumlimits_{i=1}^{n} i =frac{n(n+1)}{2}$.

$$egin{aligned} sumlimits_{i=1}^n sumlimits_{j=1}^m operatorname{lcm}(i,j)
&=sumlimits_{i=1}^n sumlimits_{j=1}^m frac{ij}{gcd(i,j)}
\&=sumlimits_{i=1}^n sumlimits_{j=1}^m sumlimits_{dmid i,dmid j} frac{ij}{d} [gcd(i,j)==d]
\&=sumlimits_{d=1}^n dsumlimits_{i=1}^{lfloor frac{n}{d} floor} sumlimits_{j=1}^{lfloor frac{m}{d} floor} ij sumlimits_{emid i,emid j} mu(e)
\&=sumlimits_{d=1}^n dsumlimits_{e=1}^{lfloor frac{n}{d} floor} e^2mu(e) sumlimits_{i=1}^{lfloor frac{n}{de} floor} sumlimits_{j=1}^{lfloor frac{m}{de} floor} ij
\&=sumlimits_{d=1}^n dsumlimits_{e=1}^{lfloor frac{n}{d} floor} e^2mu(e) S(lfloor frac{n}{de} floor) S(lfloor frac{m}{de} floor)
\&=sumlimits_{d=1}^n dsumlimits_{e=1}^{lfloor frac{n}{d} floor} e^2mu(e) S(lfloor frac{lfloor frac{n}{d} floor}{e} floor) S(lfloor frac{lfloor frac{m}{d} floor}{e} floor)
end{aligned}$$

设 $f(n,m)=sumlimits_{i=1}^n i^2 mu(i) S(lfloor frac{n}{i} floor)S(lfloor frac{m}{i} floor)$.

$$sumlimits_{d=1}^n dsumlimits_{e=1}^{lfloor frac{n}{d} floor} e^2mu(e) S(lfloor frac{lfloor frac{n}{d} floor}{e} floor) S(lfloor frac{lfloor frac{m}{d} floor}{e} floor)=sumlimits_{d=1}^n df(lfloor frac{n}{d} floor,lfloor frac{m}{d} floor)$$.

显然外层一个数论分块，内层求 $f$ 再用一个数论分块，效率 $O(n+m)$. 即可通过本题.

【代码】

#include<bits/stdc++.h>
const int mod=20101009,inv2=10050505;
const int maxn=10000000+10;
int ans,n,m,mu[maxn],sum[maxn];
std::vector<int> pr;
bool flag[maxn];
inline int S1 ( int n ) { return 1LL*n*(n+1)%mod*inv2%mod; }
inline int calc ( int n,int m )
{
    int res=0;
    for ( int i=1,j;i<=n;i=j+1 ) j=std::min(n/(n/i),m/(m/i)),res=(res+1LL*(sum[j]-sum[i-1]+mod)*S1(n/i)%mod*S1(m/i))%mod;
    return res;
}
signed main()
{
    scanf("%d%d",&n,&m);
    if ( n>m ) std::swap(n,m);
    mu[1]=1;
    for ( int i=2;i<=n;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),mu[i]=-1;
        for ( int j:pr )
        {
            if ( i*j>n ) break;
            flag[i*j]=true;mu[i*j]=-mu[i];
            if ( !(i%j) ) { mu[i*j]=0;break; }
        }
    }
    for ( int i=1;i<=n;i++ ) sum[i]=(sum[i-1]+1LL*i*i%mod*mu[i]%mod+mod)%mod;
    for ( int i=1,j;i<=n;i=j+1 ) j=std::min(n/(n/i),m/(m/i)),ans=(ans+1LL*calc(n/i,m/i)*(S1(j)-S1(i-1)+mod))%mod;
    return !printf("%d
",ans);
}

【例题】https://www.luogu.com.cn/problem/P2257

【题解】

从本篇题解起，往下均有：$n<m,mathbf{d}(n)=sumlimits_{dmid n}1,S_k(n)=sumlimits_{i=1}^n i^k$. 其中 $k=1$ 有时省略不写.

显然考虑枚举素数 $p$，考虑每个素数的贡献. 设 $pr$ 为素数集合.

$$egin{aligned} sumlimits_{pin pr} sumlimits_{i=1}^n sumlimits_{j=1}^n [gcd(i,j)==p]
&=sumlimits_{pin pr} sumlimits_{i=1}^{lfloor frac{n}{p} floor} sumlimits_{i=1}^{lfloor frac{m}{p} floor} [gcd(i,j)==1]
\&=sumlimits_{pin pr} sumlimits_{i=1}^{lfloor frac{n}{p} floor} sumlimits_{i=1}^{lfloor frac{m}{p} floor} sumlimits_{dmid i,dmid j} mu(d)
\&=sumlimits_{pint pr} sumlimits_{d=1}^{lfloor frac{n}{p} floor} mu(d) lfloor frac{n}{pd} floor lfloor frac{m}{pd} floor
end{aligned}$$

一般遇到这类整除 $pd$ 型的式子，通常会设 $T=pd$，然后转化为枚举 $T$. 此处采用同样的方式：

$$egin{aligned} sumlimits_{pint pr} sumlimits_{d=1}^{lfloor frac{n}{p} floor} mu(d) lfloor frac{n}{pd} floor lfloor frac{m}{pd} floor
&xlongequal[quad]{T=pd}sumlimits_{T=1}^n sumlimits_{pmid T,pin pr} mu(frac{T}{p}) lfloor frac{n}{T} floor lfloor frac{m}{T} floor
\&=sumlimits_{T=1}^n lfloor frac{n}{T} floor lfloor frac{m}{T} floor (sumlimits_{pmid T,pin pr} mu(frac{T}{p}))
end{aligned}$$

后面一部分可以用埃氏筛处理，即对每个素数 $p$ 枚举其倍数，分别考虑贡献即可.

【代码】

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
const int maxn=10000000+10;
int mu[maxn];long long f[maxn];
std::vector<int> pr;
bool flag[maxn];
signed main()
{
    mu[1]=1;
    for ( int i=2;i<=10000000;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),mu[i]=-1;
        for ( int j:pr )
        {
            if ( i*j>10000000 ) break;
            flag[i*j]=true;mu[i*j]=-mu[i];
            if ( !(i%j) ) { mu[i*j]=0;break; }
        }
    }
    for ( int p:pr ) for ( int i=1;i*p<=10000000;i++ ) f[i*p]+=mu[i];
    for ( int i=2;i<=10000000;i++ ) f[i]+=f[i-1];
    for ( int T=read();T--; )
    {
        int n=read(),m=read();long long ans=0;
        if ( n>m ) std::swap(n,m);
        for ( int i=1,j;i<=n;i=j+1 ) j=std::min(n/(n/i),m/(m/i)),ans+=1LL*(n/i)*(m/i)*(f[j]-f[i-1]);
        printf("%lld
",ans);
    }
    return 0;
}

【例题】http://59.61.75.5:8018/problem/1219

【题解】

显然先用二维前缀和的思想容斥. 考虑每一部分如何计算.

$$egin{aligned} sumlimits_{i=1}^n sumlimits_{j=1}^m [gcd(i,j)==k]
&= sumlimits_{i=1}^{lfloor frac{n}{k} floor} sumlimits_{j=1}^{lfloor frac{m}{k} floor} [gcd(i,j)==1]
\&=sumlimits_{i=1}^{lfloor frac{n}{k} floor} sumlimits_{j=1}^{lfloor frac{m}{k} floor} sumlimits_{dmid i,dmid j} mu(d)
\&=sumlimits_{d=1}^{lfloor frac{n}{k} floor} mu(d)  lfloor frac{lfloor frac{n}{k} floor}{d} floor lfloor frac{lfloor frac{m}{k} floor}{d} floor
end{aligned}$$

【代码】

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;
    while ( !isdigit(ch=getchar()) ) ;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return x;
}
const int maxn=50000+10;
int mu[maxn];
std::vector<int> pr;
bool flag[maxn];
inline long long solve ( int n,int m,int x )
{
    n/=x;m/=x;long long res=0;
    for ( int i=1,j;i<=n and i<=m;i=j+1 ) j=std::min(n/(n/i),m/(m/i)),res+=1LL*(mu[j]-mu[i-1])*(n/i)*(m/i);
    return res;
}
signed main()
{
    mu[1]=1;
    for ( int i=2;i<=50000;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),mu[i]=-1;
        for ( int j:pr )
        {
            if ( i*j>50000 ) break;
            flag[i*j]=true;mu[i*j]=-mu[i];
            if ( !(i%j) ) { mu[i*j]=0;break; }
        }
    }
    for ( int i=2;i<=50000;i++ ) mu[i]+=mu[i-1];
    for ( int T=read();T--; )
    {
        int a=read(),b=read(),c=read(),d=read(),k=read();
        printf("%lld
",solve(b,d,k)-solve(b,c-1,k)-solve(a-1,d,k)+solve(a-1,c-1,k));
    }
    return 0;
}

【例题】http://59.61.75.5:8018/problem/3748

【题解】

考虑枚举 $gcd=x$.

$$egin{aligned} sumlimits_{x=1}^a sumlimits_{i=1}^n sumlimits_{j=1}^m [gcd(i,j)==x] frac{ij}{gcd(i,j)}
&= sumlimits_{x=1}^a x sumlimits_{i=1}^{lfloor frac{n}{x} floor} sumlimits_{j=1}^{lfloor frac{m}{x} floor} ij[gcd(i,j)==1]
\&= sumlimits_{x=1}^a x sumlimits_{d=1}^{lfloor frac{n}{x} floor} d^2 mu(d) S(lfloor frac{n}{xd} floor) S(lfloor frac{m}{xd} floor)
\&xlongequal[quad]{T=xd} sumlimits_{T=1}^n S(lfloor frac{n}{T} floor) S(lfloor frac{m}{T} floor) sumlimits_{xmid T,xleq a} frac{T^2}{x} mu(frac{T}{x})
end{aligned}$$

后面那部分由于 $a$ 在变化无法处理. 因此考虑离线.

将所有询问按 $a$ 排序. 考虑如何插入每个 $x$ 的贡献. 显然式子是个前缀和类似的形式，可以枚举 $x$ 的倍数并用树状数组更新及查询.

【代码】

#include<bits/stdc++.h>
const int mod=1000000007;
const int maxn=100000+10;
struct Query { int n,m,a,id; } Q[maxn];
int mu[maxn],ans[maxn],t[maxn],q;
bool flag[maxn];
std::vector<int> pr;
inline void add ( int x,int y ) { for ( int i=x;i<=100000;i+=i&(-i) ) (t[i]+=y)%=mod; }
inline int query ( int x ) { int res=0;for ( int i=x;i;i-=i&(-i) ) (res+=t[i])%=mod;return res; }
inline int S1 ( int n ) { return 1LL*n*(n+1)/2%mod; }
signed main()
{
    mu[1]=1;
    for ( int i=2;i<=100000;i++ )
    {
        if ( !flag[i] ) pr.push_back(i),mu[i]=-1;
        for ( int j:pr )
        {
            if ( i*j>100000 ) break;
            flag[i*j]=true;mu[i*j]=-mu[i];
            if ( !(i%j) ) { mu[i*j]=0;break; }
        }
    }
    scanf("%d",&q);
    for ( int i=1;i<=q;i++ ) scanf("%d%d%d",&Q[i].n,&Q[i].m,&Q[i].a),Q[i].id=i;
    std::sort(Q+1,Q+q+1,[](const Query &q1,const Query &q2){return q1.a<q2.a;});
    for ( int i=1,j=1;i<=q;i++ )
    {
        while ( j<=Q[i].a )
        {
            for ( int k=j;k<=100000;k+=j ) add(k,1LL*k*(k/j*mu[k/j]%mod+mod)%mod);
            j++;
        }
        int n=std::min(Q[i].n,Q[i].m),m=std::max(Q[i].n,Q[i].m);
        for ( int l=1,r;l<=n;l=r+1 ) r=std::min(n/(n/l),m/(m/l)),ans[Q[i].id]=(ans[Q[i].id]+1LL*(query(r)-query(l-1)+mod)*S1(n/l)%mod*S1(m/l))%mod;
    }
    for ( int i=1;i<=q;i++ ) printf("%d
",ans[i]);
    return 0;
}

---

---

习题

【DTOJ4729】函数

【题解】显然 $H(d)=prodlimits_{p ext{ is prime}} g(p,d)$ 是积性函数。由二平方和定理（打表找规律）可知 $H(p^e)=3e+1\,(pequiv 1pmod 4)$，Min_25 筛（先咕着）即可.

【代码】

#include<bits/stdc++.h>
long long n,m,pr[100000],tot,Pow[100000],g[100000],S1[100000][5],S2[100000][5],sum1[100000],sum2[100000];
bool flag[100000];
inline long long f ( long long p,long long k ) { return (p%4==1) ? 3*k+1 : 1 ; }
inline long long Min_25 ( long long x,long long k )
{
    if ( x<=1 or pr[k]>x ) return 0;
    long long ans=((x<=m)?sum1[x]:sum2[n/x])-Pow[k-1];
    for ( long long i=k;i<=tot and pr[i]*pr[i]<=x;i++ ) for ( long long v=pr[i],j=1;v*pr[i]<=x;v*=pr[i],j++ ) ans+=Min_25(x/v,i+1)*f(pr[i],j)+f(pr[i],j+1);
    return ans;
}
signed main()
{
    long long T;
    for ( scanf("%lld",&T);T--; )
    {
        scanf("%lld",&n);m=(long long)sqrt(n);
        for ( long long i=1;i<=m;i++ ) flag[i]=true;
        flag[1]=false;tot=0;
        for ( long long i=2;i<=m;i++ )
        {
            if ( flag[i] ) pr[++tot]=i,Pow[tot]=Pow[tot-1]+((i%4==1)?4:1);
            for ( long long j=1;j<=tot and pr[j]*i<=m;j++ )
            {
                flag[i*pr[j]]=false;
                if ( !(i%pr[j]) ) break;
            }
        }
        for ( long long i=1;i<=n;i=n/(n/(i+1)) )
        {
            if ( i<=m ) for ( long long j=0;j<4;j++ ) S1[i][j]=(i-1)/4+((i-1)%4>=(j+3)%4);
            else for ( long long j=0;j<4;j++ ) S2[n/i][j]=(i-1)/4+((i-1)%4>=(j+3)%4);
            if ( i==n ) break;
        }
        for ( long long i=1;i<=tot;i++ )
        {
            long long p=pr[i],p2=p*p;
            for ( long long a=n;a>=p2;a=n/(n/a+1) ) for ( long long j=0;j<4;j++ )
                if ( a<=m )
                {
                    if ( a/p<=m ) S1[a][p*j%4]-=S1[a/p][j]-S1[p-1][j];
                    else S1[a][p*j%4]-=S2[n/(a/p)][j]-S1[p-1][j];
                }
                else
                {
                    if ( a/p<=m ) S2[n/a][p*j%4]-=S1[a/p][j]-S1[p-1][j];
                    else S2[n/a][p*j%4]-=S2[n/(a/p)][j]-S1[p-1][j];
                }
        }
        for ( long long i=1;i<=n;i=n/(n/(i+1)) )
        {
            if ( i<=m ) S1[i][1]--;
            else S2[n/i][1]--;
            if ( i==n ) break;
        }
        for ( long long i=1;i<=n;i=n/(n/(i+1)) )
        {
            if ( i<=m ) sum1[i]=4*S1[i][1]+S1[i][3]+1;
            else sum2[n/i]=4*S2[n/i][1]+S2[n/i][3]+1;
            if ( i==n ) break;
        }
        printf("%lld
",Min_25(n,1)+1);
    }
    return 0;
}

【DTOJ3726】旧试题

 【题解】

$$egin{aligned} sumlimits_{i=1}^A sumlimits_{j=1}^B sumlimits_{k=1}^C mathbf{d}(ijk)
&= sumlimits_{i=1}^A sumlimits_{j=1}^B sumlimits_{k=1}^C sumlimits_{amid i} sumlimits_{bmid j} sumlimits_{cmid k}[a perp b][bperp c][c perp a]
\&= sumlimits_{a=1}^A sumlimits_{b=1}^B sumlimits_{c=1}^C [a perp b][bperp c][c perp a] lfloor frac{A}{a} floor lfloor frac{B}{b} floor lfloor frac{C}{c} floor
\&= sumlimits_{a=1}^A sumlimits_{b=1}^B sumlimits_{c=1}^C lfloor frac{A}{a} floor lfloor frac{B}{b} floor lfloor frac{C}{c} floor sumlimits_{umid a,u mid b} mu(u) sumlimits_{vmid b,v mid c}mu(v) sumlimits_{wmid c,w mid a}mu(w)
\&= sumlimits_{u=1}^{min(A,B)} sumlimits_{v=1}^{min(B,C)} sumlimits_{w=1}^{min(C,A)} mu(u) mu(v) mu(w) sumlimits_{operatorname{lcm}(w,u)mid a} lfloorfrac{A}{a} floorsumlimits_{operatorname{lcm}(u,v)mid b} lfloorfrac{B}{b} floor sumlimits_{operatorname{lcm}(v,w)mid c} lfloorfrac{C}{c} floor
end{aligned}$$
设 $f(n,x)=sumlimits_{n mid d} lfloor frac{x}{d} floor$.
$$egin{aligned} sumlimits_{i=1}^A sumlimits_{j=1}^B sumlimits_{k=1}^C mathbf{d}(ijk)
&= sumlimits_{u=1}^{min(A,B)} sumlimits_{v=1}^{min(B,C)} sumlimits_{w=1}^{min(C,A)} mu(u) mu(v) mu(w) sumlimits_{operatorname{lcm}(w,u)mid a} lfloorfrac{A}{a} floorsumlimits_{operatorname{lcm}(u,v)mid b} lfloorfrac{B}{b} floorsumlimits_{operatorname{lcm}(v,w)mid c} lfloorfrac{C}{c} floor
\&=sumlimits_{u=1}^{min(A,B)} sumlimits_{v=1}^{min(B,C)} sumlimits_{w=1}^{min(C,A)} mu(u) mu(v) mu(w) f(operatorname{lcm}(w,u),A)f(operatorname{lcm}(u,v),B)f(operatorname{lcm}(v,w),C)
end{aligned}$$

反演部分结束. 考虑如何计算上面那个式子.

考虑怎样的三个数 $u,v,w$ 会对答案造成贡献，显然当且仅当 $mu(u) eq 0,mu(v) eq 0,mu(w) eq 0, operatorname{lcm}(w,u) < max(A,B,C), operatorname{lcm}(u,v) < max(A,B,C), operatorname{lcm}(v,w) < max(A,B,C)$ 时 $u,v,w$ 才会造成贡献.

把上述条件拆成三个：对于 $u,v$ 有：$mu(u) eq 0,mu(v) eq 0,operatorname{lcm}(u,v) < max(A,B,C)$. 考虑建图计算三元环个数. 对于满足条件的 $u,v$. 连边 $(u,v)$，边权 $operatorname{lcm}(u,v)$.

考虑枚举边权连边，先用埃氏筛筛出每个数的所有素因子. 考虑对于每个 $operatorname{lcm}$ 枚举子集求 $u$. 对 $v$ 枚举 $u$ 关于 $operatorname{lcm}$ 补集的超集即可. 打表证明边数为 $760741$.

$O(m sqrt m)$ 三元环计数+卡常即可.

【代码】

#include<bits/stdc++.h>
const int mod=1000000007;
const int maxn=100000+10;
const int maxm=1000000+10;
int pr[maxn],mu[maxn],a[4],d[maxn],tot;
long long f[4][maxn];
bool flag[maxn];
std::vector<int> F[maxn];
typedef struct { int v,w; } edge;
typedef struct { int u,v,w; } Edge;
std::vector<edge> G[maxn];
std::vector<Edge> E;
inline void solve()
{
    scanf("%d%d%d",&a[1],&a[2],&a[3]);std::sort(a+1,a+4);
    for ( int t=1;t<=3;t++ ) for ( int i=1;i<=a[t];i++ ) for ( int j=1;j*i<=a[t];j++ ) (f[t][i]+=a[t]/(i*j))%=mod;
    for ( int i=1;i<=a[3];i++ ) if ( mu[i] )
    {
        int n=(int)F[i].size(),U=(1<<n)-1;
        for ( int S=0;S<=U;S++ )
        {
            int u=1;
            for ( int j=0;j<n;j++ ) if ( S&(1<<j) ) u*=F[i][j];
            int C=U^S;
            for ( int T=C;;T=(T+1)|C )
            {
                int v=1;
                for ( int j=0;j<n;j++ ) if ( T&(1<<j) ) v*=F[i][j];
                d[u]++;d[v]++;
                if ( u<v ) E.push_back((Edge){u,v,i});
                if ( T==U ) break;
            }
        }
    }
    long long ans=0;
    for ( Edge e:E )
    {
        if ( d[e.u]<d[e.v] ) std::swap(e.u,e.v);
        G[e.u].push_back((edge){e.v,e.w});
    }
    for ( int u=1;u<=a[3];u++ )
    {
        for ( edge i:G[u] ) for ( edge j:G[i.v] )
        {
            int v=i.v,w=j.v,A=i.w,B=j.w,C=1LL*w*u/std::__gcd(w,u);
            if ( C>a[3] ) continue;
            (ans+=mu[u]*mu[v]*mu[w]*(f[2][B]*f[3][C]+f[2][C]*f[3][B])%mod*f[1][A])%=mod;
            (ans+=mu[u]*mu[v]*mu[w]*(f[2][C]*f[3][A]+f[2][A]*f[3][C])%mod*f[1][B])%=mod;
            (ans+=mu[u]*mu[v]*mu[w]*(f[2][A]*f[3][B]+f[2][B]*f[3][A])%mod*f[1][C])%=mod;
            (ans+=mod)%=mod;
        }
    }
    for ( Edge e:E )
    {
        int x=e.u,y=e.v,z=e.w;
        (ans+=mu[x]*(f[1][z]*f[2][z]%mod*f[3][y]+f[1][z]*f[2][y]%mod*f[3][z]+f[1][y]*f[2][z]%mod*f[3][z]))%=mod;
        (ans+=mu[y]*(f[1][z]*f[2][z]%mod*f[3][x]+f[1][z]*f[2][x]%mod*f[3][z]+f[1][x]*f[2][z]%mod*f[3][z]))%=mod;
        (ans+=mod)%=mod;
    }
    for ( int i=1;i<=a[1];i++ ) (ans+=mu[i]*f[1][i]*f[2][i]%mod*f[3][i])%=mod,(ans+=mod)%=mod;
    printf("%lld
",ans);E.clear();
    for ( int i=1;i<=a[3];i++ ) G[i].clear(),d[i]=0;
    for ( int t=1;t<=3;t++ ) for ( int i=1;i<=a[t];i++ ) f[t][i]=0;
}
signed main()
{
    mu[1]=1;
    for ( int i=2;i<=100000;i++ )
    {
        if ( !flag[i] ) pr[++tot]=i,mu[i]=-1;
        for ( int j=1;j<=tot and pr[j]*i<=100000;j++ )
        {
            flag[i*pr[j]]=true;mu[i*pr[j]]=-mu[i];
            if ( !(i%pr[j]) ) { mu[i*pr[j]]=0;break; }
        }
    }
    for ( int j=1;j<=tot;j++ ) for ( int i=1;i*pr[j]<=100000;i++ ) F[i*pr[j]].push_back(pr[j]);
    int T;
    for ( scanf("%d",&T);T--; ) solve();
    return 0;
}