zoukankan      html  css  js  c++  java
  • Solution -「LOJ #138」「模板」类欧几里得算法

    (mathcal{Description})

      Link.

      (T) 组询问,每次给出 (n,a,b,c,k_1,k_2),求

    [sum_{x=0}^nx^{k_1}leftlfloorfrac{ax+b}{c} ight floor^{k_2}mod(10^9+7) ]

      (T=1000)(n,a,b,cle10^9)(0le k_1+k_2le 10)

    (mathcal{Solution})

      类欧模板题的集大成者。

      令 (f(n,a,b,c,k_1,k_2)) 表示所求答案,尝试通过分类讨论递归到不同的情况求解。

    1. (a=0lor an+b<c)

      [f(n,a,b,c,k_1,k_2)=leftlfloorfrac{an+b}{c} ight floor^{k_2}sum_{x=0}^nx^{k_1} ]

      预处理 (k+1) 次多项式 (g_k(x)=sum_{i=0}^xi^k) 的各项系数,这个就能直接算出来。

    2. 否则若 (age c)

      [egin{aligned} f(n,a,b,c,k_1,k_2)&=sum_{x=0}^nx^{k_1} left(leftlfloorfrac{a}{c} ight floor x+leftlfloorfrac{(amod c)x+b}{c} ight floor ight)^{k_2}\ &=sum_{i=0}^{k_2} inom{k_2}{i} leftlfloorfrac{a}{c} ight floor^isum_{x=0}^n x^{k_1+i}leftlfloorfrac{(amod c)x+b}{c} ight floor^{k_2-i}\ &=sum_{i=0}^{k_2} inom{k_2}{i} leftlfloorfrac{a}{c} ight floor^i f(n,amod c,b,c,k_1+i,k_2-i) end{aligned} ]

      递归处理。

    3. 否则若 (bge c)

      类似 2.,最终得到

      [f(n,a,b,c,k_1,k_2)=sum_{i=0}^{k_2} inom{k_2}{i} leftlfloorfrac{b}{c} ight floor^i f(n,a,bmod c,c,k_1,k_2-i) ]

      递归处理。

    4. 对于其余情况

      考虑把高斯函数丢到求和指标上,注意到

      [m^k=sum_{j=0}^{m-1}[(j+1)^k-j^k] ]

      以此替代 (leftlfloorfrac{ax+b}{c} ight floor^{k_2}),并交换求和顺序,得到

      [egin{aligned} f(n,a,b,c,k_1,k_2)&=sum_{j=0}^{leftlfloorfrac{an+b}{c} ight floor-1}[(j+1)^{k_2}-j^{k_2}]sum_{x=0}^nx^{k_1}left[x>leftlfloorfrac{cj+c-b-1}{a} ight floor ight]\ &= sum_{j=0}^{leftlfloorfrac{an+b}{c} ight floor-1} [(j+1)^{k_2}-j^{k_2}]sum_{x=0}^n x^{k_1}-sum_{j=0}^{leftlfloorfrac{an+b}{c} ight floor-1} [(j+1)^{k_2}-j^{k_2}]sum_{x=0}^{leftlfloorfrac{cj+c-b-1}{a} ight floor} x^{k_1} end{aligned} ]

      前半部分直接计算,研究后半部分,发现最后一个求和符号是 (g_{k_1}left(leftlfloorfrac{cj+c-b-1}{a} ight floor ight)),枚举次数将其展开:

      [egin{aligned} &=sum_{i=0}^{k_2-1}inom{k_2}{i} sum_{j=0}^{k_1+1}([x^j]g_{k_1})sum_{x=0}^{leftlfloorfrac{an+b}{c} ight floor-1}x^ileftlfloorfrac{cx+c-b-1}{a} ight floor^j\ &=sum_{i=0}^{k_2-1}inom{k_2}{i} sum_{j=0}^{k_1+1}([x^j]g_{k_1})fleft(leftlfloorfrac{an+b}{c} ight floor-1,c,c-b-1,a,i,j ight) end{aligned} ]

      其中 (i) 枚举 ((j+1)^{k_2}) 的次数;(j) 枚举 (g_{k_1}) 的次数;(x) 枚举原来的 (j)。这样也能递归求解了。

      递归中指标变换形如欧几里得算法,设指标值域为 (V)(K=k_1+k_2),则复杂度为 (mathcal O(TK^4log V))

    (mathcal{Code})

      实现时,可以对于 (n,a,b,c),一次算出所有可能的 (f(n,a,b,c,k_1,k_2)) 的值。

    /*~Rainybunny~*/
    
    #include <cstdio>
    #include <cassert>
    #include <iostream>
    
    #define rep( i, l, r ) for ( int i = l, rep##i = r; i <= rep##i; ++i )
    #define per( i, r, l ) for ( int i = r, per##i = l; i >= per##i; --i )
    
    typedef long long LL;
    
    const int MOD = 1e9 + 7;
    int comb[15][15];
    
    inline void subeq( int& a, const int b ) { ( a -= b ) < 0 && ( a += MOD ); }
    inline int sub( int a, const int b ) { return ( a -= b ) < 0 ? a + MOD : a; }
    inline int mul( const long long a, const int b ) { return int( a * b % MOD ); }
    inline int add( int a, const int b ) { return ( a += b ) < MOD ? a : a - MOD; }
    inline void addeq( int& a, const int b ) { ( a += b ) >= MOD && ( a -= MOD ); }
    inline int mpow( int a, int b ) {
        int ret = 1;
        for ( ; b; a = mul( a, a ), b >>= 1 ) ret = mul( ret, b & 1 ? a : 1 );
        return ret;
    }
    
    struct PowerPoly {
        int n, a[12];
        PowerPoly() {}
    
        PowerPoly( const int k ): n( k ) {
            static int mat[15][15];
            for ( int i = 0, p = 0; i <= n + 1; ++i ) {
                addeq( p, mpow( i, n ) );
                mat[i][n + 2] = p;
            }
            rep ( i, 0, n + 1 ) {
                for ( int j = 0, p = 1; j <= n + 1; ++j, p = mul( p, i ) ) {
                    mat[i][j] = p;
                }
            }
    
            rep ( i, 0, n + 1 ) {
                int p = -1;
                rep ( j, i, n + 1 ) if ( mat[j][i] ) { p = j; break; }
                assert( ~p );
                if ( p != i ) std::swap( mat[i], mat[p] );
                int iv = mpow( mat[i][i], MOD - 2 );
                rep ( j, i + 1, n + 1 ) {
                    int t = mul( iv, mat[j][i] );
                    rep ( l, i, n + 2 ) subeq( mat[j][l], mul( mat[i][l], t ) );
                }
            }
            per ( i, n + 1, 0 ) {
                a[i] = mul( mat[i][n + 2], mpow( mat[i][i], MOD - 2 ) );
                rep ( j, 0, i - 1 ) subeq( mat[j][n + 2], mul( a[i], mat[j][i] ) );
            }
        }
    
        inline int operator () ( LL x ) const {
            int ret = 0; x %= MOD;
            for ( int i = 0, p = 1; i <= n + 1; ++i, p = mul( x, p ) ) {
                addeq( ret, mul( a[i], p ) );
            }
            return ret;
        }
    };
    
    const PowerPoly pwr[11]
      = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
    
    struct Atom {
        int res[11][11];
        Atom(): res{} {}
        inline int* operator [] ( const int k ) { return res[k]; }
    };
    
    /* *
     * calculate sum_{x=0}^n x^{k_1} lfloor frac{ax+b}{c} 
    floor^{k_2},
     * where k_1+k_2 <= 10.
     * */
    inline Atom solve( const LL n, const LL a, const LL b, const LL c ) {
        Atom ret; LL m = ( a * n + b ) / c;
        if ( !a || 1ll * a * n + b < c ) {
            rep ( k1, 0, 10 ) {
                int t = int( m % MOD ), p = 1;
                rep ( k2, 0, 10 - k1 ) {
                    ret[k1][k2] = mul( p, pwr[k1]( n ) );
                    p = mul( p, t );
                }
            }
            return ret;
        }
    
        if ( a >= c ) {
            Atom tmp( solve( n, a % c, b, c ) );
            rep ( k1, 0, 10 ) rep ( k2, 0, 10 - k1 ) {
                int t = int( ( a / c ) % MOD ), p = 1;
                rep ( i, 0, k2 ) {
                    addeq( ret[k1][k2], mul( comb[k2][i],
                      mul( p, tmp[k1 + i][k2 - i] ) ) );
                    p = mul( p, t );
                }
            }
            return ret;
        }
    
        if ( b >= c ) {
            Atom tmp( solve( n, a, b % c, c ) );
            rep ( k1, 0, 10 ) rep ( k2, 0, 10 ) {
                int t = int( ( b / c ) % MOD ), p = 1;
                rep ( i, 0, k2 ) {
                    addeq( ret[k1][k2], mul( comb[k2][i],
                      mul( p, tmp[k1][k2 - i] ) ) );
                    p = mul( p, t );
                }
            }
            return ret;
        }
    
        Atom tmp( solve( m - 1, c, c - b - 1, a ) );
        rep ( k1, 0, 10 ) {
            int t = int( m % MOD ), p = 1;
            rep ( k2, 0, 10 - k1 ) {
                ret[k1][k2] = mul( p, pwr[k1]( n ) );
                rep ( i, 0, k2 - 1 ) rep ( j, 0, k1 + 1 ) {
                    subeq( ret[k1][k2], mul( comb[k2][i],
                      mul( pwr[k1].a[j], tmp[i][j] ) ) );
                }
                p = mul( t, p );
            }
        }
        return ret;
    }
    
    int main() {
        comb[0][0] = 1;
        rep ( i, 1, 10 ) {
            comb[i][0] = 1;
            rep ( j, 1, i ) comb[i][j] = add( comb[i - 1][j - 1], comb[i - 1][j] );
        }
    
        int n, a, b, c, k1, k2, T;
        for ( scanf( "%d", &T ); T--; ) {
            scanf( "%d %d %d %d %d %d", &n, &a, &b, &c, &k1, &k2 );
            printf( "%d
    ", solve( n, a, b, c )[k1][k2] );
        }
        return 0;
    }
    
    
  • 相关阅读:
    eclipse + maven 环境配置
    腾讯管家去除桌面快捷小图标
    C# 在同一个项目里启动不同的类文件
    面试题-数据库篇
    面试题-编程篇
    DevExpress控件-GridControl根据条件改变单元格(Dev GridControl 单元格着色)
    Developer Express控件gridcontrol中gridView的某一个单元格是否可以自由输入
    oracle11g如何创建数据库
    通过第三方组件NPOI读取Excel的方法
    Oracle11g常用的命令
  • 原文地址:https://www.cnblogs.com/rainybunny/p/15013739.html
Copyright © 2011-2022 走看看