zoukankan      html  css  js  c++  java
  • 高斯消元

    [HNOI2013]游走 这个的方程比较经典

    hdoj 7109 n^3预处理n^2询问修改矩阵最后一列

    //#pragma GCC optimize("Ofast,no-stack-protector,unroll-loops,fast-math")
    //#pragma GCC target("sse,sse2,sse3,ssse3,sse4.1,sse4.2,avx,avx2,popcnt,tune=native")
    
    //#include <immintrin.h>
    //#include <emmintrin.h>
    #include <bits/stdc++.h>
    using namespace std;
    #define rep(i,h,t) for (int i=h;i<=t;i++)
    #define dep(i,t,h) for (int i=t;i>=h;i--)
    #define ll long long
    #define me(x) memset(x,0,sizeof(x))
    #define IL inline
    #define rint register int
    
    //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    //mt19937 myrand(seed);
    //uniform_int_distribution<ll>dist(1,10000000);
    //dist(myrand);
    inline ll rd(){
        ll x=0;char c=getchar();bool f=0;
        while(!isdigit(c)){if(c=='-')f=1;c=getchar();}
        while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
        return f?-x:x;
    }
    char ss[1<<26],*A=ss,*B=ss;
    IL char gc()
    {
        return A==B&&(B=(A=ss)+fread(ss,1,1<<26,stdin),A==B)?EOF:*A++;
    }
    template<class T>void maxa(T &x,T y)
    {
        if (y>x) x=y;
    }
    template<class T>void mina(T &x,T y)
    {
        if (y<x) x=y;
    }
    template<class T>void read(T &x)
    {
        int f=1,c; while (c=gc(),c<48||c>57) if (c=='-') f=-1; x=(c^48);
        while(c=gc(),c>47&&c<58) x=x*10+(c^48); x*=f;
    }
    char sr[1<<26],z[20]; int C=-1,Z;
    template <class T> void wer(T x)
    {
      if (x<0) sr[++C]='-',x=-x;
      while (z[++Z]=x%10+48,x/=10);
      while (sr[++C]=z[Z],--Z); sr[++C]='
    ';
    }
    const int mo=998244353;
    ll fsp(int x,int y)
    {
        if (y==1) return x;
        ll ans=fsp(x,y/2);
        ans=ans*ans%mo;
        if (y%2==1) ans=ans*x%mo;
        return ans;
    }
    struct cp {
        ll x,y;
        cp operator +(cp B)
        {
            return (cp){x+B.x,y+B.y};
        }
        cp operator -(cp B)
        {
            return (cp){x-B.x,y-B.y};
        }
        ll operator *(cp B)
        {
            return x*B.y-y*B.x;
        }
        int half() { return y < 0 || (y == 0 && x < 0); }
    };
    struct re{
        int a,b,c;
    };
    const int N=510;
    int w[N],siz[N],pp[11];
    ll h[N][10],gl[N];
    re g[N][N];
    int n,m,q;
    struct gauss{
      ll jl[N][N],f[N][N],ans[N];
      void Gauss()
      {
          n--; 
        rep(i,1,n) jl[i][i]=1;
        rep(i,1,n)
        {
          int kk=fsp(f[i][i],mo-2); 
          rep(j,1,n)
            if (i!=j)
            {
              ll t=-1ll*f[j][i]*kk%mo;
              rep(k,1,n) (f[j][k]+=t*f[i][k])%=mo;
              rep(k,1,n) (jl[j][k]+=t*jl[i][k])%=mo;
            }
        }
        rep(i,1,n)
        {
          int kk=fsp(-f[i][i],mo-2);
          rep(j,1,n) jl[i][j]=jl[i][j]*kk%mo;
        }
        n++;
      }
      inline ll js()
      {
        ll anss=0;
          rep(i,1,n-1)
          {
              ans[i]=0;
              rep(j,1,n-1) (ans[i]+=jl[i][j]*h[j][7])%=mo; 
          }
          rep(i,1,n-1)
            rep(j,1,n)
              if (g[i][j].b)
              {
                  (anss+=ans[i]*gl[i]%mo*g[i][j].b)%=mo;
              }
          return anss;
      }
    } M;
    ll ans=0,ans2=0;
    inline void cl()
    {
       ans=0;
       ll sum=0;
       rep(i,1,n) sum+=w[i]; sum=fsp(sum,mo-2); 
       rep(i,1,n) gl[i]=fsp(siz[i],mo-2);
       rep(i,1,n-1) h[i][0]=w[i]*sum%mo;
       rep(i,1,7)
         rep(j,1,n-1)
         {
           h[j][i]=0;
           rep(k,1,n)
             if (g[j][k].a) (h[j][i]+=h[k][i-1]*gl[k])%=mo;
         }
       pp[0]=1; rep(i,1,10) pp[i]=pp[i-1]*2;
       rep(i,0,6)
         rep(j,1,n-1)
           rep(k,1,n)
             if (g[j][k].a) (ans+=h[j][i]*gl[j]%mo*max(g[j][k].a/pp[i],g[j][k].b))%=mo;
    }
    inline void pr()
    {
        ll now=((ans+ans2)%mo+mo)%mo;
        wer(now);
    }
    int main()
    {
       ios::sync_with_stdio(false);
       read(n); read(m);read(q);
       rep(i,1,n-1) read(w[i]);
       rep(i,1,m)
       {
            int x,y,a,b;
            read(x); read(y); read(a); read(b); 
            siz[x]++; siz[y]++;
         g[x][y]=(re){a,b};
         g[y][x]=(re){a,b};
       }
       cl();
       rep(i,1,n-1)
       { 
         M.f[i][i]=-1;
         rep(j,1,n-1)
           if (g[i][j].a)
             M.f[i][j]=gl[j];
       }
       M.Gauss();
       ans2=M.js();
       pr();
       while (q--)
       {
            int kk,x,y,a,b,c;
            read(kk);
            if (kk==1)
            {
               read(x); read(y); read(a); read(b);
               if (x!=n) rep(i,0,6) (ans-=h[x][i]*gl[x]%mo*max(g[x][y].a/pp[i],g[x][y].b))%=mo;
               if (y!=n) rep(i,0,6) (ans-=h[y][i]*gl[y]%mo*max(g[y][x].a/pp[i],g[y][x].b))%=mo;
               if (x!=n) (ans-=M.ans[x]*gl[x]%mo*g[x][y].b)%=mo;
               if (y!=n) (ans-=M.ans[y]*gl[y]%mo*g[y][x].b)%=mo;
               g[x][y]=g[y][x]=(re){a,b};
               if (x!=n) rep(i,0,6) (ans+=h[x][i]*gl[x]%mo*max(g[x][y].a/pp[i],g[x][y].b))%=mo;
               if (y!=n) rep(i,0,6) (ans+=h[y][i]*gl[y]%mo*max(g[y][x].a/pp[i],g[y][x].b))%=mo;
               if (x!=n) (ans+=M.ans[x]*gl[x]%mo*g[x][y].b)%=mo;
               if (y!=n) (ans+=M.ans[y]*gl[y]%mo*g[y][x].b)%=mo;
            } else
            {
                read(x); read(c);
                w[x]=c; cl(); ans2=M.js();
            }
            pr();
       }
       fwrite(sr,1,C+1,stdout);
       return 0;
    }
    View Code
  • 相关阅读:
    爬虫-scrapy初试
    python-爬虫day1
    django 内存地址列表-->转换为-->字典
    django 基于 form 验证 确认密码的注册
    django 请求过程,生命周期
    django7 models 高级应用
    django6-项目练习
    Mysql之起始
    python之IO模型
    python模块之Gevent(协程)
  • 原文地址:https://www.cnblogs.com/yinwuxiao/p/15407108.html
Copyright © 2011-2022 走看看