[WC 2015复习](四)数学

都是比较简单SB的东西，求各位去WC的神犇勿喷。

1、高斯消元

(1)[BZOJ 1013][JSOI 2008]球形空间产生器sphere

http://www.lydsy.com/JudgeOnline/problem.php?id=1013

首先我们设在n维空间中，点a的坐标为(a1,a2...an),点b的坐标为(b1,b2...bn)，球心坐标(x1,x2...xn),球的半径为r

则

(a1-x1)^2+(a2-x2)^2+...+(an-xn)^2=r^2 ...①

(b1-x1)^2+(b2-x2)^2+...+(bn-xn)^2=r^2 ...②

看上去这两个式子都是二次的，不能套用一次的高斯消元，但是当我们①-②后，神奇的一幕发生了：

(a1-b1)x1+(a2-b2)x2+...+(an-bn)xn=[(a1^2-b1^2)+(a2^2-b2^2)+...+(an^2-bn^2)]/2

由于此题有n个未知数，我们需要n个方程。

以第一个点的坐标和第2~n+1个点相配对构造n个方程就ok了。常数项就是上面式子等号右边的部分

然后高斯消元搞定

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 20

using namespace std;

double pos[MAXN],a[MAXN][MAXN],ans[MAXN]; //pos[]保存第一个点的坐标，a[][]高斯消元用,a[][n+1]保存的是常数项,ans[]保存的是最终的答案、球心坐标

void Gauss(int n) //求n元一次方程组的解
{
for(int i=1;i<=n;i++) //消掉第i行之后所有行的第i列
{
int maxr=0; //找到第i行以下的第i列绝对值最大的行号maxr
for(int j=i;j<=n;j++)
if(fabs(a[j][i])>fabs(a[maxr][i]))
maxr=j;
for(int j=1;j<=n+1;j++) //交换第i行和第maxr行
swap(a[i][j],a[maxr][j]);
for(int j=i+1;j<=n;j++) //将第i行之后所有行j的第i列消去
{
double tmp=-a[j][i]/a[i][i];
for(int k=i;k<=n+1;k++)
a[j][k]+=tmp*a[i][k];
}
}
for(int i=n;i>=1;i--) //从最后一行开始回代、更新答案
{
for(int j=n;j>i;j--) //第i行的第i+1~n列中的每一列j都要消去
a[i][n+1]-=a[i][j]*ans[j];
ans[i]=a[i][n+1]/a[i][i];
}
}

int main()
{
int n;
scanf("%d",&n);
for(int i=1;i<=n;i++) //输入第一个点的坐标
scanf("%lf",&pos[i]);
for(int i=1;i<=n;i++) //输入第2~n+1号点的坐标
{
double tmp[MAXN]; //保存第i+1号点的坐标
for(int j=1;j<=n;j++)
{
scanf("%lf",&tmp[j]);
a[i][j]=pos[j]-tmp[j];
a[i][n+1]+=pos[j]*pos[j]-tmp[j]*tmp[j];
}
a[i][n+1]/=2;
}
//高斯消元
Gauss(n);
//输出答案
for(int i=1;i<=n;i++)
printf("%.3lf%c",ans[i],i==n?'\n':' ');
return 0;
}

2、线性代数

(1)[BZOJ 1297][SCOI 2009]迷路(矩阵快速幂)

http://www.lydsy.com/JudgeOnline/problem.php?id=1297

首先考虑naive的情况：任意两点之间的边权均为1，那么在n*n的矩阵map中,map[i][j]=1表示i到j有一条边，那么用T的时间从点1到点n的方案数就是自乘T次后的map的map[1]
。

但是现在任意两点之间的边权在[1,9]范围内，那么我们可以把每个点i拆成9个点，用Node i,1~Node i,9来表示它们，Node i,1~Node i,9在矩阵map中依次连边，对于边权为w的边i->j而言，我们就在矩阵map中给Node i,w向Node j,1连边，这样在矩阵map中，要从Node i,1走到Node j,1就需要经过w条边权为1的边，最终的答案就是自乘T次后的map的map[1]
。

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <algorithm>

