30
30
31
31
using PT = std::vector<size_t >;
32
32
33
+ template <typename T>
34
+ class Fraction
35
+ {
36
+ public:
37
+ Fraction (): numerator(0 ), denominator(0 ) {}
38
+ Fraction (T a): numerator(a), denominator(1 ) {}
39
+ Fraction (T a, T b): numerator(a), denominator(b) {
40
+ if (!b)
41
+ throw std::invalid_argument (" zero denominator" );
42
+ reduce_fraction ();
43
+ }
44
+ Fraction& operator = (const Fraction& rhs) {
45
+ if (&rhs == this )
46
+ return *this ;
47
+ if (!rhs.denominator )
48
+ throw std::invalid_argument (" zero denominator" );
49
+ this ->numerator = rhs.numerator ;
50
+ this ->denominator = rhs.denominator ;
51
+ return *this ;
52
+ }
53
+ const Fraction operator - () const {
54
+ if (!denominator)
55
+ throw std::invalid_argument (" zero denominator" );
56
+ return Fraction (-numerator, denominator);
57
+ }
58
+ const Fraction operator + (const Fraction& a) const {
59
+ if (!denominator || !a.denominator )
60
+ throw std::invalid_argument (" zero denominator" );
61
+ int num = numerator * a.denominator + denominator * a.numerator ;
62
+ int den = denominator * a.denominator ;
63
+ return Fraction (num, den);
64
+ }
65
+ const Fraction operator - (const Fraction& a) const {
66
+ if (!denominator || !a.denominator )
67
+ throw std::invalid_argument (" zero denominator" );
68
+ int num = numerator * a.denominator - denominator * a.numerator ;
69
+ int den = denominator * a.denominator ;
70
+ return Fraction (num, den);
71
+ }
72
+ const Fraction operator * (const Fraction& a) const {
73
+ if (!denominator || !a.denominator )
74
+ throw std::invalid_argument (" zero denominator" );
75
+ int num = numerator * a.numerator ;
76
+ int den = denominator * a.denominator ;
77
+ return Fraction (num, den);
78
+ }
79
+ const Fraction operator / (const Fraction& a) const {
80
+ if (!denominator || !a.denominator )
81
+ throw std::invalid_argument (" zero denominator" );
82
+ int num = numerator * a.denominator ;
83
+ int den = denominator * a.numerator ;
84
+ return Fraction (num, den);
85
+ }
86
+ Fraction& operator += (const Fraction& a) {
87
+ if (!denominator || !a.denominator )
88
+ throw std::invalid_argument (" zero denominator" );
89
+ int num = numerator * a.denominator + denominator * a.numerator ;
90
+ int den = denominator * a.denominator ;
91
+ numerator = num;
92
+ denominator = den;
93
+ reduce_fraction ();
94
+ return *this ;
95
+ }
96
+ Fraction& operator -= (const Fraction& a) {
97
+ if (!denominator || !a.denominator )
98
+ throw std::invalid_argument (" zero denominator" );
99
+ int num = numerator * a.denominator - denominator * a.numerator ;
100
+ int den = denominator * a.denominator ;
101
+ numerator = num;
102
+ denominator = den;
103
+ reduce_fraction ();
104
+ return *this ;
105
+ }
106
+ Fraction& operator *= (const Fraction& a) {
107
+ if (!denominator || !a.denominator )
108
+ throw std::invalid_argument (" zero denominator" );
109
+ numerator *= a.numerator ;
110
+ denominator *= a.denominator ;
111
+ reduce_fraction ();
112
+ return *this ;
113
+ }
114
+ Fraction& operator /= (const Fraction& a) {
115
+ if (!denominator || !a.denominator )
116
+ throw std::invalid_argument (" zero denominator" );
117
+ numerator *= a.denominator ;
118
+ denominator *= a.numerator ;
119
+ reduce_fraction ();
120
+ return *this ;
121
+ }
122
+ bool operator < (const Fraction& a) const {
123
+ if (!denominator || !a.denominator )
124
+ throw std::invalid_argument (" zero denominator" );
125
+ return numerator * a.denominator < denominator * a.numerator ;
126
+ }
127
+ bool operator == (const Fraction& a) const {
128
+ if (!denominator || !a.denominator )
129
+ throw std::invalid_argument (" zero denominator" );
130
+ return numerator == a.numerator && denominator == a.denominator ;
131
+ }
132
+ friend std::ostream& operator << (std::ostream& os, const Fraction& v) {
133
+ if (!v.denominator )
134
+ throw std::invalid_argument (" zero denominator" );
135
+ if (v.denominator == 1 )
136
+ return os << v.numerator ;
137
+ else
138
+ return os << v.numerator << ' /' << v.denominator ;
139
+ }
140
+ T numerator;
141
+ T denominator;
142
+ void reduce_fraction () {
143
+ // some idea from "reduce_fraction.c"
144
+ T a = numerator, b = denominator, tmp;
145
+ while (b) {
146
+ tmp = a % b;
147
+ a = b;
148
+ b = tmp;
149
+ }
150
+ // now abs(a) == abs(gcd(a, b))
151
+ if ((denominator > 0 ) ^ (a > 0 ))
152
+ a = -a;
153
+ // now a == abs(gcd(a, b)) * sgn(b)
154
+ numerator /= a;
155
+ denominator /= a;
156
+ }
157
+ };
158
+
33
159
template <typename T>
34
160
Matrix<T> LUPSolve (Matrix<T>& L, Matrix<T>& U, PT& pi, Matrix<T>& b) {
35
161
size_t n = L.rows ;
@@ -78,13 +204,13 @@ PT LUPDecomposition(Matrix<T>& A) {
78
204
T p = 0 ;
79
205
size_t kk;
80
206
for (size_t i = k; i < n; i++) {
81
- T abs = A[i][k] < 0 ? -A[i][k] : A[i][k];
82
- if (abs > p ) {
207
+ T abs = A[i][k] < (T) 0 ? -A[i][k] : A[i][k];
208
+ if (p < abs ) {
83
209
p = abs;
84
210
kk = i;
85
211
}
86
212
}
87
- if (p == 0 )
213
+ if (p == (T) 0 )
88
214
throw std::invalid_argument (" singular matrix" );
89
215
std::swap (pi[k], pi[kk]);
90
216
for (size_t i = 0 ; i < n; i++)
@@ -100,13 +226,11 @@ PT LUPDecomposition(Matrix<T>& A) {
100
226
#endif
101
227
102
228
#ifdef MAIN_LUPSolve
103
- int main (int argc, char *argv[]) {
104
- const size_t n = get_argv (argc, argv, 1 , 16 );
105
- const size_t compute = get_argv (argc, argv, 2 , 3 );
229
+ template <typename T>
230
+ void main_T (const size_t n, const size_t compute) {
106
231
std::vector<int > int_a, int_b;
107
232
random_integers (int_a, -n, n, n * n);
108
233
random_integers (int_b, -n, n, n);
109
- using T = double ;
110
234
std::vector<T> buf_a (n * n), buf_b (n);
111
235
for (size_t i = 0 ; i < int_a.size (); i++)
112
236
buf_a[i] = int_a[i];
@@ -144,6 +268,16 @@ int main(int argc, char *argv[]) {
144
268
output_integers (ans2[i], " \t " );
145
269
}
146
270
std::cout << std::endl;
271
+ }
272
+
273
+ int main (int argc, char *argv[]) {
274
+ const size_t type = get_argv (argc, argv, 1 , 0 );
275
+ const size_t n = get_argv (argc, argv, 2 , 5 );
276
+ const size_t compute = get_argv (argc, argv, 3 , 3 );
277
+ if (!type)
278
+ main_T<double >(n, compute);
279
+ else
280
+ main_T<Fraction<int >>(n, compute);
147
281
return 0 ;
148
282
}
149
283
#endif
0 commit comments