mcfx's blog - 数论

Codeforces 809E. Surprise me!

admin — Thu, 25 May 2017 20:24:07 +0800

###Description
给一棵树，每个点点权 $$a_i$$，保证 $$a_i$$各不相同，现在随机选两个点 $$u,v$$，求 $$f(u,v)=\phi(a_u\cdot a_v)\cdot dis(u,v)$$ 的期望，$$mod\ 10^9+7$$。

链接：[http://codeforces.com/contest/809/problem/E](http://codeforces.com/contest/809/problem/E "http://codeforces.com/contest/809/problem/E")

###Solution
首先，有一个结论：$$\phi(a\cdot b)=\frac{\phi(a)\cdot\phi(b)\cdot g}{\phi(g)},g=gcd(a,b)$$
证明：考虑找一个 $$x$$，使得 $$x|b,gcd(\frac{b}{x},a)=1,gcd(x,\frac{b}{x})=1$$ 且 $$a$$ 包含 $$x$$ 中所有质因子。
那么 $$\phi(a\cdot b)=\phi(a)\cdot x\cdot\phi(\frac{b}{x})=\frac{\phi(a)\cdot\phi(b)\cdot x}{\phi(x)}$$
显然，$$x$$ 是 $$g$$ 的倍数，所以 $$\frac{x}{g}=\frac{\phi(x)}{\phi(g)}$$，所以$$\phi(a\cdot b)=\frac{\phi(a)\cdot\phi(b)\cdot g}{\phi(g)}$$

有了这个结论之后可以考虑点分，那么设 $$dis(x)$$ 表示 $$x$$ 到当前根的距离。
设 $$h(x)=\frac{x}{\phi(x)}$$，那么$$f(u,v)=\phi(a_u)\cdot\phi(a_v)\cdot h(g)\cdot(dis(u)+dis(v)),g=gcd(a_u,a_v)$$
考虑当 $$u$$ 确定时，需要知道 $$\sum\phi(a_v)\cdot h(gcd(a_u,a_v))$$ 和 $$\sum\phi(a_v)\cdot h(gcd(a_u,a_v))\cdot dis(v)$$。
设 $$g(x)=h(x)-\sum[y|x,y\neq x]g(y)$$，那么 $$h(gcd(a_u,a_v))=\sum[x|a_u,x|a_v]g(x)$$
这样就可以直接枚举约数算贡献了。

###Code
```c++
#include<bits/stdc++.h>
typedef long long ll;
inline void repr(int&a,int b){if(a<b)a=b;}
#define mset(a,b) memset(a,b,sizeof(a))
#define fo0(i,n) for(int i=0,i##end=n;i<i##end;i++)
#define fo1(i,n) for(int i=1,i##end=n;i<=i##end;i++)
#define fo(i,a,b) for(int i=a,i##end=b;i<=i##end;i++)
#define foe(i,x)for(__typeof(x.end())i=x.begin();i!=x.end();++i)

const int N=200007,P=1000000007;

struct edge
{
	int to;edge*ne;
}_e[N*2],*e=_e,*p[N];

inline void add(int a,int b)
{
	*e=(edge){b,p[a]};p[a]=e++;
}

std::vector<int>v[N];
int n,pr[N],pm,phi[N],inv[N],f[N],w[N],sz[N],rsz,nro,nsz,g[N],h[N],q[N],qe,ans;
bool vis[N];

inline void dfs1(int x,int fa)
{
	sz[x]=1;
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
		{
			dfs1(i->to,x);
			sz[x]+=sz[i->to];
		}
}

inline void dfs2(int x,int fa)
{
	int t=rsz-sz[x];
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
		{
			repr(t,sz[i->to]);
			dfs2(i->to,x);
		}
	if(t<nsz)nsz=t,nro=x;
}

inline void set(int x,int y,int z)
{
	if(!~g[x])q[qe++]=x,g[x]=h[x]=0;
	g[x]=(g[x]+(ll)y*f[x])%P;
	h[x]=(h[x]+(ll)z*f[x])%P;
}

inline void get(int x,int&y,int&z)
{
	if(~g[x])y=g[x],z=h[x];
	else y=z=0;
}

inline void dfs3(int x,int fa,int dis)
{
	int a=phi[w[x]],b=(ll)a*dis%P;
	foe(i,v[w[x]])set(*i,a,b);
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
			dfs3(i->to,x,dis+1);
}

inline void dfs4(int x,int fa,int dis)
{
	int a=phi[w[x]],b=(ll)a*dis%P,A,B;
	foe(i,v[w[x]])
	{
		get(*i,A,B);
		ans=(ans+(ll)a*B+(ll)b*A)%P;
	}
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
			dfs4(i->to,x,dis+1);
}

inline void work(int x)
{
	dfs1(x,0);
	rsz=nsz=sz[x],nro=x;
	dfs2(x,0);
	vis[x=nro]=1;
	foe(i,v[w[x]])set(*i,phi[w[x]],0);
	for(edge*i=p[x];i;i=i->ne)
		if(!vis[i->to])
		{
			dfs4(i->to,x,1);
			dfs3(i->to,x,1);
		}
	fo0(i,qe)g[q[i]]=h[q[i]]=-1;
	qe=0;
	for(edge*i=p[x];i;i=i->ne)
		if(!vis[i->to])work(i->to);
}

int main()
{
	fo1(i,N-1)for(int j=i;j<N;j+=i)v[j].push_back(i);
	fo(i,2,N-1)
	{
		if(!vis[i])pr[pm++]=i,phi[i]=i-1;
		for(int j=0;i*pr[j]<N;j++)
		{
			vis[i*pr[j]]=1;
			if(i%pr[j]==0)
			{
				phi[i*pr[j]]=phi[i]*pr[j];
				break;
			}
			phi[i*pr[j]]=phi[i]*(pr[j]-1);
		}
	}
	phi[1]=1;
	inv[1]=1;
	fo(i,2,N-1)inv[i]=(ll)(P-P/i)*inv[P%i]%P;
	fo1(i,N-1)
	{
		f[i]=(ll)i*inv[phi[i]]%P;
		foe(j,v[i])if(i^*j)(f[i]+=P-f[*j])%=P;
	}
	scanf("%d",&n);
	fo1(i,n)scanf("%d",w+i);
	fo0(i,n-1)
	{
		int x,y;
		scanf("%d%d",&x,&y);
		add(x,y);
		add(y,x);
	}
	mset(vis,0);
	mset(g,0xff);
	mset(h,0xff);
	work(1);
	ans=(ll)ans*inv[n]%P*inv[n-1]%P*2%P;
	printf("%d",ans);
}
```

]]>