#define MAXN 110
#define MOD 2009

using namespace std;

struct Matrix
{
int num[MAXN][MAXN];
int n,m; //n*m大小矩阵
Matrix()
{
n=m=0;
memset(num,0,sizeof(num));
}
}map,one;

Matrix operator*(Matrix a,Matrix b)
{
Matrix c;
c.n=a.n,c.m=b.m;
for(int k=1;k<=a.m;k++)
for(int i=1;i<=c.n;i++)
for(int j=1;j<=c.m;j++)
c.num[i][j]=(c.num[i][j]+((a.num[i][k]*b.num[k][j])%MOD))%MOD;
return c;
}

Matrix fastPow(Matrix base,int pow)
{
Matrix ans=one;
while(pow)
{
if(pow&1) ans=ans*base;
base=base*base;
pow>>=1;
}
return ans;
}

int num(int x,int t) //点x的第t个时间节点
{
return (x-1)*9+t;
}

int main()
{
int n,t;
char s[MAXN];
scanf("%d%d",&n,&t);
for(int i=1;i<=9*n;i++)
one.num[i][i]=1;
map.n=map.m=one.n=one.m=n*9;
for(int i=1;i<=n;i++) //每个点的9个状态链接起来
for(int j=1;j<=8;j++)
map.num[num(i,j)][num(i,j+1)]=1;
for(int i=1;i<=n;i++)
{
scanf("%s",s+1);
for(int j=1;j<=n;j++)
if(s[j]!='0')
map.num[num(i,s[j]-'0')][num(j,1)]=1;
}
map=fastPow(map,t);
printf("%d\n",map.num[num(1,1)][num(n,1)]);
return 0;
}

3、容斥原理

(1)[BZOJ 1042][HAOI 2008]硬币购物

http://www.lydsy.com/JudgeOnline/problem.php?id=1042

下面约定c[i]=第i种硬币的面值，d[i]=第i种硬币的个数。

比较显然的思路就是用凑出面值sum的总方案数（不限制硬币个数的方案数）-凑出面值sum的不合法方案数（限制硬币个数的条件下，有些种类的硬币个数超出了限制的方案数）。

用DP求出总方案数，用f[i]表示不限制硬币个数的前提下，凑出面值i的方案数。f[i]=∑(f[j-c[k]]),1≤k≤4，为了避免重复计算方案数，可以以k为阶段，分4次统计最终的答案。

那么最终的答案=f[sum]-1号硬币超额的方案数-2号硬币超额的方案数-3号硬币超额的方案数-4号硬币超额的方案数+1号和2号都超额的方案数+...+3号和4号都超额的方案数-1号、2号、3号都超额的方案数-...-...+1、2、3、4号都超额的方案数

那么怎么求第i种硬币超额的方案数呢？不妨设当前已经用第i种硬币付了(d[i]+1)*c[i]元钱，那么还剩sum-(d[i]+1)*c[i]元钱没付，这部分剩余的钱超不超额就无所谓了，反正第i种硬币已经超额了，所以方案数为f[sum-(d[i]+1)*c[i]]

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <algorithm>

#define MAXN 100010

using namespace std;

typedef long long int LL;

LL ans,f[MAXN]; //f[i]=不考虑限制条件，得到面值i的方案数
int c[5],d[5];

void DFS(int x,int k,int sum) //在用第x种硬币，k为偶数加，k为奇数就减（也可以说k就是超额的硬币种类数），sum=还没付清的面值
{
if(sum<0) return;
if(x==5)
{
if(k&1) ans-=f[sum];
else ans+=f[sum];
return;
}
DFS(x+1,k+1,sum-(d[x]+1)*c[x]); //Case 1:第x种硬币超额了，此次多付了(d[x]+1)*c[x]元
DFS(x+1,k,sum); //Case 2:第x种硬币没超额
}

