UOJ Logo fjzzq2002的博客

博客

一个更好的多项式模板

2021-08-03 00:32:31 By fjzzq2002

一直想写个好点的多项式模板,最近终于想清楚了怎么写,于是就写了一个!代码在本文最后。


众所周知,常见的多项式模板经常有这样两个问题:速度很慢、用起来很麻烦。

速度问题是老生常谈的。NTT我们可以找一个好的板子对着抄。形式幂级数运算我们可以采用 https://negiizhao.blog.uoj.ac/blog/4671 中的优化。在线卷积我们可以使用多叉降低复杂度,循环卷积和小范围暴力降低常数。那么本文考虑的主要问题是... 如何让你的模板封装得阳间一些?

我随机选取了若干份uoj上相关题目的AC代码:

codes.png

一个字,惨。


我们先考虑形式幂级数和多项式的部分。这部分我的主要做法是,将每个多项式/形式幂级数记录精确度,即记录 $\bmod x^n$ 中的 $n$,在运算时保持该精确度不变。在牛顿迭代解 $G(f)=0$ 时,我们可以使用两个 lambda 传入 $G(f)$ 和 $G'(f)$,lambda 保持 $f$ 的精确度不变即可。我的模板中具体对应的是 shrink_len 属性,-1 表示为多项式,否则为精确次幂等于该值的集合幂级数。使用这个思路实现各运算即可。

以下各题均省略推导过程,代码均省略模板部分,模数均为 $998244353$。

loj150 挑战多项式

给定 $n$ 次多项式 $F(x)$,求 $G(x)$ 满足 $\displaystyle G(x) \equiv \left({\left({1+\ln\left({2+F(x)-F(0)-{\exp\left({\int_0^x\frac{1}{\sqrt{F(t)}}\textrm{d}t}\right)}}\right)}\right)^k}\right)^\prime \pmod {x^n}$。$n\le 10^5$

int n,m,f[123456];
int main()
{
    scanf("%d%d",&n,&m);
    default_shrink=++n;
    for(int i=0;i<n;++i) scanf("%d",f+i);
    poly p(vector<int>(f,f+n));
    p=p-gexp(gint(ginv(gsqrt(p,mod_sqrt(f[0],MOD)))));
    p=gcorner(gln(gcorner(p,{1})),{1});
    gde(gexp(mi(m)*gln(p))).print_len(n-1,'\n');
}

Quite easily done! 跑了 sum 4269ms/max 228ms,在写作本文时排在ranklist第一页。

luogu6112 直接自然溢出啥事没有 加强版 (30%)

求出多项式 $F(x)$ 使得 $\left(x^{7}+x^{6}+x^{5}-2 x^{4}+x^{2}\right) F(x)^{2}+\left(x^{5}-x^{4}-2 x^{3}+2 x^{2}+x-1\right) F(x)+\left(x^{4}-2 x^{2}+1\right)\equiv 0 \pmod {x^n}$。$n\le 10^5$

int n;
int main()
{
    cin>>n;
    poly
    p="x^7+x^6+x^5-2x^4+x^2"_p,
    q="x^5-x^4-2x^3+2x^2+x-1"_p,
    r="x^4-2x^2+1"_p,
    a=gnewton(
    [&](poly g) {
        return (g*p+q)*g+r;
    },
    [&](poly g) {
        return 2*g*p+q;
    },
    1,n+1);
    cout<<a.get(n)<<"\n";
}

C++11就是好。

uoj50 链式法则

求出多项式 $F(x)$ 使得 $F'(x) \equiv \frac{1}{2}A(x)F^2(x) + 1 \pmod {x^n}$。$n\le 200000$

模板里咋能没有一阶微分方程呢。

int n; poly a(-1);
char str[288888];
int main()
{
    scanf("%d%s",&n,str);
    a.reserve(n+1);
    for(int i=0;i<n;++i)
        if(str[i]=='1') a[i]=rfac[i];
    auto g=gnewton_d(
    [&](poly x) {
        x=mi((MOD+1)/2)*x*x*a;
        ++x[0]; return x;
    },
    [&](poly x) {
        return x*a;
    },0,n+1);
    for(int i=1;i<=n;++i) {
        mi ans=g.get(i)*fac[i];
        printf("%d\n",(int)ans);
    }
}

跑了 sum 2851ms/max 1464ms,在写作本文时(可以)排在ranklist第五页。

loj6684 有根无标号「奇树」计数

