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