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