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 }