BZOJ 2694 & 4659 Lcm

admin — Tue, 28 Feb 2017 00:16:00 +0800

###Description
![fa(1).jpg][1]

###Solution
$$\sum_i^n\sum_j^m[\mu(gcd(i,j))\neq 0]lcm(i,j)$$
$$=\sum_k^n|\mu(k)|\sum_i^{\frac{n}{k}}\sum_j^{\frac{n}{k}}[gcd(i,j)=1]i\cdot j\cdot k$$
$$=\sum_k^n|\mu(k)|\cdot k\sum_d^{\frac{n}{k}}\mu(d)\cdot d^2\sum_i^{\frac{n}{k\cdot d}}\sum_j^{\frac{m}{k\cdot d}}i\cdot j$$
$$=\sum_k^n\sum_d[d|k]\mu(d)\cdot d^2\cdot|\mu(\frac{k}{d})|\cdot\frac{p}{d}\cdot\frac{\frac{n}{k}\cdot(\frac{n}{k}+1)}{2}\cdot\frac{\frac{m}{k}\cdot(\frac{m}{k}+1)}{2}$$
最后的式子，前面是积性函数，后面根号分块

###Code
```c++
#define _GLIBCXX_IOSTREAM
#include<bits/stdc++.h>

#define N 4000001

int pri[N],pm,f[N],ps[N],pc[N];
bool np[N];

int main()
{
	for(int i=2;i<N;i++)
	{
		if(!np[i])pri[pm++]=i,f[i]=1,ps[i]=i,pc[i]=1;
		int rf;
		if(pc[i]==1)rf=f[i]*(ps[i]-ps[i]*ps[i]);
		else if(pc[i]==2)rf=-f[i]*ps[i]*ps[i]*ps[i];
		else rf=0;
		for(int j=0;i*pri[j]<N;j++)
		{
			np[i*pri[j]]=1;
			if(i%pri[j])f[i*pri[j]]=rf,ps[i*pri[j]]=pri[j],pc[i*pri[j]]=1;
			else
			{
				f[i*pri[j]]=f[i],ps[i*pri[j]]=pri[j],pc[i*pri[j]]=pc[i]+1;
				break;
			}
		}
		f[i]=rf;
	}
	f[1]=1;
	for(int i=1;i<N;i++)
		f[i]+=f[i-1];
	int T;
	for(scanf("%d",&T);T--;)
	{
		int n,m,ans=0;
		scanf("%d%d",&n,&m);
		for(int i=1,t,p,q;i<=n&&i<=m;i=t+1)
		{
			p=n/i,q=m/i,t=std::min(n/p,m/q);
			ans+=(f[t]-f[i-1])*p*(p+1)*q*(q+1)>>2;
		}
		printf("%d\n",ans&0x3fffffff);
	}
}
```

