zoukankan      html  css  js  c++  java
  • 类欧几里得入土总结 2

    类欧几里得入土总结 2

    直接说 LOJ万能欧几里得这道题

    题目要我们求这个式子

    [sum^n_{x=1} A^xB^{ax+bover c} ]

    问题转换

    考虑一件这样的事情,平面上一条直线(ax+bover c),关于这条直线,我们放到网格图上,从左到右,我们遇到一条横线执行一次(U)操作,遇到一次竖线执行一次(R)操作,经过一个整点,执行UR操作。U代表迭代一次( m curB*=B),R代表迭代一次( m curA*=A)然后累加一遍({ m sum+=}( m curA) imes ({ m curB})),这相当于在模拟这个(sum)的操作,我们称(UR)构成的序列叫做操作序列。

    若此处的(<+, imes>)有两者有结合律,(+)要对$ imes $有分配律,x,y是操作序列,xy表示直接首尾相接起来。我们容易证明以下等式

    [cases{ m curA(xy)=curA(x) imes curA(y) \ m curB(xy)=curB(x) imes curB(y) \ m sum(xy)=sum(x)+curA(x) imes sum(y) imes curB(x) } ]

    也就是说那么操作序列具有结合律

    记五元组((n,a,b,c,x,y)=s)表示我要生成一个操作序列(s),规则是对于直线(f(x)=y=lfloor {ax+bover c} floor,xin [1,n]) ,放在网格上,从左下往右上走,每经过一个竖线添加一段操作数列(y),每经过一个横线添加一段操作序列(x)(同时经过先进行(x)操作)。

    注意到这个(s)总共有(n)(y),第(i)(y)之前有({ai+bover c})(x)

    沿用【瞎讲】类欧几里得入土教程的思路,我们考虑将(a,b)都化到小于(c)的情况。

    一些铺垫


    重要补充

    (g(i))表示在第(i)(y)之前有多少个(x),我们有

    [g(i)=egin{cases}f(1)&i=1\f(i)-f(i-1) & m otherwiseend{cases} ]

    定义(f)的和定义和(g)等价,描述的是同一个情况,为了方便理解建议用差分的角度理解(i到i-1个y之间有多少个x)


    先考虑(bge c)的情况,我们要使得$b o b%c $

    (i>1)(y)之前有(f(i)-f(i-1))(x),而(b o b\% c)对这个差分毫无影响。

    然而(i=1)的时候有小问题,第一个(y)之前个(x)我们不是差分定义的,而我们使得(b o b\%c),使得在第一个(y)之前少了一些(y),具体的,少了(lfloor{a imes 1+bover c} floor -lfloor {a imes 1+b\% cover c} floor=lfloor{bover c} floor)(y)

    根据以上的分析,我们得到了

    [(n,a,b,c,x,y)=underbrace{ yydots yy}_{lfloor{bover c} floor}(n,a,b\%c,c,x,y)quad bge c ]


    再考虑(age c)的情况,我们要使得(a o a\% c)

    此时(f(x)=lfloor{(a\%c )x+bover c} floor+lfloor{aover c} floor x),设(f'(x)=lfloor{(a\%c )x+bover c} floor)(g')(f')的差分,我们发现(g(i)=lfloor {aover c} floor+g'(i)),即任何一个(y)之前都紧贴着(lfloor{aover c} floor)(x)

    因此

    [(n,a,b,c,x,y)=(n,a\%c,b,c,x,underbrace{xxx...xxx}_{lfloor{aover c} floor}space y) quad age c ]


    解决问题

    继续沿用上面那个博客的思路,到了真正解决问题的时刻了

    考虑现在我们只需要解决的问题是(a,b<c)

    边界情况:

    关于(a=0)的情况退化为一条平行于(x)的直线,就是填(n)(y)

    (b=0)无所谓

    (a,b<c)意味着一件事情,这个直线非常平缓,甚至填了若干个(y)才填一个(x)

    这也意味着另一件事情,那就是一个(x)之前有若干个(y)

    考虑 第(j>1)(x)之前有多少个(y),即考虑(maxlimits_{m} lfloor {am+bover c} floor< j)。这即(m=f'(j),j>1)

    解出来

    [m=lfloor{jc-b-1over a} floorquad j>1 ]

    也就是从第二个(x)开始到第(lfloor{an+bover c} floor)(x)是一个子问题(要记得令(j=j-1)并把(c)放出来,即)

    [f''(j)=lfloor{jc+c-b-1over a} floor quad jge 1 ]

    而第一个(x)之前还有(f''(0))(y),这一部分补上。

    还有一个问题是,最后一个(x)之后可能还有(y)没有统计到,我们手动算出来补在后面即可。

    根据以上的分析,我们得到了

    [(n,a,b,c,x,y)=underbrace{yyy...dots y}_{lfloor {c-b-1over a} floor} x (lfloor{an+bover c} floor-1,c,c-b-1,a,y,x)underbrace{yyydots yyy}_{n-lfloor{c(lfloor{an+bover c} floor-1)+(c-b-1)over a} floor} quad a,b<c ]

    总结

    先把公式全写出来

    [(n,a,b,c,x,y)= egin{cases} underbrace{ yydots yy}_{lfloor{bover c} floor}(n,a,b\%c,c,x,y)& bge c \ (n,a\%c,b,c,x,underbrace{xxx...xxx}_{lfloor{aover c} floor}space y) & age c \ underbrace{yyy...dots y}_{lfloor {c-b-1over a} floor} x (lfloor{an+bover c} floor-1,c,c-b-1,a,y,x)underbrace{yyydots yyy}_{n-lfloor{c(lfloor{an+bover c} floor-1)+(c-b-1)over a} floor} & m otherwise end{cases} ]

    然而求很多操作序列的叠加要用分治,是一个(log)的,但是发现每次分治的长度是按除法递减的,仔细分析一波发现就是(log n)。(这里地方太小我写不下)

    实现时注意每次调用递归时,要注意将第一个(y)之间的(x)补上($ m case::otherwise $ 的补(y)也可以理解为这个)。

    这种写法的一大好处是,递归部分完全不需要改动,需要改动的只有操作序列的合并方式。如果操作可以对应到矩阵那就完全不要改动了。

    这种写法的另一大好处是,任何良好定义的(<op1(+),op2( imes )>),且curB,curA的"迭代"有结合律和交换律,都可以套用这个做法。(泰拳警告:求多项式,超现实数,Floyd数组,SG函数,图连通性的sum)

    如果(sum)后面是(k)个项相乘(但是有交换律)也能解决(我愿称之为k阶-类欧几里得)

    代码实现上,由于(case1)(case3)的统一性,可以每次在递归调用前直接补上$lfloor {bover c} floor $ ( (case2)除外),那么每次进入一个递归马上可以(b o b\%c)

    此外,注意到当(f(1)=0)的时候特殊处理下。

    目前LOJ rk4 不知道为啥我跑得这么快qwq

    //@winlere
    #include<iostream>
    #include<algorithm>
    #include<cstdio>
    #include<cstring>
    #define mod 998244353
    
    typedef long long ll;
    ll qr(){
    	ll c=getchar(),ret=0,f=0;
    	while(!isdigit(c)) f|=c==45,c=getchar();
    	while( isdigit(c)) ret=ret*10+c-48,c=getchar();
    	return ret;
    }
    int MOD(const int&x){return x>=mod?x-mod:x;}
    int MOD(const int&x,const int&y){ll g=1ll*x*y;return g-g/mod*mod;}
    int N;
    
    struct EU{
    	struct MAT{
    		int v[20][20];
    		MAT(){memset(v,0,sizeof(int)*400);}
    		int*operator[](int x){return v[x];}
    		MAT operator + (const MAT&x)const{
    			MAT ret;
    			for(int t=0;t<N;++t)
    				for(int i=0;i<N;++i)
    					ret[t][i]=MOD(v[t][i]+x.v[t][i]);
    			return ret;
    		}
    		MAT operator * (const MAT&x)const{
    			MAT ret;
    			for(int k=0;k<N;++k)
    				for(int t=0;t<N;++t)
    					for(int i=0;i<N;++i)
    						ret[t][i]=MOD(ret[t][i]+MOD(v[t][k],x.v[k][i]));
    			return ret;
    		}
    	}sum,A,B;
    	EU(const int&x=1):sum(),A(),B(){for(int t=0;t<N;++t) A[t][t]=B[t][t]=x;}
    	EU operator + (const EU&x)const{
    		EU ret;
    		ret.A=A*x.A;
    		ret.B=B*x.B;
    		ret.sum=sum+A*x.sum*B;
    		return ret;
    	}
    	EU operator * (const ll&x)const{
    		EU ret,base=*this;
    		for(ll t=x;t>0;t>>=1,base=base+base)
    			if(t&1) ret=ret+base;
    		return ret;
    	}
    };
    ll f(ll x,ll a,ll b,ll c){return ((__int128)x*a+b)/c;}
    
    EU solve(ll n,ll a,ll b,ll c,const EU&x,const EU&y){
    	if(!n) return EU();
    	b%=c;
    	if(!f(n,a,b,c)) return y*n;
    	if(a>=c) return solve(n,a%c,b,c,x,x*(a/c)+y);
    	ll m=f(n,a,b,c),cnt=n-f(m-1,c,c-b-1,a);
    	return y*((c-b-1)/a)+x+solve(m-1,c,c-b-1,a,y,x)+y*cnt;
    }
    
    int main(){
    	ll a=qr(),c=qr(),b=qr(),n=qr(); N=qr();
    	EU U,R;
    	EU::MAT A,B;
    	for(int t=0;t<N;++t)
    		for(int i=0;i<N;++i)
    			A[t][i]=qr();
    	for(int t=0;t<N;++t)
    		for(int i=0;i<N;++i)
    			B[t][i]=qr();
    	U.B=B; R.A=A; R.sum=A;
    	EU::MAT ans=(U*(b/c)+solve(n,a,b,c,U,R)).sum;
    	for(int t=0;t<N;++t,putchar('
    '))
    		for(int i=0;i<N;++i)
    			printf("%d ",ans[t][i]);
    	return 0;
    }
    
    
  • 相关阅读:
    使用微软TFS代码管理工具和在金山快盘上搭建SVN的使用方法
    微软的Windows8安装体验
    软件注册码随笔
    软件注册码(算法一DES)
    PHP连接SAE平台MYSQL
    一点一滴《C++处理数据》
    BouncyCastle.Crypto的RSA算法调用源码
    一点一滴《C++学习》
    软件注册码(算法二Rijndael)
    Web 应用程序的程序常见安全防范
  • 原文地址:https://www.cnblogs.com/winlere/p/13021633.html
Copyright © 2011-2022 走看看