int main()
{
int T;
for(int i=1;i<=4;i++) scanf("%d",&c[i]);
scanf("%d",&T);
f[0]=1;
for(int i=1;i<=4;i++) //DP预处理
for(int j=c[i];j<MAXN;j++)
f[j]+=f[j-c[i]];
while(T--)
{
for(int i=1;i<=4;i++) scanf("%d",&d[i]); //输入每种硬币提供的个数
int x; //要买x元商品
ans=0;
scanf("%d",&x);
DFS(1,0,x);
printf("%lld\n",ans);
}
return 0;
}

4、数学期望

(1)[BZOJ 3566][SHOI/HBOI 2014]概率充电器(树上期望)

http://www.lydsy.com/JudgeOnline/problem.php?id=3566

首先我们要知道期望的一个性质：对于两个相互独立的事件A、B，P(A+B)=P(A)+P(B)-P(AB)。

而一个点被通上电有三种情况：1、它自己带电。2、它的儿子给它传上了点。3、它的父亲给它传上了电。

那么我们可以分两次DFS。

第一次DFS只统计一个点i自己带电和被它的儿子通上电的概率pNode[i]。很好计算，pNode[i]=pNode[i]+pNode[son]*p[edge]-pNode[i]*pNode[son]*p[edge]

第二次DFS要在第一次DFS得到每个点导电概率的基础上，再统计上一个点i被它的父亲通上电的概率，最终得到每个点i导电的概率f[i]。

第二次做得要麻烦一些，因为会出现重复计算的问题，第一次DFS时求了儿子对父亲的贡献，但是这次是父亲对儿子的贡献，对于每对父亲和儿子来说，不能把父亲对儿子的贡献算进儿子对父亲的贡献。

不妨设对于点u而言，它除了儿子v给它的贡献以外的导电概率为y，可以得出

pNode[v]*p[edge]+y-y*pNode[v]*p[edge]=f[u]

(1-pNode[v]*p[edge])y=f[u]-pNode[v]*p[edge]

得到y=(f[u]-pNode[v]*p[edge])/(1-pNode[v]*p[edge])

这样我们就可以得到，儿子v从父亲u获得的贡献就是y(也就是说y是除去了v对u的贡献的u的导电概率)，因此f[v]=pNode[v]+y*p[edge]-pNode[v]*y*p[edge]

最终答案就是所有点的导电概率之和，即∑f[i]

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 1001000
#define EPS (1e-6)

using namespace std;

struct edge
{
int u,v,next;
double p;
}edges[MAXN];

int head[MAXN],nCount=0;

void AddEdge(int U,int V,double P)
{
edges[++nCount].u=U;
edges[nCount].v=V;
edges[nCount].p=P;
edges[nCount].next=head[U];
head[U]=nCount;
}

double pNode[MAXN]; //pNode[i]=点i在子树中被冲上电的概率(不算父结点给它导电)
double f[MAXN]; //f[i]=点i充上电的概率
double ans=0;

bool dcmp(double a) //true表明a==0,false表明a!=0
{
return fabs(a-0)<EPS;
}

void DFS1(int u,int fa) //第一次DFS，累计上u的儿子给u传电增加的期望，从点u开始DFS，u的父亲是fa
{
for(int p=head[u];p!=-1;p=edges[p].next)
{
int v=edges[p].v;
if(v==fa) continue;
DFS1(v,u);
pNode[u]=pNode[u]+pNode[v]*edges[p].p-pNode[u]*edges[p].p*pNode[v];
}
}