[1]: /usr/uploads/2017/02/3245901942.jpg
  [2]: /usr/uploads/2017/02/1279320617.png

]]>

BZOJ 3552: 最右非零的数

admin — Fri, 16 Dec 2016 23:38:00 +0800

]]>

BZOJ 3113: Toy

admin — Tue, 29 Nov 2016 23:24:00 +0800

###Description
外面有一圈N个结点，中心有一个结点与N个结点都相连，总共就是2*N条边，删除N条边，使N+1个点连通，旋转相同视为等价，问有多少种情况。
![1.jpg][1]

###Input
输入N，M
3<=N<=10^9, 2<=M<=10^9

###Output
输出方案数 Mod M的结果

###Sample Input
3 10000
4 10000
4 10

###Sample Output
6
13
3

###Solution
考虑burnside引理，假设分成d块，那么每块的方案数是和bzoj1002一样的，可以用公式 f(n+1) = 3\*f(n) - f(n-1) + 2 计算
所以得到这个公式：1/n\*Sum_{d divides n} phi(n/d)*f(d)
于是直接根号n枚举d，然后根号d求phi，矩阵快速幂求f

###Code
```c++
#include<bits/stdc++.h>

typedef long long ll;
typedef long double ldb;

ll mod;

inline ll fm(ll x,ll y)
{
	ll tmp=(x*y-(ll)((ldb)x/mod*y+1e-8)*mod);
	return tmp<0?tmp+mod:tmp;
}

inline void mul(ll (&a)[3][3],ll (&b)[3][3],ll (&c)[3][3])
{
	ll t[3][3]={{0,0,0},{0,0,0},{0,0,0}};
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++)
			for(int k=0;k<3;k++)
				t[i][k]+=fm(a[i][j],b[j][k]);
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++)
			c[i][j]=t[i][j]%mod;
}

inline ll f(int x)
{
	if(x<=1)return x;
	x--;
	ll a[3][3]={{1,0,2},{0,0,mod-1},{0,1,3}},b[3][3]={{1,0,0},{0,1,0},{0,0,1}};
	for(int i=1;i<=x;i<<=1)
	{
		if(i&x)mul(b,a,b);
		mul(a,a,a);
	}
	ll ret=b[0][2]+b[2][2];
	if(ret>=mod)return ret-mod;return ret;
}

inline int phi(int x)
{
	int t=x;
	for(int i=2;i*i<=x;i++)
		if(x%i==0)
		{
			t/=i,t*=i-1;
			x/=i;
			while(x%i==0)x/=i;
		}
	if(x>1)t/=x,t*=x-1;
	return t;
}

int main()
{
	int n;
	while(~scanf("%d%lld",&n,&mod))
	{
		mod*=n;
		ll a2=0;
		for(int i=1;i*i<=n;i++)
			if(n%i==0)
			{
				if(i*i==n)
				{
					a2+=fm(phi(i),f(i));
				}
				else
				{
					a2+=fm(phi(i),f(n/i))+fm(phi(n/i),f(i));
				}
				a2%=mod;
			}
		printf("%lld\n",a2/n);
	}
}
```

[1]: /usr/uploads/2016/11/709092212.jpg

]]>

BZOJ 3884: 上帝与集合的正确用法

admin — Sun, 16 Oct 2016 00:59:00 +0800

