1 /*++
2 Copyright (c) 2013 Microsoft Corporation
3
4 Module Name:
5
6 rcf.cpp
7
8 Abstract:
9
10 Testing RCF module
11
12 Author:
13
14 Leonardo (leonardo) 2013-01-04
15
16 Notes:
17
18 --*/
19 #include "math/realclosure/realclosure.h"
20 #include "math/realclosure/mpz_matrix.h"
21 #include "util/rlimit.h"
22
tst1()23 static void tst1() {
24 unsynch_mpq_manager qm;
25 reslimit rl;
26 rcmanager m(rl, qm);
27 scoped_rcnumeral a(m);
28 #if 0
29 a = 10;
30 std::cout << sym_pp(a) << std::endl;
31 std::cout << sym_pp(eps) << std::endl;
32 std::cout << interval_pp(a) << std::endl;
33 std::cout << interval_pp(eps) << std::endl;
34 #endif
35
36 scoped_rcnumeral eps(m);
37 m.mk_infinitesimal(eps);
38 mpq aux;
39 qm.set(aux, 1, 3);
40 m.set(a, aux);
41
42 #if 0
43 std::cout << interval_pp(a) << std::endl;
44 std::cout << decimal_pp(eps, 4) << std::endl;
45 std::cout << decimal_pp(a) << std::endl;
46 std::cout << a + eps << std::endl;
47 std::cout << a * eps << std::endl;
48 std::cout << (a + eps)*eps - eps << std::endl;
49 #endif
50 std::cout << interval_pp(a - eps*2) << std::endl;
51 std::cout << interval_pp(eps + 1) << std::endl;
52 scoped_rcnumeral t(m);
53 t = (a - eps*2) / (eps + 1);
54 std::cout << t << std::endl;
55 std::cout << t * (eps + 1) << std::endl;
56 a = 10;
57 std::cout << (a + eps > a) << std::endl;
58 scoped_rcnumeral pi(m);
59 m.mk_pi(pi);
60 std::cout << pi + 1 << std::endl;
61 std::cout << decimal_pp(pi) << std::endl;
62 std::cout << decimal_pp(pi + 1) << std::endl;
63 scoped_rcnumeral e(m);
64 m.mk_e(e);
65 t = e + (pi + 1)*2;
66 std::cout << t << std::endl;
67 std::cout << decimal_pp(t, 10) << std::endl;
68 std::cout << (eps + 1 > 1) << std::endl;
69 std::cout << interval_pp((a + eps)/(a - eps)) << std::endl;
70 }
71
tst2()72 static void tst2() {
73 enable_trace("mpz_matrix");
74 unsynch_mpq_manager nm;
75 small_object_allocator allocator;
76 mpz_matrix_manager mm(nm, allocator);
77 scoped_mpz_matrix A(mm);
78 mm.mk(3, 3, A);
79 // Matrix
80 // 1 1 1
81 // 0 1 -1
82 // 0 1 1
83 A.set(0, 0, 1); A.set(0, 1, 1); A.set(0, 2, 1);
84 A.set(1, 0, 0); A.set(1, 1, 1); A.set(1, 2, -1);
85 A.set(2, 0, 0); A.set(2, 1, 1); A.set(2, 2, 1);
86 std::cout << A;
87 {
88 int b[3];
89 int c[3] = { 10, -2, 8 };
90 std::cout << "solve: " << mm.solve(A, b, c) << "\n";
91 for (unsigned i = 0; i < 3; i++) std::cout << b[i] << " ";
92 std::cout << "\n";
93 }
94 scoped_mpz_matrix A2(mm);
95 mm.tensor_product(A, A, A2);
96 std::cout << A2;
97 scoped_mpz_matrix B(mm);
98 unsigned cols[] = { 1, 3, 7, 8 };
99 mm.filter_cols(A2, 4, cols, B);
100 std::cout << B;
101 scoped_mpz_matrix C(mm);
102 unsigned perm[] = { 8, 7, 6, 5, 4, 3, 2, 1, 0 };
103 mm.permute_rows(B, perm, C);
104 std::cout << C;
105 }
106
tst_solve(unsigned n,int _A[],int _b[],int _c[],bool solved)107 static void tst_solve(unsigned n, int _A[], int _b[], int _c[], bool solved) {
108 unsynch_mpq_manager nm;
109 small_object_allocator allocator;
110 mpz_matrix_manager mm(nm, allocator);
111 scoped_mpz_matrix A(mm);
112 mm.mk(n, n, A);
113 for (unsigned i = 0; i < n; i++)
114 for (unsigned j = 0; j < n; j++)
115 A.set(i, j, _A[i*n + j]);
116 svector<int> b;
117 b.resize(n, 0);
118 if (mm.solve(A, b.data(), _c)) {
119 ENSURE(solved);
120 for (unsigned i = 0; i < n; i++) {
121 ENSURE(b[i] == _b[i]);
122 }
123 }
124 else {
125 ENSURE(!solved);
126 }
127 }
128
tst_lin_indep(unsigned m,unsigned n,int _A[],unsigned ex_sz,unsigned ex_r[])129 static void tst_lin_indep(unsigned m, unsigned n, int _A[], unsigned ex_sz, unsigned ex_r[]) {
130 unsynch_mpq_manager nm;
131 small_object_allocator allocator;
132 mpz_matrix_manager mm(nm, allocator);
133 scoped_mpz_matrix A(mm);
134 mm.mk(m, n, A);
135 for (unsigned i = 0; i < m; i++)
136 for (unsigned j = 0; j < n; j++)
137 A.set(i, j, _A[i*n + j]);
138 unsigned_vector r;
139 r.resize(A.n());
140 scoped_mpz_matrix B(mm);
141 mm.linear_independent_rows(A, r.data(), B);
142 for (unsigned i = 0; i < ex_sz; i++) {
143 ENSURE(r[i] == ex_r[i]);
144 }
145 }
146
tst_denominators()147 static void tst_denominators() {
148 reslimit rl;
149 unsynch_mpq_manager qm;
150 rcmanager m(rl, qm);
151 scoped_rcnumeral a(m);
152 scoped_rcnumeral t(m);
153 scoped_rcnumeral eps(m);
154 m.mk_pi(a);
155 m.inv(a);
156 m.mk_infinitesimal(eps);
157 t = (a - eps*2) / (a*eps + 1);
158 // t = t + a * 2;
159 scoped_rcnumeral n(m), d(m);
160 std::cout << t << "\n";
161 m.clean_denominators(t, n, d);
162 std::cout << "---->\n" << n << "\n" << d << "\n";
163 }
164
tst_rcf()165 void tst_rcf() {
166 enable_trace("rcf_clean");
167 enable_trace("rcf_clean_bug");
168 tst_denominators();
169 tst1();
170 tst2();
171 { int A[] = {0, 1, 1, 1, 0, 1, 1, 1, -1}; int c[] = {10, 4, -4}; int b[] = {-2, 4, 6}; tst_solve(3, A, b, c, true); }
172 { int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1}; int c[] = {3, 2, 2}; int b[] = {1, 1, 1}; tst_solve(3, A, b, c, false); }
173 { int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, -1}; unsigned r[] = {0, 1, 4}; tst_lin_indep(5, 3, A, 3, r); }
174 { int A[] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, -1}; unsigned r[] = {0, 4}; tst_lin_indep(5, 3, A, 2, r); }
175 { int A[] = {1, 1, 1, 1, 1, 1, 1, 0, 1, 2, 1, 2, 3, 1, 3}; unsigned r[] = {0, 2}; tst_lin_indep(5, 3, A, 2, r); }
176 }
177