2013-08-11 21:26:54 | 人氣(4,204) | 回應(0) | 上一篇 | 下一篇

[UVA][高斯消去] 12143 - Stopping Doom's Day

推薦 0 收藏 0 轉貼0 訂閱站台

 

So! The time of the universe is up and it is the dooms day after five hours :-P, and you must stop it. But to do so you have to know the value of the following expression T:

Because the secret code that will save the universe from being doomed have something to do with the value of the above expression for some value of n.

 0, 0, 120, 3200, 35160, 237120, 1164800, 4572160, 1517472

Input

The input file contains 1000 lines of inputs.

 

Each line contains a single integer n (0<n≤2000000000).

Output

For each line of input produce one line of output. This line contains the value of T.

 

 

Sample Input Output for Sample Input

12

20

1001

0

2199

803

2390


Problem setter: Shahriar Manzoor, Special Thanks: Derek Kisman

 


要推公式有點困難,完全沒有頭緒。

感謝 flere 的隊友想出的這個辦法。

\sum_{i = 1}^{n} i = n(n+1)/2
\sum_{i = 1}^{n} i*i = n(n+1)(2n+1)/6
for product for five sigma, polynomial function(x) has 10 degrees.
assume F(n) = x0 + n^{1}x1 + n^{2}x2 + ... + n^{10}x10
then, give n = 1~10 for F(n). We can solve it by Gaussian elimination.

前幾項 0, 0, 120, 3200, 35160, 237120, 1164800, 4572160, 1517472 ...

不過還是利用高斯消去與逆元分別求出方程長怎麼樣子,以及最後運算答案。

由於中間過程需要逆元,要用 float 形式的高斯消去便會有點困難,因此這裡使用分數形態。

還好沒有爆 long long


#include <stdio.h>
#include <algorithm>
#include <math.h>
using namespace std;
#define LL long long
class Frac {
    public:
        LL a, b;
        Frac() {
            a = 0, b = 1;
        }
        Frac(LL x, LL y) {
            a = x, b = y;
            reduce();
        }
        Frac operator+(const Frac &y) {
            LL ta, tb;
            tb = this->b/gcd(this->b, y.b)*y.b;
            ta = this->a*(tb/this->b) + y.a*(tb/y.b);
            Frac z(ta, tb);
            return z;
        }
        Frac operator-(const Frac &y) {
            LL ta, tb;
            tb = this->b/gcd(this->b, y.b)*y.b;
            ta = this->a*(tb/this->b) - y.a*(tb/y.b);
            Frac z(ta, tb);
            return z;
        }
        Frac operator*(const Frac &y) {
            LL tx, ty, tz, tw, g;
            tx = this->a, ty = y.b;
            g = gcd(tx, ty), tx /= g, ty /= g;
            tz = this->b, tw = y.a;
            g = gcd(tz, tw), tz /= g, tw /= g;
            Frac z(tx*tw, ty*tz);
            return z;
        }
        Frac operator/(const Frac &y) {
            LL tx, ty, tz, tw, g;
            tx = this->a, ty = y.a;
            g = gcd(tx, ty), tx /= g, ty /= g;
            tz = this->b, tw = y.b;
            g = gcd(tz, tw), tz /= g, tw /= g;
            Frac z(tx*tw, ty*tz);
            return z;
        }
    private:
        static LL gcd(LL x, LL y) {
            if(!y)  return x;
            if(x < 0)   x *= -1;
            if(y < 0)   y *= -1;
            LL t;
            while(x%y)
                t = x, x = y, y = t%y;
            return y;
        }
        void reduce() {
            LL g = gcd(a, b);
            a /= g, b /= g;
            if(b < 0)   a *= -1, b *= -1;
        }
};
long long mpow(long long x, long long y, long long mod) {
    if(y == 0)  return 1;
    if(y&1) return (x*mpow((x*x)%mod, y>>1, mod))%mod;
    return mpow((x*x)%mod, y>>1, mod);
}
Frac x[20];
void gauss() {
    int i, j, k, l, m, n;
    Frac mtx[20][20];
#define poly 11
    for(n = 1; n <= poly; n++) {
         //mtx[n][poly+1];
        long long sum = 0;
        for(i = 1; i <= n; i++)
            for(j = 1; j <= n; j++)
                for(k = 1; k <= n; k++)
                    for(l = 1; l <= n; l++)
                        for(m = 1; m <= n; m++)
                            sum += (long long)abs(i-j)*abs(j-k)*abs(k-l)*abs(l-m)*abs(m-i);
        mtx[n][poly+1] = Frac(sum, 1);
        for(i = 1; i <= poly; i++)
            mtx[n][i] = Frac(mpow(n, i-1, 1LL<<60), 1LL);
    }
    for(i = 1; i <= poly; i++) {
        k = i;
        for(j = i+1; j <= poly; j++)
            if(mtx[j][i].a > mtx[k][i].a)
                k = j;
        if(mtx[k][i].a == 0)
            continue;
        if(k != i) {
            for(j = 1; j <= poly+1; j++)
                swap(mtx[k][j], mtx[i][j]);
        }
        for(j = 1; j <= poly; j++) {
            if(i == j)  continue;
            for(k = poly+1; k >= i; k--) {
                mtx[j][k] = mtx[j][k] - mtx[j][i]/mtx[i][i]*mtx[i][k];
            }
        }
    }
    for(i = 1; i <= poly; i++)
        x[i-1] = mtx[i][poly+1]/mtx[i][i];
    /*for(i = 0; i < poly; i++) {
        printf("%lld %lld\n", x[i].a, x[i].b);
    }*/
}
int main() {
    gauss();
    int n;
    int i;
    while(scanf("%d", &n) == 1 && n) {
        long long ret = 0;
        for(i = 0; i <= 10; i++) {
            ret += x[i].a*mpow(n, i, 10007)*mpow(x[i].b, 10007-2, 10007);
            ret %= 10007;
        }
        printf("%lld\n", ret);
    }
    return 0;
}

台長: Morris
人氣(4,204) | 回應(0)| 推薦 (0)| 收藏 (0)| 轉寄
全站分類: 教育學習(進修、留學、學術研究、教育概況) | 個人分類: UVA |
此分類下一篇:[UVA] 12090 - Counting Zeroes
此分類上一篇:[UVA][博弈] 10111 - Find the Winning Move

是 (若未登入"個人新聞台帳號"則看不到回覆唷!)
* 請輸入識別碼:
(有*為必填)
TOP
詳全文