POJ 3807 - Chemist's Math

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

概要

係数が省略された化学反応式が与えられるので，正しい係数を答える．

ここで正しい係数とは，

• すべて正の整数
• 係数の gcd が1
• 左辺と右辺で元素の数が一致する

ようなものである．

制約

• 化学反応式は80文字以下
• BNF からわかるように，化学反応式中に出現する数値は100未満
• 答えの係数は 40000 以下
• うまく計算すれば途中の計算においても 32-bit の整数値に収まる
• 「正しい」係数は必ず一意に定まる (すなわち例の (8) のような入力は無い)

解法

がんばって構文解析した後，有理数でガウスの消去法で連立方程式を解く． このとき制約から，rank(係数行列) = (変数の数)-1 が保証されている．

  1 #include <iostream>
  2 #include <vector>
  3 #include <map>
  4 #include <cassert>
  5 using namespace std;
  6 
  7 long long gcd(long long a, long long b)
  8 {
  9   if (a < b) {
 10     swap(a, b);
 11   }
 12   long long r;
 13   while ((r = a % b) != 0LL) {
 14     a = b;
 15     b = r;
 16   }
 17   return b;
 18 }
 19 
 20 struct ratio
 21 {
 22   long long n, d;
 23   ratio() : n(0), d(1) {}
 24   ratio(long long a, long long b) : n(a), d(b)
 25   {
 26     normalize();
 27   }
 28 
 29   void normalize()
 30   {
 31     if (n == 0) {
 32       d = 1;
 33     } else {
 34       const long long g = gcd(n, d);
 35       n /= g;
 36       d /= g;
 37       if (d < 0) {
 38         n = -n;
 39         d = -d;
 40       }
 41     }
 42   }
 43 
 44   ratio inverse() const { return ratio(d, n); }
 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 typedef string::const_iterator Iterator;
 67 
 68 typedef map<string,int> molecule_type;
 69 molecule_type molecule(Iterator& it, const Iterator& last);
 70 
 71 void append(molecule_type& m, const molecule_type& mm)
 72 {
 73   for (molecule_type::const_iterator it = mm.begin(); it != mm.end(); ++it) {
 74     m[it->first] += it->second;
 75   }
 76 }
 77 
 78 int number(Iterator& it, const Iterator& last)
 79 {
 80   if (it == last || !isdigit(*it)) {
 81     return 1;
 82   }
 83   int n = 0;
 84   while (it != last && isdigit(*it)) {
 85     n = 10*n + *it-'0';
 86     ++it;
 87   }
 88   return n;
 89 }
 90 
 91 molecule_type group(Iterator& it, const Iterator& last)
 92 {
 93   molecule_type m;
 94   while (it != last && (*it == '(' || isupper(*it))) {
 95     molecule_type mm;
 96     if (*it == '(') {
 97       ++it;
 98       mm = molecule(it, last);
 99       assert(*it == ')');
100       ++it;
101     } else {
102       const Iterator jt = it;
103       ++it;
104       while (it != last && islower(*it)) {
105         ++it;
106       }
107       mm.insert(make_pair(string(jt, it), 1));
108     }
109     const int n = number(it, last);
110     for (molecule_type::iterator jt = mm.begin(); jt != mm.end(); ++jt) {
111       jt->second *= n;
112     }
113     append(m, mm);
114   }
115   return m;
116 }
117 
118 molecule_type molecule(Iterator& it, const Iterator& last)
119 {
120   molecule_type m = group(it, last);
121   while (it != last && (*it == '(' || isupper(*it))) {
122     const molecule_type mm = group(it, last);
123     append(m, mm);
124   }
125   return m;
126 }
127 
128 vector<molecule_type> molecule_sequence(Iterator& it, const Iterator& last)
129 {
130   vector<molecule_type> v(1, molecule(it, last));
131   while (it != last && *it == '+') {
132     ++it;
133     v.push_back(molecule(it, last));
134   }
135   return v;
136 }
137 
138 vector<long long> gaussian_elimination(vector<vector<ratio> > a)
139 {
140   const int N = a.size();
141   const int M = a[0].size();
142   const int K = min(N, M);
143   for (int i = 0; i < K; i++) {
144     for (int j = i+1; a[i][i].n == 0LL && j < N; j++) {
145       swap(a[i], a[j]);
146     }
147     if (a[i][i].n == 0LL) {
148       continue;
149     }
150     const ratio r = a[i][i].inverse();
151     for (int k = i; k < M; k++) {
152       a[i][k] = a[i][k] * r;
153     }
154     for (int j = 0; j < N; j++) {
155       if (j == i) {
156         continue;
157       }
158       const ratio t = a[j][i];
159       for (int k = i; k < M; k++) {
160         a[j][k] -= t * a[i][k];
161       }
162     }
163   }
164 
165   vector<long long> ans;
166   long long lcm = 1;
167   for (int i = 0; i < M-1; i++) {
168     const long long g = gcd(lcm, a[i][M-1].d);
169     lcm = (lcm * a[i][M-1].d)/g;
170   }
171   for (int i = 0; i < M-1; i++) {
172     ans.push_back(lcm/a[i][M-1].d * (-a[i][M-1].n));
173   }
174   ans.push_back(lcm);
175   return ans;
176 }
177 
178 vector<long long> solve(const string& str)
179 {
180   Iterator it = str.begin(), last = str.end();
181   const vector<molecule_type> lhs = molecule_sequence(it, last);
182   assert(*it == '-');
183   ++it;
184   assert(*it == '>');
185   ++it;
186   const vector<molecule_type> rhs = molecule_sequence(it, last);
187   assert(*it == '.');
188   ++it;
189   assert(it == last);
190 
191   map<string,int> dict;
192   for (vector<molecule_type>::const_iterator jt = lhs.begin(); jt != lhs.end(); ++jt) {
193     for (molecule_type::const_iterator kt = jt->begin(); kt != jt->end(); ++kt) {
194       if (!dict.count(kt->first)) {
195         const int id = dict.size();
196         dict.insert(make_pair(kt->first, id));
197       }
198     }
199   }
200   const int N1 = lhs.size(), N2 = rhs.size();
201   const int M = dict.size();
202   vector<vector<ratio> > a(M, vector<ratio>(N1+N2));
203   for (int i = 0; i < N1; i++) {
204     for (molecule_type::const_iterator kt = lhs[i].begin(); kt != lhs[i].end(); ++kt) {
205       a[dict[kt->first]][i] = ratio(kt->second, 1);
206     }
207   }
208   for (int i = 0; i < N2; i++) {
209     for (molecule_type::const_iterator kt = rhs[i].begin(); kt != rhs[i].end(); ++kt) {
210       a[dict[kt->first]][i+N1] = ratio(-kt->second, 1);
211     }
212   }
213 
214   return gaussian_elimination(a);
215 }
216 
217 int main()
218 {
219   string s;
220   while (getline(cin, s) && s != ".") {
221     const vector<long long> r = solve(s);
222     for (vector<long long>::const_iterator it = r.begin(); it != r.end(); ++it) {
223       if (it != r.begin()) {
224         cout << ' ';
225       }
226       cout << *it;
227     }
228     cout << endl;
229   }
230   return 0;
231 }