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