HDU 4686 Arc of Dream

题目大意：

        有多个测例（以EOF结束），要求输出求和值a0b0 + a1b1...+ an-1bn-1，其中a0, b0已知，并且ai = ai-1AX + AY，bi = bi-1BX + BY，其中AX、AY、BX、BY已知，其中n ≤ 10^18，a0、b0、AX、AY、BX、BY ≤ 2×10^9。

题目链接

注释代码：

/*                                                
 * Problem ID : HDU 4686 Arc of Dream  
 * Author     : Lirx.t.Una                                                
 * Language   : C                               
 * Run Time   : 250 ms                                                
 * Run Memory : 204 KB                                                
*/

#include <memory.h>
#include <stdio.h>

#define    MOD(x)        ( (x) % 1000000007 )

typedef    __int64        llg;

//思路：
//a[i+1]       =   AY(1) +   AX(a[i]) +    0(b[i]) +    0(a[i]b[i]) + 0(s[i])
//b[i+1]       =   BY(1) +    0(a[i]) +   BX(b[i]) +    0(a[i]b[i]) + 0(s[i])
//a[i+1]b[i+1] = AYBY(1) + AXBY(a[i]) + AYBX(b[i]) + AXBX(a[i]b[i]) + 0(s[i])
//s[i+1]       = AYBY(1) + AXBY(a[i]) + AYBX(b[i]) + AXBX(a[i]b[i]) + 1(s[i])
//其中s[n-1]即为问题的答案
//由上述方程组可得矩阵：
//M = 1   a0  b0 a0b0 s1=a0b0
//    0   0   0   0   0
//    0   0   0   0   0
//    0   0   0   0   0
//    0   0   0   0   0
//
//C = 1  AY  BY  AYBY  AYBY
//    0  AX  0   AXBY  AXBY
//    0  0   BX  AYBX  AYBX
//    0  0   0   AXBX  AXBX
//    0  0   0   0     1
//最终结果就是M×C^(n-1)的元素[1][5]

llg        m[6][6];//main，主矩阵矩阵
llg        c[6][6];//coefficient，系数矩阵
llg        t[6][6];//tmp，临时矩阵

void
mul( llg (*ma)[6], llg (*mb)[6] ) {//矩阵相乘，结果保存在ma中

    int        i, j, k;

    memset(t, 0, sizeof(t));

    for ( i = 1; i < 6; i++ )
        for ( k = 1; k < 6; k++ )
            if ( ma[i][k] )
                for ( j = 1; j < 6; j++ ) {
        
                    t[i][j] += MOD( ma[i][k] * mb[k][j] );
                    t[i][j]  = MOD( t[i][j] );
                }

    memcpy(ma, t, sizeof(t));
}

void
fstp(llg n) {//矩阵快速幂，结果保存在m中

    while (n) {
    
        if ( n & 1 ) mul( m, c );
        mul( c, c );
        n >>= 1;
    }
}

int
main() {

    llg        n;
    llg        a0, b0;
    llg        ax, ay;
    llg        bx, by;

    while ( ~scanf("%I64d", &n) ) {
    

        scanf("%I64d%I64d%I64d%I64d%I64d%I64d",
              &a0, &ax, &ay, &b0, &bx, &by);

        if ( !n ) {//特殊情况判断
        
            puts("0");
            continue;
        }

	if ( 1 == n ) {//特判，加速
		
		printf("%d\n", MOD( a0 * b0 ));
		continue;
	}

	//初始化
        memset(m, 0, sizeof(m));
        memset(c, 0, sizeof(c));

        m[1][1] = 1;
        m[1][2] = MOD(a0);
        m[1][3] = MOD(b0);
        m[1][4] = MOD( a0 * b0 );
        m[1][5] = m[1][4];

        c[1][1] = 1;
        c[1][2] = MOD(ay);
        c[1][3] = MOD(by);
        c[1][4] = MOD( ay * by );
        c[1][5] = c[1][4];

        c[2][2] = MOD(ax);
        c[2][4] = MOD( ax * by );
        c[2][5] = c[2][4];

        c[3][3] = MOD(bx);
        c[3][4] = MOD( ay * bx );
        c[3][5] = c[3][4];

        c[4][4] = MOD( ax * bx );
        c[4][5] = c[4][4];

        c[5][5] = 1;

        fstp(n - 1);//快速幂
        printf("%I64d\n", m[1][5]);
    }

    return 0;
}

无注释代码：

#include <memory.h>
#include <stdio.h>

#define    MOD(x)        ( (x) % 1000000007 )

typedef    __int64        llg;

llg        m[6][6];
llg        c[6][6];
llg        t[6][6];

void
mul( llg (*ma)[6], llg (*mb)[6] ) {

    int        i, j, k;

    memset(t, 0, sizeof(t));

    for ( i = 1; i < 6; i++ )
        for ( k = 1; k < 6; k++ )
            if ( ma[i][k] )
                for ( j = 1; j < 6; j++ ) {
        
                    t[i][j] += MOD( ma[i][k] * mb[k][j] );
                    t[i][j]  = MOD( t[i][j] );
                }

    memcpy(ma, t, sizeof(t));
}

void
fstp(llg n) {

    while (n) {
    
        if ( n & 1 ) mul( m, c );
        mul( c, c );
        n >>= 1;
    }
}

int
main() {

    llg        n;
    llg        a0, b0;
    llg        ax, ay;
    llg        bx, by;

    while ( ~scanf("%I64d", &n) ) {
    

        scanf("%I64d%I64d%I64d%I64d%I64d%I64d",
              &a0, &ax, &ay, &b0, &bx, &by);

        if ( !n ) {
        
            puts("0");
            continue;
        }

		if ( 1 == n ) {
		
			printf("%d\n", MOD( a0 * b0 ));
			continue;
		}

        memset(m, 0, sizeof(m));
        memset(c, 0, sizeof(c));

        m[1][1] = 1;
        m[1][2] = MOD(a0);
        m[1][3] = MOD(b0);
        m[1][4] = MOD( a0 * b0 );
        m[1][5] = m[1][4];

        c[1][1] = 1;
        c[1][2] = MOD(ay);
        c[1][3] = MOD(by);
        c[1][4] = MOD( ay * by );
        c[1][5] = c[1][4];

        c[2][2] = MOD(ax);
        c[2][4] = MOD( ax * by );
        c[2][5] = c[2][4];

        c[3][3] = MOD(bx);
        c[3][4] = MOD( ay * bx );
        c[3][5] = c[3][4];

        c[4][4] = MOD( ax * bx );
        c[4][5] = c[4][4];

        c[5][5] = 1;

        fstp(n - 1);
        printf("%I64d\n", m[1][5]);
    }

    return 0;
}