void DFS2(int u,int fa) //从点u开始DFS，u的父亲是fa
{
ans+=f[u];
for(int p=head[u];p!=-1;p=edges[p].next)
{
int v=edges[p].v;
if(v==fa) continue;
double tmp=(1-pNode[v]*edges[p].p);
if(dcmp(tmp)) f[v]=1.0; //特判，pNode[v]*edges[p].p==1的情况
else
{
double y=(f[u]-pNode[v]*edges[p].p)/(1.0-pNode[v]*edges[p].p);
f[v]=pNode[v]+y*edges[p].p-pNode[v]*y*edges[p].p;
}
DFS2(v,u);
}
}

int main()
{
memset(head,-1,sizeof(head));
int n;
scanf("%d",&n);
for(int i=1;i<n;i++)
{
int u,v,p;
scanf("%d%d%d",&u,&v,&p);
AddEdge(u,v,(double)p/100.0);
AddEdge(v,u,(double)p/100.0);
}
for(int i=1;i<=n;i++)
{
scanf("%lf",&pNode[i]);
pNode[i]/=100.0;
}
DFS1(1,-1);
f[1]=pNode[1]; //根节点1只能被儿子传上电
DFS2(1,-1);
printf("%lf\n",ans);
return 0;
}

5、欧拉函数

(1)[BZOJ 2705][SDOI 2012]Longge的问题

http://www.lydsy.com/JudgeOnline/problem.php?id=2705

我们可以把问题由Σgcd(i,N)的过程中枚举i变为枚举gcd(i,N)的值，然后加权累计答案。这样枚举次数会少很多。假设gcd(i,N)=k，则gcd(i/k,N/k)=1，因此Σgcd(i,N)可以转化为Σk*gcd(i/k,N/k) , gcd(i/k,N/k)=1，也就是i/k与N/k互质。而1<=i<=N，因此i/k<=N/k，对于给定的k而言，gcd(i/k,N/k)的个数就是在[1,N/k]范围内与N/k互质的数的个数，问题转化为了求欧拉函数，也就是gcd(i/k,N/k)=fai(N/k)。故Σgcd(i,N)=Σfai(N/k)*k，k是N的约数。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

using namespace std;

typedef long long int LL;

LL sqrtn;

LL h(LL x) //欧拉函数h(x)
{
LL ans=x;
for(LL i=2;i<=sqrtn;i++)
{
if(x%i==0)
{
ans=ans/i*(i-1);
while(x%i==0) x/=i;
}
}
if(x>1) ans=ans/x*(x-1);
return ans;
}

int main()
{
LL n,ans=0;
scanf("%lld",&n);
sqrtn=sqrt(n);
for(LL i=1;i<=sqrtn;i++) //枚举n的约数i
{
if(n%i==0) //i是n的约数
{
ans+=i*h(n/i);
if(i*i<n) ans+=(n/i)*h(i);
}
}
printf("%lld\n",ans);
return 0;
}

6、组合数学

(1)[BZOJ 3505][CQOI 2014]数三角形

http://www.lydsy.com/JudgeOnline/problem.php?id=3505

很不错的一道组合数学题。

不妨设整个网格中是n点*m点大小的，

首先我们要得到3个点放入n点*m点，总共n*m点中的方案数=C(n*m,3)

然后减去3个点在同一行的方案数=n*C(m,3)

再减去3个点在同一列的方案数=m*C(n,3)

还要减去3个点斜向共线的情况，这个情况特殊一些，

我们可以枚举斜向共线的3个点中最外面两个点之间的坐标差，假设x轴坐标差为i，y轴坐标差为j，那么假如已经确定了最外面两个点的坐标，则中间的那个点放置的方案数=gcd(i,j)-1，这个可以画图找规律发现。然后假如已知i、j，最外面两个点放置的方案数为2*(n-i)*(m-j)  (注：同一方向最外面两个点放置的方案数为(n-i)*(m-j),而最外面两个点可以按照左上-右下对角线方向放置，也可以按照左下-右上对角线方向放置，因此要乘个2)，枚举每个(i,j)，减去对应的不合法方案数即可。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>

