Note_3.31

2019/4/1 奇奇怪怪的笔记

整理了一些之前没有写过的东西，把它们拼在一起，并没有什么逻辑可言qwq

FWT快速沃尔什变换

[FWT(A)=merge(FWT(A0),FWT(A0+A1)) ]

[FWT(A)=merge(FWT(A0+A1),FWT(A1)) ]

[FWT(A)=merge(FWT(A0+A1),FWT(Ao-A1)) ]

void FWT(int *a,const int n,int type)
{
    for(int i=1;i<n;i<<=1)for(int j=0;j<n;j+=i<<1)for(k=0;k<len;++k)
    {
        int x=a[j+k],y=a[i+j+k];
        a[j+k]=x;a[i+j+k]=(type*x+y)%P;
    }
}

杜教筛

问题：求(S(n)=sum_{i=1}^{n}f(i))

考虑：

构造函数(g(x))

[egin{equation} egin{split} &sum_{i=1}^{n}(f*g)(i)\ =&sum_{d=1}^{n}sum_{i=1}^{frac{n}{d}}g(d)f(i)\ =&sum_{d=1}^{n}g(d)S(frac{n}{d})\ =&g(1)S(n)+sum_{d=2}^{n}g(d)S(frac{n}{d}) end{split} end{equation} ]

所以：

[S(n)=frac{sum_{i=1}^{n}(f*g)(i)-sum_{d=2}^{n}g(d)S(frac{n}{d})}{g(1)} ]

要找到一个使得(g)和(f*g)的前缀和都很好算的。

实现：将(S(n))存入数组(⌊frac{n}{x}⌋)位，防止冲突

为什么不会冲突呢，考虑(frac{n}{i})的(O(sqrt n))个的商，这些商的商，必然属于(O(sqrt n))个的商中

为什么呢？

假设：(n=kd+r,rleq d\k=le+s,sleq e)

那么就有：(n=(ed)l+sd+r,sd+rleq ed)

所以我们其实只访问了(O(sqrt n))个的(S(i))，把它们都存储在(⌊frac{n}{i}⌋)

这个值其实是满足([frac{n}{i}]=x)的最大(x)值，显然是不会冲突的

[mu* 1 = epsilon\ varphi* 1 = id_1\id_1 * mu = varphi ]

[∑_{i=1}^{n}id_2(i)=frac{n(n+1)(2n+1)}{6} ]

[sum id_3(i)=(sum id_1 i)^2 ]

typedef pair<ll, int> pli;
const int B = 1664510, S = 1291;
int n, T;
ll SumPhi[B + 5], ResPhi[S + 5];
int SumMu[B + 5], ResMu[S + 5]; 
namespace Pac
{
    void init()//初始化：线筛求出前B=n^{3/2}
    {
        static int cnt, vis[B + 5], prime[B + 5], phi[B + 5], mu[B + 5];
        phi[1] = mu[1] = 1;
        for (int i = 2; i <= B; i++)
        {
            if (!vis[i])
            {
                prime[cnt++] = i;
                phi[i] = i - 1;
                mu[i] = -1;
            }
            for (int j = 0, t; j < cnt && (t = i * prime[j]) <= B; j++)
            {
                vis[t] = 1;
                if (i % prime[j] == 0)
                {
                    phi[t] = phi[i] * prime[j];
                    mu[t] = 0;
                    break;
                }
                phi[t] = phi[i] * (prime[j] - 1);
                mu[t] = -mu[i];
            }
        }
        for (int i = 1; i <= B; i++)
        {
            SumPhi[i] = SumPhi[i - 1] + phi[i];
            SumMu[i] = SumMu[i - 1] + mu[i];
        }
    }

    bool vis[S + 5];

    pli GetSum(const int x) //求和部分
    {
        if (x <= B) return pli(SumPhi[x], SumMu[x]);
        int t = n / x;
        if (!vis[t])
        {
            vis[t] = 1;
            ResPhi[t] = ((ll)x + 1) * x / 2;//phi*1=id_1
            ResMu[t] = 1;//mu*1=ep
            for (uint l = 2, r; l <= x; l = r + 1)
            {
                pli res = GetSum(x / l);
                r = x / (x / l);
                ResPhi[t] -= (r - l + 1) * res.first;
                ResMu[t] -= (r - l + 1) * res.second;
            }
        }
        return pli(ResPhi[t], ResMu[t]);
    }

