2020牛客暑期多校训练营（第一场）D. Domino

考虑使用拉格朗日乘数法
$ $f\left(x_{1}, x_{2}, \cdots, x_{n}, \lambda\right)=\sum_{i=1}^{n}\left(b_ix_{i}\right)+\lambda\left(\sum_{i=1}^{n}\sum_{j=1}^{n} a_{ij} x_{i} x_{j}-t\right) \quad t \leqslant 1$ $
求偏导

$\left\{\begin{array}{l}b_{1}+2 \lambda\left(a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1n} x_{n}\right)=0 \\ b_{2}+2 \lambda\left(a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n}\right)=0 \\ b_{n}+2 \lambda\left(a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots++a_{n n} x_{n}\right)=0 \\ \sum_{i=1}^{n} \sum_{j=1}^{n}\left(a_{i j} x_{i} x_{j}\right)=t\end{array}\right.$

可得

$\left\{\begin{array}{l}B+2 \lambda A x=0 \\ x^{\top} A x=t\end{array}\right.$ $\begin{array}{l}\because B+2 \lambda A X=0 \Rightarrow 2 \lambda A X=-B \Rightarrow X=-\frac{A^{-1} B}{2 \lambda} \\ \because B^{T} x=-B^{T} \frac{A^{\top}}{2 \lambda} B=X^{T} B \Rightarrow x^{T}=-B \frac{T{A^{-1}}}{2 \lambda}\end{array}$

代入运算可得

$\begin{aligned}-B \frac{T A^{-1}}{2 \lambda} \cdot A \cdot\left(-\frac{A^{-1} B}{2 \lambda}\right) &=t \\ \frac{B^{\top} A^{-1} B}{4 \lambda^{2}} &=t \end{aligned}$ $\begin{aligned}\left(\sum_{i=1}^{n} b_{i} x_{i}\right)^{2}=\left(B^{\top} x\right)^{2} &=\left(-\frac{1}{2 x} \cdot B^{\top} \cdot A^{-1} \cdot B\right)^{2} \\ &=\frac{1}{4\lambda^{2}} \cdot\left(B^{\top} A^{-1} \cdot B\right)^{2} \\ &=t \cdot\left(B^{\top} \cdot A^{-1} \cdot B\right) \end{aligned}$ $\begin{array}{l} \because t \leq 1 \\ \therefore t 取1时最大 \\ 最大值为 B^{T} \cdot A^{-1} \cdot B\end{array}$

代码

#include <bits/stdc++.h>
using namespace std;
#define paii pair<int,int>
#define fr first
#define sc second
typedef long long ll;
const ll mod=998244353;
ll a[410][410];
ll b[410];
ll c[410];
int n,is[410],js[410];
void exgcd(int a,int b,int &x,int &y){
    if(!b)return x=1,y=0,void();
    exgcd(b,a%b,y,x);y-=x*(a/b);
}
int inv(int p){
    int x,y;exgcd(p,mod,x,y);
    return (x+mod)%mod;
}
void inv(){
    for(int k=1;k<=n;k++){
        for(int i=k;i<=n;i++) // 1
            for(int j=k;j<=n;j++)if(a[i][j]){
                is[k]=i,js[k]=j;break;
            }
        for(int i=1;i<=n;i++) // 2
            swap(a[k][i],a[is[k]][i]);
        for(int i=1;i<=n;i++)
            swap(a[i][k],a[i][js[k]]);
        if(!a[k][k]){
            puts("No Solution");
            exit(0);
        }
        a[k][k]=inv(a[k][k]); // 3
        for(int j=1;j<=n;j++)if(j!=k) // 4
            (a[k][j]*=a[k][k])%=mod;
        for(int i=1;i<=n;i++)if(i!=k) // 5
            for(int j=1;j<=n;j++)if(j!=k)
                (a[i][j]+=mod-a[i][k]*a[k][j]%mod)%=mod;
        for(int i=1;i<=n;i++)if(i!=k) // 就是这里不同
            a[i][k]=(mod-a[i][k]*a[k][k]%mod)%mod;
    }
    for(int k=n;k;k--){ // 6
        for(int i=1;i<=n;i++)
            swap(a[js[k]][i],a[k][i]);
        for(int i=1;i<=n;i++)
            swap(a[i][is[k]],a[i][k]);
    }
}
void work()
{
    while(~scanf("%d",&n)){
        memset(c,0,sizeof(c));
        memset(a,0,sizeof(a));
        memset(is,0,sizeof(is));
        memset(js,0,sizeof(js));
        memset(b,0,sizeof(b));
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                scanf("%lld",&a[i][j]);
                a[i][j]=(a[i][j]%mod+mod)%mod;
            }
        }
        for(int i=1;i<=n;i++){
            scanf("%lld",&b[i]);
            b[i]=(b[i]%mod+mod)%mod;
        }
        inv();
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                c[i]=(c[i]+a[j][i]*b[j]%mod)%mod;
            }
        }
        for(int i=1;i<=n;i++){
            c[i]=(c[i]*b[i])%mod;
        }
        ll ans=0;
        for(int i=1;i<=n;i++){
            ans=(ans+c[i])%mod;
        }
        printf("%lld\n",(ans+mod)%mod);
    }
}
int main()
{
    ios_base::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int T=1;
    //scanf("%d",&T);
    //cin>>T;
    while(T--){
        work();
    }
}