###Description
根据一些书上的记载，上帝的一次失败的创世经历是这样的：
第一天，	上帝创造了一个世界的基本元素，称做“元”。
第二天，	上帝创造了一个新的元素，称作“α”。“α”被定义为“元”构成的集合。容易发现，一共有两种不同的“α”。
第三天，	上帝又创造了一个新的元素，称作“β”。“β”被定义为“α”构成的集合。容易发现，一共有四种不同的“β”。
第四天，	上帝创造了新的元素“γ”，“γ”被定义为“β”的集合。显然，一共会有16种不同的“γ”。
如果按照这样下去，上帝创造的第四种元素将会有65536种，第五种元素将会有2^65536种。这将会是一个天文数字。
然而，上帝并没有预料到元素种类数的增长是如此的迅速。他想要让世界的元素丰富起来，因此，日复一日，年复一年，他重复地创造着新的元素……
然而不久，当上帝创造出最后一种元素“θ”时，他发现这世界的元素实在是太多了，以致于世界的容量不足，无法承受。因此在这一天，上帝毁灭了世界。
至今，上帝仍记得那次失败的创世经历，现在他想问问你，他最后一次创造的元素“θ”一共有多少种？
上帝觉得这个数字可能过于巨大而无法表示出来，因此你只需要回答这个数对p取模后的值即可。
你可以认为上帝从“α”到“θ”一共创造了10^9次元素，或10^18次，或者干脆∞次。
一句话题意：
![1.png][1]

###Input
接下来T行，每行一个正整数p，代表你需要取模的值

###Output
T行，每行一个正整数，为答案对p取模后的值

###Sample Input
3
2
3
6

###Sample Output
0
1
4

###HINT
对于100%的数据，T<=1000,p<=10^7

###Solution
首先考虑欧拉定理，如果把这个数设为S，可以发现$$2^S=2^{S\ mod\ \phi(p)}$$ （当p是2的倍数时可以先除去所有2处理）
那么问题就转化为了求$$S\ mod\ \phi(p)$$，同样考虑，最终$$\phi(p)$$会等于1，直接处理即可
由于除了第一次，每次都会除去至少一个2，所以单次时间复杂度为$$O(logp)$$

###Code
```c++
#include<cstdio>
 
#define pmax 10000000
int prime[1000000],pm,fi[pmax+1],ans[pmax+1];
bool np[pmax+1],ga[pmax+1];
 
inline int pow(int a,long long t,int mod)
{
    long long p=1,q=a,f=1;
    for(;f<=t;f<<=1)
    {
        if(f&t)p=p*q%mod;
        q=q*q%mod;
    }
    return p;
}
 
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;
}
 
inline int gen(int x,int t)
{
    if(x==1)return 0;
    int a,b,p;
    if(x>t)exgcd(x,t,a,b);else exgcd(t,x,b,a);
    p=(long long)ans[x]*t%(x*t)*b%(x*t);
    if(p<0)p+=x*t;
    return p;
}
 
int getans(int x)
{
    int t=1;
    while(~x&1)x>>=1,t<<=1;
    if(!ga[x])ga[x]=1,ans[x]=pow(2,1000ll*fi[x]+getans(fi[x]),x);
    return gen(x,t);
}
 
int main()
{
    for(int i=2;i<=pmax;i++)
    {
        if(!np[i])prime[pm++]=i,fi[i]=i-1;
        for(int j=0;j<pm&&i*prime[j]<=pmax;j++)
        {
            np[i*prime[j]]=1;
            if(i%prime[j])
            {
                fi[i*prime[j]]=fi[i]*(prime[j]-1);
            }
            else
            {
                fi[i*prime[j]]=fi[i]*prime[j];
                break;
            }
        }
    }
    ga[1]=1,ans[1]=0;
    int T,n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        printf("%d\n",getans(n));
    }
}
```

[1]: /usr/uploads/2016/10/774614809.png

]]>

BZOJ1041: [HAOI2008]圆上的整点

admin — Sun, 18 Sep 2016 23:41:00 +0800

###Description
求一个给定的圆(x^2+y^2=r^2)，在圆周上有多少个点的坐标是整数。

###Input
只有一个正整数n,n<=2000 000 000

###Output
整点个数

###Sample Input
4

###Sample Output
4

###Solution
只考虑x,y均大于0的情况
首先，暴力枚举肯定是O(n)的，T
对于勾股数有一个结论，就是如果$$x^2+y^2=n^2$$，那么$$x=t(a^2-b^2),y=2tab,n=t(a^2+b^2)$$
那么我们可以枚举t，显然，只需要枚举n的每个质因数是否在t中就可以了（t中不需要有重复质因子）
枚举t后，需要将n/t写成a^2+b^2的形式，枚举a从1到sqrt(n/t)，时间复杂度O(sqrt(n/t))
注意去掉重复的a,b，用map实现即可
总时间复杂度应该是根号乘log吧。。
<del>这显然不是正解，不过跑的还是挺快的</del>

