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) \bmod n$. 其中 $c$ 为一个固定的随机数.
先任意取 $x_1$，然后令 $x_i=f(x_{i-1})(i \geq 2)$ 即可.
观察生成的随机数列，发现分布较为均匀.
但注意到一个问题：生成的数列在一定长度后会陷入循环，形成 $\rho$ 形数列.
因此每次判断 $d=\gcd(\lvert x_i-x_1 \rvert,n)$ 是否符合要求. 当 $x_i=x_1$ 时退出循环. 由引理知这样一个环的期望长度为 $O(\sqrt n)$.
考虑对求解的过程进行优化. 由于计算 $\gcd$ 的时间较慢，考虑每隔一段时间统一计算. 采用倍增的思想，统计 $res=(\prod \lvert x_i - x_1 \rvert) \bmod n$.
当 $res=0$ 时出现环. 每隔 $2^k$ 次计算一次 $\gcd$，重新初始化即可. 代码如下：
1 inline long long work ( long long x ) 2 { 3 long long s=0,t=0,c=1LL*rand()*rand()*rand()%(x-1)+1; 4 for ( int step=1;;step<<=1,s=t ) 5 { 6 long long res=1; 7 for ( int i=1;i<=step;i++ ) 8 { 9 t=(mul(t,t,x)+c)%x; 10 res=mul(res,std::abs(t-s),x); 11 if ( !(i%127) ) 12 { 13 long long d=std::__gcd(res,x); 14 if ( d!=1 ) return d; 15 } 16 } 17 long long d=std::__gcd(res,x); 18 if ( d!=1 ) return d; 19 } 20 return -1; 21 } 22 inline void Pollard_Rho ( long long x ) 23 { 24 if ( x==1 or x<=ans ) return; 25 if ( Miller_Rabbin(x) ) { ans=std::max(ans,x);return; } 26 long long y=x; 27 while ( y>=x ) y=work(x); 28 while ( !(x%y) ) x/=y; 29 Pollard_Rho(x);Pollard_Rho(y); 30 }
四、BSGS/exBSGS
1. BSGS
求解方程：$y^x \equiv z \pmod p$（$p$ 为素数）.
令 $x=am-b(0 \leq b \leq m-1)$，则有：$y^{am} \equiv z \times y^b \pmod p$.
等式右边有 $m$ 种取值，用 hash_table 或 unordered_map 记录.
等式左边枚举 $a$，查询右边是否出现过该值即可.
当 $m=\sqrt p$ 时效率最优，效率为 $O(\sqrt p)$.
实际写代码时由于数据原因，常将 $m$ 与一个常量（例如 $1000$，视数据而定）取 $\min$.
代码如下：
1 std::unordered_map<int,int> map; 2 inline int BSGS ( int y,int z,int p ) 3 { 4 int m=(int)sqrt(p)+1; 5 if ( z%p==1 ) return 0; 6 if ( !(y%p) ) return -1; 7 map.clear(); 8 for ( int i=0,t=z;i<m;i++,t=1LL*t*y%p ) map[t]=i; 9 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]; 10 return -1; 11 }
2. exBSGS
求解方程：$y^x \equiv z \pmod p$（$p$ 为任意正整数）.
将等式展开有：$y \times y^{x-1} + k \times p = z(k \in \mathbb{Z})$.
由裴属定理（Day2 内容）知，当且仅当 $d=gcd(y,k) \mid z$ 时有解.
若有解则有：$\frac{y}{d} \times y^{x-1} + k \times \frac{p}{d} = z$.
不断递归求解，最终得到 $\gcd(y,p)=1$ 时即可用 BSGS 求解.
设递归 $c$ 次，$g=\prod\limits_{i=1}^{c} d_i$.
则有：$y^{x-c} \times \frac{y^c}{g} \equiv \frac{z}{g} \pmod {\frac{p}{g}}$.
代码如下：
1 std::unordered_map<int,int> map; 2 inline int exBSGS ( int y,int p,int z ) 3 { 4 if ( z==1 ) return 0; 5 int c=0,g=1; 6 for ( int d=std::__gcd(y,p);d!=1;d=std::__gcd(y,p) ) 7 { 8 if ( z%d ) return -1; 9 c++;z/=d;p/=d;g=1LL*g*y/d%p; 10 if ( z==g ) return c; 11 } 12 int m=sqrt(p)+1;map.clear(); 13 for ( int i=0,t=z;i<m;i++,t=1LL*t*y%p ) map[t]=i; 14 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; 15 return -1; 16 }
Day2 20.02.12
一、裴蜀定理
裴蜀定理：方程 $\sum\limits_{i=1}^{n} a_i \times 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$.
$\therefore d \mid \sum\limits_{i=1}^{n} a_i \times x_i$，即 $\gcd(a_1,a_2,\dots,a_n) \mid b$.
②充分性：
考虑数学归纳法.
Ⅰ. $n=2$ 时成立（详见 exgcd 部分）.
Ⅱ. 假设 $n=k$ 时成立，考虑 $n=k+1$ 时，
设 $S=\sum\limits_{i=1}{k} a_i,d=\gcd(a_1,a_2,\dots,a_k)$.
则方程 $S \times X + a_{k+1} \times y = \gcd(d,a_{k+1})$ 有整数解.
即存在 $(X,Y)$ 使得 $\sum\limits_{i=1}^{k+1} a_i \times x_i = \gcd(a_1,a_2,\dots,a_{k+1})$ 一组解为 $x_i=X(i \in [1,k]),x_{k+1}=Y$.
$\therefore 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 \bmod m,p_2 \bmod m,\dots,p_{\varphi(m)} \bmod m\}$ 与 $\{x_1,x_2,\dots,x_{\varphi(m)}$ 一一对应.
累乘得：$\prod\limits_{i=1}^{\varphi(m)} p_i =\prod\limits_{i=1}^{\varphi(m)} x_i $，则 $a^{\varphi(m)} \equiv 1 \pmod m$.
扩展欧拉定理：当 $b \geq \varphi(m)$ 时，$a^b \equiv a^{b \bmod \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)}\times x$ 即可，$y$ 同理.
设：$\begin{cases}ax_1+by_1=\gcd(a,b)\\bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)\end{cases}$.
又 $\gcd(a,b)=\gcd(b,a \bmod b)$ 得，$ax_1+by_1=bx_2+(a \bmod b)y_2$.
将 $a\bmod b$ 展开并整理得：$ax_1+by_1=ay_2+b(x_2-\lfloor \frac{a}{b} \rfloor y_2)$.
不断递归处理即可. 边界条件为当 $b=0$ 时，$x=1,y=0$. 代码如下：
1 inline void exgcd ( int a,int b,int &x,int &y ) 2 { 3 if ( !b ) { x=1;y=0;return; } 4 exgcd(b,a%b,y,x);y-=a/b*x; 5 }
四、CRT/exCRT
求解方程组：$\begin{cases}x\equiv a_1 \pmod {n_1}\\x\equiv a_2 \pmod {n_2}\\ \dots \\x\equiv a_k \pmod {n_k}\end{cases}$. 其中 $\forall i \neq j, n_i \perp n_j$.
CRT（中国剩余定理）：
① 令 $n=\prod\limits_{i=1}^{k} n_i$.
② 对于第 $i$ 个方程：
Ⅰ. 计算 $m_i=\frac{n}{n_i}$.
Ⅱ. 计算 $m_i$ 在 $\bmod n_i$ 意义下的逆元 $m_i^{-1}$. 即方程 $am_i+bn_i=1$ 的最小正整数解 $a$，可使用 exgcd 求解.
Ⅲ. 令 $c_i=m_im_i^{-1}$（不对 $n_i$ 取模）.
③ 则唯一解为 $x \equiv \sum\limits_{i=1}^{k} a_ic_i \pmod n$.
证明直接代入即可. 代码如下：
1 inline long long CRT ( int n,long long *a,long long *m ) 2 { 3 long long M=1,res=0; 4 for ( int i=1;i<=n;i++ ) M*=m[i]; 5 for ( int i=1;i<=n;i++ ) 6 { 7 long long w=M/m[i],x,y;exgcd(w,m[i],x,y);x=(x%m[i]+m[i])%m[i]; 8 res=(res+mul(mul(a[i],w,M),x,M)+M)%M; 9 } 10 return res; 11 }
求解方程组：$\begin{cases}x\equiv a_1 \pmod {n_1}\\x\equiv a_2 \pmod {n_2}\\ \dots \\x\equiv a_k \pmod {n_k}\end{cases}$. 其中 $n_i$ 为任意正整数.
exCRT：
考虑 $\begin{cases} x\equiv a_1 \pmod {n_1}\\x\equiv a_2 \pmod {n_2} \end{cases}$. 有：$\begin{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)}$. 两两合并即可.
代码如下：
1 inline long long exCRT ( void ) 2 { 3 long long m=a[1],ans=r[1],x=0,y=0; 4 for ( int i=1;i<=n;i++ ) 5 { 6 long long B=((r[i]-ans)%a[i]+a[i])%a[i],gcd=exgcd(m,a[i],x,y); 7 x=mul(x,B/gcd,a[i]);ans+=m*x;m*=a[i]/gcd;ans=(ans+m)%m; 8 } 9 return (ans%m+m)%m; 10 }
五、Lucas/exLucas
Lucas 定理：对于素数 $p$ 有：$C_n^m \bmod p = C_{\lfloor n/p \rfloor}^{\lfloor m/p \rfloor} \times C_{n \bmod p}^{m \bmod p} \bmod 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 \bmod p$，其中 $p$ 为任意正整数.
考虑对 $p$ 进行分解. 转化为每个在模素因数的幂意义下的组合数取余. 分别解出后用 CRT 合并.
现在考虑 $C_n^m \bmod p^k$. 将组合数展开为阶乘形式，分别计算.
对于 $n!$，考虑提取所有 $p$ 的倍数的因子，以 $n=19,p=3$ 为例有：
$19!=1\times 2\times \dots \times 19=1\times 2 \times 4\times 5\times \dots \times 16 \times 17 \times 19 \times (3 \times 6 \times \dots \times 18)$.
前面一部分整段的直接算，最后一段每个单独计算即可. 后面部分提取出一个 $p$ 后转化为更小的阶乘，递归求解.
代码如下：
1 inline long long FCT ( long long n,long long x,long long p ) 2 { 3 if ( !n ) return 1; 4 long long s=1; 5 for ( long long i=1;i<=p;i++ ) if ( i%x ) s=mul(s,i,p); 6 s=power(s,n/p,p); 7 for ( long long i=n/p*p+1;i<=n;i++ ) if ( i%x ) s=mul(s,i,p); 8 return mul(s,FCT(n/x,x,p),p); 9 } 10 inline long long INV ( long long a,long long b ) { long long x,y;exgcd(a,b,x,y);return (x%b+b)%b; } 11 inline long long multiLucas ( long long n,long long m,long long x,long long p ) 12 { 13 int cnt=0; 14 for ( long long i=n;i;i/=x ) cnt+=i/x; 15 for ( long long i=m;i;i/=x ) cnt-=i/x; 16 for ( long long i=n-m;i;i/=x ) cnt-=i/x; 17 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); 18 } 19 inline long long exLucas ( long long n,long long m,long long p ) 20 { 21 int cnt=0;long long mod[20]={0},a[20]={0}; 22 for ( long long i=2;i*i<=p;i++ ) if ( !(p%i) ) 23 { 24 mod[++cnt]=1; 25 while ( !(p%i) ) mod[cnt]*=i,p/=i; 26 a[cnt]=multiLucas(n,m,i,mod[cnt]); 27 } 28 if ( p>1 ) mod[++cnt]=p,a[cnt]=multiLucas(n,m,p,p); 29 return CRT(cnt,a,mod); 30 }
Day3 20.02.13
一、常系数一阶线性递推
递推式：$a_i=ka_{i-1}+b(k \neq 0)$.
通项公式：$a_n=\begin{cases}a_0+nb&k=1\\a_0k^n+\frac{b}{k-1}k^n-\frac{b}{k-1}&k\neq 1\end{cases}$.
证明：构造 $a_i+x_0=k(a_{i-1}+x_0)$ 后运用等比数列相关知识即可.
【例题】DTOJ 1378 随机数生成器
【题目链接】http://59.61.75.5:8018/problem/1378
【题解】
直接分类讨论即可. 解方程时用 BSGS 求解.
【代码】
1 #include<bits/stdc++.h> 2 inline int read ( void ) 3 { 4 int x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 inline int power ( int x,int y,int mod ) 10 { 11 int z=1; 12 for ( ;y;y>>=1,x=1LL*x*x%mod ) if ( y&1 ) z=1LL*z*x%mod; 13 return z; 14 } 15 std::unordered_map<int,int> map; 16 inline int BSGS ( int y,int z,int p ) 17 { 18 int m=(int)sqrt(p)+1;map.clear(); 19 if ( z%p==1 ) return 1; 20 for ( int i=0,t=z;i<m;i++,t=1LL*t*y%p ) map[t]=i; 21 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; 22 return -1; 23 } 24 signed main() 25 { 26 for ( int T=read();T--; ) 27 { 28 int p=read(),a=read(),b=read(),X1=read(),t=read(),ans=-1; 29 if ( X1==t ) ans=1; 30 else if ( a==1 ) 31 { 32 if ( b ) ans=(1LL*(t-X1+p)*power(b,p-2,p)+1)%p; 33 } 34 else if ( a==0 ) 35 { 36 if ( b==t ) ans=2; 37 } 38 else 39 { 40 int res=1LL*b*power((a-1+p)%p,p-2,p)%p; 41 ans=BSGS(a,1LL*(t+res)%p*power((X1+res)%p,p-2,p)%p,p); 42 } 43 if ( !ans ) ans+=p; 44 printf("%d\n",ans); 45 } 46 return 0; 47 }
二、常系数二阶齐次线性递推
递推式：$a_i=pa_{i-1}+qa_{i-2}(p,q \neq 0)$.
通项公式：设 $\alpha,\beta$ 为方程 $x^2-px-q=0$ 两根，则 $a_n=A\alpha^n+B\beta^n$. 其中 $A,B$ 为常数，可用待定系数法求解.
证明：构造 $a_n - \alpha a_{n-1}= \beta(a_{n-1}-\alpha a_{n-2})$ 即可.
注意到 $\alpha,\beta$ 常为无理数，可用二次剩余或矩阵维护.
【例题】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$ 较小，对其讨论即可.
【代码】
1 #include<bits/stdc++.h> 2 const unsigned long long mod=7528443412579576937llu; 3 unsigned long long b,d,n,ans; 4 inline unsigned long long mul ( unsigned long long a,unsigned long long b ) 5 { 6 unsigned long long res=0; 7 while ( b ) 8 { 9 if ( b&1 ) res=(res+a)%mod; 10 a=(a+a)%mod;b>>=1; 11 } 12 return res; 13 } 14 struct matrix 15 { 16 unsigned long long a[2][2]; 17 inline friend matrix operator * ( const matrix &A,const matrix &B ) 18 { 19 matrix C;memset(C.a,0,sizeof(C.a)); 20 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; 21 return C; 22 } 23 inline friend matrix operator ^ ( matrix A,unsigned long long n ) 24 { 25 matrix B;B.a[0][0]=B.a[1][1]=1;B.a[0][1]=B.a[1][0]=0; 26 while ( n ) 27 { 28 if ( n&1 ) B=B*A; 29 A=A*A;n>>=1; 30 } 31 return B; 32 } 33 }A; 34 signed main() 35 { 36 scanf("%llu%llu%llu",&b,&d,&n); 37 if ( n==0 ) ans=2; 38 else if ( n==1 ) ans=b; 39 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; 40 if ( b*b!=d and !(n&1) ) ans=(ans-1+mod)%mod; 41 return !printf("%llu\n",ans); 42 }
三、Fibonacci 循环节
显然 Fibonacci 数列在 $\bmod p$ 意义下有循环节，且循环节长度不大于 $p^2$.
由于证明较为复杂，此处仅给出结论：
设 $G(p)$ 为 $\bmod p$ 意义下循环节. 对 $p$ 分解有： $p=\prod\limits_{i=1}^{k} p_i^{a_i}$.
则 $G(p)=\prod\limits_{i=1}^{k} G(p_i^{a_i})=\prod\limits_{i=1}^{k} G(p_i) \times 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$ 次，每次求出循环节即可.
【代码】
1 #include<bits/stdc++.h> 2 inline long long read ( void ) 3 { 4 long long x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 long long a,b,pr[1000000],cnt[1000000],tot; 10 std::unordered_map<long long,long long> map; 11 inline long long mul ( long long x,long long y,long long p ) { return x*y-(long long)((long double)x*y/p)*p; } 12 std::pair<long long,long long> calc ( long long n,long long p ) 13 { 14 if ( n==-1 ) return std::make_pair(1,0); 15 if ( n==0 ) return std::make_pair(0,1); 16 if ( n%2==1 ) 17 { 18 std::pair<long long,long long> res=calc(n-1,p); 19 return std::make_pair(res.second,(res.first+res.second)%p); 20 } 21 else 22 { 23 std::pair<long long,long long> res=calc(n/2,p); 24 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); 25 } 26 } 27 inline long long get ( long long n,long long p ) 28 { 29 std::pair<long long,long long> res=calc(2*n-1,p); 30 return (mul(a,res.first,p)+mul((b-a),res.second,p))%p; 31 } 32 inline long long G ( long long p ) 33 { 34 if ( p==2 ) return 3; 35 if ( p==3 ) return 8; 36 if ( p==5 ) return 20; 37 return ( p%5==4 or p%5==1 ) ? p-1 : 2*(p+1) ; 38 } 39 inline long long find_loop ( long long n ) 40 { 41 tot=0; 42 for ( long long i=2;i*i<=n;i++ ) if ( !(n%i) ) 43 { 44 pr[++tot]=i;cnt[tot]=0; 45 while ( !(n%i) ) cnt[tot]++,n/=i; 46 } 47 if ( n>1 ) pr[++tot]=n,cnt[tot]=1; 48 long long ans=1; 49 for ( int i=1;i<=tot;i++ ) 50 { 51 long long res=G(pr[i]); 52 for ( int j=2;j<=cnt[i];j++ ) res*=pr[i]; 53 ans=ans/std::__gcd(ans,res)*res; 54 } 55 return ans; 56 } 57 long long n,k,Mod[150]; 58 signed main() 59 { 60 for ( int T=read();T--; ) 61 { 62 a=read();b=read();n=read();k=read();Mod[1]=read(); 63 for ( int i=2;i<=k;i++ ) Mod[i]=find_loop(Mod[i-1]); 64 n%=find_loop(Mod[k]); 65 for ( int i=k;i;i-- ) n=get(n,Mod[i]); 66 printf("%lld\n",n); 67 } 68 return 0; 69 }
【例题】DTOJ 3098 斐波那契数
【题目链接】http://59.61.75.5:8018/problem/3098
【题解】
预处理两个循环节，每次迭代即可.
【代码】
1 #include<bits/stdc++.h> 2 inline long long read ( void ) 3 { 4 long long x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 inline long long read ( long long mod ) 10 { 11 long long x=0;char ch; 12 while ( !isdigit(ch=getchar()) ) ; 13 for ( x=ch^48;isdigit(ch=getchar()); ) x=((x<<1)+(x<<3)+(ch^48))%mod; 14 return x; 15 } 16 long long a,b,pr[1000000],cnt[1000000],tot; 17 std::unordered_map<long long,long long> map; 18 inline long long mul ( long long x,long long y,long long p ) 19 { 20 long long z=0; 21 for ( ;y;y>>=1,x=(x+x)%p ) if ( y&1 ) z=(z+x)%p; 22 return z; 23 } 24 std::pair<long long,long long> calc ( long long n,long long p ) 25 { 26 if ( n==-1 ) return std::make_pair(1,0); 27 if ( n==0 ) return std::make_pair(0,1); 28 if ( n%2==1 ) 29 { 30 std::pair<long long,long long> res=calc(n-1,p); 31 return std::make_pair(res.second,(res.first+res.second)%p); 32 } 33 else 34 { 35 std::pair<long long,long long> res=calc(n/2,p); 36 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); 37 } 38 } 39 inline long long get ( long long n,long long p ) 40 { 41 std::pair<long long,long long> res=calc(n,p); 42 return res.first; 43 } 44 inline long long G ( long long p ) 45 { 46 if ( p==2 ) return 3; 47 if ( p==3 ) return 8; 48 if ( p==5 ) return 20; 49 return ( p%5==4 or p%5==1 ) ? p-1 : 2*(p+1) ; 50 } 51 inline long long find_loop ( long long n ) 52 { 53 tot=0; 54 for ( long long i=2;i*i<=n;i++ ) if ( !(n%i) ) 55 { 56 pr[++tot]=i;cnt[tot]=0; 57 while ( !(n%i) ) cnt[tot]++,n/=i; 58 } 59 if ( n>1 ) pr[++tot]=n,cnt[tot]=1; 60 long long ans=1; 61 for ( int i=1;i<=tot;i++ ) 62 { 63 long long res=G(pr[i]); 64 for ( int j=2;j<=cnt[i];j++ ) res*=pr[i]; 65 ans=ans/std::__gcd(ans,res)*res; 66 } 67 return ans; 68 } 69 long long n,Mod[150]; 70 signed main() 71 { 72 Mod[1]=1000000007;Mod[2]=find_loop(Mod[1]);Mod[3]=find_loop(Mod[2]); 73 for ( int T=read();T--; ) printf("%lld\n",get(get(read(Mod[3]),Mod[2]),Mod[1])); 74 return 0; 75 }
Day4 20.02.14
一、狄利克雷卷积
定义：数论函数 $f,g$ 的狄利克雷卷积为：$f*g(n)=\sum\limits_{d \mid n} f(d) \times 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)(n\perp 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)=\sum\limits_{d\mid n} f(d)$，则 $f(n)=\sum\limits_{d\mid n} g(d)\mu(\frac{n}{d})$.
关于 $\mu$ 有：$\mu(n)=\begin{cases}1&n=1\\(-1)^k&n=\prod\limits_{i=1}^k p_i\\0&otherwise\end{cases}$.
另一方向，即枚举的 $d$ 改为 $n$ 的倍数，也有类似结论的反演定理.
具体为：若 $g(n)=\sum\limits_{n\mid d} f(d)$，则 $f(n)=\sum\limits_{n\mid d} g(d)\mu(\frac{d}{n})$.
这个式子常常作为容斥的思路出现，往往答案需要求的为 $f(1)$，而 $f$ 本身不容易求出，$g$ 容易求出，则有 $f(1)=\sum\limits_{i=1}^{n}g(i)\mu(i)$.
【例题】DTOJ 4698 亚特兰大
【题目链接】http://59.61.75.5:8018/problem/4690
【题解】
显然 $\gcd$ 为 $x$ 的路径条数不好求，但 $\gcd$ 为 $x$ 的倍数的可以求出.
考虑反演，用第二个公式可得：$f(1)=\sum\limits_{i=1}^{n}g(i)\mu(i)$.
现在考虑如何求 $g$. 显然对于 $g(d)$，将所有边权为 $d$ 的边提出处理即可. 具体处理方式可以用可持久化并查集实现.
【代码】
1 #include<bits/stdc++.h> 2 inline int read ( void ) 3 { 4 int x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 const int maxn=100000+10; 10 const int maxm=1000000+10; 11 int n,m=1000000,Q,mu[maxm];bool flag[maxm],vis[200]; 12 struct edge { int u,v,w,t; } E[maxn]; 13 std::vector<edge> G[maxm]; 14 std::pair<int,int> q[200]; 15 std::vector<int> pr,B; 16 int f[maxn],siz[maxn],W[200][200],map[maxn],tot; 17 long long Ans[200],res; 18 inline void R ( int x ) { f[x]=x;siz[x]=1; } 19 inline int find ( int x ) { return f[x]==x ? x : find(f[x]); } 20 std::stack<std::pair<int,int>> st; 21 std::set<int> A; 22 inline void merge ( int u,int v ) 23 { 24 u=find(u);v=find(v); 25 res-=1LL*siz[u]*(siz[u]-1)/2; 26 res-=1LL*siz[v]*(siz[v]-1)/2; 27 siz[u]+=siz[v];f[v]=u; 28 res+=1LL*siz[u]*(siz[u]-1)/2; 29 st.push(std::make_pair(u,v)); 30 } 31 inline void split ( void ) 32 { 33 auto p=st.top();st.pop(); 34 int u=p.first,v=p.second; 35 res-=1LL*siz[u]*(siz[u]-1)/2; 36 siz[u]-=siz[v];f[v]=v; 37 res+=1LL*siz[v]*(siz[v]-1)/2; 38 res+=1LL*siz[u]*(siz[u]-1)/2; 39 } 40 signed main() 41 { 42 mu[1]=1; 43 for ( int i=2;i<=m;i++ ) 44 { 45 if ( !flag[i] ) pr.push_back(i),mu[i]=-1; 46 for ( int j:pr ) 47 { 48 if ( i*j>m ) break; 49 flag[i*j]=true;mu[i*j]=-mu[i]; 50 if ( !(i%j) ) { mu[i*j]=0;break; } 51 } 52 } 53 n=read(); 54 for ( int i=1;i<n;i++ ) E[i].u=read(),E[i].v=read(),E[i].w=read(),E[i].t=-1; 55 Q=read(); 56 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); 57 for ( int k:A ) B.push_back(k),map[k]=tot++; 58 for ( int k=0;k<tot;k++ ) for ( int i=0;i<=Q;i++ ) W[k][i]=E[B[k]].w; 59 for ( int i=1;i<=Q;i++ ) for ( int j=i;j<=Q;j++ ) W[map[q[i].first]][j]=q[i].second; 60 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) ) 61 { 62 if ( mu[j] ) G[j].push_back(E[i]); 63 if ( j*j!=E[i].w and mu[E[i].w/j] ) G[E[i].w/j].push_back(E[i]); 64 } 65 for ( int j=0;j<=Q;j++ ) for ( int i=0;i<tot;i++ ) 66 { 67 int u=E[B[i]].u,v=E[B[i]].v; 68 for ( int k=1;k*k<=W[i][j];k++ ) if ( !(W[i][j]%k) ) 69 { 70 if ( mu[k] ) G[k].push_back((edge){u,v,W[i][j],j}); 71 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}); 72 } 73 } 74 for ( int x=1;x<=m;x++ ) if ( mu[x] ) 75 { 76 res=0; 77 for ( edge e:G[x] ) R(e.u),R(e.v); 78 for ( int i=0,j;i<(int)G[x].size();i=j ) 79 { 80 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); 81 if ( ~G[x][i].t ) 82 { 83 Ans[G[x][i].t]+=mu[x]*res; 84 for ( int k=j-1;k>=i;k-- ) split(); 85 } 86 else 87 { 88 for ( int k=0;k<=Q;k++ ) vis[k]=false; 89 for ( int k=j;k<(int)G[x].size();k++ ) vis[G[x][k].t]=true; 90 for ( int k=0;k<=Q;k++ ) if ( !vis[k] ) Ans[k]+=mu[x]*res; 91 } 92 } 93 } 94 for ( int i=0;i<=Q;i++ ) printf("%lld\n",Ans[i]); 95 return 0; 96 }
【例题】求 $\sum\limits_{i=1}^n \sum\limits_{j=1}^m \gcd(i,j)$.
【题解一】
设 $N=\min(m,n)$，考虑每个数的贡献：
$$\begin{aligned} \sum\limits_{i=1}^n \sum\limits_{j=1}^m \gcd(i,j) &=\sum\limits_{i=1}^N \sum\limits_{j=1}^N \gcd(i,j)
\\ &= \sum\limits_{d=1}^N \sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{d} \rfloor} [\gcd(i,j)==1]
\\ &=\sum\limits_{d=1}^N d \sum\limits_{i=1}^{\lfloor \frac{N}{d} \rfloor} \mu(i) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor
\\&= \sum\limits_{d=1}^N d \sum\limits_{d \mid T} \mu(\frac{T}{d}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor
\\&=\sum\limits_{T=1}^{N} \sum\limits_{d \mid T} d\mu(\frac{T}{d}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor
\\&=\sum\limits_{T=1}^{N} \varphi(T) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \end{aligned}$$
这里应用到了一个重要结论：$\varphi=\mathbf{id}*\mu$，证明由狄利克雷卷积部分最后给出的一个式子得到.
【题解二】
$$\begin{aligned} \sum\limits_{i=1}^n \sum\limits_{j=1}^m \gcd(i,j) &= \sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{d\mid i, d\mid j} \varphi(d)
\\&=\sum\limits_{d=1}^{\min(n,m)} \varphi(d) \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \end{aligned}$$
Day5 20.02.15
一、整除分块（数论分块）
在数论的题目中，常常出现 $\sum\limits_{i=1}^{n} f(i) \lfloor \frac{n}{i} \rfloor$ 形式的式子（如 Day4 莫比乌斯反演例题）。
考虑到 $\lfloor \frac{n}{i} \rfloor$ 的取值只有 $\sqrt n$ 种（易证），若 $f$ 的前缀和能够求出，则可以 $O(\sqrt n)$ 计算上式.
代码如下：
1 inline long long solve ( int n ) 2 { 3 long long res=0; 4 for ( int i=1,j;i<=n;i=j+1 ) j=n/(n/i),res+=(calc(j)-calc(i-1))*(n/i); 5 return res; 6 }
二、杜教筛
令 $S(n)=\sum\limits_{i=1}^{n} f(i)$.
则有：$\sum\limits_{i=1}^{n} f*g(i)=\sum\limits_{i=1}^{n} \sum\limits_{xy=i}f(x)g(y)=\sum\limits_{y=1}^{n}g(y)\sum\limits_{xy \leq n} f(x)=\sum\limits_{y=1}^{n}g(y)S(\lfloor \frac{n}{y} \rfloor)$.
因此可得 $g(1)S(n)=\sum\limits_{i=1}^{n} f*g(i)-\sum\limits_{y=2}^{n}g(y)S(\lfloor \frac{n}{y} \rfloor)$.
因此若 $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}$. 杜教筛.
【代码】
1 #include<bits/stdc++.h> 2 #define print(n) std::cerr<<#n<<'='<<n<<'\n' 3 inline int read ( void ) 4 { 5 int n=0;char ch;bool flag=true; 6 while ( !isdigit(ch=getchar()) ) if ( ch=='-' ) flag=false; 7 for ( n=ch^48;isdigit(ch=getchar()); ) n=(n<<1)+(n<<3)+(ch^48); 8 return flag ? n : -n ; 9 } 10 const int m=1700000; 11 int pr[m+10],tot; 12 long long mu[m+10],phi[m+10]; 13 bool flag[m+10]; 14 std::unordered_map<int,long long> summu,sumphi; 15 inline int sum_mu ( int n ) 16 { 17 if ( n<=m ) return mu[n]; 18 if ( summu.count(n) ) return summu[n]; 19 long long ans=1; 20 for ( register int i=2,j;i<=n;i=j+1 ) j=n/(n/i),ans-=(j-i+1)*sum_mu(n/i); 21 return summu[n]=ans; 22 } 23 inline long long sum_phi ( int n ) 24 { 25 if ( n<=m ) return phi[n]; 26 if ( sumphi.count(n) ) return sumphi[n]; 27 long long ans=1LL*n*(n+1)/2; 28 for ( register int i=2,j;i<=n;i=j+1 ) j=n/(n/i),ans-=(j-i+1)*sum_phi(n/i); 29 return sumphi[n]=ans; 30 } 31 int main() 32 { 33 mu[1]=phi[1]=1; 34 for ( register int i=2;i<=m;i++ ) 35 { 36 if ( !flag[i] ) pr[++tot]=i,phi[i]=i-1,mu[i]=-1; 37 for ( register int j=1;j<=tot && i*pr[j]<=m;j++ ) 38 { 39 flag[i*pr[j]]=true; 40 if ( !(i%pr[j]) ) { mu[i*pr[j]]=0;phi[i*pr[j]]=phi[i]*pr[j];break; } 41 mu[i*pr[j]]=-mu[i];phi[i*pr[j]]=phi[i]*(pr[j]-1); 42 } 43 } 44 for ( register int i=2;i<=m;i++ ) mu[i]+=mu[i-1],phi[i]+=phi[i-1]; 45 for ( int T=read(),n;T;T-- ) n=read(),printf("%lld %d\n",sum_phi(n),sum_mu(n)); 46 return 0; 47 }
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
【题解】
$$\begin{aligned}\sum\limits_{i=1}^n \sum\limits_{j=1}^n ij\gcd(i,j) &= \sum\limits_{i=1}^n \sum\limits_{j=1}^n ij \sum\limits_{d \mid i,d \mid j} \varphi(d)
\\&=\sum\limits_{d=1}^n d^2\varphi(d) (\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor} i) (\sum\limits_{j=1}^{\lfloor \frac{n}{d} \rfloor} j)
\\&=\sum\limits_{d=1}^n d^2\varphi(d) (S(\lfloor \frac{n}{d} \rfloor))^2 \end{aligned}$$
其中 $S(n)=\sum\limits_{i=1}^{n} i=\frac{n(n+1)}{2}$. 设 $f(n)=n^2\varphi(n),g(n)=n^2,h(n)=\sum\limits_{d=1}^n d^2\varphi(d)$. 则有：
$$h(n)=\sum\limits_{i=1}^{n} f*g(i)-\sum\limits_{y=2}^{n}g(y)S(\lfloor \frac{n}{y} \rfloor)$$
$$f*g(i)=\sum\limits_{d\mid i} f(d)g(\lfloor \frac{n}{d} \rfloor)=i^2 \sum\limits_{d\mid i} \varphi(d)=i^3$$
$$h(n)=\sum\limits_{i=1}^{n} i^3-\sum\limits_{y=2}^{n}y^2S(\lfloor \frac{n}{y} \rfloor)$$
杜教筛.
【代码】
1 #include<bits/stdc++.h> 2 const int maxm=10000000+10; 3 int m=10000000,p,val[maxm],inv2,inv6; 4 long long n,phi[maxm]; 5 bool flag[maxm]; 6 std::vector<int> pr; 7 std::map<long long,int> map; 8 inline int power ( int x,int y ) 9 { 10 int z=1; 11 for ( ;y;y>>=1,x=1LL*x*x%p ) if ( y&1 ) z=1LL*z*x%p; 12 return z; 13 } 14 inline int S1 ( long long n ) { return n%p*((n+1)%p)%p*inv2%p; } 15 inline int S2 ( long long n ) { return n%p*((n+1)%p)%p*((2*n+1)%p)%p*inv6%p; } 16 inline int S3 ( long long n ) { return 1LL*S1(n)*S1(n)%p; } 17 inline int h ( long long n ) 18 { 19 if ( n<=m ) return val[n]; 20 if ( map.count(n) ) return map[n]; 21 int res=S3(n); 22 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; 23 return map[n]=res; 24 } 25 signed main() 26 { 27 scanf("%d%lld",&p,&n);int ans=0; 28 inv2=power(2,p-2);inv6=power(6,p-2); 29 phi[1]=1; 30 for ( int i=2;i<=m;i++ ) 31 { 32 if ( !flag[i] ) pr.push_back(i),phi[i]=i-1; 33 for ( int j:pr ) 34 { 35 if ( i*j>m ) break; 36 flag[i*j]=true;phi[i*j]=phi[i]*(j-1); 37 if ( !(i%j) ) { phi[i*j]=phi[i]*j;break; } 38 } 39 } 40 for ( int i=1;i<=m;i++ ) val[i]=(val[i-1]+1LL*i*i%p*phi[i])%p; 41 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; 42 return !printf("%d\n",ans); 43 }
【例题】http://59.61.75.5:8018/problem/1551
【题解】
显然 $O(abc)$ 枚举. 考虑线性求 $\mathbf{d}(n)$.
考虑约数个数定理，有约数个数为素因数幂的指数加一之积.
记 $num(i)$ 为 $i$ 包含的最小素因数的幂的指数. 线性筛处理即可. 细节参见代码.
【代码】
1 #include<bits/stdc++.h> 2 const unsigned int mod=1073741824; 3 const int maxn=8000000+10; 4 int a,b,c,n,num[maxn]; 5 unsigned d[maxn],ans; 6 bool flag[maxn]; 7 std::vector<int> pr; 8 signed main() 9 { 10 scanf("%d%d%d",&a,&b,&c);n=a*b*c; 11 d[1]=1; 12 for ( int i=2;i<=n;i++ ) 13 { 14 if ( !flag[i] ) pr.push_back(i),num[i]=1,d[i]=2; 15 for ( int j:pr ) 16 { 17 if ( i*j>n ) break; 18 flag[i*j]=true;num[i*j]=1;d[i*j]=d[i]*d[j]; 19 if ( !(i%j) ) { num[i*j]=num[i]+1;d[i*j]=d[i]/(num[i]+1)*(num[i*j]+1);break; } 20 } 21 } 22 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; 23 return !printf("%u\n",ans); 24 }
【例题】http://59.61.75.5:8018/problem/2304
【题解】
设 $f(n)=\sum\limits_{i=1}^n \lfloor \frac{n}{i} \rfloor=\sum\limits_{i=1}^n \mathbf{d}(i)$.
$$\begin{aligned}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \mathbf{d}(ij)
&= \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \sum\limits_{a\mid i} \sum\limits_{b \mid j} [gcd(a,b)==1]
\\&= \sum\limits_{a=1}^n \sum\limits_{b=1}^m [gcd(a,b)==1] \lfloor \frac{n}{a} \rfloor \lfloor \frac{m}{b} \rfloor
\\&= \sum\limits_{d=1}^{\min(n,m)} \mu(d) (\sum\limits_{a=1}^{\lfloor \frac{n}{d} \rfloor} \lfloor \frac{n}{ad} \rfloor)(\sum\limits_{b=1}^{\lfloor \frac{m}{d} \rfloor} \lfloor \frac{m}{bd} \rfloor)
\\&= \sum\limits_{d=1}^{\min(n,m)} \mu(d)f(\lfloor \frac{n}{d} \rfloor)f(\lfloor \frac{m}{d} \rfloor)\end{aligned}$$
约数个数可以线性筛处理（同上题），数论分块即可.
【代码】
1 #include<bits/stdc++.h> 2 inline int read ( void ) 3 { 4 int x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 const int maxn=50000+10; 10 int a,b,c,mu[maxn],num[maxn]; 11 long long d[maxn];bool flag[maxn]; 12 std::vector<int> pr; 13 signed main() 14 { 15 mu[1]=d[1]=1; 16 for ( int i=2;i<=50000;i++ ) 17 { 18 if ( !flag[i] ) pr.push_back(i),num[i]=1,d[i]=2,mu[i]=-1; 19 for ( int j:pr ) 20 { 21 if ( i*j>50000 ) break; 22 flag[i*j]=true;num[i*j]=1;d[i*j]=d[i]*d[j];mu[i*j]=-mu[i]; 23 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; } 24 } 25 } 26 for ( int i=2;i<=50000;i++ ) mu[i]+=mu[i-1],d[i]+=d[i-1]; 27 for ( int T=read();T--; ) 28 { 29 int n=read(),m=read();long long ans=0; 30 for ( int i=1,j;i<=std::min(n,m);i=j+1 ) 31 { 32 j=std::min(n/(n/i),m/(m/i)); 33 ans+=1LL*d[n/i]*d[m/i]*(mu[j]-mu[i-1]); 34 } 35 printf("%lld\n",ans); 36 } 37 return 0; 38 }
Day7 20.02.17
【例题】https://www.luogu.com.cn/problem/P1829
【题解】
设 $n < m, S(n) = \sum\limits_{i=1}^{n} i =\frac{n(n+1)}{2}$.
$$\begin{aligned} \sum\limits_{i=1}^n \sum\limits_{j=1}^m \operatorname{lcm}(i,j)
&=\sum\limits_{i=1}^n \sum\limits_{j=1}^m \frac{ij}{\gcd(i,j)}
\\&=\sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{d\mid i,d\mid j} \frac{ij}{d} [\gcd(i,j)==d]
\\&=\sum\limits_{d=1}^n d\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{d} \rfloor} ij \sum\limits_{e\mid i,e\mid j} \mu(e)
\\&=\sum\limits_{d=1}^n d\sum\limits_{e=1}^{\lfloor \frac{n}{d} \rfloor} e^2\mu(e) \sum\limits_{i=1}^{\lfloor \frac{n}{de} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{de} \rfloor} ij
\\&=\sum\limits_{d=1}^n d\sum\limits_{e=1}^{\lfloor \frac{n}{d} \rfloor} e^2\mu(e) S(\lfloor \frac{n}{de} \rfloor) S(\lfloor \frac{m}{de} \rfloor)
\\&=\sum\limits_{d=1}^n d\sum\limits_{e=1}^{\lfloor \frac{n}{d} \rfloor} e^2\mu(e) S(\lfloor \frac{\lfloor \frac{n}{d} \rfloor}{e} \rfloor) S(\lfloor \frac{\lfloor \frac{m}{d} \rfloor}{e} \rfloor)
\end{aligned}$$
设 $f(n,m)=\sum\limits_{i=1}^n i^2 \mu(i) S(\lfloor \frac{n}{i} \rfloor)S(\lfloor \frac{m}{i} \rfloor)$.
$$\sum\limits_{d=1}^n d\sum\limits_{e=1}^{\lfloor \frac{n}{d} \rfloor} e^2\mu(e) S(\lfloor \frac{\lfloor \frac{n}{d} \rfloor}{e} \rfloor) S(\lfloor \frac{\lfloor \frac{m}{d} \rfloor}{e} \rfloor)=\sum\limits_{d=1}^n df(\lfloor \frac{n}{d} \rfloor,\lfloor \frac{m}{d} \rfloor)$$.
显然外层一个数论分块，内层求 $f$ 再用一个数论分块，效率 $O(n+m)$. 即可通过本题.
【代码】
1 #include<bits/stdc++.h> 2 const int mod=20101009,inv2=10050505; 3 const int maxn=10000000+10; 4 int ans,n,m,mu[maxn],sum[maxn]; 5 std::vector<int> pr; 6 bool flag[maxn]; 7 inline int S1 ( int n ) { return 1LL*n*(n+1)%mod*inv2%mod; } 8 inline int calc ( int n,int m ) 9 { 10 int res=0; 11 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; 12 return res; 13 } 14 signed main() 15 { 16 scanf("%d%d",&n,&m); 17 if ( n>m ) std::swap(n,m); 18 mu[1]=1; 19 for ( int i=2;i<=n;i++ ) 20 { 21 if ( !flag[i] ) pr.push_back(i),mu[i]=-1; 22 for ( int j:pr ) 23 { 24 if ( i*j>n ) break; 25 flag[i*j]=true;mu[i*j]=-mu[i]; 26 if ( !(i%j) ) { mu[i*j]=0;break; } 27 } 28 } 29 for ( int i=1;i<=n;i++ ) sum[i]=(sum[i-1]+1LL*i*i%mod*mu[i]%mod+mod)%mod; 30 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; 31 return !printf("%d\n",ans); 32 }
【例题】https://www.luogu.com.cn/problem/P2257
【题解】
从本篇题解起，往下均有：$n<m,\mathbf{d}(n)=\sum\limits_{d\mid n}1,S_k(n)=\sum\limits_{i=1}^n i^k$. 其中 $k=1$ 有时省略不写.
显然考虑枚举素数 $p$，考虑每个素数的贡献. 设 $pr$ 为素数集合.
$$\begin{aligned} \sum\limits_{p\in pr} \sum\limits_{i=1}^n \sum\limits_{j=1}^n [\gcd(i,j)==p]
&=\sum\limits_{p\in pr} \sum\limits_{i=1}^{\lfloor \frac{n}{p} \rfloor} \sum\limits_{i=1}^{\lfloor \frac{m}{p} \rfloor} [\gcd(i,j)==1]
\\&=\sum\limits_{p\in pr} \sum\limits_{i=1}^{\lfloor \frac{n}{p} \rfloor} \sum\limits_{i=1}^{\lfloor \frac{m}{p} \rfloor} \sum\limits_{d\mid i,d\mid j} \mu(d)
\\&=\sum\limits_{p\int pr} \sum\limits_{d=1}^{\lfloor \frac{n}{p} \rfloor} \mu(d) \lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor
\end{aligned}$$
一般遇到这类整除 $pd$ 型的式子，通常会设 $T=pd$，然后转化为枚举 $T$. 此处采用同样的方式：
$$\begin{aligned} \sum\limits_{p\int pr} \sum\limits_{d=1}^{\lfloor \frac{n}{p} \rfloor} \mu(d) \lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor
&\xlongequal[\quad]{T=pd}\sum\limits_{T=1}^n \sum\limits_{p\mid T,p\in pr} \mu(\frac{T}{p}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor
\\&=\sum\limits_{T=1}^n \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor (\sum\limits_{p\mid T,p\in pr} \mu(\frac{T}{p}))
\end{aligned}$$
后面一部分可以用埃氏筛处理，即对每个素数 $p$ 枚举其倍数，分别考虑贡献即可.
【代码】
1 #include<bits/stdc++.h> 2 inline int read ( void ) 3 { 4 int x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 const int maxn=10000000+10; 10 int mu[maxn];long long f[maxn]; 11 std::vector<int> pr; 12 bool flag[maxn]; 13 signed main() 14 { 15 mu[1]=1; 16 for ( int i=2;i<=10000000;i++ ) 17 { 18 if ( !flag[i] ) pr.push_back(i),mu[i]=-1; 19 for ( int j:pr ) 20 { 21 if ( i*j>10000000 ) break; 22 flag[i*j]=true;mu[i*j]=-mu[i]; 23 if ( !(i%j) ) { mu[i*j]=0;break; } 24 } 25 } 26 for ( int p:pr ) for ( int i=1;i*p<=10000000;i++ ) f[i*p]+=mu[i]; 27 for ( int i=2;i<=10000000;i++ ) f[i]+=f[i-1]; 28 for ( int T=read();T--; ) 29 { 30 int n=read(),m=read();long long ans=0; 31 if ( n>m ) std::swap(n,m); 32 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]); 33 printf("%lld\n",ans); 34 } 35 return 0; 36 }
【例题】http://59.61.75.5:8018/problem/1219
【题解】
显然先用二维前缀和的思想容斥. 考虑每一部分如何计算.
$$\begin{aligned} \sum\limits_{i=1}^n \sum\limits_{j=1}^m [\gcd(i,j)==k]
&= \sum\limits_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{k} \rfloor} [\gcd(i,j)==1]
\\&=\sum\limits_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{k} \rfloor} \sum\limits_{d\mid i,d\mid j} \mu(d)
\\&=\sum\limits_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d) \lfloor \frac{\lfloor \frac{n}{k} \rfloor}{d} \rfloor \lfloor \frac{\lfloor \frac{m}{k} \rfloor}{d} \rfloor
\end{aligned}$$
【代码】
1 #include<bits/stdc++.h> 2 inline int read ( void ) 3 { 4 int x=0;char ch; 5 while ( !isdigit(ch=getchar()) ) ; 6 for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48); 7 return x; 8 } 9 const int maxn=50000+10; 10 int mu[maxn]; 11 std::vector<int> pr; 12 bool flag[maxn]; 13 inline long long solve ( int n,int m,int x ) 14 { 15 n/=x;m/=x;long long res=0; 16 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); 17 return res; 18 } 19 signed main() 20 { 21 mu[1]=1; 22 for ( int i=2;i<=50000;i++ ) 23 { 24 if ( !flag[i] ) pr.push_back(i),mu[i]=-1; 25 for ( int j:pr ) 26 { 27 if ( i*j>50000 ) break; 28 flag[i*j]=true;mu[i*j]=-mu[i]; 29 if ( !(i%j) ) { mu[i*j]=0;break; } 30 } 31 } 32 for ( int i=2;i<=50000;i++ ) mu[i]+=mu[i-1]; 33 for ( int T=read();T--; ) 34 { 35 int a=read(),b=read(),c=read(),d=read(),k=read(); 36 printf("%lld\n",solve(b,d,k)-solve(b,c-1,k)-solve(a-1,d,k)+solve(a-1,c-1,k)); 37 } 38 return 0; 39 }
【例题】http://59.61.75.5:8018/problem/3748
【题解】
考虑枚举 $\gcd=x$.
$$\begin{aligned} \sum\limits_{x=1}^a \sum\limits_{i=1}^n \sum\limits_{j=1}^m [\gcd(i,j)==x] \frac{ij}{\gcd(i,j)}
&= \sum\limits_{x=1}^a x \sum\limits_{i=1}^{\lfloor \frac{n}{x} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{x} \rfloor} ij[\gcd(i,j)==1]
\\&= \sum\limits_{x=1}^a x \sum\limits_{d=1}^{\lfloor \frac{n}{x} \rfloor} d^2 \mu(d) S(\lfloor \frac{n}{xd} \rfloor) S(\lfloor \frac{m}{xd} \rfloor)
\\&\xlongequal[\quad]{T=xd} \sum\limits_{T=1}^n S(\lfloor \frac{n}{T} \rfloor) S(\lfloor \frac{m}{T} \rfloor) \sum\limits_{x\mid T,x\leq a} \frac{T^2}{x} \mu(\frac{T}{x})
\end{aligned}$$
后面那部分由于 $a$ 在变化无法处理. 因此考虑离线.
将所有询问按 $a$ 排序. 考虑如何插入每个 $x$ 的贡献. 显然式子是个前缀和类似的形式，可以枚举 $x$ 的倍数并用树状数组更新及查询.
【代码】
1 #include<bits/stdc++.h> 2 const int mod=1000000007; 3 const int maxn=100000+10; 4 struct Query { int n,m,a,id; } Q[maxn]; 5 int mu[maxn],ans[maxn],t[maxn],q; 6 bool flag[maxn]; 7 std::vector<int> pr; 8 inline void add ( int x,int y ) { for ( int i=x;i<=100000;i+=i&(-i) ) (t[i]+=y)%=mod; } 9 inline int query ( int x ) { int res=0;for ( int i=x;i;i-=i&(-i) ) (res+=t[i])%=mod;return res; } 10 inline int S1 ( int n ) { return 1LL*n*(n+1)/2%mod; } 11 signed main() 12 { 13 mu[1]=1; 14 for ( int i=2;i<=100000;i++ ) 15 { 16 if ( !flag[i] ) pr.push_back(i),mu[i]=-1; 17 for ( int j:pr ) 18 { 19 if ( i*j>100000 ) break; 20 flag[i*j]=true;mu[i*j]=-mu[i]; 21 if ( !(i%j) ) { mu[i*j]=0;break; } 22 } 23 } 24 scanf("%d",&q); 25 for ( int i=1;i<=q;i++ ) scanf("%d%d%d",&Q[i].n,&Q[i].m,&Q[i].a),Q[i].id=i; 26 std::sort(Q+1,Q+q+1,[](const Query &q1,const Query &q2){return q1.a<q2.a;}); 27 for ( int i=1,j=1;i<=q;i++ ) 28 { 29 while ( j<=Q[i].a ) 30 { 31 for ( int k=j;k<=100000;k+=j ) add(k,1LL*k*(k/j*mu[k/j]%mod+mod)%mod); 32 j++; 33 } 34 int n=std::min(Q[i].n,Q[i].m),m=std::max(Q[i].n,Q[i].m); 35 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; 36 } 37 for ( int i=1;i<=q;i++ ) printf("%d\n",ans[i]); 38 return 0; 39 }
习题
【DTOJ4729】函数
【题解】显然 $H(d)=\prod\limits_{p\text{ is prime}} g(p,d)$ 是积性函数。由二平方和定理（打表找规律）可知 $H(p^e)=3e+1\,(p\equiv 1\pmod 4)$，Min_25 筛（先咕着）即可.
【代码】
1 #include<bits/stdc++.h> 2 long long n,m,pr[100000],tot,Pow[100000],g[100000],S1[100000][5],S2[100000][5],sum1[100000],sum2[100000]; 3 bool flag[100000]; 4 inline long long f ( long long p,long long k ) { return (p%4==1) ? 3*k+1 : 1 ; } 5 inline long long Min_25 ( long long x,long long k ) 6 { 7 if ( x<=1 or pr[k]>x ) return 0; 8 long long ans=((x<=m)?sum1[x]:sum2[n/x])-Pow[k-1]; 9 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); 10 return ans; 11 } 12 signed main() 13 { 14 long long T; 15 for ( scanf("%lld",&T);T--; ) 16 { 17 scanf("%lld",&n);m=(long long)sqrt(n); 18 for ( long long i=1;i<=m;i++ ) flag[i]=true; 19 flag[1]=false;tot=0; 20 for ( long long i=2;i<=m;i++ ) 21 { 22 if ( flag[i] ) pr[++tot]=i,Pow[tot]=Pow[tot-1]+((i%4==1)?4:1); 23 for ( long long j=1;j<=tot and pr[j]*i<=m;j++ ) 24 { 25 flag[i*pr[j]]=false; 26 if ( !(i%pr[j]) ) break; 27 } 28 } 29 for ( long long i=1;i<=n;i=n/(n/(i+1)) ) 30 { 31 if ( i<=m ) for ( long long j=0;j<4;j++ ) S1[i][j]=(i-1)/4+((i-1)%4>=(j+3)%4); 32 else for ( long long j=0;j<4;j++ ) S2[n/i][j]=(i-1)/4+((i-1)%4>=(j+3)%4); 33 if ( i==n ) break; 34 } 35 for ( long long i=1;i<=tot;i++ ) 36 { 37 long long p=pr[i],p2=p*p; 38 for ( long long a=n;a>=p2;a=n/(n/a+1) ) for ( long long j=0;j<4;j++ ) 39 if ( a<=m ) 40 { 41 if ( a/p<=m ) S1[a][p*j%4]-=S1[a/p][j]-S1[p-1][j]; 42 else S1[a][p*j%4]-=S2[n/(a/p)][j]-S1[p-1][j]; 43 } 44 else 45 { 46 if ( a/p<=m ) S2[n/a][p*j%4]-=S1[a/p][j]-S1[p-1][j]; 47 else S2[n/a][p*j%4]-=S2[n/(a/p)][j]-S1[p-1][j]; 48 } 49 } 50 for ( long long i=1;i<=n;i=n/(n/(i+1)) ) 51 { 52 if ( i<=m ) S1[i][1]--; 53 else S2[n/i][1]--; 54 if ( i==n ) break; 55 } 56 for ( long long i=1;i<=n;i=n/(n/(i+1)) ) 57 { 58 if ( i<=m ) sum1[i]=4*S1[i][1]+S1[i][3]+1; 59 else sum2[n/i]=4*S2[n/i][1]+S2[n/i][3]+1; 60 if ( i==n ) break; 61 } 62 printf("%lld\n",Min_25(n,1)+1); 63 } 64 return 0; 65 }
【DTOJ3726】旧试题
【题解】
$$\begin{aligned} \sum\limits_{i=1}^A \sum\limits_{j=1}^B \sum\limits_{k=1}^C \mathbf{d}(ijk)
&= \sum\limits_{i=1}^A \sum\limits_{j=1}^B \sum\limits_{k=1}^C \sum\limits_{a\mid i} \sum\limits_{b\mid j} \sum\limits_{c\mid k}[a \perp b][b\perp c][c \perp a]
\\&= \sum\limits_{a=1}^A \sum\limits_{b=1}^B \sum\limits_{c=1}^C [a \perp b][b\perp c][c \perp a] \lfloor \frac{A}{a} \rfloor \lfloor \frac{B}{b} \rfloor \lfloor \frac{C}{c} \rfloor
\\&= \sum\limits_{a=1}^A \sum\limits_{b=1}^B \sum\limits_{c=1}^C \lfloor \frac{A}{a} \rfloor \lfloor \frac{B}{b} \rfloor \lfloor \frac{C}{c} \rfloor \sum\limits_{u\mid a,u \mid b} \mu(u) \sum\limits_{v\mid b,v \mid c}\mu(v) \sum\limits_{w\mid c,w \mid a}\mu(w)
\\&= \sum\limits_{u=1}^{\min(A,B)} \sum\limits_{v=1}^{\min(B,C)} \sum\limits_{w=1}^{\min(C,A)} \mu(u) \mu(v) \mu(w) \sum\limits_{\operatorname{lcm}(w,u)\mid a} \lfloor\frac{A}{a} \rfloor\sum\limits_{\operatorname{lcm}(u,v)\mid b} \lfloor\frac{B}{b}\rfloor \sum\limits_{\operatorname{lcm}(v,w)\mid c} \lfloor\frac{C}{c}\rfloor
\end{aligned}$$
设 $f(n,x)=\sum\limits_{n \mid d} \lfloor \frac{x}{d} \rfloor$.
$$\begin{aligned} \sum\limits_{i=1}^A \sum\limits_{j=1}^B \sum\limits_{k=1}^C \mathbf{d}(ijk)
&= \sum\limits_{u=1}^{\min(A,B)} \sum\limits_{v=1}^{\min(B,C)} \sum\limits_{w=1}^{\min(C,A)} \mu(u) \mu(v) \mu(w) \sum\limits_{\operatorname{lcm}(w,u)\mid a} \lfloor\frac{A}{a} \rfloor\sum\limits_{\operatorname{lcm}(u,v)\mid b} \lfloor\frac{B}{b} \rfloor\sum\limits_{\operatorname{lcm}(v,w)\mid c} \lfloor\frac{C}{c}\rfloor
\\&=\sum\limits_{u=1}^{\min(A,B)} \sum\limits_{v=1}^{\min(B,C)} \sum\limits_{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)\neq 0,\mu(v) \neq 0,\mu(w) \neq 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)\neq 0,\mu(v) \neq 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)$ 三元环计数+卡常即可.
【代码】
1 #include<bits/stdc++.h> 2 const int mod=1000000007; 3 const int maxn=100000+10; 4 const int maxm=1000000+10; 5 int pr[maxn],mu[maxn],a[4],d[maxn],tot; 6 long long f[4][maxn]; 7 bool flag[maxn]; 8 std::vector<int> F[maxn]; 9 typedef struct { int v,w; } edge; 10 typedef struct { int u,v,w; } Edge; 11 std::vector<edge> G[maxn]; 12 std::vector<Edge> E; 13 inline void solve() 14 { 15 scanf("%d%d%d",&a[1],&a[2],&a[3]);std::sort(a+1,a+4); 16 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; 17 for ( int i=1;i<=a[3];i++ ) if ( mu[i] ) 18 { 19 int n=(int)F[i].size(),U=(1<<n)-1; 20 for ( int S=0;S<=U;S++ ) 21 { 22 int u=1; 23 for ( int j=0;j<n;j++ ) if ( S&(1<<j) ) u*=F[i][j]; 24 int C=U^S; 25 for ( int T=C;;T=(T+1)|C ) 26 { 27 int v=1; 28 for ( int j=0;j<n;j++ ) if ( T&(1<<j) ) v*=F[i][j]; 29 d[u]++;d[v]++; 30 if ( u<v ) E.push_back((Edge){u,v,i}); 31 if ( T==U ) break; 32 } 33 } 34 } 35 long long ans=0; 36 for ( Edge e:E ) 37 { 38 if ( d[e.u]<d[e.v] ) std::swap(e.u,e.v); 39 G[e.u].push_back((edge){e.v,e.w}); 40 } 41 for ( int u=1;u<=a[3];u++ ) 42 { 43 for ( edge i:G[u] ) for ( edge j:G[i.v] ) 44 { 45 int v=i.v,w=j.v,A=i.w,B=j.w,C=1LL*w*u/std::__gcd(w,u); 46 if ( C>a[3] ) continue; 47 (ans+=mu[u]*mu[v]*mu[w]*(f[2][B]*f[3][C]+f[2][C]*f[3][B])%mod*f[1][A])%=mod; 48 (ans+=mu[u]*mu[v]*mu[w]*(f[2][C]*f[3][A]+f[2][A]*f[3][C])%mod*f[1][B])%=mod; 49 (ans+=mu[u]*mu[v]*mu[w]*(f[2][A]*f[3][B]+f[2][B]*f[3][A])%mod*f[1][C])%=mod; 50 (ans+=mod)%=mod; 51 } 52 } 53 for ( Edge e:E ) 54 { 55 int x=e.u,y=e.v,z=e.w; 56 (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; 57 (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; 58 (ans+=mod)%=mod; 59 } 60 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; 61 printf("%lld\n",ans);E.clear(); 62 for ( int i=1;i<=a[3];i++ ) G[i].clear(),d[i]=0; 63 for ( int t=1;t<=3;t++ ) for ( int i=1;i<=a[t];i++ ) f[t][i]=0; 64 } 65 signed main() 66 { 67 mu[1]=1; 68 for ( int i=2;i<=100000;i++ ) 69 { 70 if ( !flag[i] ) pr[++tot]=i,mu[i]=-1; 71 for ( int j=1;j<=tot and pr[j]*i<=100000;j++ ) 72 { 73 flag[i*pr[j]]=true;mu[i*pr[j]]=-mu[i]; 74 if ( !(i%pr[j]) ) { mu[i*pr[j]]=0;break; } 75 } 76 } 77 for ( int j=1;j<=tot;j++ ) for ( int i=1;i*pr[j]<=100000;i++ ) F[i*pr[j]].push_back(pr[j]); 78 int T; 79 for ( scanf("%d",&T);T--; ) solve(); 80 return 0; 81 }
