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 }