#define MAXN 1005010

using namespace std;

typedef long long int LL;

LL C[MAXN][5];
int n,m;

void GetCombination() //打表求组合数C[][]
{
C[0][0]=1;
for(int i=1;i<=n*m;i++)
for(int j=0;j<=4;j++)
{
if(i==j||j==0) C[i][j]=1;
else C[i][j]=C[i-1][j]+C[i-1][j-1];
}
}

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

int main()
{
LL ans=0;
scanf("%d%d",&n,&m);
n++,m++; //n、m变成了点的个数
GetCombination();
ans+=C[n*m][3]; //答案初始化为n*m个点里任选3个点的方案数
ans-=n*C[m][3]; //去掉行上共线的方案数
ans-=m*C
[3]; //去掉列上共线的方案数
//最后去掉斜线上共线的方案数
for(int i=1;i<n;i++) //枚举斜线上共线的最外面两点之间的x轴坐标差i、y轴坐标差j
for(int j=1;j<m;j++)
{
int tmp=gcd(i,j)-1; //tmp=中间的点的位置个数
if(tmp>=1) //中间的点有位置放
ans-=tmp*2*(n-i)*(m-j); //ans-=(固定了外边两点的位置，中间的点放置的方案数)*(已知外边两点坐标差i、j，外边两点在n*m个点中放置的方案数)
}
printf("%lld\n",ans);
return 0;
}

7、莫比乌斯函数、莫比乌斯反演

(1)[BZOJ 2440][中山市选2011]完全平方数

http://www.lydsy.com/JudgeOnline/problem.php?id=2440

题目大意：求区间[1,+∞)中第k小的不是完全平方数倍数的数。

整体思路：容斥原理、莫比乌斯函数。

不妨设集合S包含了所有不是完全平方数的数。观察到随着区间[1,n]的右端点数字的增大，集合S中的数字个数也会增多，也就是说答案具备单调性，因此很自然的想法便是二分答案。

于是问题转化为了给定n，求[1,n]中有多少数字在集合S内。我们可以用容斥原理求出这个问题的反问题的答案（求[1,n]中有多少数字不在集合S内）。

这里约定一个概念：平方质数，就是是一个质数的二次方的数。

最终答案=n-( n/4+n/9+n/25+……{平方数因数是一个平方质数})+( n/36+n/100+……{含两个平方质数的因数})-{含三个平方质数的因数} ……

但是这样做很麻烦，不妨考虑用莫比乌斯函数来解决此题。



而我们知道在区间[1,n]中是x倍数的数字个数有个，不妨用i2来表示一个完全平方数，则我们可以枚举i，答案是：



公式的解释：很显然i不会超过sqrt(n)，因为有条件i2≤n，另外实际上这个公式简化过了，题目说1不是完全平方数，那么i=1的情况实际上并不能算进去，答案实际上应该是

而i=1时，

，故可以把上式简化。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>

#define MAXN 51000
#define INF 2000000000

using namespace std;

typedef long long int LL;

int miu[MAXN];

void getMiu() //打表求莫比乌斯函数
{
for(LL i=2;i<MAXN;i++) //初始时miu[i]=i表示i是质数，否则为i最大的素因子
{
if(!miu[i]) //i是质数
{
miu[i]=i;
for(LL j=i*i;j<MAXN;j+=i) //所有的j都是质数i的倍数
miu[j]=i;
}
}
//之后miu[i]就是莫比乌斯函数，初始化μ(1)=1
miu[1]=1;
for(LL i=2;i<MAXN;i++)
{
if(i%(miu[i]*miu[i])==0) //数字i有大于1的平方素因子
miu[i]=0;
else //否则数字i无大于1的平方素因子，不妨设p为i最大的质因数，miu[i]=-miu[i/p]
miu[i]=-miu[i/miu[i]];
}
}

