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-2011 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  22  */
  23
  24 #include "ginac.h"
  25 using namespace GiNaC;
  26
  27 #include <iostream>
  28 #include <sstream>
  29 using namespace std;
  30
  31 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
  32  * conversions. */
  33 static unsigned exam_numeric1()
  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         if (test_crat.info(info_flags::nonnegative)) {
  72                 clog << test_crat
  73                      << " erroneously recognized as non-negative number" << endl;
  74                 ++result;
  75         }
  76
  77         int i = numeric(1984).to_int();
  78         if (i-1984) {
  79                 clog << "conversion of " << i
  80                      << " from numeric to int failed" << endl;
  81                 ++result;
  82         }
  83
  84         e1 = test_int1;
  85         if (!e1.info(info_flags::posint)) {
  86                 clog << "expression " << e1
  87                      << " erroneously not recognized as positive integer" << endl;
  88                 ++result;
  89         }
  90
  91         e2 = test_int1 + a;
  92         if (e2.info(info_flags::integer)) {
  93                 clog << "expression " << e2
  94                      << " erroneously recognized as integer" << endl;
  95                 ++result;
  96         }
  97
  98         // The next two were two actual bugs in CLN till June, 12, 1999:
  99         test_rat1 = numeric(3)/numeric(2);
 100         test_rat1 += test_rat1;
 101         if (!test_rat1.is_integer()) {
 102                 clog << "3/2 + 3/2 erroneously not integer 3 but instead "
 103                      << test_rat1 << endl;
 104                 ++result;
 105         }
 106         test_rat1 = numeric(3)/numeric(2);
 107         numeric test_rat2 = test_rat1 + numeric(1);  // 5/2
 108         test_rat2 -= test_rat1;  // 1
 109         if (!test_rat2.is_integer()) {
 110                 clog << "5/2 - 3/2 erroneously not integer 1 but instead "
 111                      << test_rat2 << 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 original bug in CLN was finally killed on September 2nd. */
 123 static unsigned exam_numeric2()
 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.is_zero()) {
 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 exam_numeric3()
 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 exam_numeric4()
 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         if (!passed) {
 291                 clog << "One or more square roots of squares of integers did not return exact integers" << endl;
 292                 ++result;
 293         }
 294
 295         // square roots of squares of rationals:
 296         passed = true;
 297         for (int num=0; num<41; ++num)
 298                 for (int den=1; den<42; ++den)
 299                         if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational())
 300                                 passed = false;
 301         if (!passed) {
 302                 clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
 303                 ++result;
 304         }
 305
 306         return result;
 307 }
 308
 309 /* This test examines that simplifications of the form 5^(3/2) -> 5*5^(1/2)
 310  * are carried out properly. */
 311 static unsigned exam_numeric5()
 312 {
 313         unsigned result = 0;
 314
 315         // A variation of one of Ramanujan's wonderful identities must be
 316         // verifiable with very primitive means:
 317         ex e1 = pow(1 + pow(3,numeric(1,5)) - pow(3,numeric(2,5)),3);
 318         ex e2 = expand(e1 - 10 + 5*pow(3,numeric(3,5)));
 319         if (!e2.is_zero()) {
 320                 clog << "expand((1+3^(1/5)-3^(2/5))^3-10+5*3^(3/5)) returned "
 321                      << e2 << " instead of 0." << endl;
 322                 ++result;
 323         }
 324
 325         return result;
 326 }
 327
 328 /* This test checks whether the numeric output/parsing routines are
 329    consistent. */
 330 static unsigned exam_numeric6()
 331 {
 332         unsigned result = 0;
 333
 334         symbol sym("sym");
 335         vector<ex> test_numbers;
 336         test_numbers.push_back(numeric(0));                     // zero
 337         test_numbers.push_back(numeric(1));                     // one
 338         test_numbers.push_back(numeric(-1));            // minus one
 339         test_numbers.push_back(numeric(42));            // positive integer
 340         test_numbers.push_back(numeric(-42));           // negative integer
 341         test_numbers.push_back(numeric(14,3));          // positive rational
 342         test_numbers.push_back(numeric(-14,3));         // negative rational
 343         test_numbers.push_back(numeric(3.141));         // positive decimal
 344         test_numbers.push_back(numeric(-3.141));        // negative decimal
 345         test_numbers.push_back(numeric(0.1974));        // positive decimal, leading zero
 346         test_numbers.push_back(numeric(-0.1974));       // negative decimal, leading zero
 347         test_numbers.push_back(sym);                            // symbol
 348
 349         for (vector<ex>::const_iterator br=test_numbers.begin(); br<test_numbers.end(); ++br) {
 350                 for (vector<ex>::const_iterator bi=test_numbers.begin(); bi<test_numbers.end(); ++bi) {
 351
 352                         for (vector<ex>::const_iterator er=test_numbers.begin(); er<test_numbers.end(); ++er) {
 353                                 for (vector<ex>::const_iterator ei=test_numbers.begin(); ei<test_numbers.end(); ++ei) {
 354
 355                                         // Construct expression, don't test invalid ones
 356                                         ex base = (*br) + (*bi)*I, exponent = (*er) + (*ei)*I, x;
 357                                         try {
 358                                                 x = pow(base, exponent);
 359                                         } catch (...) {
 360                                                 continue;
 361                                         }
 362
 363                                         // Print to string
 364                                         std::ostringstream s;
 365                                         s << x;
 366
 367                                         // Read back expression from string
 368                                         string x_as_string = s.str();
 369                                         ex x_again(x_as_string, lst(sym));
 370
 371                                         // They should be equal
 372                                         if (!x_again.is_equal(x)) {
 373                                                 clog << x << " was read back as " << x_again << endl;
 374                                                 ++result;
 375                                         }
 376                                 }
 377                         }
 378                 }
 379         }
 380
 381         return result;
 382 }
 383
 384 unsigned exam_numeric()
 385 {
 386         unsigned result = 0;
 387
 388         cout << "examining consistency of numeric types" << flush;
 389
 390         result += exam_numeric1();  cout << '.' << flush;
 391         result += exam_numeric2();  cout << '.' << flush;
 392         result += exam_numeric3();  cout << '.' << flush;
 393         result += exam_numeric4();  cout << '.' << flush;
 394         result += exam_numeric5();  cout << '.' << flush;
 395         result += exam_numeric6();  cout << '.' << flush;
 396
 397         return result;
 398 }
 399
 400 int main(int argc, char** argv)
 401 {
 402         return exam_numeric();
 403 }