check/numeric_consist.cpp

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