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