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