求出多项式 $F(x)$ 使得 $F(x) \equiv x\text{MSet}(\text{MSet}(F(x))-1) \pmod {x^n}$($\text{MSet}$ 的定义可参见 此处)。$n\le 100000$

这里我们进行牛顿迭代,注意到 $\text{MSet} (f(z))=\exp\left(\sum_{k\ge 1}f(z^k)/k\right)$,而倍增时 $\exp\left(\sum_{k\ge 2}f(z^k)/k\right)$ 有足够的精确度,因此我们可以将该部分视为一个常多项式,认为 $\text{MSet} (f(z))=\exp f(z)~\underbrace{\exp\left(\sum_{k\ge 2}f(z^k)/k\right)}_{\text{常多项式}}$。

由于模板中显然没有提供计算 $\sum_{k\ge 2}f(z^k)/k$ 的接口,我们需要手(kao)写(bei)。

poly MSet_(poly p,int st) {
    p.fit_shrink();
    assert(p.get(0)==0);
    poly q(p.shrink_len);
    q.fit_shrink();
    for(int i=st;i<p.size();++i) {
        mi iv=rfac[i]*fac[i-1];
        for(int j=1;i*j<p.size();++j)
            q[i*j]=q[i*j]+p[j]*iv;
    }
    return q;
}
int n;
int main()
{
    scanf("%d",&n);
    poly o=gnewton(
    [&](poly u) {
        return gshl(MSet(gcorner(MSet(u),{0})),1)-u;
    },
    [&](poly u) {
        poly P=MSet_(u,2),
            z=gcorner(gexp(P+u),{0}),
            Q=MSet_(z,2);
        //x*exp(Pexp(F)-1)Q-F=0
        return gcorner(gshl(gexp(z+u+P+Q),1),{-1});
    },
    0,n+1);
    for(int i=1;i<=n;++i)
        printf("%d\n",(int)o.get(i));
}

