FWT快速沃尔什变换——基于朴素数学原理的卷积算法
阅读原文时间:2023年07月09日阅读:1

这是我的第一篇学习笔记,如有差错,请海涵…


目录

引子

卷积形式

算法流程

OR卷积

AND卷积

XOR卷积

模板


引子

首先,考虑这是兔子

数一数,会发现你有一只兔子,现在,我再给你一只兔子

再数一数,会发现什么?没错,你有两只兔子,也就是说,1+1=2

这就是算数的基本原理了,聪明的你懂了吗?

好,我们可以学FWT了..


卷积形式

我们回忆一下多项式乘法的式子:

这个可以用FFT或NTT优化到O(nlogn)求出每一个Ci,但不是本章的重点,只是引出卷积的概念:

而FWT主要是解决以下三种卷积形式:


算法流程

卷积的算法原理就是把一个数列快速转换成另一种数列,然后每一位元素之间就可以直接单独相乘计算,最后再把答案数列快速转换回来。

FFT体现这个原理的方式就是把多项式转换成点值表达式,然后由于每个点的横坐标相同,纵坐标直接乘起来就得到最终的点值表达式,最后把答案的多项式表达通过点值表达式解出来。

那FWT怎么做呢?

首先就是数列长度的问题,我们知道,多项式乘法最终会得到一个长为lenA+lenB-1的多项式,而考虑位运算的卷积——很容易想出,最终的数列长度一定是,n是A、B大小转换为二进制后的数的最大位数。

我们设数列A的转换数列是DWT(A),转换后的数列A的原数列是IDWT(A)

既然它是位运算,那么我们就按位分治

我们从二进制最高位考虑起,每次把当前位为0或1的元素分开成两个数列,很显然,由于数列长度为,直接每次从中间分开就好了,

那么

这里的“{  ,  }”是把两个数列前后拼一起,A+B是把两个数列排头对齐,然后每一位相加。

具体的系数a,b,c,d是怎么样,or , and 和 xor 的情况是不一样的。

因为是按位或,所以当前位为1的对0没有影响,而0的元素都要对1有影响(0可以 | 1变成1,但是1怎么 | 都不会变成0),于是它的DWT就是这样

这样DWT(A)[i]就相当于下标按位或 i 后等于 i 的元素和,转换回去刚好就把当前位为1的减去为0的就行,即

这就是DWT的逆运算形式吧。

ps:巧合的是,这个玩意其实也是快速莫比乌斯变换FMT,两个是一样的,完全没有区别,也就是说DWT(A)[i]其实也是i的所有子集元素和。

举个栗子

,

解决了!

和or很相像

因为是按位与,所以当前位为0的对1没有影响,而1的元素都要对0有影响(1可以&0变成0,但是0怎么&都不会变成1),于是它的DWT就是这样

这样DWT(A)[i]就相当于下标按位与 i 后等于 i 的元素和,转换回去刚好就把当前位为0的减去为1的就行,即

这又刚好是DWT的逆运算了。

再举个栗子

,

这个就比较特殊了

我们从栗子里会发现,对于异或,我们最后其实要把 a0b0+a1b1 和 a1b0+a0b1 单独刨出来。(这不是废话!)

那么在DWT(C)中,a0b0的系数要和a1b1一样,a1b0的系数要和a0b1一样

……

于是它的DWT就是这样!:

这样DWT(C)就符合条件了,它的IDWT是

这个得看栗子才明白

,


模板

下面是非递归版本的DWT以及IDWT,m为数列长度(

inline void DWTOR(int *s,int m) {
    for(int k = m;k > 1;k >>= 1) {
        for(int i = 0;i < m;i += k) {
            for(int j = i+(k>>1);j < i+k;j ++) {
                int s0 = s[j-(k>>1)],s1 = s[j];
                s[j] = qm((s0 +0ll+ s1) , zxy);
            }
        }
    }
    return ;
}
inline void IDWTOR(int *s,int m) {
    for(int k = 2;k <= m;k <<= 1) {
        for(int i = 0;i < m;i += k) {
            for(int j = i+(k>>1);j < i+k;j ++) {
                int s0 = s[j-(k>>1)],s1 = s[j];
                s[j] = qm((s1 +0ll+ zxy - s0) , zxy);
            }
        }
    }
    return ;
}

inline void DWTAND(int *s,int m) {
    for(int k = m;k > 1;k >>= 1) {
        for(int i = 0;i < m;i += k) {
            for(int j = i+(k>>1);j < i+k;j ++) {
                LL s0 = s[j-(k>>1)],s1 = s[j];
                s[j-(k>>1)] = qm((s0 +0ll+ s1) , zxy);
            }
        }
    }
    return ;
}
inline void IDWTAND(int *s,int m) {
    for(int k = 2;k <= m;k <<= 1) {
        for(int i = 0;i < m;i += k) {
            for(int j = i+(k>>1);j < i+k;j ++) {
                int s0 = s[j-(k>>1)],s1 = s[j];
                s[j-(k>>1)] = qm((s0 +0ll+ zxy - s1) , zxy);
            }
        }
    }
    return ;
}

inline void DWTXOR(int *s,int m) {
    for(int k = m;k > 1;k >>= 1) {
        for(int i = 0;i < m;i += k) {
            for(int j = i+(k>>1);j < i+k;j ++) {
                int s0 = s[j-(k>>1)],s1 = s[j];
                s[j] = qm((s0 +0ll+ zxy - s1) , zxy);
                s[j-(k>>1)] = qm((s0 +0ll+ s1) , zxy);
            }
        }
    }
    return ;
}
inline void IDWTXOR(int *s,int m) {
    for(int k = 2;k <= m;k <<= 1) {
        for(int i = 0;i < m;i += k) {
            for(int j = i+(k>>1);j < i+k;j ++) {
                int s0 = s[j-(k>>1)],s1 = s[j];
                s[j-(k>>1)] = qm((s0 +0ll+ s1) , zxy) *1ll* inv2 % zxy;
                s[j] = qm((s0 +0ll+ zxy - s1) , zxy) *1ll* inv2 % zxy;
            }
        }
    }
    return ;
}

FWT就到这里了,大家都懂了吧

手机扫一扫

移动阅读更方便

阿里云服务器
腾讯云服务器
七牛云服务器

你可能感兴趣的文章