LL cal(LL n) //求区间[1,n]中有多少个不是完全平方数
{
LL ans=0;
for(LL i=1;i*i<=n;i++)
if(miu[i])
ans+=miu[i]*n/i/i;
return ans;
}

int main()
{
int TestCase;
getMiu();
scanf("%d",&TestCase);
while(TestCase--)
{
LL n;
scanf("%lld",&n);
LL ans,lowerBound=1,upperBound=2*n+1;
while(lowerBound<=upperBound)
{
LL mid=(lowerBound+upperBound)>>1;
if(cal(mid)>=n)
{
ans=mid;
upperBound=mid-1;
}
else lowerBound=mid+1;
}
printf("%lld\n",ans);
}
return 0;
}

(2)[BZOJ 2440][HAOI 2011]Problem b

http://www.lydsy.com/JudgeOnline/problem.php?id=2440

比较显然的是gcd(i,j)=k等价于gcd(i/k,j/k)=1。不妨设solve(n,m)表示求i属于[1,n],j属于[1,m]的互质的有序的(i,j)对数。考虑到最终要求的是有序的gcd(i,j)=k的(i,j)对数，因此最终的答案很显然为solve(b/k,d/k)-solve((a-1)/k,d/k)-solve(b/k,(c-1)/k)+solve((a-1)/k,(c-1)/k)。

那么考虑solve(n,m)的过程，不妨设n<m，显然我们可以枚举一个i，得到所有的有序对(x,y)，其中x是i的倍数，y也是i的倍数，那么这样的x的个数为n/i，y的个数为m/i。有序对(x,y)的个数为(n/i)*(m/i)。

一个很朴素的写法如下：

LL _G(int a,int b)          //朴素算法
{
int tx=min(a,b),ty=max(a,b);
LL ans = 0;
for(int i = 1; i <= tx; i++)
ans += (LL)mu[i]*(tx/i)*(ty/i);
return ans;
}

但是遗憾的是，由于题目范围很大，这种朴素做法会TLE。

注意到这里的(tx/i)和(ty/i)在连续的一段i中是相同的，那么我们可以将这一部分连续的i合并成一个块，这一块的miu[i]之和乘上(tx/i)*(ty/i)即可。

这时候发现虽然经过了分块，但是由于有对miu[i]求和的过程，问题的复杂度仍然没有降低。

但是显然这里的miu[i]之和是可以用前缀和优化，于是成功降维！

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>

#define MAXN 100050

using namespace std;

typedef long long int LL;

int miu[MAXN];

void getMiu()
{
memset(miu,0,sizeof(miu));
for(LL i=2;i<MAXN;i++)
{
if(!miu[i])
{
miu[i]=i;
for(LL j=i*i;j<MAXN;j+=i)
miu[j]=i;
}
}
//求莫比乌斯函数
miu[1]=1;
for(LL i=2;i<MAXN;i++)
{
if(i%(miu[i]*miu[i])==0)
miu[i]=0;
else
miu[i]=-miu[i/miu[i]];
}
}

int sum[MAXN];

LL solve(int n,int m) //分块优化求i属于[1,n],j属于[1,m]互质的(i,j)对数
{
LL ans=0;
if(n>m) swap(n,m);
for(int i=1,p=0;i<=n;i=p+1) //p=分块的编号
{
p=min(n/(n/i),m/(m/i));
ans+=(LL)(sum[p]-sum[i-1])*(n/i)*(m/i);
}
return ans;
}

int main()
{
getMiu(); //求莫比乌斯函数
for(int i=1;i<MAXN;i++)
sum[i]=sum[i-1]+miu[i];
int a,b,c,d,k;
int T;
scanf("%d",&T);
while(T--)
{
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
LL ans=solve(b/k,d/k)-solve((a-1)/k,d/k)-solve(b/k,(c-1)/k)+solve((a-1)/k,(c-1)/k);
printf("%lld\n",ans);
}
return 0;
}