[NOI2013]向量内积

目录

• 壹、题目描述 ¶
• 贰、题解 ¶
  • ◆ 前言...
  • ◆ 从简单想起。
  • ◆ 归一！
  • ◆ 额外的问题？
• 叁、参考代码 ¶
• 肆、关键之处 ¶

壹、题目描述 ¶

传送门 to LOJ.

贰、题解 ¶

◆ 前言...

没想到这辈子还会这样大力分类讨论。（然后打了一个极长极冗杂的代码）

---

◆ 从简单想起。

我们进行分类讨论，先看看 \(k=2\) 的情况。

如果我们要大力判断两个向量相乘，记 \(X\) 为原向量的堆叠矩阵，那么向量两两相乘的结果就是 \(X\times X^T\). 但是暴力计算这个东西，时间复杂度达到惊人的 \(\mathcal O(n^2)\).

但是，从另外一个方面想，如果不存在任意两个向量相乘在 \(\bmod 2\) 下为 \(0\)，那么最后得到的 \(XX^T\) 势必等于 \(11^T\)，那么，我们可以从这方面入手。

我们可以随机向量 \(\alpha\)，判断 \(\alpha XX^T\) 与 \(\alpha 11^T\) 是否相同，如果存在某一位 \(i\) 不同，那么一定是 \(v_i\) 与某个向量相乘不为 \(1\)，再暴力判断即可。

计算 \(\alpha XX^T\) 复杂度为 \(\mathcal O(nd^2)\)，而 \(\alpha 11^T\) 怎么算不必我多说了吧，而暴力判断 \(i\) 的部分复杂度为 \(\mathcal O(nd)\)，这样是可行的。

---

◆ 归一！

但是，使用相同的方法在 \(k=3\) 似乎就走不动了，因为余数不只是 \(1\)，还可能是 \(2\).

但是，我们发现一个东西 —— \(1^2\equiv2^2\equiv 1\pmod 3\)，那么，继承相似的思路，我们得到 \(XX^T\)，然后将每一个位置上的数字取平方（注意，并不是 \((XX^T)^2\)，是每个位置上的数自己的平方），然后再使用类似的思路计算。

---

◆ 额外的问题？

