在任意代数结构上的多项式乘法学习笔记

前言

Stop learning useless algorithms, go and solve some problems, learn how to use binary search.

以下内容大多是作者看完《如何在任意代数结构上做多项式乘法》^[1] 后口胡的，所以可能和原文章不太一样。如果错了或者有更好的做法请告诉我。

分圆多项式

定义为 $Φ_{n} (x) = \prod_{1 \leq k < n, gcd (k, n) = 1} (x - ω_{n}^{k})$ .

也可以感性理解为 $x = ω_{n}$ 时 $x^{n} - 1 = 0$ ，约掉一些“显然”不为 $0$ 的因式后剩下的素多项式。

分圆多项式都是整系数素多项式，且 $Φ_{n} (x)$ 最高次数为 $φ (n)$ 。

结论：多项式 $f (x)$ 代入 $x = ω_{n}$ 后做运算得到的结果（用 $ω_{n}$ 表示，最后再把 $ω_{n}$ 换成 $x$ ）等于先做运算再 $mod Φ_{n} (x)$ 得到的结果。感性理解就是 $Φ_{n} (ω_{n}) = 0$ 。

当 $n = p^{m}$ （ $p$ 为素数， $m \geq 1$ ）时， $Φ_{n} (x) = \sum_{i = 0}^{φ (n)} x^{i}$ 。

算法原理

要求：三个群 $(A, +_{A}), (B, +_{B}), (C, +_{C})$ ，乘法运算 $\cdot : A \times B \to C$ 具有分配律 $(a_{1} +_{A} a_{2}) \cdot (b_{1} +_{B} b_{2}) = a_{1} \cdot b_{1} +_{C} a_{1} \cdot b_{2} +_{C} a_{2} \cdot b_{1} +_{C} a_{2} \cdot b_{2}$ 。

此时必有 $\forall b \in B, e_{A} \times b = e_{c}$ 和 $\forall a \in A, a \times e_{b} = e_{c}$ ，其中 $e_{A}, e_{B}, e_{C}$ 分别为 $A, B, C$ 中的单位元。于是将 $A$ ， $B$ 的高位填对应的单位元即可。

证明：

$a \cdot b = (a +_{A} e_{A}) \cdot b = a \cdot b +_{C} e_{A} \cdot b$ ，两边加上 $a \cdot b$ 在 $C$ 中的逆元即可。

$e_{B}$ 是类似的。

$C$ 最好还能有较快（ $O (1)$ ）的自然数乘，定义为多个相同的元素加在一起。

Part 1. 解决除法

IDFT 最后要除以长度，而 $C$ 中没有定义自然数乘的逆。

一个解决方法是，分别做长为 $2$ 的幂的 DFT 和长为 $3$ 的幂的 DFT，这样每个元素的 $2^{c_{2}}$ 倍和 $3^{c_{3}}$ 倍都已知（ $c_{2}$ 和 $c_{3}$ 取决于长度），类似辗转相除做即可。

Part 2. 解决单位根

这是一个很神仙的做法。

考虑把一部分 $x$ 代入 $ω_{m}$ 满足 $φ (m) > \deg A (x) + \deg B (x)$ ，然后将 $m$ 拆成 $m = p q$ 。取 $p = q = \sqrt{m}$ 能保证最优复杂度，读者自（wo）证（bu）不（hui）难（zheng）。具体实现可以参考代码。

具体地，

\begin{aligned} A (x) & = \sum_{j = 0}^{q - 1} (\sum_{i = 0}^{φ (p) - 1} a_{i q + j} x^{i q}) x^{j} \\ = \sum_{j = 0}^{q - 1} (\sum_{i = 0}^{φ (p) - 1} a_{i q + j} ω_{p q}^{i q}) x^{j} \\ = \sum_{j = 0}^{q - 1} (\sum_{i = 0}^{φ (p) - 1} a_{i q + j} ω_{p}^{i}) x^{j} \end{aligned}

然后将内层带 $ω_{p}$ 的东西看成系数对外层做 DFT。实现时可以做成指针套数组的形式。这个部分可能不太好理解，可以看代码。

做完 DFT 要进行内层元素相乘，可以递归。

最后对分圆多项式取模即可。实现时可以暴力将高位减到低位。

应用

好像没啥用...

目前想到的就是做 $c_{k} = \prod_{i + j = k} a_{i}^{b_{j}}$ 之类的卷积？

实现

给出一份大常数的实现。期待有大佬能优化。

题目是 lgP3803.