###Code
```c++
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<map>

std::map<int,bool> f;
  
#define pmax 65536
int prime[10000],pm,x[10],y[10],cc,ans=0,n;
bool np[pmax+1];

inline int solve(int nn)
{
	int ans=0;
	for(int i=sqrt(nn);i>0;i--)
	{
		int a=nn-i*i,b=sqrt(nn-i*i);
		if(b*b==a)
		{
			int c=2*b*i,d=std::abs(b*b-i*i);
			if(!(c&&d))continue;
			if(!(f[c*(n/nn)]||f[d*(n/nn)]))
			{
				f[c*(n/nn)]=f[d*(n/nn)]=1;
				ans+=8;
			}
		}
	}
	return ans;
}
  
void dfs(int id,int nn)
{
    if(id==cc)
    {
        ans+=solve(n/nn);
        return;
    }
    dfs(id+1,nn);
    dfs(id+1,nn*x[id]);
}

int main()
{
    for(int i=2;i<=pmax;i++)
    {
        if(!np[i])prime[pm++]=i;
        for(int j=0;j<pm&&i*prime[j]<=pmax;j++)
        {
            np[i*prime[j]]=1;
            if(i%prime[j]==0)break;
        }
    }
    int m;
    scanf("%d",&n);
    if(n==0)
    {
    	printf("1\n");
    	return 0;
	}
    m=n;
    for(int i=0;i<pm;i++)
        if(n%prime[i]==0)
        {
            x[cc]=prime[i];
            while(n%prime[i]==0)n/=prime[i],y[cc]++;
            cc++;
        }
    if(n>1)x[cc++]=n;
    n=m;
    dfs(0,1);
    printf("%d\n",ans+4);
    return 0;
}
```

]]>

BZOJ 1407: [Noi2002]Savage

admin — Wed, 14 Sep 2016 19:14:00 +0800

###Description
![1407.jpg][1]

###Input
第1行为一个整数N(1<=N<=15)，即野人的数目。
第2行到第N+1每行为三个整数Ci, Pi, Li表示每个野人所住的初始洞穴编号，每年走过的洞穴数及寿命值。
(1<=Ci,Pi<=100, 0<=Li<=10^6)

###Output
仅包含一个数M，即最少可能的山洞数。输入数据保证有解，且M不大于10^6。

###Sample Input
3
1 3 4
2 7 3
3 2 1

###Sample Output
6

###Solution
枚举m，每次枚举两个野人，将他们的ci,pi相减，记为c,p，那么如果有一个k在他们的年龄范围内，且c+pk模m等于0，则这个m不能取，k的最小值可以用扩展欧几里得求出
预处理每两个野人的ci,pi之差，li的最小值，可以优化时间

###Code
```c++
#include<cstdio>

int n,c[15],p[15],l[15],c2[105],p2[105],l2[105],n2=0;

inline int min(int a,int b)
{
	if(a<b)return a;else return b;
}

inline int max(int a,int b)
{
	if(a>b)return a;else return b;
}

int sol,gcd;

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

inline int safeexgcd(int a,int b,int &x,int &y)
{
	if(a<b)exgcd(b,a,y,x);else exgcd(a,b,x,y);
}

int main()
{
	int maxc=0;
	scanf("%d",&n);
	for(int i=0;i<n;i++)scanf("%d%d%d",c+i,p+i,l+i),maxc=max(c[i],maxc);
	for(int i=0;i<n;i++)
		for(int j=0;j<i;j++)
		{
			c2[n2]=c[i]-c[j];
			p2[n2]=p[i]-p[j];
			l2[n2]=min(l[i],l[j]);
			n2++;
		}
	for(int i=maxc;1;i++)
	{
		bool flag=true;
		int x,y;
		for(int j=0;j<n2;j++)
		{
			sol=-c2[j];
			safeexgcd(p2[j],i,x,y);
			if(p2[j]*x+i*y!=sol)continue;
			int tmp=i/gcd;
			x%=tmp;
			if(x<0)x+=tmp;
			if(x<=l2[j])
			{
				flag=false;
				break;
			}
		}
		if(flag)
		{
			printf("%d\n",i);
			return 0;
		}
	}
}
```
另外随便写的代码就第4了
![puh.png][2]

[1]: /usr/uploads/2016/09/2556788862.jpg
  [2]: /usr/uploads/2016/09/909201119.png

]]>