    void work()
    {
        init();
        scanf("%d", &T);
        while (T--)
        {
            scanf("%d", &n);
            memset(vis + 1, 0, sizeof(bool) * n / B);
            pli ans = GetSum(n);
            printf("%lld %d
", ans.first, ans.second);
        }
    }
}

最小圆覆盖[随机增量法]

随机大法好

只需要确定三个点，答案即为它们的外接圆

每次加点，如果不属于之前的圆，再枚举另外两个点，确定圆即可

看上去是(O(n^3))的，实际期望是(O(n))

#include<bits/stdc++.h>
#define ll long long
#define ld long double
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)>(b)?(b):(a))
#define reg register
using namespace std;
const int MN=100005;
const ld eps=1e-12;
int N;ld R;
struct P
{
	ld x,y;
	P(ld _x=0,ld _y=0):x(_x),y(_y){}
}a[MN],O;
inline P operator + (const P&o,const P&oo){return P(o.x+oo.x,o.y+oo.y);}
inline P operator - (const P&o,const P&oo){return P(o.x-oo.x,o.y-oo.y);}
inline P operator * (const P&o,const ld k){return P(o.x*k,o.y*k);}
inline ld dot(const P&o,const P&oo){return o.x*oo.x+o.y*oo.y;}
inline ld len(const P&o){return sqrt(dot(o,o));}
inline ld dis(const P&o,const P&oo){return sqrt((o.x-oo.x)*(o.x-oo.x)+(o.y-oo.y)*(o.y-oo.y));}
inline ld cross(const P&o,const P&oo){return o.x*oo.y-o.y*oo.x;}
inline P Mid(const P&A,const P&B){return P((A.x+B.x)/2,(A.y+B.y)/2);}
//中垂线 
inline void Vertical(const P&A,const P&B,P&p,P&v){p=Mid(A,B);P tmp=B-A;v.x=-tmp.y,v.y=tmp.x;}
//外接圆 
/*
P=A+(CA x CD)/(CD x AB) AB
*/ 
inline void Circle(const P&A,const P&B,const P&C)
{
    P p,q,v,w;Vertical(A,B,p,v),Vertical(A,C,q,w);
    P u=p-q;
    ld t=cross(w,u)/cross(v,w);
    O=p+v*t,R=dis(O,A); 
}
int main()
{
    scanf("%d",&N);
    reg int i,j,k;
    for(i=1;i<=N;++i)scanf("%Lf%Lf",&a[i].x,&a[i].y);
    random_shuffle(a+1,a+N);
    O=a[1],R=0;
    for(i=2;i<=N;++i)if(dis(a[i],O)>R)
    {
        O=a[i],R=0;
        for(j=1;j<i;++j)if(dis(a[j],O)>R)
        {
            O=Mid(a[i],a[j]),R=dis(a[i],O);
            for(k=1;k<j;++k)if(dis(a[k],O)>R)Circle(a[i],a[j],a[k]);
        }
    }
    printf("%.10f
%.10f %.10f",(double)R,(double)O.x,(double)O.y);
    return 0;
}

MTT(任意模数NTT)

求(F(x)*G(x) mod a)

(1leq n,mleq 10^5,0leq f_i,g_ileq10^9,2leq pleq 10^9+9)

法一：

卷积后每一位不超过(10^{24})

考虑对三个模数：(469762049,998244353,1004535809)，分别作(NTT)，然后采用(CRT)合并

这(3)个质数的原根均为(3)

考虑(CRT)合并怎么做

先合并前两个数，再合并后面两个

但是要用快速乘：

inline ll mul(ll a, ll b,const ll p)
{
    static const long double eps = 1e-8;
    a = (a % p + p) % p;
    b = (b % p + p) % p;
    return ((a * b - ll((long double)a / p * b + eps) * p) % p + p) % p;
}

完整代码

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define reg register
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}
const int N=1<<20|5,mod[]={469762049,998244353,1004535809},g=3;
#define Mul(a,b,p) (1ll*(a)*(b)%(p))
inline int fpow(int x,int m,int p){int r=1;for(;m;m>>=1,x=Mul(x,x,p))if(m&1)r=Mul(r,x,p);return r;}
int pos[N],MOD,len,di,F[N],G[N];
struct Polynomial
{
	int P,invg,r[N];
	void NTT(int *a,int n,int type)
	{
		reg int i,j,p,k,w,wn,X,Y;
		for(i=0;i<n;++i) if(i<pos[i]) std::swap(a[i],a[pos[i]]);
		for(i=1;i<n;i<<=1)
		{
			wn=fpow(type>0?g:invg,(P-1)/(i<<1),P);
			for(p=i<<1,j=0;j<n;j+=p)
			for(w=1,k=0;k<i;++k,w=Mul(w,wn,P))
			{
				X=a[j+k];Y=Mul(w,a[j+k+i],P);
				a[j+k]=(X+Y)%P;a[j+k+i]=(X-Y+P)%P;
			}
		}
		reg int invn=fpow(n,P-2,P);
		if(type==-1)for(i=0;i<n;++i)a[i]=Mul(a[i],invn,P);
	}
	void combine(int *A,int *B,int n)
	{
		memcpy(F,A,sizeof(int[n]));
        memcpy(G,B,sizeof(int[n]));
        NTT(F,n,1),NTT(G,n,1);
        for(reg int i=0;i<n;++i)r[i]=Mul(F[i],G[i],P);
        NTT(r,n,-1);
	}
}poly[3];
int a[N],b[N];
inline ll mul(ll a,ll b,ll p)
{
	static const long double eps=1e-8;
	a=(a%p+p)%p;b=(b%p+p)%p;
    return ((a*b-(ll)((long double)a/p*b+eps)*p)%p+p)%p;
}
//t0=a0p1 inv(p1,p0),a1p0 inv(p0,p1) mod p0p1
//t1=(a2-t0)p0^{-1}p1^{-1}
//ans=t1p0p1+t0
void CRT(int n)
{
	ll t0,t1;
	int p0=mod[0],p1=mod[1],p2=mod[2];
	int *u=poly[0].r,*v=poly[1].r,*w=poly[2].r;
	ll p0p1=1ll*p0*p1;
	ll _p0p1=mul(fpow(p0,p2-2,p2),fpow(p1,p2-2,p2),p2);
	ll _0=mul(p1,fpow(p1,p0-2,p0),p0p1);
	ll _1=mul(p0,fpow(p0,p1-2,p1),p0p1);
	for(reg int i=0;i<n;++i)
	{
		t0=(mul(u[i],_0,p0p1)+mul(v[i],_1,p0p1))%p0p1;
		t1=mul((w[i]-t0%p2+p2)%p2,_p0p1,p2);
		printf("%d ",int((mul(t1,p0p1,MOD)+t0%MOD)%MOD));
	}
}
int main()
{
	reg int n,m,i;
	n=read()+1;m=read()+1;MOD=read();
	for(i=0;i<n;++i) a[i]=read();
	for(i=0;i<m;++i) b[i]=read();
	for(len=1;len<n+m;len<<=1,di++);
    for(i=0;i<len;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(di-1));
    for(i=0;i<3;++i)
	{
		poly[i].P=mod[i];
		poly[i].invg=fpow(g,mod[i]-2,mod[i]);
		poly[i].combine(a,b,len);
	}
    CRT(n+m-1);
    return 0;
}

法二：

设(p=sqrt p)，向上取整

把每个数表示成(X=a_xp+b_x)

那么就有(ans_i=sum_{x+y=i}a_xa_yp^2+(a_xb_y+b_xa_y)p+b_xb_y)

分别用(FFT)即可

大概就这些？

接下几天会陆续补齐多项式全家桶（其实只是为了更熟练地打(FFT) ）吧，不做一只鸽子，嗯呐

---

Blog来自PaperCloud，未经允许，请勿转载，TKS！