#define DEBUG 0
#include <iostream>
#include <algorithm>
#include <cmath>
#define UP(i,s,e) for(auto i=s; i<e; ++i)
#define DOWN(i,e,s) for(auto i=e; i-->s;)
using std::cin; using std::cout;
namespace Poly{ // }{{{
template<int BASE, typename T>
void change(T* arr, int len){
    int *rev = new int[len];
    rev[0] = 0;
    UP(i, 1, len){
        rev[i] = rev[i/BASE]/BASE;
        rev[i] += i%BASE*(len/BASE);
    }
    UP(i, 0, len) if(rev[i] > i) std::swap(arr[i], arr[rev[i]]);
    delete[] rev;
}
template<int BASE, class A>
void fft(A **a, int len, int siz, bool idft){ // siz == len(a[0])
    static A *tmp[BASE];
    UP(i, 0, BASE){
        tmp[i] = new A[siz];
        //UP(j, siz, siz*BASE){
        //    tmp[i][j].unit();
        //}
    }
    change<BASE>(a, len);
    int wn = siz/BASE;
    for(int h=BASE; h<=len; h*=BASE){
        for(int st=0; st<len; st+=h){
            int w=0;
            UP(i, st, st+h/BASE){
                UP(j, 0, BASE) std::swap(a[i+h/BASE*j], tmp[j]);
                UP(j, 0, BASE){
                    auto &now = a[i+h/BASE*j];
                    std::copy(tmp[0], tmp[0]+siz, now);
                    UP(k, 1, BASE){
                        UP(l, 0, siz){
                            int idx = l-(idft?-1:1)*(w+siz/BASE*j)*k;
                            idx %= siz; idx = idx < 0 ? idx + siz : idx;
                            now[l] += tmp[k][idx];
                        }
                    }
                }
                w += wn;
            }
        }
        wn /= BASE;
    }
    UP(i, 0, BASE) delete[] tmp[i];
    //delete[] tmp;
}
// mod Phi_len(x)
// len = BASE**n
template<int BASE, class A, class B, class C>
int polymul_base(A *a, B *b, C *ret, int len
#if DEBUG
        , int test=0
#endif
        ){
    UP(i, 0, len/BASE*(BASE-1)) ret[i].unit();
    if(len < 100
#if DEBUG
            && !test
#endif
            ){
        int phi_len = len / BASE * (BASE-1);
        UP(i, 0, len) UP(j, 0, len){
            if((i+j)%len >= phi_len) UP(k, 1, BASE){
                ret[(i+j)%len-len/BASE*k] += (a[i]*b[j]).inv();
            } else {
                ret[(i+j)%len] += a[i]*b[j];
            }
        }
        return 1;
    }
    int tim = std::round(std::log(len)/std::log(BASE));
    int p = std::round(std::pow(BASE, tim/2+1));
    int q = std::round(std::pow(BASE, (tim-1)/2));
    A **aa = new A*[BASE*q];
    B **bb = new B*[BASE*q];
    C **cc = new C*[BASE*q];
    UP(i, 0, BASE*q){
        aa[i] = new A[p];
        bb[i] = new B[p];
        cc[i] = new C[p];
    }
    UP(i, 0, q*BASE) UP(j, 0, p){
        aa[i][j].unit(); bb[i][j].unit();// cc[i][j].unit();
    }
    UP(i, 0, q*BASE) UP(j, p/BASE*(BASE-1), p){
        cc[i][j].unit();
    }
    UP(i, 0, q){
        UP(j, 0, p){
            if(j*q+i >= len){ 
                break;
                //aa[i][j].unit(); bb[i][j].unit();
            }
            else {
                aa[i][j] = a[j*q+i];
                bb[i][j] = b[j*q+i];
            }
        }
        //UP(j, p/BASE*(BASE-1), p){ aa[i][j].unit(); bb[i][j].unit(); }
    }
    //UP(i, q, BASE*q){
        //UP(j, 0, p){ aa[i][j].unit(); bb[i][j].unit(); }
    //}
    fft<BASE>(aa, BASE*q, p, false);
    fft<BASE>(bb, BASE*q, p, false);
    int scale;
    UP(i, 0, BASE*q){
        scale = polymul_base<BASE>(aa[i], bb[i], cc[i], p
#if DEBUG
                , test ? test-1 : 0
#endif
                );
    }
    UP(i, 0, BASE*q){
        delete[] aa[i];
        delete[] bb[i];
    }
    delete[] aa;
    delete[] bb;
    fft<BASE>(cc, BASE*q, p, true);
    int pq = p*q;
    int phi_pq = pq/BASE*(BASE-1);
    UP(i, 0, BASE*q) UP(j, 0, p){
        int pl = (i+j*q)%pq;
        if(pl >= phi_pq) UP(k, 1, BASE) ret[(pl-pq/BASE*k)%len] += cc[i][j].inv();
        else ret[pl%len] += cc[i][j];
    }
    UP(i, 0, BASE*q) delete[] cc[i];
    delete[] cc;
    return scale * BASE * q;
}
template<class A, class B, class C>
void polymul(A *a, B *b, C *ret, int len){
    bool swapped = false;
    C *tmp = new C[len*2];
    int l2 = std::round(std::pow(2, std::ceil(std::log(len*2) / std::log(2))));
    int l3 = std::round(std::pow(3, std::ceil(std::log(len*3/2) / std::log(3))));
    int tim2 = polymul_base<2>(a, b, tmp, l2);
    int tim3 = polymul_base<3>(a, b, ret, l3);
    while(tim3 != 1){
        if(tim2 > tim3){
            int scale = tim2 / tim3;
            UP(i, 0, len) tmp[i] += ret[i].inv() * scale;
            tim2 %= tim3;
        }
        std::swap(tim2, tim3);
        std::swap(ret, tmp);
        swapped ^= 1;
    }
    if(swapped){ std::swap(ret, tmp); std::copy(tmp, tmp+len, ret); }
    delete[] tmp;
}
} // {}}}
namespace m{ // }{{{
constexpr int N = 5e6+2;
struct u32{
    unsigned val;
    u32(){}
    u32(unsigned v):val(v){}
    void unit(){val = 0;}
    u32 inv(){ return -val; }
    u32 &operator+=(u32 b){ val += b.val; return *this; }
    u32 &operator*=(u32 b){ val *= b.val; return *this; }
    u32 operator*(u32 b){ return b *= *this;}
    u32 operator*(unsigned x){ return val*x; }
} ia[N], ib[N], ic[N];
int in, im;
void work(){
    cin >> in >> im;
    UP(i, 0, in+1){
        cin >> ia[i].val;
    }
    UP(i, 0, im+1){
        cin >> ib[i].val;
    }
    Poly::polymul(ia, ib, ic, in+im+1);
    UP(i, 0, in+im+1){
        cout << ic[i].val << ' ';
    }
}
} // {}}}
int main(){cin.tie(0)->sync_with_stdio(0); m::work(); return 0;}

