- changed behaviour of numeric::is_rational() and added numeric::is_cinteger()
[ginac.git] / check / numeric_consist.cpp
1 /** @file numeric_consist.cpp
2  *
3  *  This test routine creates some numbers and check the result of several
4  *  boolean tests on these numbers like is_integer() etc... */
5
6 /*
7  *  GiNaC Copyright (C) 1999 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 <stdlib.h>
25 #include <ginac/ginac.h>
26
27 #ifndef NO_GINAC_NAMESPACE
28 using namespace GiNaC;
29 #endif // ndef NO_GINAC_NAMESPACE
30
31 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
32  * conversions. */
33 static unsigned numeric_consist1(void)
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 (ex_to_numeric(e2).is_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     // Check some numerator and denominator calculations:
111     for (int i=0; i<10; ++i) {
112         int re_q, im_q;
113         do { re_q = rand(); } while (re_q == 0);
114         do { im_q = rand(); } while (im_q == 0);
115         numeric r(rand()-RAND_MAX/2, re_q);
116         numeric i(rand()-RAND_MAX/2, im_q);
117         numeric z = r + I*i;
118         numeric p = numer(z);
119         numeric q = denom(z);
120         numeric res = p/q;
121         if (res != z) {
122             clog << z << " erroneously transformed into " 
123                  << p << "/" << q << " by numer() and denom()" << endl;
124             ++result;
125         }
126     }    
127     return result;
128 }
129
130 /* We had some fun with a bug in CLN that caused it to loop forever when
131  * calculating expt(a,b) if b is a rational and a a nonnegative integer.
132  * Implementing a workaround sadly introduced another bug on May 28th 1999
133  * that was fixed on May 31st.  The workaround turned out to be stupid and
134  * the original bug in CLN was finally killed on September 2nd. */
135 static unsigned numeric_consist2(void)
136 {
137     unsigned result = 0;
138     
139     ex zero = numeric(0);
140     ex two = numeric(2);
141     ex three = numeric(3);
142     
143     // The hang in this code was the reason for the original workaround
144     if (pow(two,two/three)==42) {
145         clog << "pow(2,2/3) erroneously returned 42" << endl;
146         ++result;  // cannot happen
147     }
148     
149     // Actually, this used to raise a FPE after introducing the workaround
150     if (two*zero!=zero) {
151         clog << "2*0 erroneously returned " << two*zero << endl;
152         ++result;
153     }
154     
155     // And this returned a cl_F due to the implicit call of numeric::power()
156     ex six = two*three;
157     if (!six.info(info_flags::integer)) {
158         clog << "2*3 erroneously returned the non-integer " << six << endl;
159         ++result;
160     }
161     
162     // The fix in the workaround left a whole which was fixed hours later...
163     ex another_zero = pow(zero,numeric(1)/numeric(2));
164     if (another_zero.compare(exZERO())) {
165         clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
166         ++result;
167     }
168     
169     return result;
170 }
171
172 /* Assorted tests to ensure some crucial functions behave exactly as specified
173  * in the documentation. */
174 static unsigned numeric_consist3(void)
175 {
176     unsigned result = 0;
177     numeric calc_rem, calc_quo;
178     numeric a, b;
179     
180     // check if irem(a, b) and irem(a, b, q) really behave like Maple's 
181     // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
182     // These overloaded routines indeed need to be checked separately since
183     // internally they might be doing something completely different:
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     a = -23; b = 4; calc_rem = irem(a, b);
197     if (calc_rem != -3) {
198         clog << "irem(" << a << "," << b << ") erroneously returned "
199              << calc_rem << endl;
200         ++result;
201     }
202     a = -23; b = -4; calc_rem = irem(a, b);
203     if (calc_rem != -3) {
204         clog << "irem(" << a << "," << b << ") erroneously returned "
205              << calc_rem << endl;
206         ++result;
207     }
208     // and now the overloaded irem(a,b,q):
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     a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
222     if (calc_rem != -3 || calc_quo != -5) {
223         clog << "irem(" << a << "," << b << ",q) erroneously returned "
224              << calc_rem << " with q=" << calc_quo << endl;
225         ++result;
226     }
227     a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
228     if (calc_rem != -3 || calc_quo != 5) {
229         clog << "irem(" << a << "," << b << ",q) erroneously returned "
230              << calc_rem << " with q=" << calc_quo << endl;
231         ++result;
232     }
233     // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's 
234     // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
235     // These overloaded routines indeed need to be checked separately since
236     // internally they might be doing something completely different:
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     a = -23; b = 4; calc_quo = iquo(a, b);
250     if (calc_quo != -5) {
251         clog << "iquo(" << a << "," << b << ") erroneously returned "
252              << calc_quo << endl;
253         ++result;
254     }
255     a = -23; b = -4; calc_quo = iquo(a, b);
256     if (calc_quo != 5) {
257         clog << "iquo(" << a << "," << b << ") erroneously returned "
258              << calc_quo << endl;
259         ++result;
260     }
261     // and now the overloaded iquo(a,b,r):
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     a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
275     if (calc_quo != -5 || calc_rem != -3) {
276         clog << "iquo(" << a << "," << b << ",r) erroneously returned "
277              << calc_quo << " with r=" << calc_rem << endl;
278         ++result;
279     }
280     a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
281     if (calc_quo != 5 || calc_rem != -3) {
282         clog << "iquo(" << a << "," << b << ",r) erroneously returned "
283              << calc_quo << " with r=" << calc_rem << endl;
284         ++result;
285     }
286     
287     return result;
288 }
289
290 /* Now we perform some less trivial checks about several functions which should
291  * return exact numbers if possible. */
292 static unsigned numeric_consist4(void)
293 {
294     unsigned result = 0;
295     bool passed;
296     
297     // square roots of squares of integers:
298     passed = true;
299     for (int i=0; i<42; ++i) {
300         if ( !sqrt(numeric(i*i)).is_integer() ) {
301             passed = false;
302         }
303     }
304     if (!passed) {
305         clog << "One or more square roots of squares of integers did not return exact integers" << endl;
306         ++result;
307     }
308     
309     // square roots of squares of rationals:
310     passed = true;
311     for (int num=0; num<41; ++num) {
312         for (int den=1; den<42; ++den) {
313             if ( !sqrt(numeric(num*num)/numeric(den*den)).is_rational() ) {
314                 passed = false;
315             }
316         }
317     }
318     if (!passed) {
319         clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
320         ++result;
321     }
322     
323     return result;
324 }
325
326 unsigned numeric_consist(void)
327 {
328     unsigned result = 0;
329
330     cout << "checking consistency of numeric types..." << flush;
331     clog << "---------consistency of numeric types:" << endl;
332     
333     result += numeric_consist1();
334     result += numeric_consist2();
335     result += numeric_consist3();
336     result += numeric_consist4();
337
338     if (!result) {
339         cout << " passed ";
340         clog << "(no output)" << endl;
341     } else {
342         cout << " failed ";
343     }
344     
345     return result;
346 }