- enforced GiNaC coding standards :-)
[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  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include <stdlib.h>
24 #include <ginac/ginac.h>
25
26 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
27  * conversions. */
28 static unsigned numeric_consist1(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     symbol a("a");
36     ex e1, e2;
37     
38     if ( !test_int1.is_integer() ) {
39         clog << test_int1
40              << " erroneously not recognized as integer" << endl;
41         ++result;
42     }
43     if ( !test_int1.is_rational() ) {
44         clog << test_int1
45              << " erroneously not recognized as rational" << endl;
46         ++result;
47     }
48     
49     if ( !test_rat1.is_rational() ) {
50         clog << test_rat1
51              << " erroneously not recognized as rational" << endl;
52         ++result;
53     }
54     if ( test_rat1.is_integer() ) {
55         clog << test_rat1
56              << " erroneously recognized as integer" << endl;
57         ++result;
58     }
59     
60     int i = numeric(1984).to_int();
61     if ( i-1984 ) {
62         clog << "conversion of " << i
63              << " from numeric to int failed" << endl;
64         ++result;
65     }
66     
67     e1 = test_int1;
68     if ( !e1.info(info_flags::posint) ) {
69         clog << "expression " << e1
70              << " erroneously not recognized as positive integer" << endl;
71         ++result;
72     }
73     
74     e2 = test_int1 + a;
75     if ( ex_to_numeric(e2).is_integer() ) {
76         clog << "expression " << e2
77              << " erroneously recognized as integer" << endl;
78         ++result;
79     }
80     
81     // The next two were two actual bugs in CLN till June, 12, 1999:
82     test_rat1 = numeric(3)/numeric(2);
83     test_rat1 += test_rat1;
84     if ( !test_rat1.is_integer() ) {
85         clog << "3/2 + 3/2 erroneously not integer 3 but instead "
86              << test_rat1 << endl;
87         ++result;
88     }
89     test_rat1 = numeric(3)/numeric(2);
90     numeric test_rat2 = test_rat1 + numeric(1);  // 5/2
91     test_rat2 -= test_rat1;             // 1
92     if ( !test_rat2.is_integer() ) {
93         clog << "5/2 - 3/2 erroneously not integer 1 but instead "
94              << test_rat2 << endl;
95         ++result;
96     }
97     
98     // Check some numerator and denominator calculations:
99     for (int i=0; i<10; ++i) {
100         int re_q, im_q;
101         do { re_q = rand(); } while (re_q == 0);
102         do { im_q = rand(); } while (im_q == 0);
103         numeric r(rand()-RAND_MAX/2, re_q);
104         numeric i(rand()-RAND_MAX/2, im_q);
105         numeric z = r + I*i;
106         numeric p = numer(z);
107         numeric q = denom(z);
108         numeric res = p/q;
109         if (res != z) {
110             clog << z << " erroneously transformed into " 
111                  << p << "/" << q << " by numer() and denom()" << 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 bug was finally killed in CLN on September 2nd. */
123 static unsigned numeric_consist2(void)
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.compare(exZERO()) ) {
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 numeric_consist3(void)
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 numeric_consist4(void)
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         }
291     }
292     if ( !passed ) {
293         clog << "One or more square roots of squares of integers did not return exact integers" << endl;
294         ++result;
295     }
296     // square roots of squares of rationals:
297     passed = true;
298     for (int num=0; num<41; ++num) {
299         for (int den=1; den<42; ++den) {
300             if ( !sqrt(numeric(num*num)/numeric(den*den)).is_rational() ) {
301                 passed = false;
302             }
303         }
304     }
305     if ( !passed ) {
306         clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
307         ++result;
308     }
309     
310     return result;
311 }
312
313 unsigned numeric_consist(void)
314 {
315     unsigned result = 0;
316
317     cout << "checking consistency of numeric types..." << flush;
318     clog << "---------consistency of numeric types:" << endl;
319     
320     result += numeric_consist1();
321     result += numeric_consist2();
322     result += numeric_consist3();
323     result += numeric_consist4();
324
325     if ( !result ) {
326         cout << " passed ";
327         clog << "(no output)" << endl;
328     } else {
329         cout << " failed ";
330     }
331     
332     return result;
333 }