update：自己造了道题，欢迎爆踩标算（

附代码（仅给出不同部分）：

namespace m{ // }{{{
constexpr int N = 5e6+2;
using u32 = unsigned;
using ll = long long;
u32 qpow(u32 x, u32 t){
    u32 ans = 1;
    while(t){
        if(t&1){ ans = ans * x ; }
        x = x * x;
        t >>= 1;
    }
    return ans;
}
struct Modbase{
    u32 val;
    Modbase(){}
    Modbase(u32 x):val(x){}
    Modbase &operator+=(Modbase x){ val = x.val*val; return *this; }
    Modbase operator*(int x){ return qpow(val, x); }
    Modbase inv(){ return qpow(val, (~0u)>>1); }
    void unit(){ val = 1;}
} ia[N], ic[N];
struct Modexp{
    u32 val;
    Modexp &operator+=(Modexp x){ val += x.val; return *this; }
    void unit(){ val = 0; }
} ib[N];
Modbase operator*(Modbase x, Modexp y){
    return qpow(x.val, y.val&((~0u)>>1));
}
int in, im;
void work(){
    cin >> in >> im;
    UP(i, 0, in){
        cin >> ia[i].val;
    }
    UP(i, 0, im){
        cin >> ib[i].val;
    }
    int len = in+im-1;
    int l2 = std::round(std::pow(2, std::ceil(std::log(len*2) / std::log(2))));
    int l3 = std::round(std::pow(3, std::ceil(std::log(len*3/2) / std::log(3))));
    len = std::max(l2, l3);
    UP(i, in, len) ia[i].unit();
    Poly::polymul(ia, ib, ic, in+im-1);
    UP(i, 0, in+im-1){
        cout << ic[i].val << ' ';
    }
    cout << '\n';
}
} // {}}}

https://www.cnblogs.com/whx1003/p/16214952.html ↩︎

posted @ 2023-12-26 13:47 383494 阅读(95) 评论(3) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 多项式全家桶

· 128MTT 学习笔记

· 如何在任意代数结构上做多项式乘法

· 多项式乘法

· 任意类型多项式乘法

公告

昵称： 383494
园龄： 2年3个月
粉丝： 8
关注： 11

+加关注

2025年3月

日

一

二

三

四

五

六

x383494

在任意代数结构上的多项式乘法学习笔记

前言

分圆多项式

算法原理

Part 1. 解决除法

Part 2. 解决单位根

应用

实现

公告

搜索

常用链接

合集

随笔分类

随笔档案

阅读排行榜

评论排行榜

推荐排行榜

最新评论

x383494

在任意代数结构上的多项式乘法 学习笔记

前言

分圆多项式

算法原理

Part 1. 解决除法

Part 2. 解决单位根

应用

实现

公告

搜索

常用链接

合集

随笔分类

随笔档案

阅读排行榜

评论排行榜

推荐排行榜

最新评论

在任意代数结构上的多项式乘法学习笔记