[bzoj1129][数论][扩展欧几里得]Per

原创于 2018-12-15 09:19:05 发布 · 497 阅读

1 ·

CC 4.0 BY-SA版权

bzoj 同时被 3 个专栏收录

482 篇文章

订阅专栏

数论

63 篇文章

订阅专栏

拓展欧几里得

7 篇文章

订阅专栏

本文介绍了一种计算序列在所有可能排列中字典序排名的算法，通过康托展开方法，结合树状数组和质因数分解，解决了序列长度可达30万的大规模问题。算法核心在于对序列进行预处理，利用树状数组维护每种元素出现次数，并通过质因数分解处理模运算。

Description

给你一个序列s，你把这个序列的所有不同排列按字典序排列后，求s的排名mod m

Input

序列的长度n<300 000,m n个数，代表序列s

Output

排名mod m

Sample Input

4 1000

2 1 10 2

Sample Output

5

题解

题挺不错
大概也想到了一些
可是就是没写出来…
容易发现，我们可以用类似康托展开的方法来做
枚举前 $i$ 个是相同的，设后面第 $i$ 大的数出现了 $v i s [i]$ 次，第 $i$ 位的数大小是 $c a l [i]$
这里的贡献是
$∑u=1cal[i+1]−1(n−i)!(vis[u]−1)!Πk!=uvis[k]!\sum_{u=1}^{cal[i+1]-1}\frac{(n-i)!}{(vis[u]-1)!\Pi_{k!=u}vis[k]!}$
上下同时乘一个 $v i s [u]!$ 可以得到
$∑u=1cal[i+1]−1(n−i)!vis[u]Πvis[k]!\sum_{u=1}^{cal[i+1]-1}\frac{(n-i)!vis[u]}{\Pi vis[k]!}$
提出来 $v i s [u]$ 可以得到
$(n−i)!Πvis[k]!∗∑u=1cal[i+1]−1vis[u]\frac{(n-i)!}{\Pi vis[k]!}*\sum_{u=1}^{cal[i+1]-1}vis[u]$
于是就可以 $i$ 从前往后扫，后面用一个树状数组，前面直接维护就行了
如果 $m$ 是一个质数这里就做完了…
但是他不是，所以你得分解一下质因数
对于每个 $p_i^{k_i}$ 做一遍上面的过程，除法的时候要把 $p_i$ 提出来单独减掉然后乘上去或者乘个逆元
最后 $e x g c d$ 合并一下就可以了

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<vector>
#include<ctime>
#include<map>
#include<bitset>
#define LL long long
#define mp(x,y) make_pair(x,y)
#define lc now<<1
#define rc now<<1|1
using namespace std;
inline int read()
{
	int f=1,x=0;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
int stack[20];
inline void write(int x)
{
	if(x<0){putchar('-');x=-x;}
    if(!x){putchar('0');return;}
    int top=0;
    while(x)stack[++top]=x%10,x/=10;
    while(top)putchar(stack[top--]+'0');
}
inline void pr1(int x){write(x);putchar(' ');}
inline void pr2(int x){write(x);putchar('\n');}
struct LSnode{int y,p;}w[310000];
bool cmp(LSnode n1,LSnode n2){return n1.y<n2.y;}
int s[310000],tt;
int lowbit(int x){return x&-x;}
void change(int x,int c){for(;x<=tt;x+=lowbit(x))s[x]+=c;}
int getsum(int x){int ret=0;for(;x>=1;x-=lowbit(x))ret+=s[x];return ret;}

int cal[310000],vis[310000],n;
LL fla[310000],MD1;
LL mod,fac,u1[310000],u2[310000],p1[310000],p2[310000];


LL ok[310000*4],oo[310000];
void buildtree(int now,int l,int r)
{
	if(l==r){ok[now]=oo[l];return ;}
	int mid=(l+r)/2;
	buildtree(lc,l,mid);buildtree(rc,mid+1,r);
	ok[now]=ok[lc]*ok[rc]%mod;
}
void modify(int now,int l,int r,int p,int c)
{
	if(l==r){ok[now]=c;return ;}
	int mid=(l+r)/2;
	if(p<=mid)modify(lc,l,mid,p,c);
	else modify(rc,mid+1,r,p,c);
	ok[now]=ok[lc]*ok[rc]%mod;
}

LL inv[310000];
LL exgcd(LL a,LL b,LL &x,LL &y)
{
	if(a==0)
	{
		x=0;y=1;
		return b;
	}
	else
	{
		LL tx,ty;
		LL d=exgcd(b%a,a,tx,ty);
		x=ty-(b/a)*tx;
		y=tx;
		return d;
	}
}
LL getinv(LL u)
{
	LL A=u,B=mod,x,y,K=1;
	LL d=exgcd(A,B,x,y);
	x=(x*(K/d)%(B/d)+(B/d))%(B/d);
	return x;
}
LL pow_mod(LL a,int b)
{
	LL ret=1;
	while(b)
	{
		if(b&1)ret=ret*a%mod;
		a=a*a%mod;b>>=1;
	}
	return ret;
}
void init()
{
	for(int i=1;i<=n;i++)
	{
		p1[i]=i;p2[i]=0;
		while(!(p1[i]%fac))p1[i]/=fac,p2[i]++;
	}
	u1[0]=1;inv[0]=getinv(1);
	u2[0]=0;
	for(int i=1;i<=n;i++)
	{
		u1[i]=u1[i-1]*p1[i]%mod;
		inv[i]=getinv(u1[i]);
		u2[i]=u2[i-1]+p2[i];
	}
	memset(s,0,sizeof(s));
	memset(vis,0,sizeof(vis));
	for(int i=1;i<=n;i++)vis[cal[i]]++;
	for(int i=1;i<=tt;i++)change(i,vis[i]);
}
//u1 阶乘去fac的乘积
//u2 阶乘有多少fac
//p1 单点去fac的数 
//p2 单点有多少fac
LL query()
{
	init();
	LL ret=0;
	LL g1=1,g2=0;
//	for(int i=1;i<=tt;i++)oo[i]=u1[vis[i]],g2=g2+u2[vis[i]];
//	buildtree(1,1,tt);
	for(int i=1;i<=tt;i++)g1=g1*inv[vis[i]]%mod,g2=g2+u2[vis[i]];
	for(int i=1;i<=n;i++)
	{
		LL s1=u1[n-i],s2=u2[n-i];
//		g1=ok[1];
		LL sum;
		if(s2-g2>=0)sum=s1*g1%mod*pow_mod(fac,s2-g2)%mod;
		else sum=s1*g1%mod*getinv(pow_mod(fac,abs(s2-g2)))%mod;
		ret=(ret+sum*getsum(cal[i]-1)%mod)%mod;
		
		g2=g2-u2[vis[cal[i]]];
		g1=g1*u1[vis[cal[i]]]%mod;vis[cal[i]]--; 
//		modify(1,1,tt,cal[i],u1[--vis[cal[i]]]);
		change(cal[i],-1);
		g2=g2+u2[vis[cal[i]]];
		g1=g1*inv[vis[cal[i]]]%mod;
	}
	return ret;
}
LL a1[310000],mm1[310000],mm2[310000],m3[310000],tot;
int main()
{
//	freopen("a.in","r",stdin);
//	freopen("a.out","w",stdout);
	n=read();MD1=read();
	for(int i=1;i<=n;i++)w[i].y=read(),w[i].p=i;
	sort(w+1,w+1+n,cmp);
	for(int i=1;i<=n;i++)
	{
		if(w[i].y!=w[i-1].y)tt++;
		cal[w[i].p]=tt;fla[tt]=w[i].y;
	}
	LL temp=MD1;
	for(int i=2;i*i<=temp;i++)
		if(!(temp%i))
		{
			tot++;
			mm1[tot]=i;mm2[tot]=1;
			while(!(temp%i))temp/=i,mm2[tot]*=i;
		}
	if(temp!=1)mm1[++tot]=temp,mm2[tot]=temp;
	for(int i=1;i<=tot;i++)
	{
		mod=mm2[i];fac=mm1[i];
		a1[i]=query();
	}
	LL m1=mm2[1],b1=a1[1],m2,b2;
	for(int i=2;i<=tot;i++)
	{
		m2=mm2[i];b2=a1[i];
		LL A=m1,B=m2,x,y,K=b2-b1;
		LL d=exgcd(A,B,x,y);
		x=(x*(K/d)%(B/d)+(B/d))%(B/d);
		b1=m1*x+b1;m1=m1/d*m2;
	}
	pr2((b1+1)%MD1);
	return 0;
}