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 }