Add remark about complex conjugation.
[ginac.git] / check / exam_numeric.cpp
index d3c1264a3aa1fdde148631f3f9f18b4fcfd19dfa..90e435e8f53f03243ee9723325f848f269012b64 100644 (file)
@@ -4,7 +4,7 @@
  *  tests on these numbers like is_integer() etc... */
 
 /*
- *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
 
 #include "exams.h"
 
+#include <sstream>
+
 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
  * conversions. */
-static unsigned exam_numeric1(void)
+static unsigned exam_numeric1()
 {
-    unsigned result = 0;
-    numeric test_int1(42);
-    numeric test_int2(5);
-    numeric test_rat1 = test_int1; test_rat1 /= test_int2;
-    test_rat1 = -test_rat1;         // -42/5
-    numeric test_crat = test_rat1+I*test_int2;  // 5*I-42/5
-    symbol a("a");
-    ex e1, e2;
-    
-    if (!test_int1.is_integer()) {
-        clog << test_int1
-             << " erroneously not recognized as integer" << endl;
-        ++result;
-    }
-    if (!test_int1.is_rational()) {
-        clog << test_int1
-             << " erroneously not recognized as rational" << endl;
-        ++result;
-    }
-    
-    if (!test_rat1.is_rational()) {
-        clog << test_rat1
-             << " erroneously not recognized as rational" << endl;
-        ++result;
-    }
-    if (test_rat1.is_integer()) {
-        clog << test_rat1
-             << " erroneously recognized as integer" << endl;
-        ++result;
-    }
-    
-    if (!test_crat.is_crational()) {
-        clog << test_crat
-             << " erroneously not recognized as complex rational" << endl;
-        ++result;
-    }
-    
-    int i = numeric(1984).to_int();
-    if (i-1984) {
-        clog << "conversion of " << i
-             << " from numeric to int failed" << endl;
-        ++result;
-    }
-    
-    e1 = test_int1;
-    if (!e1.info(info_flags::posint)) {
-        clog << "expression " << e1
-             << " erroneously not recognized as positive integer" << endl;
-        ++result;
-    }
-    
-    e2 = test_int1 + a;
-    if (ex_to_numeric(e2).is_integer()) {
-        clog << "expression " << e2
-             << " erroneously recognized as integer" << endl;
-        ++result;
-    }
-    
-    // The next two were two actual bugs in CLN till June, 12, 1999:
-    test_rat1 = numeric(3)/numeric(2);
-    test_rat1 += test_rat1;
-    if (!test_rat1.is_integer()) {
-        clog << "3/2 + 3/2 erroneously not integer 3 but instead "
-             << test_rat1 << endl;
-        ++result;
-    }
-    test_rat1 = numeric(3)/numeric(2);
-    numeric test_rat2 = test_rat1 + numeric(1);  // 5/2
-    test_rat2 -= test_rat1;  // 1
-    if (!test_rat2.is_integer()) {
-        clog << "5/2 - 3/2 erroneously not integer 1 but instead "
-             << test_rat2 << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       numeric test_int1(42);
+       numeric test_int2(5);
+       numeric test_rat1 = test_int1; test_rat1 /= test_int2;
+       test_rat1 = -test_rat1;          // -42/5
+       numeric test_crat = test_rat1+I*test_int2;  // 5*I-42/5
+       symbol a("a");
+       ex e1, e2;
+       
+       if (!test_int1.is_integer()) {
+               clog << test_int1
+                    << " erroneously not recognized as integer" << endl;
+               ++result;
+       }
+       if (!test_int1.is_rational()) {
+               clog << test_int1
+                    << " erroneously not recognized as rational" << endl;
+               ++result;
+       }
+       
+       if (!test_rat1.is_rational()) {
+               clog << test_rat1
+                    << " erroneously not recognized as rational" << endl;
+               ++result;
+       }
+       if (test_rat1.is_integer()) {
+               clog << test_rat1
+                        << " erroneously recognized as integer" << endl;
+               ++result;
+       }
+       
+       if (!test_crat.is_crational()) {
+               clog << test_crat
+                    << " erroneously not recognized as complex rational" << endl;
+               ++result;
+       }
+       
+       int i = numeric(1984).to_int();
+       if (i-1984) {
+               clog << "conversion of " << i
+                    << " from numeric to int failed" << endl;
+               ++result;
+       }
+       
+       e1 = test_int1;
+       if (!e1.info(info_flags::posint)) {
+               clog << "expression " << e1
+                    << " erroneously not recognized as positive integer" << endl;
+               ++result;
+       }
+       
+       e2 = test_int1 + a;
+       if (e2.info(info_flags::integer)) {
+               clog << "expression " << e2
+                    << " erroneously recognized as integer" << endl;
+               ++result;
+       }
+       
+       // The next two were two actual bugs in CLN till June, 12, 1999:
+       test_rat1 = numeric(3)/numeric(2);
+       test_rat1 += test_rat1;
+       if (!test_rat1.is_integer()) {
+               clog << "3/2 + 3/2 erroneously not integer 3 but instead "
+                    << test_rat1 << endl;
+               ++result;
+       }
+       test_rat1 = numeric(3)/numeric(2);
+       numeric test_rat2 = test_rat1 + numeric(1);  // 5/2
+       test_rat2 -= test_rat1;  // 1
+       if (!test_rat2.is_integer()) {
+               clog << "5/2 - 3/2 erroneously not integer 1 but instead "
+                    << test_rat2 << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 /* We had some fun with a bug in CLN that caused it to loop forever when
@@ -110,215 +112,287 @@ static unsigned exam_numeric1(void)
  * Implementing a workaround sadly introduced another bug on May 28th 1999
  * that was fixed on May 31st.  The workaround turned out to be stupid and
  * the original bug in CLN was finally killed on September 2nd. */
-static unsigned exam_numeric2(void)
+static unsigned exam_numeric2()
 {
-    unsigned result = 0;
-    
-    ex zero = numeric(0);
-    ex two = numeric(2);
-    ex three = numeric(3);
-    
-    // The hang in this code was the reason for the original workaround
-    if (pow(two,two/three)==42) {
-        clog << "pow(2,2/3) erroneously returned 42" << endl;
-        ++result;  // cannot happen
-    }
-    
-    // Actually, this used to raise a FPE after introducing the workaround
-    if (two*zero!=zero) {
-        clog << "2*0 erroneously returned " << two*zero << endl;
-        ++result;
-    }
-    
-    // And this returned a cl_F due to the implicit call of numeric::power()
-    ex six = two*three;
-    if (!six.info(info_flags::integer)) {
-        clog << "2*3 erroneously returned the non-integer " << six << endl;
-        ++result;
-    }
-    
-    // The fix in the workaround left a whole which was fixed hours later...
-    ex another_zero = pow(zero,numeric(1)/numeric(2));
-    if (!another_zero.is_zero()) {
-        clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       
+       ex zero = numeric(0);
+       ex two = numeric(2);
+       ex three = numeric(3);
+       
+       // The hang in this code was the reason for the original workaround
+       if (pow(two,two/three)==42) {
+               clog << "pow(2,2/3) erroneously returned 42" << endl;
+               ++result;  // cannot happen
+       }
+       
+       // Actually, this used to raise a FPE after introducing the workaround
+       if (two*zero!=zero) {
+               clog << "2*0 erroneously returned " << two*zero << endl;
+               ++result;
+       }
+       
+       // And this returned a cl_F due to the implicit call of numeric::power()
+       ex six = two*three;
+       if (!six.info(info_flags::integer)) {
+               clog << "2*3 erroneously returned the non-integer " << six << endl;
+               ++result;
+       }
+       
+       // The fix in the workaround left a whole which was fixed hours later...
+       ex another_zero = pow(zero,numeric(1)/numeric(2));
+       if (!another_zero.is_zero()) {
+               clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 /* Assorted tests to ensure some crucial functions behave exactly as specified
  * in the documentation. */
-static unsigned exam_numeric3(void)
+static unsigned exam_numeric3()
 {
-    unsigned result = 0;
-    numeric calc_rem, calc_quo;
-    numeric a, b;
-    
-    // check if irem(a, b) and irem(a, b, q) really behave like Maple's 
-    // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
-    // These overloaded routines indeed need to be checked separately since
-    // internally they might be doing something completely different:
-    a = 23; b = 4; calc_rem = irem(a, b);
-    if (calc_rem != 3) {
-        clog << "irem(" << a << "," << b << ") erroneously returned "
-             << calc_rem << endl;
-        ++result;
-    }
-    a = 23; b = -4; calc_rem = irem(a, b);
-    if (calc_rem != 3) {
-        clog << "irem(" << a << "," << b << ") erroneously returned "
-             << calc_rem << endl;
-        ++result;
-    }
-    a = -23; b = 4; calc_rem = irem(a, b);
-    if (calc_rem != -3) {
-        clog << "irem(" << a << "," << b << ") erroneously returned "
-             << calc_rem << endl;
-        ++result;
-    }
-    a = -23; b = -4; calc_rem = irem(a, b);
-    if (calc_rem != -3) {
-        clog << "irem(" << a << "," << b << ") erroneously returned "
-             << calc_rem << endl;
-        ++result;
-    }
-    // and now the overloaded irem(a,b,q):
-    a = 23; b = 4; calc_rem = irem(a, b, calc_quo);
-    if (calc_rem != 3 || calc_quo != 5) {
-        clog << "irem(" << a << "," << b << ",q) erroneously returned "
-             << calc_rem << " with q=" << calc_quo << endl;
-        ++result;
-    }
-    a = 23; b = -4; calc_rem = irem(a, b, calc_quo);
-    if (calc_rem != 3 || calc_quo != -5) {
-        clog << "irem(" << a << "," << b << ",q) erroneously returned "
-             << calc_rem << " with q=" << calc_quo << endl;
-        ++result;
-    }
-    a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
-    if (calc_rem != -3 || calc_quo != -5) {
-        clog << "irem(" << a << "," << b << ",q) erroneously returned "
-             << calc_rem << " with q=" << calc_quo << endl;
-        ++result;
-    }
-    a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
-    if (calc_rem != -3 || calc_quo != 5) {
-        clog << "irem(" << a << "," << b << ",q) erroneously returned "
-             << calc_rem << " with q=" << calc_quo << endl;
-        ++result;
-    }
-    // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's 
-    // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
-    // These overloaded routines indeed need to be checked separately since
-    // internally they might be doing something completely different:
-    a = 23; b = 4; calc_quo = iquo(a, b);
-    if (calc_quo != 5) {
-        clog << "iquo(" << a << "," << b << ") erroneously returned "
-             << calc_quo << endl;
-        ++result;
-    }
-    a = 23; b = -4; calc_quo = iquo(a, b);
-    if (calc_quo != -5) {
-        clog << "iquo(" << a << "," << b << ") erroneously returned "
-             << calc_quo << endl;
-        ++result;
-    }
-    a = -23; b = 4; calc_quo = iquo(a, b);
-    if (calc_quo != -5) {
-        clog << "iquo(" << a << "," << b << ") erroneously returned "
-             << calc_quo << endl;
-        ++result;
-    }
-    a = -23; b = -4; calc_quo = iquo(a, b);
-    if (calc_quo != 5) {
-        clog << "iquo(" << a << "," << b << ") erroneously returned "
-             << calc_quo << endl;
-        ++result;
-    }
-    // and now the overloaded iquo(a,b,r):
-    a = 23; b = 4; calc_quo = iquo(a, b, calc_rem);
-    if (calc_quo != 5 || calc_rem != 3) {
-        clog << "iquo(" << a << "," << b << ",r) erroneously returned "
-             << calc_quo << " with r=" << calc_rem << endl;
-        ++result;
-    }
-    a = 23; b = -4; calc_quo = iquo(a, b, calc_rem);
-    if (calc_quo != -5 || calc_rem != 3) {
-        clog << "iquo(" << a << "," << b << ",r) erroneously returned "
-             << calc_quo << " with r=" << calc_rem << endl;
-        ++result;
-    }
-    a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
-    if (calc_quo != -5 || calc_rem != -3) {
-        clog << "iquo(" << a << "," << b << ",r) erroneously returned "
-             << calc_quo << " with r=" << calc_rem << endl;
-        ++result;
-    }
-    a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
-    if (calc_quo != 5 || calc_rem != -3) {
-        clog << "iquo(" << a << "," << b << ",r) erroneously returned "
-             << calc_quo << " with r=" << calc_rem << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       numeric calc_rem, calc_quo;
+       numeric a, b;
+       
+       // check if irem(a, b) and irem(a, b, q) really behave like Maple's 
+       // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
+       // These overloaded routines indeed need to be checked separately since
+       // internally they might be doing something completely different:
+       a = 23; b = 4; calc_rem = irem(a, b);
+       if (calc_rem != 3) {
+               clog << "irem(" << a << "," << b << ") erroneously returned "
+                    << calc_rem << endl;
+               ++result;
+       }
+       a = 23; b = -4; calc_rem = irem(a, b);
+       if (calc_rem != 3) {
+               clog << "irem(" << a << "," << b << ") erroneously returned "
+                        << calc_rem << endl;
+               ++result;
+       }
+       a = -23; b = 4; calc_rem = irem(a, b);
+       if (calc_rem != -3) {
+               clog << "irem(" << a << "," << b << ") erroneously returned "
+                    << calc_rem << endl;
+               ++result;
+       }
+       a = -23; b = -4; calc_rem = irem(a, b);
+       if (calc_rem != -3) {
+               clog << "irem(" << a << "," << b << ") erroneously returned "
+                    << calc_rem << endl;
+               ++result;
+       }
+       // and now the overloaded irem(a,b,q):
+       a = 23; b = 4; calc_rem = irem(a, b, calc_quo);
+       if (calc_rem != 3 || calc_quo != 5) {
+               clog << "irem(" << a << "," << b << ",q) erroneously returned "
+                    << calc_rem << " with q=" << calc_quo << endl;
+               ++result;
+       }
+       a = 23; b = -4; calc_rem = irem(a, b, calc_quo);
+       if (calc_rem != 3 || calc_quo != -5) {
+               clog << "irem(" << a << "," << b << ",q) erroneously returned "
+                    << calc_rem << " with q=" << calc_quo << endl;
+               ++result;
+       }
+       a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
+       if (calc_rem != -3 || calc_quo != -5) {
+               clog << "irem(" << a << "," << b << ",q) erroneously returned "
+                    << calc_rem << " with q=" << calc_quo << endl;
+               ++result;
+       }
+       a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
+       if (calc_rem != -3 || calc_quo != 5) {
+               clog << "irem(" << a << "," << b << ",q) erroneously returned "
+                    << calc_rem << " with q=" << calc_quo << endl;
+               ++result;
+       }
+       // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's 
+       // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
+       // These overloaded routines indeed need to be checked separately since
+       // internally they might be doing something completely different:
+       a = 23; b = 4; calc_quo = iquo(a, b);
+       if (calc_quo != 5) {
+               clog << "iquo(" << a << "," << b << ") erroneously returned "
+                    << calc_quo << endl;
+               ++result;
+       }
+       a = 23; b = -4; calc_quo = iquo(a, b);
+       if (calc_quo != -5) {
+               clog << "iquo(" << a << "," << b << ") erroneously returned "
+                    << calc_quo << endl;
+               ++result;
+       }
+       a = -23; b = 4; calc_quo = iquo(a, b);
+       if (calc_quo != -5) {
+               clog << "iquo(" << a << "," << b << ") erroneously returned "
+                    << calc_quo << endl;
+               ++result;
+       }
+       a = -23; b = -4; calc_quo = iquo(a, b);
+       if (calc_quo != 5) {
+               clog << "iquo(" << a << "," << b << ") erroneously returned "
+                    << calc_quo << endl;
+               ++result;
+       }
+       // and now the overloaded iquo(a,b,r):
+       a = 23; b = 4; calc_quo = iquo(a, b, calc_rem);
+       if (calc_quo != 5 || calc_rem != 3) {
+               clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+                    << calc_quo << " with r=" << calc_rem << endl;
+               ++result;
+       }
+       a = 23; b = -4; calc_quo = iquo(a, b, calc_rem);
+       if (calc_quo != -5 || calc_rem != 3) {
+               clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+                    << calc_quo << " with r=" << calc_rem << endl;
+               ++result;
+       }
+       a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
+       if (calc_quo != -5 || calc_rem != -3) {
+               clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+                    << calc_quo << " with r=" << calc_rem << endl;
+               ++result;
+       }
+       a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
+       if (calc_quo != 5 || calc_rem != -3) {
+               clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+                    << calc_quo << " with r=" << calc_rem << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 /* Now we perform some less trivial checks about several functions which should
  * return exact numbers if possible. */
-static unsigned exam_numeric4(void)
+static unsigned exam_numeric4()
+{
+       unsigned result = 0;
+       bool passed;
+       
+       // square roots of squares of integers:
+       passed = true;
+       for (int i=0; i<42; ++i)
+               if (!sqrt(numeric(i*i)).is_integer())
+                       passed = false;
+       if (!passed) {
+               clog << "One or more square roots of squares of integers did not return exact integers" << endl;
+               ++result;
+       }
+       
+       // square roots of squares of rationals:
+       passed = true;
+       for (int num=0; num<41; ++num)
+               for (int den=1; den<42; ++den)
+                       if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational())
+                               passed = false;
+       if (!passed) {
+               clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
+               ++result;
+       }
+       
+       return result;
+}
+
+/* This test examines that simplifications of the form 5^(3/2) -> 5*5^(1/2)
+ * are carried out properly. */
+static unsigned exam_numeric5()
 {
-    unsigned result = 0;
-    bool passed;
-    
-    // square roots of squares of integers:
-    passed = true;
-    for (int i=0; i<42; ++i) {
-        if (!sqrt(numeric(i*i)).is_integer()) {
-            passed = false;
-        }
-    }
-    if (!passed) {
-        clog << "One or more square roots of squares of integers did not return exact integers" << endl;
-        ++result;
-    }
-    
-    // square roots of squares of rationals:
-    passed = true;
-    for (int num=0; num<41; ++num) {
-        for (int den=1; den<42; ++den) {
-            if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational()) {
-                passed = false;
-            }
-        }
-    }
-    if (!passed) {
-        clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       
+       // A variation of one of Ramanujan's wonderful identities must be
+       // verifiable with very primitive means:
+       ex e1 = pow(1 + pow(3,numeric(1,5)) - pow(3,numeric(2,5)),3);
+       ex e2 = expand(e1 - 10 + 5*pow(3,numeric(3,5)));
+       if (!e2.is_zero()) {
+               clog << "expand((1+3^(1/5)-3^(2/5))^3-10+5*3^(3/5)) returned "
+                    << e2 << " instead of 0." << endl;
+               ++result;
+       }
+       
+       return result;
+}
+
+/* This test checks whether the numeric output/parsing routines are
+   consistent. */
+static unsigned exam_numeric6()
+{
+       unsigned result = 0;
+
+       symbol sym("sym");
+       vector<ex> test_numbers;
+       test_numbers.push_back(numeric(0));                     // zero
+       test_numbers.push_back(numeric(1));                     // one
+       test_numbers.push_back(numeric(-1));            // minus one
+       test_numbers.push_back(numeric(42));            // positive integer
+       test_numbers.push_back(numeric(-42));           // negative integer
+       test_numbers.push_back(numeric(14,3));          // positive rational
+       test_numbers.push_back(numeric(-14,3));         // negative rational
+       test_numbers.push_back(numeric(3.141));         // positive decimal
+       test_numbers.push_back(numeric(-3.141));        // negative decimal
+       test_numbers.push_back(numeric(0.1974));        // positive decimal, leading zero
+       test_numbers.push_back(numeric(-0.1974));       // negative decimal, leading zero
+       test_numbers.push_back(sym);                            // symbol
+
+       for (vector<ex>::const_iterator br=test_numbers.begin(); br<test_numbers.end(); ++br) {
+               for (vector<ex>::const_iterator bi=test_numbers.begin(); bi<test_numbers.end(); ++bi) {
+
+                       for (vector<ex>::const_iterator er=test_numbers.begin(); er<test_numbers.end(); ++er) {
+                               for (vector<ex>::const_iterator ei=test_numbers.begin(); ei<test_numbers.end(); ++ei) {
+
+                                       // Construct expression, don't test invalid ones
+                                       ex base = (*br) + (*bi)*I, exponent = (*er) + (*ei)*I, x;
+                                       try {
+                                               x = pow(base, exponent);
+                                       } catch (...) {
+                                               continue;
+                                       }
+
+                                       // Print to string
+                                       std::ostringstream s;
+                                       s << x;
+
+                                       // Read back expression from string
+                                       string x_as_string = s.str();
+                                       ex x_again(x_as_string, lst(sym));
+
+                                       // They should be equal
+                                       if (!x_again.is_equal(x)) {
+                                               clog << x << " was read back as " << x_again << endl;
+                                               ++result;
+                                       }
+                               }
+                       }
+               }
+       }
+
+       return result;
 }
 
-unsigned exam_numeric(void)
+unsigned exam_numeric()
 {
-    unsigned result = 0;
-    
-    cout << "examining consistency of numeric types" << flush;
-    clog << "----------consistency of numeric types:" << endl;
-    
-    result += exam_numeric1();  cout << '.' << flush;
-    result += exam_numeric2();  cout << '.' << flush;
-    result += exam_numeric3();  cout << '.' << flush;
-    result += exam_numeric4();  cout << '.' << flush;
-    
-    if (!result) {
-        cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
-        cout << " failed " << endl;
-    }
-    
-    return result;
+       unsigned result = 0;
+       
+       cout << "examining consistency of numeric types" << flush;
+       clog << "----------consistency of numeric types:" << endl;
+       
+       result += exam_numeric1();  cout << '.' << flush;
+       result += exam_numeric2();  cout << '.' << flush;
+       result += exam_numeric3();  cout << '.' << flush;
+       result += exam_numeric4();  cout << '.' << flush;
+       result += exam_numeric5();  cout << '.' << flush;
+       result += exam_numeric6();  cout << '.' << flush;
+       
+       if (!result) {
+               cout << " passed " << endl;
+               clog << "(no output)" << endl;
+       } else {
+               cout << " failed " << endl;
+       }
+       
+       return result;
 }