check/exam_numeric.cpp

   1 /** @file exam_numeric.cpp
   2  *
   3  *  These exams creates some numbers and check the result of several boolean
   4  *  tests on these numbers like is_integer() etc... */
   5
   6 /*
   7  *  GiNaC Copyright (C) 1999-2003 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 "exams.h"
  25
  26 #include <sstream>
  27
  28 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
  29  * conversions. */
  30 static unsigned exam_numeric1(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         numeric test_crat = test_rat1+I*test_int2;  // 5*I-42/5
  38         symbol a("a");
  39         ex e1, e2;
  40
  41         if (!test_int1.is_integer()) {
  42                 clog << test_int1
  43                      << " erroneously not recognized as integer" << endl;
  44                 ++result;
  45         }
  46         if (!test_int1.is_rational()) {
  47                 clog << test_int1
  48                      << " erroneously not recognized as rational" << endl;
  49                 ++result;
  50         }
  51
  52         if (!test_rat1.is_rational()) {
  53                 clog << test_rat1
  54                      << " erroneously not recognized as rational" << endl;
  55                 ++result;
  56         }
  57         if (test_rat1.is_integer()) {
  58                 clog << test_rat1
  59                          << " erroneously recognized as integer" << endl;
  60                 ++result;
  61         }
  62
  63         if (!test_crat.is_crational()) {
  64                 clog << test_crat
  65                      << " erroneously not recognized as complex rational" << endl;
  66                 ++result;
  67         }
  68
  69         int i = numeric(1984).to_int();
  70         if (i-1984) {
  71                 clog << "conversion of " << i
  72                      << " from numeric to int failed" << endl;
  73                 ++result;
  74         }
  75
  76         e1 = test_int1;
  77         if (!e1.info(info_flags::posint)) {
  78                 clog << "expression " << e1
  79                      << " erroneously not recognized as positive integer" << endl;
  80                 ++result;
  81         }
  82
  83         e2 = test_int1 + a;
  84         if (e2.info(info_flags::integer)) {
  85                 clog << "expression " << e2
  86                      << " erroneously recognized as integer" << endl;
  87                 ++result;
  88         }
  89
  90         // The next two were two actual bugs in CLN till June, 12, 1999:
  91         test_rat1 = numeric(3)/numeric(2);
  92         test_rat1 += test_rat1;
  93         if (!test_rat1.is_integer()) {
  94                 clog << "3/2 + 3/2 erroneously not integer 3 but instead "
  95                      << test_rat1 << endl;
  96                 ++result;
  97         }
  98         test_rat1 = numeric(3)/numeric(2);
  99         numeric test_rat2 = test_rat1 + numeric(1);  // 5/2
 100         test_rat2 -= test_rat1;  // 1
 101         if (!test_rat2.is_integer()) {
 102                 clog << "5/2 - 3/2 erroneously not integer 1 but instead "
 103                      << test_rat2 << endl;
 104                 ++result;
 105         }
 106
 107         return result;
 108 }
 109
 110 /* We had some fun with a bug in CLN that caused it to loop forever when
 111  * calculating expt(a,b) if b is a rational and a a nonnegative integer.
 112  * Implementing a workaround sadly introduced another bug on May 28th 1999
 113  * that was fixed on May 31st.  The workaround turned out to be stupid and
 114  * the original bug in CLN was finally killed on September 2nd. */
 115 static unsigned exam_numeric2(void)
 116 {
 117         unsigned result = 0;
 118
 119         ex zero = numeric(0);
 120         ex two = numeric(2);
 121         ex three = numeric(3);
 122
 123         // The hang in this code was the reason for the original workaround
 124         if (pow(two,two/three)==42) {
 125                 clog << "pow(2,2/3) erroneously returned 42" << endl;
 126                 ++result;  // cannot happen
 127         }
 128
 129         // Actually, this used to raise a FPE after introducing the workaround
 130         if (two*zero!=zero) {
 131                 clog << "2*0 erroneously returned " << two*zero << endl;
 132                 ++result;
 133         }
 134
 135         // And this returned a cl_F due to the implicit call of numeric::power()
 136         ex six = two*three;
 137         if (!six.info(info_flags::integer)) {
 138                 clog << "2*3 erroneously returned the non-integer " << six << endl;
 139                 ++result;
 140         }
 141
 142         // The fix in the workaround left a whole which was fixed hours later...
 143         ex another_zero = pow(zero,numeric(1)/numeric(2));
 144         if (!another_zero.is_zero()) {
 145                 clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
 146                 ++result;
 147         }
 148
 149         return result;
 150 }
 151
 152 /* Assorted tests to ensure some crucial functions behave exactly as specified
 153  * in the documentation. */
 154 static unsigned exam_numeric3(void)
 155 {
 156         unsigned result = 0;
 157         numeric calc_rem, calc_quo;
 158         numeric a, b;
 159
 160         // check if irem(a, b) and irem(a, b, q) really behave like Maple's
 161         // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
 162         // These overloaded routines indeed need to be checked separately since
 163         // internally they might be doing something completely different:
 164         a = 23; b = 4; calc_rem = irem(a, b);
 165         if (calc_rem != 3) {
 166                 clog << "irem(" << a << "," << b << ") erroneously returned "
 167                      << calc_rem << endl;
 168                 ++result;
 169         }
 170         a = 23; b = -4; calc_rem = irem(a, b);
 171         if (calc_rem != 3) {
 172                 clog << "irem(" << a << "," << b << ") erroneously returned "
 173                          << calc_rem << endl;
 174                 ++result;
 175         }
 176         a = -23; b = 4; calc_rem = irem(a, b);
 177         if (calc_rem != -3) {
 178                 clog << "irem(" << a << "," << b << ") erroneously returned "
 179                      << calc_rem << endl;
 180                 ++result;
 181         }
 182         a = -23; b = -4; calc_rem = irem(a, b);
 183         if (calc_rem != -3) {
 184                 clog << "irem(" << a << "," << b << ") erroneously returned "
 185                      << calc_rem << endl;
 186                 ++result;
 187         }
 188         // and now the overloaded irem(a,b,q):
 189         a = 23; b = 4; calc_rem = irem(a, b, calc_quo);
 190         if (calc_rem != 3 || calc_quo != 5) {
 191                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
 192                      << calc_rem << " with q=" << calc_quo << endl;
 193                 ++result;
 194         }
 195         a = 23; b = -4; calc_rem = irem(a, b, calc_quo);
 196         if (calc_rem != 3 || calc_quo != -5) {
 197                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
 198                      << calc_rem << " with q=" << calc_quo << endl;
 199                 ++result;
 200         }
 201         a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
 202         if (calc_rem != -3 || calc_quo != -5) {
 203                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
 204                      << calc_rem << " with q=" << calc_quo << endl;
 205                 ++result;
 206         }
 207         a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
 208         if (calc_rem != -3 || calc_quo != 5) {
 209                 clog << "irem(" << a << "," << b << ",q) erroneously returned "
 210                      << calc_rem << " with q=" << calc_quo << endl;
 211                 ++result;
 212         }
 213         // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's
 214         // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
 215         // These overloaded routines indeed need to be checked separately since
 216         // internally they might be doing something completely different:
 217         a = 23; b = 4; calc_quo = iquo(a, b);
 218         if (calc_quo != 5) {
 219                 clog << "iquo(" << a << "," << b << ") erroneously returned "
 220                      << calc_quo << endl;
 221                 ++result;
 222         }
 223         a = 23; b = -4; calc_quo = iquo(a, b);
 224         if (calc_quo != -5) {
 225                 clog << "iquo(" << a << "," << b << ") erroneously returned "
 226                      << calc_quo << endl;
 227                 ++result;
 228         }
 229         a = -23; b = 4; calc_quo = iquo(a, b);
 230         if (calc_quo != -5) {
 231                 clog << "iquo(" << a << "," << b << ") erroneously returned "
 232                      << calc_quo << endl;
 233                 ++result;
 234         }
 235         a = -23; b = -4; calc_quo = iquo(a, b);
 236         if (calc_quo != 5) {
 237                 clog << "iquo(" << a << "," << b << ") erroneously returned "
 238                      << calc_quo << endl;
 239                 ++result;
 240         }
 241         // and now the overloaded iquo(a,b,r):
 242         a = 23; b = 4; calc_quo = iquo(a, b, calc_rem);
 243         if (calc_quo != 5 || calc_rem != 3) {
 244                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
 245                      << calc_quo << " with r=" << calc_rem << endl;
 246                 ++result;
 247         }
 248         a = 23; b = -4; calc_quo = iquo(a, b, calc_rem);
 249         if (calc_quo != -5 || calc_rem != 3) {
 250                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
 251                      << calc_quo << " with r=" << calc_rem << endl;
 252                 ++result;
 253         }
 254         a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
 255         if (calc_quo != -5 || calc_rem != -3) {
 256                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
 257                      << calc_quo << " with r=" << calc_rem << endl;
 258                 ++result;
 259         }
 260         a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
 261         if (calc_quo != 5 || calc_rem != -3) {
 262                 clog << "iquo(" << a << "," << b << ",r) erroneously returned "
 263                      << calc_quo << " with r=" << calc_rem << endl;
 264                 ++result;
 265         }
 266
 267         return result;
 268 }
 269
 270 /* Now we perform some less trivial checks about several functions which should
 271  * return exact numbers if possible. */
 272 static unsigned exam_numeric4(void)
 273 {
 274         unsigned result = 0;
 275         bool passed;
 276
 277         // square roots of squares of integers:
 278         passed = true;
 279         for (int i=0; i<42; ++i)
 280                 if (!sqrt(numeric(i*i)).is_integer())
 281                         passed = false;
 282         if (!passed) {
 283                 clog << "One or more square roots of squares of integers did not return exact integers" << endl;
 284                 ++result;
 285         }
 286
 287         // square roots of squares of rationals:
 288         passed = true;
 289         for (int num=0; num<41; ++num)
 290                 for (int den=1; den<42; ++den)
 291                         if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational())
 292                                 passed = false;
 293         if (!passed) {
 294                 clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
 295                 ++result;
 296         }
 297
 298         return result;
 299 }
 300
 301 /* This test examines that simplifications of the form 5^(3/2) -> 5*5^(1/2)
 302  * are carried out properly. */
 303 static unsigned exam_numeric5(void)
 304 {
 305         unsigned result = 0;
 306
 307         // A variation of one of Ramanujan's wonderful identities must be
 308         // verifiable with very primitive means:
 309         ex e1 = pow(1 + pow(3,numeric(1,5)) - pow(3,numeric(2,5)),3);
 310         ex e2 = expand(e1 - 10 + 5*pow(3,numeric(3,5)));
 311         if (!e2.is_zero()) {
 312                 clog << "expand((1+3^(1/5)-3^(2/5))^3-10+5*3^(3/5)) returned "
 313                      << e2 << " instead of 0." << endl;
 314                 ++result;
 315         }
 316
 317         return result;
 318 }
 319
 320 /* This test checks whether the numeric output/parsing routines are
 321    consistent. */
 322 static unsigned exam_numeric6(void)
 323 {
 324         unsigned result = 0;
 325
 326         symbol sym("sym");
 327         vector<ex> test_numbers;
 328         test_numbers.push_back(numeric(0));                     // zero
 329         test_numbers.push_back(numeric(1));                     // one
 330         test_numbers.push_back(numeric(-1));            // minus one
 331         test_numbers.push_back(numeric(42));            // positive integer
 332         test_numbers.push_back(numeric(-42));           // negative integer
 333         test_numbers.push_back(numeric(14,3));          // positive rational
 334         test_numbers.push_back(numeric(-14,3));         // negative rational
 335         test_numbers.push_back(numeric(3.141));         // positive decimal
 336         test_numbers.push_back(numeric(-3.141));        // negative decimal
 337         test_numbers.push_back(numeric(0.1974));        // positive decimal, leading zero
 338         test_numbers.push_back(numeric(-0.1974));       // negative decimal, leading zero
 339         test_numbers.push_back(sym);                            // symbol
 340
 341         for (vector<ex>::const_iterator br=test_numbers.begin(); br<test_numbers.end(); ++br) {
 342                 for (vector<ex>::const_iterator bi=test_numbers.begin(); bi<test_numbers.end(); ++bi) {
 343
 344                         for (vector<ex>::const_iterator er=test_numbers.begin(); er<test_numbers.end(); ++er) {
 345                                 for (vector<ex>::const_iterator ei=test_numbers.begin(); ei<test_numbers.end(); ++ei) {
 346
 347                                         // Construct expression, don't test invalid ones
 348                                         ex base = (*br) + (*bi)*I, exponent = (*er) + (*ei)*I, x;
 349                                         try {
 350                                                 x = pow(base, exponent);
 351                                         } catch (...) {
 352                                                 continue;
 353                                         }
 354
 355                                         // Print to string
 356                                         std::ostringstream s;
 357                                         s << x;
 358
 359                                         // Read back expression from string
 360                                         string x_as_string = s.str();
 361                                         ex x_again(x_as_string, lst(sym));
 362
 363                                         // They should be equal
 364                                         if (!x_again.is_equal(x)) {
 365                                                 clog << x << " was read back as " << x_again << endl;
 366                                                 ++result;
 367                                         }
 368                                 }
 369                         }
 370                 }
 371         }
 372
 373         return result;
 374 }
 375
 376 unsigned exam_numeric(void)
 377 {
 378         unsigned result = 0;
 379
 380         cout << "examining consistency of numeric types" << flush;
 381         clog << "----------consistency of numeric types:" << endl;
 382
 383         result += exam_numeric1();  cout << '.' << flush;
 384         result += exam_numeric2();  cout << '.' << flush;
 385         result += exam_numeric3();  cout << '.' << flush;
 386         result += exam_numeric4();  cout << '.' << flush;
 387         result += exam_numeric5();  cout << '.' << flush;
 388         result += exam_numeric6();  cout << '.' << flush;
 389
 390         if (!result) {
 391                 cout << " passed " << endl;
 392                 clog << "(no output)" << endl;
 393         } else {
 394                 cout << " failed " << endl;
 395         }
 396
 397         return result;
 398 }