POJ 1707 - Sum of powers

http://poj.org/problem?id=1707

概要

\(S_k(n) = \Sigma_{i=1} ^ n i ^ k\)

と \(S_k(n)\) を定義する。これを

\(S_k(n) = \frac{1}{M}(a_{k+1}n ^ {k+1} + a_{k}n ^ k + \cdots + a_1 n + a_0\)

と表したとき、\(M, a_{k+1}, a_k, \cdots, a_1, a_0\) を答える。 ただし、\(a_{k+1}, a_k, \cdots, a_1, a_0\) は整数であり、\(M\) は最小のものをとる。

制約

解法

定義から \(S_k(n) - S_k(n-1) = n ^ k\) であり、 \(S_k(n)\) を多項式で表したもので同様に \(S_k(n) - S_k(n-1)\) を求めて、 この恒等式を係数に関する連立方程式で解いて \(a_{k+1}, a_k, \cdots, a_1\) を有理数で求める。 \(a_0\) は常に 0。

  1 #include <iostream>
  2 #include <vector>
  3 using namespace std;
  4 
  5 long long gcd(long long a, long long b)
  6 {
  7   if (a < b) {
  8     swap(a, b);
  9   }
 10   long long r;
 11   while ((r = a % b) != 0LL) {
 12     a = b;
 13     b = r;
 14   }
 15   return b;
 16 }
 17 
 18 struct ratio
 19 {
 20   long long n, d;
 21   ratio() : n(0), d(1) {}
 22   ratio(long long a) : n(a), d(1) {}
 23   ratio(long long a, long long b) : n(a), d(b) { normalize(); }
 24 
 25   void normalize()
 26   {
 27     if (n == 0) {
 28       d = 1;
 29     } else {
 30       const long long g = gcd(n, d);
 31       n /= g;
 32       d /= g;
 33       if (d < 0) {
 34         n = -n;
 35         d = -d;
 36       }
 37     }
 38   }
 39 
 40   ratio inverse() const { return ratio(d, n); }
 41 
 42   bool operator==(const ratio& r) const { return n == r.n && d == r.d; }
 43   bool operator==(int x) const { return d == 1 && n == x; }
 44   ratio operator-() const { return ratio(-n, d); }
 45 
 46   ratio& operator-=(const ratio& r)
 47   {
 48     const long long dd = d * r.d;
 49     const long long nn = n*r.d - r.n*d;
 50     d = dd;
 51     n = nn;
 52     normalize();
 53     return *this;
 54   }
 55   ratio operator*(const ratio& r) const { return ratio(n * r.n, d * r.d); }
 56 };
 57 ostream& operator<<(ostream& os, const ratio& r)
 58 {
 59   if (r.d == 1) {
 60     return os << r.n;
 61   } else {
 62     return os << r.n << '/' << r.d;
 63   }
 64 }
 65 
 66 void gaussian_elimination(vector<vector<ratio> >& a, vector<ratio>& b)
 67 {
 68   const int N = a.size();
 69   const int M = a[0].size();
 70   // assert(N <= M)
 71 
 72   for (int i = 0, p = 0; i < M; i++, p++) {
 73     int q;
 74     for (q = p; q < N && a[q][i] == 0; q++);
 75     if (q == N) {
 76       --p;
 77       continue;
 78     }
 79     swap(a[i], a[q]);
 80     swap(b[i], b[q]);
 81 
 82     const ratio r = a[i][i].inverse();
 83     for (int k = i; k < M; k++) {
 84       a[i][k] = a[i][k] * r;
 85     }
 86     b[i] = b[i] * r;
 87 
 88     for (int j = 0; j < N; j++) {
 89       if (j == i) {
 90         continue;
 91       }
 92       const ratio t = a[j][i];
 93       for (int k = i; k < M; k++) {
 94         a[j][k] -= t * a[i][k];
 95       }
 96       b[j] -= t * b[i];
 97     }
 98   }
 99 }
100 
101 int main()
102 {
103   static long long comb[22][22];
104   for (int i = 0; i <= 21; i++) {
105     comb[i][0] = comb[i][i] = 1;
106   }
107   for (int i = 1; i <= 21; i++) {
108     for (int j = 1; j <= 21; j++) {
109       comb[i][j] = comb[i-1][j] + comb[i-1][j-1];
110     }
111   }
112 
113   int K;
114   cin >> K;
115   const int N = K+1;
116   vector<vector<ratio> > a(N, vector<ratio>(N));
117   for (int i = 0; i < N; i++) {
118     for (int j = i; j < N; j++) {
119       a[i][j] = comb[j+1][i];
120       if ((j-i) % 2 != 0) {
121         a[i][j] = -a[i][j];
122       }
123     }
124   }
125   vector<ratio> b(N);
126   b[N-1] = 1;
127   gaussian_elimination(a, b);
128   vector<long long> ans;
129   long long lcm = 1;
130   for (int i = 0; i < N; i++) {
131     const long long g = gcd(lcm, b[i].d);
132     lcm = (lcm * b[i].d)/g;
133   }
134   cout << lcm;
135   for (int i = N-1; i >= 0; i--) {
136     cout << " " << (b[i].n * lcm / b[i].d);
137   }
138   cout << " 0" << endl;
139   return 0;
140 }
poj/1707.cc