理解[bzoj 3243][Noi2013]向量内积

传送门
题目大意：给定n个d维向量，输出任意一对向量，满足他俩的内积为k的倍数。
我想了一个裸的随机算法：随机选两个向量求内积，然后。。。后10个点只过了两个。。。看到标准解法也是用的随机，一开始感觉非常不服他的随机到底比我的强在哪。。。（后来发现是我太弱）
值得一提的是，此题必须使用随机算法来确定答案。只取全1向量的做法是错误的。
一组很简单的数据就能卡掉（不信？试一试）

meow.in
4 6 2
1 1 0 0 0 0
1 1 1 1 0 0
1 1 1 0 0 1
1 0 0 1 1 1

meow.out(样例)
1 2 

题解吧……先放在后面

#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;
#define rep(I, S, T) for (int I = S; I <= T; I ++)
#define rst(ARR) memset(ARR, 0, sizeof(ARR))
int read(){
    int ret = 0; char ch;
    do ch = getchar(); while (ch<'0' || ch>'9');
    do ret = ret*10+ch-'0', ch = getchar(); while (ch>='0' && ch<='9');
    return ret;
}
const int MAXN = 100005;
const int MAXK = 105;
int n, d, k;
int seq[MAXN][MAXK], diag[MAXN], X[MAXN], res[MAXN], Y[MAXN];
void check(int i, int j) {
    int sum = 0;
    rep(p, 1, d) sum += seq[i][p]*seq[j][p];
    sum %=k;
    if (sum==0 ) {
        if (i>j) swap(i, j);
        printf("%d %d\n", i, j);
        exit(0);
    }
}
void solve2() {
    rep(T, 1, 10) {
        int s = 0;
        rst(res);rst(Y);
        rep(i, 1, n) X[i] = rand()%k, s+=X[i];
        s%=k;
        rep(i, 1, n) rep(j, 1, d) res[j] += seq[i][j]*X[i];
        rep(i, 1, d) res[i]%=k;
        rep(i, 1, n) {
            rep(j, 1, d) Y[i] += res[j]*seq[i][j];
            Y[i] += k-diag[i]*X[i];
            Y[i]%=k;
        }
        rep(i, 1, n) if ((Y[i] + X[i])%k!=s) 
            rep(j, 1, n) if (j!=i)check(i, j);
    }
}
void solve3() {
    rep(i, 1, n) diag[i]=(bool)diag[i];
    rep(T, 1, 7) {
        int s = 0;
        rep(i, 1, d*d) res[i] = 0;
        rep(i, 1, n) X[i] = rand()%k, s+=X[i];
        s%=k; 
        rep(i, 1, n) {
            int *s = seq[i], pt = 0;
            rep(j, 1, d) rep(p, 1, d) res[++pt] += s[j]*s[p]*X[i];
        }
        rep(i, 1, d*d) res[i]%=k;
        rep(i, 1, n) {
            int *sq = seq[i], pt = 0;
            Y[i] = 0;
            rep(j, 1, d) rep(p, 1, d) Y[i] += res[++pt]*sq[j]*sq[p];
            Y[i] += k-diag[i]*X[i]; Y[i]%=k;
            if ((Y[i] + X[i])%k!=s) 
                rep(j, 1, n) if (j!=i) check(i, j);
        }
    }
}
int main()
{
    srand(223333333);
    scanf("%d%d%d", &n, &d, &k);
    rep(i, 1, n) rep(j, 1, d) seq[i][j] = read()%k;
    rep(i, 1, n) {
        rep(j, 1, d) diag[i] += seq[i][j]*seq[i][j];
        diag[i]%=k;
    }
    if (k==2) solve2(); else solve3();
    puts("-1 -1");
    return 0;
}