在经过一些简单的常数优化(注意到 $G$ 和 $G'$ 都需要求出 $\text{MSet}$,可以进行记忆化)后,跑了 max 267ms,在写作本文时(可以)排在ranklist第一(忽略一位打表的同学)。

loj6538 烷基计数 加强版 加强版

求出多项式 $F(x)$ 使得 $F(z)\equiv z \frac{F(z)^{3}+3 F(z) F\left(z^{2}\right)+2 F\left(z^{3}\right)}{6}+1\pmod {x^n}$。$n\le 100000$

类似上题进行牛顿迭代。

int n;
int main() {
    cin>>n;
    cout<<gnewton(
    [&](const poly& s) {
        poly t=s;
        t=inv(6)*(t*(t*t+3*gamp(t,2))+2*gamp(t,3));
        t.insert(0); t=t+1-s; return t;
    },
    [&](poly t) {
        t=inv(2)*(t*t+gamp(t,2));
        t.insert(0); t=t-1; return t;
    },
    1,n+1).get(n)<<"\n";
}

跑了 sum 298ms/max 98ms,在写作本文时排在ranklist第一页。

此外,模板中还实现了快速线性递推、快速多点求值插值、幂和、简单函数复合逆等算法。


2021年了,牛顿法等朴素迭代方法已经不足以满足我们。接下来我们来考虑在线卷积算法。

常见(wo'jian'guo'de)的半在线卷积算法实现方法为,传入 $i$,更新到 $f[i]$ 位置时调用 $relax(i)$,但这种方法实现扩展性较差。

于是我们可以考虑把序列真的作为一个流来处理,主要定义如下:

struct ocpoly {
    function<mi(int)> get_handle;
    vector<int> vis; vector<mi> val;
    mi get(int x) {
        if((int)vis.size()>x&&vis[x]) return val[x];
        if((int)vis.size()<=x)
            vis.resize(x+1),val.resize(x+1);
        val[x]=get_handle(x); vis[x]=1; return val[x];
    }
    poly await(int n) {
        poly s(-1);
        for(int i=0;i<n;++i)
            s[i]=get(i);
        return s;
    }
    void copy(ocpoly*b) {
        get_handle=[=](int c) {
            return b->get(c);
        };
    }
    ocpoly() {}
    ocpoly(ocpoly const&)=delete;
    ocpoly& operator = (ocpoly const&)=delete;
}

即我们对外的接口为 get 函数,该函数实现了对lambda get_handle 返回值的记忆化(我们希望 $\text{get_handle}(i)$ 返回 $x^i$ 位置上的系数)。这种写法的好处在于扩展性较强,我们只需要改变lambda就可以派生出各种流。例如我们可以定义在一个流上进行积分:

struct ocint: public ocpoly {
    ocint(ocpoly*p,mi c=0) {
        get_handle=[=](int u) {
            if(u==0) return c;
            return p->get(u-1)*rfac[u]*fac[u-1];
        };
    }
}

我们可以定义更多的组合子,例如如果我们定义 $\text{ocde}$ 为流上的求导,$\text{ocquo}(A,B)$ 为 $A(x)B(x)^{-1}$ 形成的形式幂级数,我们可以组合成流上的 $\ln$:

struct ocln: public ocpoly {
    ocln(ocpoly*a) {
        copy(new ocint(new ocquo(new ocde(a),a)));
    }
};

(在内存管理的角度上这种写法不是很优美,但是懒得用智能指针)

当然,实现最困难的部分在于全在线多项式卷积 $\text{ocmul}$。首先全在线转半在线只需要倍增最大长度,即我们对于当前考虑的位置 $i$,我们保证所有 $i=a+b$ 有 $\min(a,b)\le u$,即要求 $2(u+1)>i$,并在不满足要求时修改 $u$。这样我们就可以将在线卷积转换为 $A[0..u],B[u+1..]$ 和 $B[0..u],A[u+1..]$ 的半在线卷积。

为了实现非递归的 $4$-叉 半在线卷积,我们进行形如二项式分组的过程(也可以理解为记录递归栈),对于每一组我们记录 DFT 后的点值,每次维护分组并计算贡献。在小范围上我们进行循环展开的暴力。

illus.png

这里DFT的总长是 $O(n\log n)$ 级别的,因此需要进行内存回收。这里的内存管理需要格外小心。

利用 ocmul 全在线地写完 $\exp$, 根号, 求逆, 快速幂... 大功告成!

举个例子,实现快速幂 $F^k$:注意到 $(F^k)'F=kF^kF'$,在线地计算 $\int kF^kF'F^{-1}$ 即可($F^k$ 即本身)。

如果你和我一样认为形如 copy(new ocint(new ocquo(new ocde( 的一大串不好看,我们可以使用一个包裹类包起来:

struct pipeline {
    ocpoly*p;
    pipeline() {p=new ocpoly();}
    pipeline(ocpoly *q) {assert(q);p=q;}
    pipeline(poly s) {p=new ocfixed(s);}
    mi get(int n) {return p->get(n);}
    poly await(int n) {return p->await(n);}
    void set(pipeline q) {p->copy(q.p);}
};
pipeline operator + (pipeline a,pipeline b)
{return pipeline(new ocadd(a.p,b.p));}

这样就可以把这个包裹类 pipeline 当做对象愉快地进行操作了!我们甚至可以保留与多项式部分相同的接口。

loj150 挑战多项式

没错,我们可以全在线地实现本题!代码也长得差不多。

int n,m,f[123456];
int main()
{
    scanf("%d%d",&n,&m); ++n;
    for(int i=0;i<n;++i) scanf("%d",f+i);
    pipeline p(vector<int>(f,f+n));
    p=p-gexp(gint(ginv(gsqrt(p,mod_sqrt(f[0],MOD)))));
    p=gcorner(gln(gcorner(p,{1})),{1});
    gde(gexp(gln(p)*m)).await(n).print_len(n-1,'\n');
}

跑了 sum 8958ms/max 478ms,也就慢了一倍左右。

uoj50 链式法则

$F(x) \equiv \int\left(\frac{1}{2}A(x)F^2(x) + 1 \right)dx$,直接照抄式子!

int n; poly a(-1);
char str[288888];
int main()
{
    scanf("%d%s",&n,str);
    a.reserve(n+1);
    for(int i=0;i<n;++i)
        if(str[i]=='1') a[i]=rfac[i];
    pipeline o;
    o.set(gint(a*o*o/2+1_p)); //!
    poly g=o.await(n+1);
    for(int i=1;i<=n;++i) {
        mi ans=g.get(i)*fac[i];
        printf("%d\n",(int)ans);
    }
}

主要就一行代码,sum 1042ms/max 432ms,(在写作本文时)喜提ranklist第三。

loj6684 有根无标号「奇树」计数

$F(x) \equiv x\text{MSet}(\text{MSet}(F(x))-1)$,直接照抄。

int n;
int main() {
    scanf("%d",&n);
    pipeline p;
    p.set(gshl(MSet(MSet(p)-1_p),1));
    poly o=p.await(n+1);
    for(int i=1;i<=n;++i)
        printf("%d\n",(int)o.get(i));
}

同样是一行代码,跑了 max 200ms,在写作本文时喜提ranklist第一(忽略一位打表的同学)。

loj6538 烷基计数 加强版 加强版

也是照抄一行式子。

int n;
int main() {
    scanf("%d",&n);
    pipeline T;
    T.set(gcorner(gshl((gamp(T,3)*2+T*(gamp(T,2)*3+T*T))/6,1),{1}));
    poly p=T.await(n+1);
    printf("%d\n",p.get(n));
}

跑了 sum 582ms/max 167ms,在写作本文时(可以)排在ranklist第六页(虽然这次有点惨但好歹能过)。

我们也可以方便地实现更为复杂的全在线过程,例如 $F\equiv GH(1+x)+3x^2$,$G\equiv FH(1+x^2)+5x$,$H\equiv FG(1-x^3)+3x^4$。

int main() {
    pipeline f,g,h;
    f.set(gcorner(g*h*"1+x"_p+"3x^2"_p,{0}));
    g.set(gcorner(f*h*"1+x^2"_p+"5x"_p,{0}));
    h.set(gcorner(f*g*"1-x^3"_p+"3x^4"_p,{0}));
    poly F=f.await(100000);
    poly G=g.await(100000);
    poly H=h.await(100000);
    F.dump('\n');G.dump('\n');H.dump('\n');
}

可以看到全在线方法写起来比较优美,适用性广泛,效率也比较快,是我们熟悉的迭代方法的有力补充。


最后我们给出以上模板的实现,如果有更好的实现 / 更好的函数类型名 / 发现bug(这个比较重要)欢迎在评论区讨论。由于代码较长,建议自行截取需要的部分进行使用。

当然,发这份代码的本意并不是希望大家做题的时候上来就贴这份板子,主要还是希望给大家提供写板子的思路,以及普及一下在线卷积的作用?

代码更新:2021.08.07 (paste, gist) 修复了一处内存泄漏;修改了一些代码风格;修改了幂和的实现;换了个快点的ntt、细节上的常数优化(减少bitrev的调用、优化ocmul的小范围暴力)、加入了半在线卷积;加入了多项式欧几里得算法和有理分式重建

如何hack统一省选d2t3(详细揭秘)

2020-06-24 23:47:37 By fjzzq2002

有许多同学写了统一省选d2t3,由于种种原因,当生成树个数在模意义下为0但存在生成树时会输出错误的答案。那么下面就让小编带大家一起了解一下如何生成一个在模意义下生成树个数为0的连通图。

我们考虑这样的一个构造,我们生成两个子图 $G_1$ 和 $G_2$,$G_1$ 和 $G_2$ 中前两个点都有一条边,我们考虑这两个图包含这条边的生成树个数和总生成树个数,记为 $(f_1,g_1),(f_2,g_2)$,那么我们考虑生成这样的一张图:

aa.png

容易算出该图恰有 $f_1g_2+f_2g_1+g_1g_2=(f_1+g_1)(f_2+g_2)-f_1f_2$ 个生成树。$(f_1+g_1)(f_2+g_2)-f_1f_2\equiv 0$即$\frac{f_1+g_1}{f_1}\equiv \frac{f_2}{f_2+g_2}$。我们生成若干($O(\sqrt{mod})$)个随机图并meet in the middle即可。

下面附一组生成的数据,边权是在 $\{1,2\}$ 中随机的。

17 44
1 3 1
1 4 2
1 7 1
1 8 2
1 9 1
2 5 1
2 7 1
2 9 1
3 5 2
3 6 1
3 7 1
3 8 2
3 9 1
4 7 2
4 9 1
5 6 1
5 8 2
5 9 1
6 7 1
6 8 2
6 9 1
8 9 2
8 10 2
9 10 2
9 11 2
9 12 2
9 13 1
9 14 2
9 15 2
9 16 1
9 17 1
10 11 1
10 12 2
10 15 2
10 17 2
12 13 2
12 14 2
12 15 1
12 16 1
12 17 2
13 15 1
13 16 2
14 17 1
16 17 1

这就是关于如何hack统一省选d2t3的事情了,大家有什么想法呢,欢迎在评论区告诉小编一起讨论哦!

fjzzq2002 Avatar