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