首先！模型转换：$n$个$d$维向量构成了一个$n*d$的矩阵$A$，现在计算$A*A^T$，得到的$n$阶方阵$B$中
$B_{i,j}$表示向量i与向量j的内积
好啊！非常优美，但是直接计算复杂度仍然是$O(n^2d)$
暂时只考虑$k=2$
我们更关心是否存在0元素对吧
可以拿$B$和全一方阵$C$比较一下，若有不同，说明存在0元素
为了降低复杂度，我们不能对$n*n$矩阵进行一一计算，所以不妨取一个$1*n$的随机向量$X$
根据结合律，计算$X*A*A^T$和$X*C$，比较结果。而向量乘矩阵的复杂度是$O(nd)$的，比较资磁
（似乎这是判断矩阵是否相等的经典办法？）
只要有一个元素不同=>对应列上存在一个0向量，只需要暴力枚举寻找位置
没有元素不同=>X选得不好或者本来就没有
根据定理**，正确概率至少有$1/2$，那么重复十余次就可以得到结果了
（2014年胡泽聪的集训队论文可以参考一下）

然后问题来了：
Q1.$B$ 中 对角线上可能存在0丫（$a_i*a_i=0$）
A1:这个比较重要，要将干扰排除。幸运的是，只有n个内积需要计算，这样计算$X*A*A^T$时候需要把这n个内积的结果减去（消除影响），然后在对角线上加一进行比较。复杂度还是$O(n*d)$

Q2.为什么要随机。。直接取全一向量不行么？
A2:设$D=X*A*A^T$，这里解释一下$D$代表的意义：
把$A*A^T$得到的$n$阶方阵$B$中，把第$j$列拿出来，把$B$中每一行的元素与$X$中元素对应相乘后相加，得到的结果就是$D_j$。如果$X$取全一向量，那么$D_j$就是这一列上所有元素的和。如果不全是一，那么$D_j$是其中一个子集的和。
看起来取全一向量非常靠谱，可是不要忘记——我们是在$mod2$意义下操作。如果这一列1的取值个数正好是2的倍数，那么这一个位置$D_j=0$。如果所有位置都是0，会被认为是无解，尽管可能有不止一对向量的内积是0。
理想的情况当然应该是先把$B$对2取模，然后对于转换后的矩阵进行整数系中的运算。但是……这样就不能应用乘法分配律了……
网上有一部分没有使用随机化的做法，随便找出一个，就能用开头的办法卡掉。
所以靠谱的做法是：在$k=2$下随机，计算结果第$i$列的值为从$B$矩阵中找出第$i$列一个随机子集的点积之和，还是有挺大的概率得到一个非0数字的

Q3.同样是随机，为什么裸随机不如套用矩阵进行随机？
A3:别忘了，在标准做法中，每次可以对一整列进行判定，只要有一个出错，整个就会出错，而n列的判定问题又可以在$O(nd)$时间内解决（随机化）。相当于是把多个内积打包在一起进行判定。

Q4.光顾着说$k=2$去了，$k=3$该怎么做？和$k=2$有什么区别？
A4:这一次，矩阵B中不全是1了（有可能是2），不能把它和全一矩阵比较了。
怎么办呢？把每一个内积的值平方！因为$2\equiv-1(mod3)$，平方以后还是1，就转化成了全1矩阵
把内积平方拆开：

$$(\sum_{i=1}^d a_ib_i) * (\sum_{j=1}^d a_jb_j) = \sum_{i=1}^d{\sum_{j=1}^d{a_{i}a_{j}b_{i}b_{j}}}$$
可以看成一个d^2维向量 $\vec c$，其中 $$c_{(i-1)*d+j}=a_i*a_j,\quad1\leq i,j\leq n$$
$b$与$a$同理
这样就完成了问题转化，只是需要用 $O(d^2)$时间计算了

$BTW$
当初上uoj想找一份标程扒一下（因为太弱了，并不会），发现好多人用了一样的代码汗，而且真正用随机向量的人不多啊。
构造了反例，想Hack一下，可惜那道题是一道spj，没法Hack
于是出于业（xian）界（zhe）良（mei）心（shi），我想vfk大大发了邮件，希望加强一下数据
结果第二天发现vfk凌晨回复了邮件！并且添加了Extra Test！这种敬业精神必须要赞！！
我还想加强一下BZOJ数据。。不过TA1111爷似乎也注意到了。。比我早几天。。那么我就不把事情做绝了