[UVA][高斯消去] 12143 - Stopping Doom's Day
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).
OutputFor 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;
}