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