715acff9d15783f57c3db1ccc52936ff5ad89113
[ginac.git] / check / exam_numeric.cpp
1 /** @file exam_numeric.cpp
2  *
3  *  These exams creates some numbers and check the result of several boolean
4  *  tests on these numbers like is_integer() etc... */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  */
23
24 #include "ginac.h"
25 using namespace GiNaC;
26
27 #include <iostream>
28 #include <sstream>
29 using namespace std;
30
31 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
32  * conversions. */
33 static unsigned exam_numeric1()
34 {
35         unsigned result = 0;
36         numeric test_int1(42);
37         numeric test_int2(5);
38         numeric test_rat1 = test_int1; test_rat1 /= test_int2;
39         test_rat1 = -test_rat1;          // -42/5
40         numeric test_crat = test_rat1+I*test_int2;  // 5*I-42/5
41         symbol a("a");
42         ex e1, e2;
43         
44         if (!test_int1.is_integer()) {
45                 clog << test_int1
46                      << " erroneously not recognized as integer" << endl;
47                 ++result;
48         }
49         if (!test_int1.is_rational()) {
50                 clog << test_int1
51                      << " erroneously not recognized as rational" << endl;
52                 ++result;
53         }
54         
55         if (!test_rat1.is_rational()) {
56                 clog << test_rat1
57                      << " erroneously not recognized as rational" << endl;
58                 ++result;
59         }
60         if (test_rat1.is_integer()) {
61                 clog << test_rat1
62                          << " erroneously recognized as integer" << endl;
63                 ++result;
64         }
65         
66         if (!test_crat.is_crational()) {
67                 clog << test_crat
68                      << " erroneously not recognized as complex rational" << endl;
69                 ++result;
70         }
71         
72         int i = numeric(1984).to_int();
73         if (i-1984) {
74                 clog << "conversion of " << i
75                      << " from numeric to int failed" << endl;
76                 ++result;
77         }
78         
79         e1 = test_int1;
80         if (!e1.info(info_flags::posint)) {
81                 clog << "expression " << e1
82                      << " erroneously not recognized as positive integer" << endl;
83                 ++result;
84         }
85         
86         e2 = test_int1 + a;
87         if (e2.info(info_flags::integer)) {
88                 clog << "expression " << e2
89                      << " erroneously recognized as integer" << endl;
90                 ++result;
91         }
92         
93         // The next two were two actual bugs in CLN till June, 12, 1999:
94         test_rat1 = numeric(3)/numeric(2);
95         test_rat1 += test_rat1;
96         if (!test_rat1.is_integer()) {
97                 clog << "3/2 + 3/2 erroneously not integer 3 but instead "
98                      << test_rat1 << endl;
99                 ++result;
100         }
101         test_rat1 = numeric(3)/numeric(2);
102         numeric test_rat2 = test_rat1 + numeric(1);  // 5/2
103         test_rat2 -= test_rat1;  // 1
104         if (!test_rat2.is_integer()) {
105                 clog << "5/2 - 3/2 erroneously not integer 1 but instead "
106                      << test_rat2 << endl;
107                 ++result;
108         }
109         
110         return result;
111 }
112
113 /* We had some fun with a bug in CLN that caused it to loop forever when
114  * calculating expt(a,b) if b is a rational and a a nonnegative integer.
115  * Implementing a workaround sadly introduced another bug on May 28th 1999
116  * that was fixed on May 31st.  The workaround turned out to be stupid and
117  * the original bug in CLN was finally killed on September 2nd. */
118 static unsigned exam_numeric2()
119 {
120         unsigned result = 0;
121         
122         ex zero = numeric(0);
123         ex two = numeric(2);
124         ex three = numeric(3);
125         
126         // The hang in this code was the reason for the original workaround
127         if (pow(two,two/three)==42) {
128                 clog << "pow(2,2/3) erroneously returned 42" << endl;
129                 ++result;  // cannot happen
130         }
131         
132         // Actually, this used to raise a FPE after introducing the workaround
133         if (two*zero!=zero) {
134                 clog << "2*0 erroneously returned " << two*zero << endl;
135                 ++result;
136         }
137         
138         // And this returned a cl_F due to the implicit call of numeric::power()
139         ex six = two*three;
140         if (!six.info(info_flags::integer)) {
141                 clog << "2*3 erroneously returned the non-integer " << six << endl;
142                 ++result;
143         }
144         
145         // The fix in the workaround left a whole which was fixed hours later...
146         ex another_zero = pow(zero,numeric(1)/numeric(2));
147         if (!another_zero.is_zero()) {
148                 clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
149                 ++result;
150         }
151         
152         return result;
153 }
154
155 /* Assorted tests to ensure some crucial functions behave exactly as specified
156  * in the documentation. */
157 static unsigned exam_numeric3()
158 {
159         unsigned result = 0;
160         numeric calc_rem, calc_quo;
161         numeric a, b;
162         
163         // check if irem(a, b) and irem(a, b, q) really behave like Maple's 
164         // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
165         // These overloaded routines indeed need to be checked separately since
166         // internally they might be doing something completely different:
167         a = 23; b = 4; calc_rem = irem(a, b);
168         if (calc_rem != 3) {
169                 clog << "irem(" << a << "," << b << ") erroneously returned "
170                      << calc_rem << endl;
171                 ++result;
172         }
173         a = 23; b = -4; calc_rem = irem(a, b);
174         if (calc_rem != 3) {
175                 clog << "irem(" << a << "," << b << ") erroneously returned "
176                          << calc_rem << endl;
177                 ++result;
178         }
179         a = -23; b = 4; calc_rem = irem(a, b);
180         if (calc_rem != -3) {
181                 clog << "irem(" << a << "," << b << ") erroneously returned "
182                      << calc_rem << endl;
183                 ++result;
184         }
185         a = -23; b = -4; calc_rem = irem(a, b);
186         if (calc_rem != -3) {
187                 clog << "irem(" << a << "," << b << ") erroneously returned "
188                      << calc_rem << endl;
189                 ++result;
190         }
191         // and now the overloaded irem(a,b,q):
192         a = 23; b = 4; calc_rem = irem(a, b, calc_quo);
193         if (calc_rem != 3 || calc_quo != 5) {
194                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
195                      << calc_rem << " with q=" << calc_quo << endl;
196                 ++result;
197         }
198         a = 23; b = -4; calc_rem = irem(a, b, calc_quo);
199         if (calc_rem != 3 || calc_quo != -5) {
200                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
201                      << calc_rem << " with q=" << calc_quo << endl;
202                 ++result;
203         }
204         a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
205         if (calc_rem != -3 || calc_quo != -5) {
206                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
207                      << calc_rem << " with q=" << calc_quo << endl;
208                 ++result;
209         }
210         a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
211         if (calc_rem != -3 || calc_quo != 5) {
212                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
213                      << calc_rem << " with q=" << calc_quo << endl;
214                 ++result;
215         }
216         // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's 
217         // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
218         // These overloaded routines indeed need to be checked separately since
219         // internally they might be doing something completely different:
220         a = 23; b = 4; calc_quo = iquo(a, b);
221         if (calc_quo != 5) {
222                 clog << "iquo(" << a << "," << b << ") erroneously returned "
223                      << calc_quo << endl;
224                 ++result;
225         }
226         a = 23; b = -4; calc_quo = iquo(a, b);
227         if (calc_quo != -5) {
228                 clog << "iquo(" << a << "," << b << ") erroneously returned "
229                      << calc_quo << endl;
230                 ++result;
231         }
232         a = -23; b = 4; calc_quo = iquo(a, b);
233         if (calc_quo != -5) {
234                 clog << "iquo(" << a << "," << b << ") erroneously returned "
235                      << calc_quo << endl;
236                 ++result;
237         }
238         a = -23; b = -4; calc_quo = iquo(a, b);
239         if (calc_quo != 5) {
240                 clog << "iquo(" << a << "," << b << ") erroneously returned "
241                      << calc_quo << endl;
242                 ++result;
243         }
244         // and now the overloaded iquo(a,b,r):
245         a = 23; b = 4; calc_quo = iquo(a, b, calc_rem);
246         if (calc_quo != 5 || calc_rem != 3) {
247                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
248                      << calc_quo << " with r=" << calc_rem << endl;
249                 ++result;
250         }
251         a = 23; b = -4; calc_quo = iquo(a, b, calc_rem);
252         if (calc_quo != -5 || calc_rem != 3) {
253                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
254                      << calc_quo << " with r=" << calc_rem << endl;
255                 ++result;
256         }
257         a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
258         if (calc_quo != -5 || calc_rem != -3) {
259                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
260                      << calc_quo << " with r=" << calc_rem << endl;
261                 ++result;
262         }
263         a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
264         if (calc_quo != 5 || calc_rem != -3) {
265                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
266                      << calc_quo << " with r=" << calc_rem << endl;
267                 ++result;
268         }
269         
270         return result;
271 }
272
273 /* Now we perform some less trivial checks about several functions which should
274  * return exact numbers if possible. */
275 static unsigned exam_numeric4()
276 {
277         unsigned result = 0;
278         bool passed;
279         
280         // square roots of squares of integers:
281         passed = true;
282         for (int i=0; i<42; ++i)
283                 if (!sqrt(numeric(i*i)).is_integer())
284                         passed = false;
285         if (!passed) {
286                 clog << "One or more square roots of squares of integers did not return exact integers" << endl;
287                 ++result;
288         }
289         
290         // square roots of squares of rationals:
291         passed = true;
292         for (int num=0; num<41; ++num)
293                 for (int den=1; den<42; ++den)
294                         if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational())
295                                 passed = false;
296         if (!passed) {
297                 clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
298                 ++result;
299         }
300         
301         return result;
302 }
303
304 /* This test examines that simplifications of the form 5^(3/2) -> 5*5^(1/2)
305  * are carried out properly. */
306 static unsigned exam_numeric5()
307 {
308         unsigned result = 0;
309         
310         // A variation of one of Ramanujan's wonderful identities must be
311         // verifiable with very primitive means:
312         ex e1 = pow(1 + pow(3,numeric(1,5)) - pow(3,numeric(2,5)),3);
313         ex e2 = expand(e1 - 10 + 5*pow(3,numeric(3,5)));
314         if (!e2.is_zero()) {
315                 clog << "expand((1+3^(1/5)-3^(2/5))^3-10+5*3^(3/5)) returned "
316                      << e2 << " instead of 0." << endl;
317                 ++result;
318         }
319         
320         return result;
321 }
322
323 /* This test checks whether the numeric output/parsing routines are
324    consistent. */
325 static unsigned exam_numeric6()
326 {
327         unsigned result = 0;
328
329         symbol sym("sym");
330         vector<ex> test_numbers;
331         test_numbers.push_back(numeric(0));                     // zero
332         test_numbers.push_back(numeric(1));                     // one
333         test_numbers.push_back(numeric(-1));            // minus one
334         test_numbers.push_back(numeric(42));            // positive integer
335         test_numbers.push_back(numeric(-42));           // negative integer
336         test_numbers.push_back(numeric(14,3));          // positive rational
337         test_numbers.push_back(numeric(-14,3));         // negative rational
338         test_numbers.push_back(numeric(3.141));         // positive decimal
339         test_numbers.push_back(numeric(-3.141));        // negative decimal
340         test_numbers.push_back(numeric(0.1974));        // positive decimal, leading zero
341         test_numbers.push_back(numeric(-0.1974));       // negative decimal, leading zero
342         test_numbers.push_back(sym);                            // symbol
343
344         for (vector<ex>::const_iterator br=test_numbers.begin(); br<test_numbers.end(); ++br) {
345                 for (vector<ex>::const_iterator bi=test_numbers.begin(); bi<test_numbers.end(); ++bi) {
346
347                         for (vector<ex>::const_iterator er=test_numbers.begin(); er<test_numbers.end(); ++er) {
348                                 for (vector<ex>::const_iterator ei=test_numbers.begin(); ei<test_numbers.end(); ++ei) {
349
350                                         // Construct expression, don't test invalid ones
351                                         ex base = (*br) + (*bi)*I, exponent = (*er) + (*ei)*I, x;
352                                         try {
353                                                 x = pow(base, exponent);
354                                         } catch (...) {
355                                                 continue;
356                                         }
357
358                                         // Print to string
359                                         std::ostringstream s;
360                                         s << x;
361
362                                         // Read back expression from string
363                                         string x_as_string = s.str();
364                                         ex x_again(x_as_string, lst(sym));
365
366                                         // They should be equal
367                                         if (!x_again.is_equal(x)) {
368                                                 clog << x << " was read back as " << x_again << endl;
369                                                 ++result;
370                                         }
371                                 }
372                         }
373                 }
374         }
375
376         return result;
377 }
378
379 unsigned exam_numeric()
380 {
381         unsigned result = 0;
382         
383         cout << "examining consistency of numeric types" << flush;
384         
385         result += exam_numeric1();  cout << '.' << flush;
386         result += exam_numeric2();  cout << '.' << flush;
387         result += exam_numeric3();  cout << '.' << flush;
388         result += exam_numeric4();  cout << '.' << flush;
389         result += exam_numeric5();  cout << '.' << flush;
390         result += exam_numeric6();  cout << '.' << flush;
391         
392         return result;
393 }
394
395 int main(int argc, char** argv)
396 {
397         return exam_numeric();
398 }