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