但是，如何计算 \(\alpha (XX^T)'\) 呢？我们考虑得到的是 \(\beta\) 向量，那么

\[\begin{aligned} \beta_{j}&=\sum_i\alpha_i\times \left(\sum_{k}X_{i,k}X^T_{k,j}\right)^2 \\ &=\sum_i\alpha_i\times \left(\sum_{k}X_{i,k}X^T_{k,j}\right)\left(\sum_{k}X_{i,k}X^T_{k,j}\right) \\ &=\sum_i\alpha_i\times \sum_{k_1}\sum_{k_2}X_{i,k_1}X^T_{k_1,j}X_{i,k_2}X^T_{k_2,j} \\ &=\sum_i\sum_{k_1}\sum_{k_2}\alpha_iX_{i,k_1}X^T_{k_1,j}X_{i,k_2}X^T_{k_2,j} \\ &=\sum_{k_1}\sum_{k_2}X^T_{k_1,j}X^T_{k_2,j}\sum_i \alpha_iX_{i,k_1}X_{i,k_2} \\ \textbf{if }f(k_1,k_2)&\overset\Delta=\sum_i \alpha_iX_{i,k_1}X_{i,k_2}\quad\textbf{then } \\ \beta_j&=\sum_{k_1}\sum_{k_2}X^T_{k_1,j}X^T_{k_2,j}f(k_1,k_2) \end{aligned} \]

预处理 \(f(k_1,k_2)\) 复杂度为 \(\mathcal O(nd^2)\) 的，那么计算 \(\beta_j\) 就是 \(\mathcal O(d^2)\) 的，计算整个 \(\beta\) 就是 \(\mathcal O(nd^2)\) 的了。

叁、参考代码 ¶

#pragma GCC optimize(2)
#include<cstdio>
#include<vector>
#include<cstring>
#include<algorithm>
#include<ctime>
using namespace std;

// #define NDEBUG
#include<cassert>

namespace Elaina{
    #define rep(i, l, r) for(int i=(l), i##_end_=(r); i<=i##_end_; ++i)
    #define drep(i, l, r) for(int i=(l), i##_end_=(r); i>=i##_end_; --i)
    #define fi first
    #define se second
    #define mp(a, b) make_pair(a, b)
    #define Endl putchar('\n')
    #define mmset(a, b) memset(a, b, sizeof a)
    // #define int long long
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int, int> pii;
    typedef pair<ll, ll> pll;
    template<class T>inline T fab(T x){ return x<0? -x: x; }
    template<class T>inline void getmin(T& x, const T rhs){ x=min(x, rhs); }
    template<class T>inline void getmax(T& x, const T rhs){ x=max(x, rhs); }
    template<class T>inline T readin(T x){
        x=0; int f=0; char c;
        while((c=getchar())<'0' || '9'<c) if(c=='-') f=1;
        for(x=(c^48); '0'<=(c=getchar()) && c<='9'; x=(x<<1)+(x<<3)+(c^48));
        return f? -x: x;
    }
    template<class T>inline void writc(T x, char s='\n'){
        static int fwri_sta[1005], fwri_ed=0;
        if(x<0) putchar('-'), x=-x;
        do fwri_sta[++fwri_ed]=x%10, x/=10; while(x);
        while(putchar(fwri_sta[fwri_ed--]^48), fwri_ed);
        putchar(s);
    }
}
using namespace Elaina;

const int maxn=1e5;
const int maxd=1e2;

template<class T> struct matrix{
    int n, m;
    vector< vector<T> >a;
    inline matrix(){ n=m=0; a.clear(); }
    inline matrix(const int N, const int M): n(N), m(M){
        a=vector< vector<T> >(n);
        for(int i=0; i<n; ++i)
            a[i]=vector<T>(m, 0);
    }
    inline matrix epsilon(const int N){
        (*this)=matrix(N, N);
        for(int i=0; i<n; ++i) a[i][i]=1;
        return (*this);
    }
    inline matrix I(const int N){
        n=m=N;
        a=vector< vector<T> >(n);
        for(int i=0; i<n; ++i) a[i]=vector<T>(m, 1);
        return (*this);
    }
    inline matrix operator *(const matrix rhs) const{
        assert(m==rhs.n); matrix ret(n, rhs.m);
        for(int i=0; i<n; ++i) for(int j=0; j<m; ++j) if(a[i][j])
            for(int k=0; k<rhs.m; ++k)
                ret.a[i][k]+=a[i][j]*rhs.a[j][k];
        return ret;
    }
    inline bool operator ==(const matrix rhs) const{
        if(n!=rhs.n || m!=rhs.m) return false;
        for(int i=0; i<n; ++i) for(int j=0; j<m; ++j)
            if(a[i][j]!=rhs.a[i][j])
                return false;
        return true;
    }
    inline bool operator !=(const matrix rhs) const{
        return !((*this)==rhs);
    }
    inline void print() const{
        printf("N == %d, M == %d\n", n, m);
        for(int i=0; i<n; ++i){
            for(int j=0; j<m; ++j)
                writc(a[i][j], ' ');
            Endl;
        }
    }
};

matrix<int>A, B;
vector<int>x[maxn+5];

int n, d, k;

inline void input(){
    n=readin(1), d=readin(1), k=readin(1);
    for(int i=0; i<n; ++i){
        x[i].resize(d);
        for(int j=0; j<d; ++j)
            x[i][j]=readin(1)%k;
    }
    A=matrix<int>(n, d), B=matrix<int>(d, n);
    for(int i=0; i<n; ++i) A.a[i]=x[i];
    for(int j=0; j<n; ++j) for(int i=0; i<d; ++i)
        B.a[i][j]=x[j][i];
}

inline int getrnd(const int Mod){
    return (1ll*rand()*rand())%Mod;
}

inline matrix<int> randVec(){
    matrix<int>vec=matrix<int>(1, n);
    for(int i=0; i<n; ++i)
        vec.a[0][i]=getrnd(100);
    return vec;
}

namespace moduleEqualTwo{
    matrix<int>vec, H, G;
    inline void launch(){
        int p=-1;
        auto check=[&](){
            for(int i=0; i<n; ++i)
                if(H.a[0][i]%2!=G.a[0][i]%2){ p=i; return false; }
            return true;
        };
        auto mulI=[&](){
            int sum=0;
            for(int i=0; i<n; ++i)
                sum+=vec.a[0][i];
            matrix<int>ret(1, n);
            for(int i=0; i<n; ++i)
                ret.a[0][i]=sum;
            return ret;
        };
        rep(_, 1, 20){
            vec=randVec();
            H=vec*A*B, G=mulI();
            if(!check()) break;
        }
        if(p==-1) printf("-1 -1\n");
        else{
            auto vecMul=[](const int a, const int b){
                int ret=0;
                for(int i=0; i<d; ++i) ret+=x[a][i]*x[b][i]%k;
                return ret;
            };
            for(int i=0; i<n; ++i) if(i!=p && vecMul(p, i)%2==0){
                printf("%d %d\n", min(p+1, i+1), max(p+1, i+1)); return;
            }
        }
    }
}

namespace moduleEqualThree{
    matrix<int>vec, H, G;
    int F[maxd+5][maxd+5];
    inline void launch(){
        auto mulI=[&](){
            int sum=0;
            for(int i=0; i<n; ++i)
                sum=(sum+vec.a[0][i])%k;
            matrix<int>ret(1, n);
            for(int i=0; i<n; ++i)
                ret.a[0][i]=sum;
            return ret;
        };
        int p=-1;
        auto check=[&](){
            for(int i=0; i<n; ++i)
                if(H.a[0][i]%3!=G.a[0][i]%3){ p=i; return false; }
            return true;
        };
        rep(_, 1, 20){
            vec=randVec();
            G=mulI();
            for(int k1=0; k1<d; ++k1) for(int k2=0; k2<d; ++k2)
                for(int j=0; j<n; ++j)
                    F[k1][k2]=(F[k1][k2]+1ll*vec.a[0][j]*x[j][k1]%k*x[j][k2]%k)%k;
            H=matrix<int>(1, n);
            for(int i=0; i<n; ++i){
                for(int k1=0; k1<d; ++k1) for(int k2=0; k2<d; ++k2)
                    H.a[0][i]=(H.a[0][i]+1ll*x[i][k1]*x[i][k2]%k*F[k1][k2]%k)%k;
            }
            if(!check()) break;
        }
        if(p==-1) printf("-1 -1\n");
        else{
            auto vecMul=[](const int a, const int b){
                int ret=0;
                for(int i=0; i<d; ++i) ret=(ret+1ll*x[a][i]*x[b][i])%k;
                return ret;
            };
            // printf("p == %d\n", p);
            for(int i=0; i<n; ++i) if(vecMul(p, i)%3==0){
                printf("%d %d\n", min(p+1, i+1), max(p+1, i+1)); return;
            }
        }
    }
}

signed main(){
    srand(114514+1919810);
    input();
    if(k==2) moduleEqualTwo::launch();
    else moduleEqualThree::launch();
    return 0;
}

肆、关键之处 ¶

在 \(\bmod 3\) 下的一个特性 \(1^2\equiv2^2\equiv 1\pmod 3\)；

另外，从 \(k=2\) 的简单想起也很重要，还有就是考察无解时的情形，对求解或许有必要的帮助。