From: Richard Kreckel Date: Sat, 12 Aug 2000 21:59:10 +0000 (+0000) Subject: - Same shit as Christian did yesterday in ginac/. X-Git-Tag: release_0-7-0~41 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=af922d5eb36ed70e4a9e3ffaf4c24492cf89a1a6 - Same shit as Christian did yesterday in ginac/. --- diff --git a/check/check_inifcns.cpp b/check/check_inifcns.cpp index 9cebceeb..bec2ed97 100644 --- a/check/check_inifcns.cpp +++ b/check/check_inifcns.cpp @@ -26,188 +26,188 @@ /* Some tests on the sine trigonometric function. */ static unsigned inifcns_check_sin(void) { - unsigned result = 0; - bool errorflag = false; - - // sin(n*Pi) == 0? - errorflag = false; - for (int n=-10; n<=10; ++n) { - if (sin(n*Pi).eval() != numeric(0) || - !sin(n*Pi).eval().info(info_flags::integer)) - errorflag = true; - } - if (errorflag) { - // we don't count each of those errors - clog << "sin(n*Pi) with integer n does not always return exact 0" - << endl; - ++result; - } - - // sin((n+1/2)*Pi) == {+|-}1? - errorflag = false; - for (int n=-10; n<=10; ++n) { - if (!sin((n+numeric(1,2))*Pi).eval().info(info_flags::integer) || - !(sin((n+numeric(1,2))*Pi).eval() == numeric(1) || - sin((n+numeric(1,2))*Pi).eval() == numeric(-1))) - errorflag = true; - } - if (errorflag) { - clog << "sin((n+1/2)*Pi) with integer n does not always return exact {+|-}1" - << endl; - ++result; - } - - // compare sin((q*Pi).evalf()) with sin(q*Pi).eval().evalf() at various - // points. E.g. if sin(Pi/10) returns something symbolic this should be - // equal to sqrt(5)/4-1/4. This routine will spot programming mistakes - // of this kind: - errorflag = false; - ex argument; - numeric epsilon(double(1e-8)); - for (int n=-340; n<=340; ++n) { - argument = n*Pi/60; - if (abs(sin(evalf(argument))-evalf(sin(argument)))>epsilon) { - clog << "sin(" << argument << ") returns " - << sin(argument) << endl; - errorflag = true; - } - } - if (errorflag) - ++result; - - return result; + unsigned result = 0; + bool errorflag = false; + + // sin(n*Pi) == 0? + errorflag = false; + for (int n=-10; n<=10; ++n) { + if (sin(n*Pi).eval() != numeric(0) || + !sin(n*Pi).eval().info(info_flags::integer)) + errorflag = true; + } + if (errorflag) { + // we don't count each of those errors + clog << "sin(n*Pi) with integer n does not always return exact 0" + << endl; + ++result; + } + + // sin((n+1/2)*Pi) == {+|-}1? + errorflag = false; + for (int n=-10; n<=10; ++n) { + if (!sin((n+numeric(1,2))*Pi).eval().info(info_flags::integer) || + !(sin((n+numeric(1,2))*Pi).eval() == numeric(1) || + sin((n+numeric(1,2))*Pi).eval() == numeric(-1))) + errorflag = true; + } + if (errorflag) { + clog << "sin((n+1/2)*Pi) with integer n does not always return exact {+|-}1" + << endl; + ++result; + } + + // compare sin((q*Pi).evalf()) with sin(q*Pi).eval().evalf() at various + // points. E.g. if sin(Pi/10) returns something symbolic this should be + // equal to sqrt(5)/4-1/4. This routine will spot programming mistakes + // of this kind: + errorflag = false; + ex argument; + numeric epsilon(double(1e-8)); + for (int n=-340; n<=340; ++n) { + argument = n*Pi/60; + if (abs(sin(evalf(argument))-evalf(sin(argument)))>epsilon) { + clog << "sin(" << argument << ") returns " + << sin(argument) << endl; + errorflag = true; + } + } + if (errorflag) + ++result; + + return result; } /* Simple tests on the cosine trigonometric function. */ static unsigned inifcns_check_cos(void) { - unsigned result = 0; - bool errorflag; - - // cos((n+1/2)*Pi) == 0? - errorflag = false; - for (int n=-10; n<=10; ++n) { - if (cos((n+numeric(1,2))*Pi).eval() != numeric(0) || - !cos((n+numeric(1,2))*Pi).eval().info(info_flags::integer)) - errorflag = true; - } - if (errorflag) { - clog << "cos((n+1/2)*Pi) with integer n does not always return exact 0" - << endl; - ++result; - } - - // cos(n*Pi) == 0? - errorflag = false; - for (int n=-10; n<=10; ++n) { - if (!cos(n*Pi).eval().info(info_flags::integer) || - !(cos(n*Pi).eval() == numeric(1) || - cos(n*Pi).eval() == numeric(-1))) - errorflag = true; - } - if (errorflag) { - clog << "cos(n*Pi) with integer n does not always return exact {+|-}1" - << endl; - ++result; - } - - // compare cos((q*Pi).evalf()) with cos(q*Pi).eval().evalf() at various - // points. E.g. if cos(Pi/12) returns something symbolic this should be - // equal to 1/4*(1+1/3*sqrt(3))*sqrt(6). This routine will spot - // programming mistakes of this kind: - errorflag = false; - ex argument; - numeric epsilon(double(1e-8)); - for (int n=-340; n<=340; ++n) { - argument = n*Pi/60; - if (abs(cos(evalf(argument))-evalf(cos(argument)))>epsilon) { - clog << "cos(" << argument << ") returns " - << cos(argument) << endl; - errorflag = true; - } - } - if (errorflag) - ++result; - - return result; + unsigned result = 0; + bool errorflag; + + // cos((n+1/2)*Pi) == 0? + errorflag = false; + for (int n=-10; n<=10; ++n) { + if (cos((n+numeric(1,2))*Pi).eval() != numeric(0) || + !cos((n+numeric(1,2))*Pi).eval().info(info_flags::integer)) + errorflag = true; + } + if (errorflag) { + clog << "cos((n+1/2)*Pi) with integer n does not always return exact 0" + << endl; + ++result; + } + + // cos(n*Pi) == 0? + errorflag = false; + for (int n=-10; n<=10; ++n) { + if (!cos(n*Pi).eval().info(info_flags::integer) || + !(cos(n*Pi).eval() == numeric(1) || + cos(n*Pi).eval() == numeric(-1))) + errorflag = true; + } + if (errorflag) { + clog << "cos(n*Pi) with integer n does not always return exact {+|-}1" + << endl; + ++result; + } + + // compare cos((q*Pi).evalf()) with cos(q*Pi).eval().evalf() at various + // points. E.g. if cos(Pi/12) returns something symbolic this should be + // equal to 1/4*(1+1/3*sqrt(3))*sqrt(6). This routine will spot + // programming mistakes of this kind: + errorflag = false; + ex argument; + numeric epsilon(double(1e-8)); + for (int n=-340; n<=340; ++n) { + argument = n*Pi/60; + if (abs(cos(evalf(argument))-evalf(cos(argument)))>epsilon) { + clog << "cos(" << argument << ") returns " + << cos(argument) << endl; + errorflag = true; + } + } + if (errorflag) + ++result; + + return result; } /* Simple tests on the tangent trigonometric function. */ static unsigned inifcns_check_tan(void) { - unsigned result = 0; - bool errorflag; - - // compare tan((q*Pi).evalf()) with tan(q*Pi).eval().evalf() at various - // points. E.g. if tan(Pi/12) returns something symbolic this should be - // equal to 2-sqrt(3). This routine will spot programming mistakes of - // this kind: - errorflag = false; - ex argument; - numeric epsilon(double(1e-8)); - for (int n=-340; n<=340; ++n) { - if (!(n%30) && (n%60)) // skip poles - ++n; - argument = n*Pi/60; - if (abs(tan(evalf(argument))-evalf(tan(argument)))>epsilon) { - clog << "tan(" << argument << ") returns " - << tan(argument) << endl; - errorflag = true; - } - } - if (errorflag) - ++result; - - return result; + unsigned result = 0; + bool errorflag; + + // compare tan((q*Pi).evalf()) with tan(q*Pi).eval().evalf() at various + // points. E.g. if tan(Pi/12) returns something symbolic this should be + // equal to 2-sqrt(3). This routine will spot programming mistakes of + // this kind: + errorflag = false; + ex argument; + numeric epsilon(double(1e-8)); + for (int n=-340; n<=340; ++n) { + if (!(n%30) && (n%60)) // skip poles + ++n; + argument = n*Pi/60; + if (abs(tan(evalf(argument))-evalf(tan(argument)))>epsilon) { + clog << "tan(" << argument << ") returns " + << tan(argument) << endl; + errorflag = true; + } + } + if (errorflag) + ++result; + + return result; } /* Simple tests on the dilogarithm function. */ static unsigned inifcns_check_Li2(void) { - // NOTE: this can safely be removed once CLN supports dilogarithms and - // checks them itself. - unsigned result = 0; - bool errorflag; - - // check the relation Li2(z^2) == 2 * (Li2(z) + Li2(-z)) numerically, which - // should hold in the entire complex plane: - errorflag = false; - ex argument; - numeric epsilon(double(1e-16)); - for (int n=0; n<200; ++n) { - argument = numeric(20.0*rand()/(RAND_MAX+1.0)-10.0) - + numeric(20.0*rand()/(RAND_MAX+1.0)-10.0)*I; - if (abs(Li2(pow(argument,2))-2*Li2(argument)-2*Li2(-argument)) > epsilon) { - cout << "Li2(z) at z==" << argument - << " failed to satisfy Li2(z^2)==2*(Li2(z)+Li2(-z))" << endl; - errorflag = true; - } - } - - if (errorflag) - ++result; - - return result; + // NOTE: this can safely be removed once CLN supports dilogarithms and + // checks them itself. + unsigned result = 0; + bool errorflag; + + // check the relation Li2(z^2) == 2 * (Li2(z) + Li2(-z)) numerically, which + // should hold in the entire complex plane: + errorflag = false; + ex argument; + numeric epsilon(double(1e-16)); + for (int n=0; n<200; ++n) { + argument = numeric(20.0*rand()/(RAND_MAX+1.0)-10.0) + + numeric(20.0*rand()/(RAND_MAX+1.0)-10.0)*I; + if (abs(Li2(pow(argument,2))-2*Li2(argument)-2*Li2(-argument)) > epsilon) { + cout << "Li2(z) at z==" << argument + << " failed to satisfy Li2(z^2)==2*(Li2(z)+Li2(-z))" << endl; + errorflag = true; + } + } + + if (errorflag) + ++result; + + return result; } unsigned check_inifcns(void) { - unsigned result = 0; + unsigned result = 0; - cout << "checking consistency of symbolic functions" << flush; - clog << "---------consistency of symbolic functions:" << endl; - - result += inifcns_check_sin(); cout << '.' << flush; - result += inifcns_check_cos(); cout << '.' << flush; - result += inifcns_check_tan(); cout << '.' << flush; - result += inifcns_check_Li2(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + cout << "checking consistency of symbolic functions" << flush; + clog << "---------consistency of symbolic functions:" << endl; + + result += inifcns_check_sin(); cout << '.' << flush; + result += inifcns_check_cos(); cout << '.' << flush; + result += inifcns_check_tan(); cout << '.' << flush; + result += inifcns_check_Li2(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/check_lsolve.cpp b/check/check_lsolve.cpp index 1873bb2d..4a84a13c 100644 --- a/check/check_lsolve.cpp +++ b/check/check_lsolve.cpp @@ -24,177 +24,177 @@ #include "checks.h" static unsigned check_matrix_solve(unsigned m, unsigned n, unsigned p, - unsigned degree) + unsigned degree) { - const symbol a("a"); - matrix A(m,n); - matrix B(m,p); - // set the first min(m,n) rows of A and B - for (unsigned ro=0; (ro x; - matrix X(n,p); - for (unsigned i=0; i x; + matrix X(n,p); + for (unsigned i=0; i a; - vector x; - for (unsigned i=0; i a; + vector x; + for (unsigned i=0; i1.0001 && ratio<0.9999) { - clog << "(" << num << "/" << den << ")^(" << nm - << ") erroneously evaluated to " << radical; - errorflag = true; - } - } - } - if (errorflag) - ++result; - - return result; + unsigned result = 0; + bool errorflag = false; + int i_num, i_den; + + // Check non-nested radicals (n/d)^(m/n) in ex wrapper class: + for (int i=0; i<200; ++i) { // FIXME: run to ~200 + for (int j=2; j<13; ++j) { + // construct an exponent 1/j... + numeric nm(1,j); + nm += numeric(int(20.0*rand()/(RAND_MAX+1.0))-10); + // ...a numerator... + do { i_num = rand(); } while (i_num == 0); + numeric num(i_num); + // ...and a denominator. + do { i_den = (rand())/100; } while (i_den == 0); + numeric den(i_den); + // construct the radicals: + ex radical = pow(ex(num)/ex(den),ex(nm)); + numeric floating = pow(num/den,nm); + // test the result: + if (is_ex_of_type(radical,numeric)) { + clog << "(" << num << "/" << den << ")^(" << nm + << ") should have been a product, instead it's " + << radical << endl; + errorflag = true; + } + numeric ratio = ex_to_numeric(evalf(radical))/floating; + if (ratio>1.0001 && ratio<0.9999) { + clog << "(" << num << "/" << den << ")^(" << nm + << ") erroneously evaluated to " << radical; + errorflag = true; + } + } + } + if (errorflag) + ++result; + + return result; } unsigned check_numeric(void) { - unsigned result = 0; - - cout << "checking consistency of numeric types" << flush; - clog << "---------consistency of numeric types:" << endl; - - result += check_numeric1(); cout << '.' << flush; - result += check_numeric2(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + unsigned result = 0; + + cout << "checking consistency of numeric types" << flush; + clog << "---------consistency of numeric types:" << endl; + + result += check_numeric1(); cout << '.' << flush; + result += check_numeric2(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/checks.cpp b/check/checks.cpp index b7827381..87782fcc 100644 --- a/check/checks.cpp +++ b/check/checks.cpp @@ -28,52 +28,52 @@ int main() { - unsigned result = 0; - - srand((unsigned)time(NULL)); - - try { - for (int i=0; i<1; ++i) - result += check_numeric(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - for (int i=0; i<1; ++i) - result += check_inifcns(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - for (int i=0; i<1; ++i) - result += check_matrices(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - for (int i=0; i<1; ++i) - result += check_lsolve(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - if (result) { - cout << "Error: something went wrong. "; - if (result == 1) { - cout << "(one failure)" << endl; - } else { - cout << "(" << result << " individual failures)" << endl; - } - cout << "please check checks.out against check.ref for more details." - << endl << "happy debugging!" << endl; - } - - return result; + unsigned result = 0; + + srand((unsigned)time(NULL)); + + try { + for (int i=0; i<1; ++i) + result += check_numeric(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + for (int i=0; i<1; ++i) + result += check_inifcns(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + for (int i=0; i<1; ++i) + result += check_matrices(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + for (int i=0; i<1; ++i) + result += check_lsolve(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + if (result) { + cout << "Error: something went wrong. "; + if (result == 1) { + cout << "(one failure)" << endl; + } else { + cout << "(" << result << " individual failures)" << endl; + } + cout << "please check checks.out against check.ref for more details." + << endl << "happy debugging!" << endl; + } + + return result; } diff --git a/check/checks.h b/check/checks.h index 5827f223..e585eab1 100644 --- a/check/checks.h +++ b/check/checks.h @@ -37,8 +37,8 @@ using namespace GiNaC; const ex dense_univariate_poly(const symbol & x, unsigned degree); const ex dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree); const ex sparse_tree(const symbol & x, const symbol & y, const symbol & z, - int level, - bool trig = false, bool rational = true, bool complex = false); + int level, + bool trig = false, bool rational = true, bool complex = false); // prototypes for all individual checks should be unsigned fcn(): unsigned check_numeric(); diff --git a/check/exam_differentiation.cpp b/check/exam_differentiation.cpp index 1a5aabb1..13e54711 100644 --- a/check/exam_differentiation.cpp +++ b/check/exam_differentiation.cpp @@ -23,141 +23,141 @@ #include "exams.h" static unsigned check_diff(const ex &e, const symbol &x, - const ex &d, unsigned nth=1) + const ex &d, unsigned nth=1) { - ex ed = e.diff(x, nth); - if ((ed - d).compare(ex(0)) != 0) { - switch (nth) { - case 0: - clog << "zeroth "; - break; - case 1: - break; - case 2: - clog << "second "; - break; - case 3: - clog << "third "; - break; - default: - clog << nth << "th "; - } - clog << "derivative of " << e << " by " << x << " returned " - << ed << " instead of " << d << endl; - clog << "returned:" << endl; - ed.printtree(clog); - clog << endl << "instead of" << endl; - d.printtree(clog); + ex ed = e.diff(x, nth); + if ((ed - d).compare(ex(0)) != 0) { + switch (nth) { + case 0: + clog << "zeroth "; + break; + case 1: + break; + case 2: + clog << "second "; + break; + case 3: + clog << "third "; + break; + default: + clog << nth << "th "; + } + clog << "derivative of " << e << " by " << x << " returned " + << ed << " instead of " << d << endl; + clog << "returned:" << endl; + ed.printtree(clog); + clog << endl << "instead of" << endl; + d.printtree(clog); - return 1; - } - return 0; + return 1; + } + return 0; } // Simple (expanded) polynomials static unsigned exam_differentiation1(void) { - unsigned result = 0; - symbol x("x"), y("y"); - ex e1, e2, e, d; - - // construct bivariate polynomial e to be diff'ed: - e1 = pow(x, -2) * 3 + pow(x, -1) * 5 + 7 + x * 11 + pow(x, 2) * 13; - e2 = pow(y, -2) * 5 + pow(y, -1) * 7 + 11 + y * 13 + pow(y, 2) * 17; - e = (e1 * e2).expand(); - - // d e / dx: - d = ex("121-55/x^2-66/x^3-30/x^3/y^2-42/x^3/y-78/x^3*y-102/x^3*y^2-25/x^2/y^2-35/x^2/y-65/x^2*y-85/x^2*y^2+77/y+143*y+187*y^2+130*x/y^2+182/y*x+338*x*y+442*x*y^2+55/y^2+286*x",lst(x,y)); - result += check_diff(e, x, d); - - // d e / dy: - d = ex("91-30/x^2/y^3-21/x^2/y^2+39/x^2+102/x^2*y-50/x/y^3-35/x/y^2+65/x+170/x*y-77*x/y^2+143*x+374*x*y-130/y^3*x^2-91/y^2*x^2+169*x^2+442*x^2*y-110/y^3*x-70/y^3+238*y-49/y^2",lst(x,y)); - result += check_diff(e, y, d); - - // d^2 e / dx^2: - d = ex("286+90/x^4/y^2+126/x^4/y+234/x^4*y+306/x^4*y^2+50/x^3/y^2+70/x^3/y+130/x^3*y+170/x^3*y^2+130/y^2+182/y+338*y+442*y^2+198/x^4+110/x^3",lst(x,y)); - result += check_diff(e, x, d, 2); - - // d^2 e / dy^2: - d = ex("238+90/x^2/y^4+42/x^2/y^3+102/x^2+150/x/y^4+70/x/y^3+170/x+330*x/y^4+154*x/y^3+374*x+390*x^2/y^4+182*x^2/y^3+442*x^2+210/y^4+98/y^3",lst(x,y)); - result += check_diff(e, y, d, 2); - - return result; + unsigned result = 0; + symbol x("x"), y("y"); + ex e1, e2, e, d; + + // construct bivariate polynomial e to be diff'ed: + e1 = pow(x, -2) * 3 + pow(x, -1) * 5 + 7 + x * 11 + pow(x, 2) * 13; + e2 = pow(y, -2) * 5 + pow(y, -1) * 7 + 11 + y * 13 + pow(y, 2) * 17; + e = (e1 * e2).expand(); + + // d e / dx: + d = ex("121-55/x^2-66/x^3-30/x^3/y^2-42/x^3/y-78/x^3*y-102/x^3*y^2-25/x^2/y^2-35/x^2/y-65/x^2*y-85/x^2*y^2+77/y+143*y+187*y^2+130*x/y^2+182/y*x+338*x*y+442*x*y^2+55/y^2+286*x",lst(x,y)); + result += check_diff(e, x, d); + + // d e / dy: + d = ex("91-30/x^2/y^3-21/x^2/y^2+39/x^2+102/x^2*y-50/x/y^3-35/x/y^2+65/x+170/x*y-77*x/y^2+143*x+374*x*y-130/y^3*x^2-91/y^2*x^2+169*x^2+442*x^2*y-110/y^3*x-70/y^3+238*y-49/y^2",lst(x,y)); + result += check_diff(e, y, d); + + // d^2 e / dx^2: + d = ex("286+90/x^4/y^2+126/x^4/y+234/x^4*y+306/x^4*y^2+50/x^3/y^2+70/x^3/y+130/x^3*y+170/x^3*y^2+130/y^2+182/y+338*y+442*y^2+198/x^4+110/x^3",lst(x,y)); + result += check_diff(e, x, d, 2); + + // d^2 e / dy^2: + d = ex("238+90/x^2/y^4+42/x^2/y^3+102/x^2+150/x/y^4+70/x/y^3+170/x+330*x/y^4+154*x/y^3+374*x+390*x^2/y^4+182*x^2/y^3+442*x^2+210/y^4+98/y^3",lst(x,y)); + result += check_diff(e, y, d, 2); + + return result; } // Trigonometric functions static unsigned exam_differentiation2(void) { - unsigned result = 0; - symbol x("x"), y("y"), a("a"), b("b"); - ex e1, e2, e, d; - - // construct expression e to be diff'ed: - e1 = y*pow(x, 2) + a*x + b; - e2 = sin(e1); - e = b*pow(e2, 2) + y*e2 + a; - - d = 2*b*e2*cos(e1)*(2*x*y + a) + y*cos(e1)*(2*x*y + a); - result += check_diff(e, x, d); - - d = 2*b*pow(cos(e1),2)*pow(2*x*y + a, 2) + 4*b*y*e2*cos(e1) - - 2*b*pow(e2,2)*pow(2*x*y + a, 2) - y*e2*pow(2*x*y + a, 2) - + 2*pow(y,2)*cos(e1); - result += check_diff(e, x, d, 2); - - d = 2*b*e2*cos(e1)*pow(x, 2) + e2 + y*cos(e1)*pow(x, 2); - result += check_diff(e, y, d); + unsigned result = 0; + symbol x("x"), y("y"), a("a"), b("b"); + ex e1, e2, e, d; + + // construct expression e to be diff'ed: + e1 = y*pow(x, 2) + a*x + b; + e2 = sin(e1); + e = b*pow(e2, 2) + y*e2 + a; + + d = 2*b*e2*cos(e1)*(2*x*y + a) + y*cos(e1)*(2*x*y + a); + result += check_diff(e, x, d); + + d = 2*b*pow(cos(e1),2)*pow(2*x*y + a, 2) + 4*b*y*e2*cos(e1) + - 2*b*pow(e2,2)*pow(2*x*y + a, 2) - y*e2*pow(2*x*y + a, 2) + + 2*pow(y,2)*cos(e1); + result += check_diff(e, x, d, 2); + + d = 2*b*e2*cos(e1)*pow(x, 2) + e2 + y*cos(e1)*pow(x, 2); + result += check_diff(e, y, d); - d = 2*b*pow(cos(e1),2)*pow(x,4) - 2*b*pow(e2,2)*pow(x,4) - + 2*cos(e1)*pow(x,2) - y*e2*pow(x,4); - result += check_diff(e, y, d, 2); - - // construct expression e to be diff'ed: - e2 = cos(e1); - e = b*pow(e2, 2) + y*e2 + a; - - d = -2*b*e2*sin(e1)*(2*x*y + a) - y*sin(e1)*(2*x*y + a); - result += check_diff(e, x, d); - - d = 2*b*pow(sin(e1),2)*pow(2*y*x + a,2) - 4*b*e2*sin(e1)*y - - 2*b*pow(e2,2)*pow(2*y*x + a,2) - y*e2*pow(2*y*x + a,2) - - 2*pow(y,2)*sin(e1); - result += check_diff(e, x, d, 2); - - d = -2*b*e2*sin(e1)*pow(x,2) + e2 - y*sin(e1)*pow(x, 2); - result += check_diff(e, y, d); - - d = -2*b*pow(e2,2)*pow(x,4) + 2*b*pow(sin(e1),2)*pow(x,4) - - 2*sin(e1)*pow(x,2) - y*e2*pow(x,4); - result += check_diff(e, y, d, 2); + d = 2*b*pow(cos(e1),2)*pow(x,4) - 2*b*pow(e2,2)*pow(x,4) + + 2*cos(e1)*pow(x,2) - y*e2*pow(x,4); + result += check_diff(e, y, d, 2); + + // construct expression e to be diff'ed: + e2 = cos(e1); + e = b*pow(e2, 2) + y*e2 + a; + + d = -2*b*e2*sin(e1)*(2*x*y + a) - y*sin(e1)*(2*x*y + a); + result += check_diff(e, x, d); + + d = 2*b*pow(sin(e1),2)*pow(2*y*x + a,2) - 4*b*e2*sin(e1)*y + - 2*b*pow(e2,2)*pow(2*y*x + a,2) - y*e2*pow(2*y*x + a,2) + - 2*pow(y,2)*sin(e1); + result += check_diff(e, x, d, 2); + + d = -2*b*e2*sin(e1)*pow(x,2) + e2 - y*sin(e1)*pow(x, 2); + result += check_diff(e, y, d); + + d = -2*b*pow(e2,2)*pow(x,4) + 2*b*pow(sin(e1),2)*pow(x,4) + - 2*sin(e1)*pow(x,2) - y*e2*pow(x,4); + result += check_diff(e, y, d, 2); return result; } - + // exp function static unsigned exam_differentiation3(void) { - unsigned result = 0; - symbol x("x"), y("y"), a("a"), b("b"); - ex e1, e2, e, d; + unsigned result = 0; + symbol x("x"), y("y"), a("a"), b("b"); + ex e1, e2, e, d; - // construct expression e to be diff'ed: - e1 = y*pow(x, 2) + a*x + b; - e2 = exp(e1); - e = b*pow(e2, 2) + y*e2 + a; - - d = 2*b*pow(e2, 2)*(2*x*y + a) + y*e2*(2*x*y + a); - result += check_diff(e, x, d); - - d = 4*b*pow(e2,2)*pow(2*y*x + a,2) + 4*b*pow(e2,2)*y - + 2*pow(y,2)*e2 + y*e2*pow(2*y*x + a,2); - result += check_diff(e, x, d, 2); - - d = 2*b*pow(e2,2)*pow(x,2) + e2 + y*e2*pow(x,2); - result += check_diff(e, y, d); - - d = 4*b*pow(e2,2)*pow(x,4) + 2*e2*pow(x,2) + y*e2*pow(x,4); - result += check_diff(e, y, d, 2); + // construct expression e to be diff'ed: + e1 = y*pow(x, 2) + a*x + b; + e2 = exp(e1); + e = b*pow(e2, 2) + y*e2 + a; + + d = 2*b*pow(e2, 2)*(2*x*y + a) + y*e2*(2*x*y + a); + result += check_diff(e, x, d); + + d = 4*b*pow(e2,2)*pow(2*y*x + a,2) + 4*b*pow(e2,2)*y + + 2*pow(y,2)*e2 + y*e2*pow(2*y*x + a,2); + result += check_diff(e, x, d, 2); + + d = 2*b*pow(e2,2)*pow(x,2) + e2 + y*e2*pow(x,2); + result += check_diff(e, y, d); + + d = 4*b*pow(e2,2)*pow(x,4) + 2*e2*pow(x,2) + y*e2*pow(x,4); + result += check_diff(e, y, d, 2); return result; } @@ -165,29 +165,29 @@ static unsigned exam_differentiation3(void) // log functions static unsigned exam_differentiation4(void) { - unsigned result = 0; - symbol x("x"), y("y"), a("a"), b("b"); - ex e1, e2, e, d; - - // construct expression e to be diff'ed: - e1 = y*pow(x, 2) + a*x + b; - e2 = log(e1); - e = b*pow(e2, 2) + y*e2 + a; - - d = 2*b*e2*(2*x*y + a)/e1 + y*(2*x*y + a)/e1; - result += check_diff(e, x, d); - - d = 2*b*pow((2*x*y + a),2)*pow(e1,-2) + 4*b*y*e2/e1 - - 2*b*e2*pow(2*x*y + a,2)*pow(e1,-2) + 2*pow(y,2)/e1 - - y*pow(2*x*y + a,2)*pow(e1,-2); - result += check_diff(e, x, d, 2); - - d = 2*b*e2*pow(x,2)/e1 + e2 + y*pow(x,2)/e1; - result += check_diff(e, y, d); - - d = 2*b*pow(x,4)*pow(e1,-2) - 2*b*e2*pow(e1,-2)*pow(x,4) - + 2*pow(x,2)/e1 - y*pow(x,4)*pow(e1,-2); - result += check_diff(e, y, d, 2); + unsigned result = 0; + symbol x("x"), y("y"), a("a"), b("b"); + ex e1, e2, e, d; + + // construct expression e to be diff'ed: + e1 = y*pow(x, 2) + a*x + b; + e2 = log(e1); + e = b*pow(e2, 2) + y*e2 + a; + + d = 2*b*e2*(2*x*y + a)/e1 + y*(2*x*y + a)/e1; + result += check_diff(e, x, d); + + d = 2*b*pow((2*x*y + a),2)*pow(e1,-2) + 4*b*y*e2/e1 + - 2*b*e2*pow(2*x*y + a,2)*pow(e1,-2) + 2*pow(y,2)/e1 + - y*pow(2*x*y + a,2)*pow(e1,-2); + result += check_diff(e, x, d, 2); + + d = 2*b*e2*pow(x,2)/e1 + e2 + y*pow(x,2)/e1; + result += check_diff(e, y, d); + + d = 2*b*pow(x,4)*pow(e1,-2) - 2*b*e2*pow(e1,-2)*pow(x,4) + + 2*pow(x,2)/e1 - y*pow(x,4)*pow(e1,-2); + result += check_diff(e, y, d, 2); return result; } @@ -195,102 +195,102 @@ static unsigned exam_differentiation4(void) // Functions with two variables static unsigned exam_differentiation5(void) { - unsigned result = 0; - symbol x("x"), y("y"), a("a"), b("b"); - ex e1, e2, e, d; - - // test atan2 - e1 = y*pow(x, 2) + a*x + b; - e2 = x*pow(y, 2) + b*y + a; - e = atan2(e1,e2); - /* - d = pow(y,2)*(-b-y*pow(x,2)-a*x)/(pow(b+y*pow(x,2)+a*x,2)+pow(x*pow(y,2)+b*y+a,2)) - +(2*y*x+a)/((x*pow(y,2)+b*y+a)*(1+pow(b*y*pow(x,2)+a*x,2)/pow(x*pow(y,2)+b*y+a,2))); - */ - /* - d = ((a+2*y*x)*pow(y*b+pow(y,2)*x+a,-1)-(a*x+b+y*pow(x,2))* - pow(y*b+pow(y,2)*x+a,-2)*pow(y,2))* - pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1); - */ - /* - d = pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1) - *pow(y*b+pow(y,2)*x+a,-1)*(a+2*y*x) - +pow(y,2)*(-a*x-b-y*pow(x,2))* - pow(pow(y*b+pow(y,2)*x+a,2)+pow(a*x+b+y*pow(x,2),2),-1); - */ - d = pow(y,2)*pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)* - (-b-y*pow(x,2)-x*a)+ - pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)* - (y*b+pow(y,2)*x+a)*(2*y*x+a); - result += check_diff(e, x, d); - - return result; + unsigned result = 0; + symbol x("x"), y("y"), a("a"), b("b"); + ex e1, e2, e, d; + + // test atan2 + e1 = y*pow(x, 2) + a*x + b; + e2 = x*pow(y, 2) + b*y + a; + e = atan2(e1,e2); + /* + d = pow(y,2)*(-b-y*pow(x,2)-a*x)/(pow(b+y*pow(x,2)+a*x,2)+pow(x*pow(y,2)+b*y+a,2)) + +(2*y*x+a)/((x*pow(y,2)+b*y+a)*(1+pow(b*y*pow(x,2)+a*x,2)/pow(x*pow(y,2)+b*y+a,2))); + */ + /* + d = ((a+2*y*x)*pow(y*b+pow(y,2)*x+a,-1)-(a*x+b+y*pow(x,2))* + pow(y*b+pow(y,2)*x+a,-2)*pow(y,2))* + pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1); + */ + /* + d = pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1) + *pow(y*b+pow(y,2)*x+a,-1)*(a+2*y*x) + +pow(y,2)*(-a*x-b-y*pow(x,2))* + pow(pow(y*b+pow(y,2)*x+a,2)+pow(a*x+b+y*pow(x,2),2),-1); + */ + d = pow(y,2)*pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)* + (-b-y*pow(x,2)-x*a)+ + pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)* + (y*b+pow(y,2)*x+a)*(2*y*x+a); + result += check_diff(e, x, d); + + return result; } // Series static unsigned exam_differentiation6(void) { - symbol x("x"); - ex e, d, ed; - - e = sin(x).series(x==0, 8); - d = cos(x).series(x==0, 7); - ed = e.diff(x); - ed = series_to_poly(ed); - d = series_to_poly(d); - - if ((ed - d).compare(ex(0)) != 0) { - clog << "derivative of " << e << " by " << x << " returned " - << ed << " instead of " << d << ")" << endl; - return 1; - } - return 0; + symbol x("x"); + ex e, d, ed; + + e = sin(x).series(x==0, 8); + d = cos(x).series(x==0, 7); + ed = e.diff(x); + ed = series_to_poly(ed); + d = series_to_poly(d); + + if ((ed - d).compare(ex(0)) != 0) { + clog << "derivative of " << e << " by " << x << " returned " + << ed << " instead of " << d << ")" << endl; + return 1; + } + return 0; } // Hashing can help a lot, if differentiation is done cleverly static unsigned exam_differentiation7(void) { - symbol x("x"); - ex P = x + pow(x,3); - ex e = (P.diff(x) / P).diff(x, 2); - ex d = 6/P - 18*x/pow(P,2) - 54*pow(x,3)/pow(P,2) + 2/pow(P,3) - +18*pow(x,2)/pow(P,3) + 54*pow(x,4)/pow(P,3) + 54*pow(x,6)/pow(P,3); - - if (!(e-d).expand().is_zero()) { - clog << "expanded second derivative of " << (P.diff(x) / P) << " by " << x - << " returned " << e.expand() << " instead of " << d << endl; - return 1; - } - if (e.nops() > 3) { - clog << "second derivative of " << (P.diff(x) / P) << " by " << x - << " has " << e.nops() << " operands. " - << "The result is still correct but not optimal: 3 are enough! " - << "(Hint: maybe the product rule for objects of class mul should be more careful about assembling the result?)" << endl; - return 1; - } - return 0; + symbol x("x"); + ex P = x + pow(x,3); + ex e = (P.diff(x) / P).diff(x, 2); + ex d = 6/P - 18*x/pow(P,2) - 54*pow(x,3)/pow(P,2) + 2/pow(P,3) + +18*pow(x,2)/pow(P,3) + 54*pow(x,4)/pow(P,3) + 54*pow(x,6)/pow(P,3); + + if (!(e-d).expand().is_zero()) { + clog << "expanded second derivative of " << (P.diff(x) / P) << " by " << x + << " returned " << e.expand() << " instead of " << d << endl; + return 1; + } + if (e.nops() > 3) { + clog << "second derivative of " << (P.diff(x) / P) << " by " << x + << " has " << e.nops() << " operands. " + << "The result is still correct but not optimal: 3 are enough! " + << "(Hint: maybe the product rule for objects of class mul should be more careful about assembling the result?)" << endl; + return 1; + } + return 0; } unsigned exam_differentiation(void) { - unsigned result = 0; - - cout << "examining symbolic differentiation" << flush; - clog << "----------symbolic differentiation:" << endl; - - result += exam_differentiation1(); cout << '.' << flush; - result += exam_differentiation2(); cout << '.' << flush; - result += exam_differentiation3(); cout << '.' << flush; - result += exam_differentiation4(); cout << '.' << flush; - result += exam_differentiation5(); cout << '.' << flush; - result += exam_differentiation6(); cout << '.' << flush; - result += exam_differentiation7(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - return result; + unsigned result = 0; + + cout << "examining symbolic differentiation" << flush; + clog << "----------symbolic differentiation:" << endl; + + result += exam_differentiation1(); cout << '.' << flush; + result += exam_differentiation2(); cout << '.' << flush; + result += exam_differentiation3(); cout << '.' << flush; + result += exam_differentiation4(); cout << '.' << flush; + result += exam_differentiation5(); cout << '.' << flush; + result += exam_differentiation6(); cout << '.' << flush; + result += exam_differentiation7(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + return result; } diff --git a/check/exam_inifcns.cpp b/check/exam_inifcns.cpp index 823e25ec..baeacfd4 100644 --- a/check/exam_inifcns.cpp +++ b/check/exam_inifcns.cpp @@ -26,137 +26,137 @@ /* Assorted tests on other transcendental functions. */ static unsigned inifcns_consist_trans(void) { - unsigned result = 0; - symbol x("x"); - ex chk; - - chk = asin(1)-acos(0); - if (!chk.is_zero()) { - clog << "asin(1)-acos(0) erroneously returned " << chk - << " instead of 0" << endl; - ++result; - } - - // arbitrary check of type sin(f(x)): - chk = pow(sin(acos(x)),2) + pow(sin(asin(x)),2) - - (1+pow(x,2))*pow(sin(atan(x)),2); - if (chk != 1-pow(x,2)) { - clog << "sin(acos(x))^2 + sin(asin(x))^2 - (1+x^2)*sin(atan(x))^2 " - << "erroneously returned " << chk << " instead of 1-x^2" << endl; - ++result; - } - - // arbitrary check of type cos(f(x)): - chk = pow(cos(acos(x)),2) + pow(cos(asin(x)),2) - - (1+pow(x,2))*pow(cos(atan(x)),2); - if (!chk.is_zero()) { - clog << "cos(acos(x))^2 + cos(asin(x))^2 - (1+x^2)*cos(atan(x))^2 " - << "erroneously returned " << chk << " instead of 0" << endl; - ++result; - } - - // arbitrary check of type tan(f(x)): - chk = tan(acos(x))*tan(asin(x)) - tan(atan(x)); - if (chk != 1-x) { - clog << "tan(acos(x))*tan(asin(x)) - tan(atan(x)) " - << "erroneously returned " << chk << " instead of -x+1" << endl; - ++result; - } - - // arbitrary check of type sinh(f(x)): - chk = -pow(sinh(acosh(x)),2).expand()*pow(sinh(atanh(x)),2) - - pow(sinh(asinh(x)),2); - if (!chk.is_zero()) { - clog << "expand(-(sinh(acosh(x)))^2)*(sinh(atanh(x))^2) - sinh(asinh(x))^2 " - << "erroneously returned " << chk << " instead of 0" << endl; - ++result; - } - - // arbitrary check of type cosh(f(x)): - chk = (pow(cosh(asinh(x)),2) - 2*pow(cosh(acosh(x)),2)) - * pow(cosh(atanh(x)),2); - if (chk != 1) { - clog << "(cosh(asinh(x))^2 - 2*cosh(acosh(x))^2) * cosh(atanh(x))^2 " - << "erroneously returned " << chk << " instead of 1" << endl; - ++result; - } - - // arbitrary check of type tanh(f(x)): - chk = (pow(tanh(asinh(x)),-2) - pow(tanh(acosh(x)),2)).expand() - * pow(tanh(atanh(x)),2); - if (chk != 2) { - clog << "expand(tanh(acosh(x))^2 - tanh(asinh(x))^(-2)) * tanh(atanh(x))^2 " - << "erroneously returned " << chk << " instead of 2" << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x("x"); + ex chk; + + chk = asin(1)-acos(0); + if (!chk.is_zero()) { + clog << "asin(1)-acos(0) erroneously returned " << chk + << " instead of 0" << endl; + ++result; + } + + // arbitrary check of type sin(f(x)): + chk = pow(sin(acos(x)),2) + pow(sin(asin(x)),2) + - (1+pow(x,2))*pow(sin(atan(x)),2); + if (chk != 1-pow(x,2)) { + clog << "sin(acos(x))^2 + sin(asin(x))^2 - (1+x^2)*sin(atan(x))^2 " + << "erroneously returned " << chk << " instead of 1-x^2" << endl; + ++result; + } + + // arbitrary check of type cos(f(x)): + chk = pow(cos(acos(x)),2) + pow(cos(asin(x)),2) + - (1+pow(x,2))*pow(cos(atan(x)),2); + if (!chk.is_zero()) { + clog << "cos(acos(x))^2 + cos(asin(x))^2 - (1+x^2)*cos(atan(x))^2 " + << "erroneously returned " << chk << " instead of 0" << endl; + ++result; + } + + // arbitrary check of type tan(f(x)): + chk = tan(acos(x))*tan(asin(x)) - tan(atan(x)); + if (chk != 1-x) { + clog << "tan(acos(x))*tan(asin(x)) - tan(atan(x)) " + << "erroneously returned " << chk << " instead of -x+1" << endl; + ++result; + } + + // arbitrary check of type sinh(f(x)): + chk = -pow(sinh(acosh(x)),2).expand()*pow(sinh(atanh(x)),2) + - pow(sinh(asinh(x)),2); + if (!chk.is_zero()) { + clog << "expand(-(sinh(acosh(x)))^2)*(sinh(atanh(x))^2) - sinh(asinh(x))^2 " + << "erroneously returned " << chk << " instead of 0" << endl; + ++result; + } + + // arbitrary check of type cosh(f(x)): + chk = (pow(cosh(asinh(x)),2) - 2*pow(cosh(acosh(x)),2)) + * pow(cosh(atanh(x)),2); + if (chk != 1) { + clog << "(cosh(asinh(x))^2 - 2*cosh(acosh(x))^2) * cosh(atanh(x))^2 " + << "erroneously returned " << chk << " instead of 1" << endl; + ++result; + } + + // arbitrary check of type tanh(f(x)): + chk = (pow(tanh(asinh(x)),-2) - pow(tanh(acosh(x)),2)).expand() + * pow(tanh(atanh(x)),2); + if (chk != 2) { + clog << "expand(tanh(acosh(x))^2 - tanh(asinh(x))^(-2)) * tanh(atanh(x))^2 " + << "erroneously returned " << chk << " instead of 2" << endl; + ++result; + } + + return result; } /* Simple tests on the tgamma function. We stuff in arguments where the results * exists in closed form and check if it's ok. */ static unsigned inifcns_consist_gamma(void) { - unsigned result = 0; - ex e; - - e = tgamma(ex(1)); - for (int i=2; i<8; ++i) - e += tgamma(ex(i)); - if (e != numeric(874)) { - clog << "tgamma(1)+...+tgamma(7) erroneously returned " - << e << " instead of 874" << endl; - ++result; - } - - e = tgamma(ex(1)); - for (int i=2; i<8; ++i) - e *= tgamma(ex(i)); - if (e != numeric(24883200)) { - clog << "tgamma(1)*...*tgamma(7) erroneously returned " - << e << " instead of 24883200" << endl; - ++result; - } - - e = tgamma(ex(numeric(5, 2)))*tgamma(ex(numeric(9, 2)))*64; - if (e != 315*Pi) { - clog << "64*tgamma(5/2)*tgamma(9/2) erroneously returned " - << e << " instead of 315*Pi" << endl; - ++result; - } - - e = tgamma(ex(numeric(-13, 2))); - for (int i=-13; i<7; i=i+2) - e += tgamma(ex(numeric(i, 2))); - e = (e*tgamma(ex(numeric(15, 2)))*numeric(512)); - if (e != numeric(633935)*Pi) { - clog << "512*(tgamma(-13/2)+...+tgamma(5/2))*tgamma(15/2) erroneously returned " - << e << " instead of 633935*Pi" << endl; - ++result; - } - - return result; + unsigned result = 0; + ex e; + + e = tgamma(ex(1)); + for (int i=2; i<8; ++i) + e += tgamma(ex(i)); + if (e != numeric(874)) { + clog << "tgamma(1)+...+tgamma(7) erroneously returned " + << e << " instead of 874" << endl; + ++result; + } + + e = tgamma(ex(1)); + for (int i=2; i<8; ++i) + e *= tgamma(ex(i)); + if (e != numeric(24883200)) { + clog << "tgamma(1)*...*tgamma(7) erroneously returned " + << e << " instead of 24883200" << endl; + ++result; + } + + e = tgamma(ex(numeric(5, 2)))*tgamma(ex(numeric(9, 2)))*64; + if (e != 315*Pi) { + clog << "64*tgamma(5/2)*tgamma(9/2) erroneously returned " + << e << " instead of 315*Pi" << endl; + ++result; + } + + e = tgamma(ex(numeric(-13, 2))); + for (int i=-13; i<7; i=i+2) + e += tgamma(ex(numeric(i, 2))); + e = (e*tgamma(ex(numeric(15, 2)))*numeric(512)); + if (e != numeric(633935)*Pi) { + clog << "512*(tgamma(-13/2)+...+tgamma(5/2))*tgamma(15/2) erroneously returned " + << e << " instead of 633935*Pi" << endl; + ++result; + } + + return result; } /* Simple tests on the Psi-function (aka polygamma-function). We stuff in arguments where the result exists in closed form and check if it's ok. */ static unsigned inifcns_consist_psi(void) { - unsigned result = 0; - symbol x; - ex e, f; - - // We check psi(1) and psi(1/2) implicitly by calculating the curious - // little identity tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) == 2*log(2). - e += (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1)); - e -= (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1,2)); - if (e!=2*log(2)) { - clog << "tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) erroneously returned " - << e << " instead of 2*log(2)" << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x; + ex e, f; + + // We check psi(1) and psi(1/2) implicitly by calculating the curious + // little identity tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) == 2*log(2). + e += (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1)); + e -= (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1,2)); + if (e!=2*log(2)) { + clog << "tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) erroneously returned " + << e << " instead of 2*log(2)" << endl; + ++result; + } + + return result; } /* Simple tests on the Riemann Zeta function. We stuff in arguments where the @@ -164,47 +164,47 @@ static unsigned inifcns_consist_psi(void) * the Bernoulli numbers as a side effect. */ static unsigned inifcns_consist_zeta(void) { - unsigned result = 0; - ex e; - - for (int i=0; i<13; i+=2) - e += zeta(i)/pow(Pi,i); - if (e!=numeric(-204992279,638512875)) { - clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned " - << e << " instead of -204992279/638512875" << endl; - ++result; - } - - e = 0; - for (int i=-1; i>-16; i--) - e += zeta(i); - if (e!=numeric(487871,1633632)) { - clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned " - << e << " instead of 487871/1633632" << endl; - ++result; - } - - return result; + unsigned result = 0; + ex e; + + for (int i=0; i<13; i+=2) + e += zeta(i)/pow(Pi,i); + if (e!=numeric(-204992279,638512875)) { + clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned " + << e << " instead of -204992279/638512875" << endl; + ++result; + } + + e = 0; + for (int i=-1; i>-16; i--) + e += zeta(i); + if (e!=numeric(487871,1633632)) { + clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned " + << e << " instead of 487871/1633632" << endl; + ++result; + } + + return result; } unsigned exam_inifcns(void) { - unsigned result = 0; - - cout << "examining consistency of symbolic functions" << flush; - clog << "----------consistency of symbolic functions:" << endl; - - result += inifcns_consist_trans(); cout << '.' << flush; - result += inifcns_consist_gamma(); cout << '.' << flush; - result += inifcns_consist_psi(); cout << '.' << flush; - result += inifcns_consist_zeta(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + unsigned result = 0; + + cout << "examining consistency of symbolic functions" << flush; + clog << "----------consistency of symbolic functions:" << endl; + + result += inifcns_consist_trans(); cout << '.' << flush; + result += inifcns_consist_gamma(); cout << '.' << flush; + result += inifcns_consist_psi(); cout << '.' << flush; + result += inifcns_consist_zeta(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/exam_lsolve.cpp b/check/exam_lsolve.cpp index 04eb438b..cafda13a 100644 --- a/check/exam_lsolve.cpp +++ b/check/exam_lsolve.cpp @@ -24,185 +24,185 @@ static unsigned exam_lsolve1(void) { - // A trivial example. - unsigned result = 0; - symbol x("x"); - ex eq, aux; - - eq = (3*x+5 == numeric(8)); - aux = lsolve(eq, x); - if (aux != 1) { - ++result; - clog << "solution of 3*x+5==8 erroneously returned " - << aux << endl; - } - - return result; + // A trivial example. + unsigned result = 0; + symbol x("x"); + ex eq, aux; + + eq = (3*x+5 == numeric(8)); + aux = lsolve(eq, x); + if (aux != 1) { + ++result; + clog << "solution of 3*x+5==8 erroneously returned " + << aux << endl; + } + + return result; } static unsigned exam_lsolve2a(void) { - // An example from the Maple online help. - unsigned result = 0; - symbol a("a"), b("b"), x("x"), y("y"); - lst eqns, vars; - ex sol; - - // Create the linear system [a*x+b*y==3,x-y==b]... - eqns.append(a*x+b*y==3).append(x-y==b); - // ...to be solved for [x,y]... - vars.append(x).append(y); - // ...and solve it: - sol = lsolve(eqns, vars); - ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) - ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) - - // It should have returned [x==(3+b^2)/(a+b),y==(3-a*b)/(a+b)] - if (!normal(sol_x - (3+pow(b,2))/(a+b)).is_zero() || - !normal(sol_y - (3-a*b)/(a+b)).is_zero()) { - ++result; - clog << "solution of the system " << eqns << " for " << vars - << " erroneously returned " << sol << endl; - } - - return result; + // An example from the Maple online help. + unsigned result = 0; + symbol a("a"), b("b"), x("x"), y("y"); + lst eqns, vars; + ex sol; + + // Create the linear system [a*x+b*y==3,x-y==b]... + eqns.append(a*x+b*y==3).append(x-y==b); + // ...to be solved for [x,y]... + vars.append(x).append(y); + // ...and solve it: + sol = lsolve(eqns, vars); + ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) + ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) + + // It should have returned [x==(3+b^2)/(a+b),y==(3-a*b)/(a+b)] + if (!normal(sol_x - (3+pow(b,2))/(a+b)).is_zero() || + !normal(sol_y - (3-a*b)/(a+b)).is_zero()) { + ++result; + clog << "solution of the system " << eqns << " for " << vars + << " erroneously returned " << sol << endl; + } + + return result; } static unsigned exam_lsolve2b(void) { - // A boring example from Mathematica's online help. - unsigned result = 0; - symbol x("x"), y("y"); - lst eqns, vars; - ex sol; - - // Create the linear system [3*x+y==7,2*x-5*y==8]... - eqns.append(3*x+y==7).append(2*x-5*y==8); - // ...to be solved for [x,y]... - vars.append(x).append(y); - // ...and solve it: - sol = lsolve(eqns, vars); - ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) - ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) - - // It should have returned [x==43/17,y==-10/17] - if ((sol_x != numeric(43,17)) || - (sol_y != numeric(-10,17))) { - ++result; - clog << "solution of the system " << eqns << " for " << vars - << " erroneously returned " << sol << endl; - } - - return result; + // A boring example from Mathematica's online help. + unsigned result = 0; + symbol x("x"), y("y"); + lst eqns, vars; + ex sol; + + // Create the linear system [3*x+y==7,2*x-5*y==8]... + eqns.append(3*x+y==7).append(2*x-5*y==8); + // ...to be solved for [x,y]... + vars.append(x).append(y); + // ...and solve it: + sol = lsolve(eqns, vars); + ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) + ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) + + // It should have returned [x==43/17,y==-10/17] + if ((sol_x != numeric(43,17)) || + (sol_y != numeric(-10,17))) { + ++result; + clog << "solution of the system " << eqns << " for " << vars + << " erroneously returned " << sol << endl; + } + + return result; } static unsigned exam_lsolve2c(void) { - // A more interesting example from the Maple online help. - unsigned result = 0; - symbol x("x"), y("y"); - lst eqns, vars; - ex sol; - - // Create the linear system [I*x+y==1,I*x-y==2]... - eqns.append(I*x+y==1).append(I*x-y==2); - // ...to be solved for [x,y]... - vars.append(x).append(y); - // ...and solve it: - sol = lsolve(eqns, vars); - ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) - ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) - - // It should have returned [x==-3/2*I,y==-1/2] - if ((sol_x != numeric(-3,2)*I) || - (sol_y != numeric(-1,2))) { - ++result; - clog << "solution of the system " << eqns << " for " << vars - << " erroneously returned " << sol << endl; - } - - return result; + // A more interesting example from the Maple online help. + unsigned result = 0; + symbol x("x"), y("y"); + lst eqns, vars; + ex sol; + + // Create the linear system [I*x+y==1,I*x-y==2]... + eqns.append(I*x+y==1).append(I*x-y==2); + // ...to be solved for [x,y]... + vars.append(x).append(y); + // ...and solve it: + sol = lsolve(eqns, vars); + ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) + ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) + + // It should have returned [x==-3/2*I,y==-1/2] + if ((sol_x != numeric(-3,2)*I) || + (sol_y != numeric(-1,2))) { + ++result; + clog << "solution of the system " << eqns << " for " << vars + << " erroneously returned " << sol << endl; + } + + return result; } static unsigned exam_lsolve2S(void) { - // A degenerate example that went wrong in GiNaC 0.6.2. - unsigned result = 0; - symbol x("x"), y("y"), t("t"); - lst eqns, vars; - ex sol; - - // Create the linear system [0*x+0*y==0,0*x+1*y==t]... - eqns.append(0*x+0*y==0).append(0*x+1*y==t); - // ...to be solved for [x,y]... - vars.append(x).append(y); - // ...and solve it: - sol = lsolve(eqns, vars); - ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) - ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) - - // It should have returned [x==x,y==t] - if ((sol_x != x) || - (sol_y != t)) { - ++result; - clog << "solution of the system " << eqns << " for " << vars - << " erroneously returned " << sol << endl; - } - - return result; + // A degenerate example that went wrong in GiNaC 0.6.2. + unsigned result = 0; + symbol x("x"), y("y"), t("t"); + lst eqns, vars; + ex sol; + + // Create the linear system [0*x+0*y==0,0*x+1*y==t]... + eqns.append(0*x+0*y==0).append(0*x+1*y==t); + // ...to be solved for [x,y]... + vars.append(x).append(y); + // ...and solve it: + sol = lsolve(eqns, vars); + ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) + ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) + + // It should have returned [x==x,y==t] + if ((sol_x != x) || + (sol_y != t)) { + ++result; + clog << "solution of the system " << eqns << " for " << vars + << " erroneously returned " << sol << endl; + } + + return result; } static unsigned exam_lsolve3S(void) { - // A degenerate example that went wrong while trying to improve elimination - unsigned result = 0; - symbol b("b"), c("c"); - symbol x("x"), y("y"), z("z"); - lst eqns, vars; - ex sol; - - // Create the linear system [y+z==b,-y+z==c] with one additional row... - eqns.append(ex(0)==ex(0)).append(b==z+y).append(c==z-y); - // ...to be solved for [x,y,z]... - vars.append(x).append(y).append(z); - // ...and solve it: - sol = lsolve(eqns, vars); - ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) - ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) - ex sol_z = sol.op(2).rhs(); // rhs of solution for third variable (z) - - // It should have returned [x==x,y==t,] - if ((sol_x != x) || - (sol_y != (b-c)/2) || - (sol_z != (b+c)/2)) { - ++result; - clog << "solution of the system " << eqns << " for " << vars - << " erroneously returned " << sol << endl; - } - - return result; + // A degenerate example that went wrong while trying to improve elimination + unsigned result = 0; + symbol b("b"), c("c"); + symbol x("x"), y("y"), z("z"); + lst eqns, vars; + ex sol; + + // Create the linear system [y+z==b,-y+z==c] with one additional row... + eqns.append(ex(0)==ex(0)).append(b==z+y).append(c==z-y); + // ...to be solved for [x,y,z]... + vars.append(x).append(y).append(z); + // ...and solve it: + sol = lsolve(eqns, vars); + ex sol_x = sol.op(0).rhs(); // rhs of solution for first variable (x) + ex sol_y = sol.op(1).rhs(); // rhs of solution for second variable (y) + ex sol_z = sol.op(2).rhs(); // rhs of solution for third variable (z) + + // It should have returned [x==x,y==t,] + if ((sol_x != x) || + (sol_y != (b-c)/2) || + (sol_z != (b+c)/2)) { + ++result; + clog << "solution of the system " << eqns << " for " << vars + << " erroneously returned " << sol << endl; + } + + return result; } unsigned exam_lsolve(void) { - unsigned result = 0; - - cout << "examining linear solve" << flush; - clog << "----------linear solve:" << endl; - - result += exam_lsolve1(); cout << '.' << flush; - result += exam_lsolve2a(); cout << '.' << flush; - result += exam_lsolve2b(); cout << '.' << flush; - result += exam_lsolve2c(); cout << '.' << flush; - result += exam_lsolve2S(); cout << '.' << flush; - result += exam_lsolve3S(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + unsigned result = 0; + + cout << "examining linear solve" << flush; + clog << "----------linear solve:" << endl; + + result += exam_lsolve1(); cout << '.' << flush; + result += exam_lsolve2a(); cout << '.' << flush; + result += exam_lsolve2b(); cout << '.' << flush; + result += exam_lsolve2c(); cout << '.' << flush; + result += exam_lsolve2S(); cout << '.' << flush; + result += exam_lsolve3S(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp index 9f17bb3d..a032bfe8 100644 --- a/check/exam_matrices.cpp +++ b/check/exam_matrices.cpp @@ -25,268 +25,268 @@ static unsigned matrix_determinants(void) { - unsigned result = 0; - ex det; - matrix m1(1,1), m2(2,2), m3(3,3), m4(4,4); - symbol a("a"), b("b"), c("c"); - symbol d("d"), e("e"), f("f"); - symbol g("g"), h("h"), i("i"); - - // check symbolic trivial matrix determinant - m1.set(0,0,a); - det = m1.determinant(); - if (det != a) { - clog << "determinant of 1x1 matrix " << m1 - << " erroneously returned " << det << endl; - ++result; - } - - // check generic dense symbolic 2x2 matrix determinant - m2.set(0,0,a).set(0,1,b); - m2.set(1,0,c).set(1,1,d); - det = m2.determinant(); - if (det != (a*d-b*c)) { - clog << "determinant of 2x2 matrix " << m2 - << " erroneously returned " << det << endl; - ++result; - } - - // check generic dense symbolic 3x3 matrix determinant - m3.set(0,0,a).set(0,1,b).set(0,2,c); - m3.set(1,0,d).set(1,1,e).set(1,2,f); - m3.set(2,0,g).set(2,1,h).set(2,2,i); - det = m3.determinant(); - if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) { - clog << "determinant of 3x3 matrix " << m3 - << " erroneously returned " << det << endl; - ++result; - } - - // check dense numeric 3x3 matrix determinant - m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3)); - m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2)); - m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2)); - det = m3.determinant(); - if (det != 42) { - clog << "determinant of 3x3 matrix " << m3 - << " erroneously returned " << det << endl; - ++result; - } - - // check dense symbolic 2x2 matrix determinant - m2.set(0,0,a/(a-b)).set(0,1,1); - m2.set(1,0,b/(a-b)).set(1,1,1); - det = m2.determinant(); - if (det != 1) { - if (det.normal() == 1) // only half wrong - clog << "determinant of 2x2 matrix " << m2 - << " was returned unnormalized as " << det << endl; - else // totally wrong - clog << "determinant of 2x2 matrix " << m2 - << " erroneously returned " << det << endl; - ++result; - } - - // check sparse symbolic 4x4 matrix determinant - m4.set(0,1,a).set(1,0,b).set(3,2,c).set(2,3,d); - det = m4.determinant(); - if (det != a*b*c*d) { - clog << "determinant of 4x4 matrix " << m4 - << " erroneously returned " << det << endl; - ++result; - } - - // check characteristic polynomial - m3.set(0,0,a).set(0,1,-2).set(0,2,2); - m3.set(1,0,3).set(1,1,a-1).set(1,2,2); - m3.set(2,0,3).set(2,1,4).set(2,2,a-3); - ex p = m3.charpoly(a); - if (p != 0) { - clog << "charpoly of 3x3 matrix " << m3 - << " erroneously returned " << p << endl; - ++result; - } - - return result; + unsigned result = 0; + ex det; + matrix m1(1,1), m2(2,2), m3(3,3), m4(4,4); + symbol a("a"), b("b"), c("c"); + symbol d("d"), e("e"), f("f"); + symbol g("g"), h("h"), i("i"); + + // check symbolic trivial matrix determinant + m1.set(0,0,a); + det = m1.determinant(); + if (det != a) { + clog << "determinant of 1x1 matrix " << m1 + << " erroneously returned " << det << endl; + ++result; + } + + // check generic dense symbolic 2x2 matrix determinant + m2.set(0,0,a).set(0,1,b); + m2.set(1,0,c).set(1,1,d); + det = m2.determinant(); + if (det != (a*d-b*c)) { + clog << "determinant of 2x2 matrix " << m2 + << " erroneously returned " << det << endl; + ++result; + } + + // check generic dense symbolic 3x3 matrix determinant + m3.set(0,0,a).set(0,1,b).set(0,2,c); + m3.set(1,0,d).set(1,1,e).set(1,2,f); + m3.set(2,0,g).set(2,1,h).set(2,2,i); + det = m3.determinant(); + if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) { + clog << "determinant of 3x3 matrix " << m3 + << " erroneously returned " << det << endl; + ++result; + } + + // check dense numeric 3x3 matrix determinant + m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3)); + m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2)); + m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2)); + det = m3.determinant(); + if (det != 42) { + clog << "determinant of 3x3 matrix " << m3 + << " erroneously returned " << det << endl; + ++result; + } + + // check dense symbolic 2x2 matrix determinant + m2.set(0,0,a/(a-b)).set(0,1,1); + m2.set(1,0,b/(a-b)).set(1,1,1); + det = m2.determinant(); + if (det != 1) { + if (det.normal() == 1) // only half wrong + clog << "determinant of 2x2 matrix " << m2 + << " was returned unnormalized as " << det << endl; + else // totally wrong + clog << "determinant of 2x2 matrix " << m2 + << " erroneously returned " << det << endl; + ++result; + } + + // check sparse symbolic 4x4 matrix determinant + m4.set(0,1,a).set(1,0,b).set(3,2,c).set(2,3,d); + det = m4.determinant(); + if (det != a*b*c*d) { + clog << "determinant of 4x4 matrix " << m4 + << " erroneously returned " << det << endl; + ++result; + } + + // check characteristic polynomial + m3.set(0,0,a).set(0,1,-2).set(0,2,2); + m3.set(1,0,3).set(1,1,a-1).set(1,2,2); + m3.set(2,0,3).set(2,1,4).set(2,2,a-3); + ex p = m3.charpoly(a); + if (p != 0) { + clog << "charpoly of 3x3 matrix " << m3 + << " erroneously returned " << p << endl; + ++result; + } + + return result; } static unsigned matrix_invert1(void) { - unsigned result = 0; - matrix m(1,1); - symbol a("a"); - - m.set(0,0,a); - matrix m_i = m.inverse(); - - if (m_i(0,0) != pow(a,-1)) { - clog << "inversion of 1x1 matrix " << m - << " erroneously returned " << m_i << endl; - ++result; - } - - return result; + unsigned result = 0; + matrix m(1,1); + symbol a("a"); + + m.set(0,0,a); + matrix m_i = m.inverse(); + + if (m_i(0,0) != pow(a,-1)) { + clog << "inversion of 1x1 matrix " << m + << " erroneously returned " << m_i << endl; + ++result; + } + + return result; } static unsigned matrix_invert2(void) { - unsigned result = 0; - matrix m(2,2); - symbol a("a"), b("b"), c("c"), d("d"); - m.set(0,0,a).set(0,1,b); - m.set(1,0,c).set(1,1,d); - matrix m_i = m.inverse(); - ex det = m.determinant(); - - if ((normal(m_i(0,0)*det) != d) || - (normal(m_i(0,1)*det) != -b) || - (normal(m_i(1,0)*det) != -c) || - (normal(m_i(1,1)*det) != a)) { - clog << "inversion of 2x2 matrix " << m - << " erroneously returned " << m_i << endl; - ++result; - } - - return result; + unsigned result = 0; + matrix m(2,2); + symbol a("a"), b("b"), c("c"), d("d"); + m.set(0,0,a).set(0,1,b); + m.set(1,0,c).set(1,1,d); + matrix m_i = m.inverse(); + ex det = m.determinant(); + + if ((normal(m_i(0,0)*det) != d) || + (normal(m_i(0,1)*det) != -b) || + (normal(m_i(1,0)*det) != -c) || + (normal(m_i(1,1)*det) != a)) { + clog << "inversion of 2x2 matrix " << m + << " erroneously returned " << m_i << endl; + ++result; + } + + return result; } static unsigned matrix_invert3(void) { - unsigned result = 0; - matrix m(3,3); - symbol a("a"), b("b"), c("c"); - symbol d("d"), e("e"), f("f"); - symbol g("g"), h("h"), i("i"); - m.set(0,0,a).set(0,1,b).set(0,2,c); - m.set(1,0,d).set(1,1,e).set(1,2,f); - m.set(2,0,g).set(2,1,h).set(2,2,i); - matrix m_i = m.inverse(); - ex det = m.determinant(); - - if ((normal(m_i(0,0)*det) != (e*i-f*h)) || - (normal(m_i(0,1)*det) != (c*h-b*i)) || - (normal(m_i(0,2)*det) != (b*f-c*e)) || - (normal(m_i(1,0)*det) != (f*g-d*i)) || - (normal(m_i(1,1)*det) != (a*i-c*g)) || - (normal(m_i(1,2)*det) != (c*d-a*f)) || - (normal(m_i(2,0)*det) != (d*h-e*g)) || - (normal(m_i(2,1)*det) != (b*g-a*h)) || - (normal(m_i(2,2)*det) != (a*e-b*d))) { - clog << "inversion of 3x3 matrix " << m - << " erroneously returned " << m_i << endl; - ++result; - } - - return result; + unsigned result = 0; + matrix m(3,3); + symbol a("a"), b("b"), c("c"); + symbol d("d"), e("e"), f("f"); + symbol g("g"), h("h"), i("i"); + m.set(0,0,a).set(0,1,b).set(0,2,c); + m.set(1,0,d).set(1,1,e).set(1,2,f); + m.set(2,0,g).set(2,1,h).set(2,2,i); + matrix m_i = m.inverse(); + ex det = m.determinant(); + + if ((normal(m_i(0,0)*det) != (e*i-f*h)) || + (normal(m_i(0,1)*det) != (c*h-b*i)) || + (normal(m_i(0,2)*det) != (b*f-c*e)) || + (normal(m_i(1,0)*det) != (f*g-d*i)) || + (normal(m_i(1,1)*det) != (a*i-c*g)) || + (normal(m_i(1,2)*det) != (c*d-a*f)) || + (normal(m_i(2,0)*det) != (d*h-e*g)) || + (normal(m_i(2,1)*det) != (b*g-a*h)) || + (normal(m_i(2,2)*det) != (a*e-b*d))) { + clog << "inversion of 3x3 matrix " << m + << " erroneously returned " << m_i << endl; + ++result; + } + + return result; } static unsigned matrix_solve2(void) { - // check the solution of the multiple system A*X = B: - // [ 1 2 -1 ] [ x0 y0 ] [ 4 0 ] - // [ 1 4 -2 ]*[ x1 y1 ] = [ 7 0 ] - // [ a -2 2 ] [ x2 y2 ] [ a 4 ] - unsigned result = 0; - symbol a("a"); - symbol x0("x0"), x1("x1"), x2("x2"); - symbol y0("y0"), y1("y1"), y2("y2"); - matrix A(3,3); - A.set(0,0,1).set(0,1,2).set(0,2,-1); - A.set(1,0,1).set(1,1,4).set(1,2,-2); - A.set(2,0,a).set(2,1,-2).set(2,2,2); - matrix B(3,2); - B.set(0,0,4).set(1,0,7).set(2,0,a); - B.set(0,1,0).set(1,1,0).set(2,1,4); - matrix X(3,2); - X.set(0,0,x0).set(1,0,x1).set(2,0,x2); - X.set(0,1,y0).set(1,1,y1).set(2,1,y2); - matrix cmp(3,2); - cmp.set(0,0,1).set(1,0,3).set(2,0,3); - cmp.set(0,1,0).set(1,1,2).set(2,1,4); - matrix sol(A.solve(X, B)); - for (unsigned ro=0; ro<3; ++ro) - for (unsigned co=0; co<2; ++co) - if (cmp(ro,co) != sol(ro,co)) - result = 1; - if (result) { - clog << "Solving " << A << " * " << X << " == " << B << endl - << "erroneously returned " << sol << endl; - } - - return result; + // check the solution of the multiple system A*X = B: + // [ 1 2 -1 ] [ x0 y0 ] [ 4 0 ] + // [ 1 4 -2 ]*[ x1 y1 ] = [ 7 0 ] + // [ a -2 2 ] [ x2 y2 ] [ a 4 ] + unsigned result = 0; + symbol a("a"); + symbol x0("x0"), x1("x1"), x2("x2"); + symbol y0("y0"), y1("y1"), y2("y2"); + matrix A(3,3); + A.set(0,0,1).set(0,1,2).set(0,2,-1); + A.set(1,0,1).set(1,1,4).set(1,2,-2); + A.set(2,0,a).set(2,1,-2).set(2,2,2); + matrix B(3,2); + B.set(0,0,4).set(1,0,7).set(2,0,a); + B.set(0,1,0).set(1,1,0).set(2,1,4); + matrix X(3,2); + X.set(0,0,x0).set(1,0,x1).set(2,0,x2); + X.set(0,1,y0).set(1,1,y1).set(2,1,y2); + matrix cmp(3,2); + cmp.set(0,0,1).set(1,0,3).set(2,0,3); + cmp.set(0,1,0).set(1,1,2).set(2,1,4); + matrix sol(A.solve(X, B)); + for (unsigned ro=0; ro<3; ++ro) + for (unsigned co=0; co<2; ++co) + if (cmp(ro,co) != sol(ro,co)) + result = 1; + if (result) { + clog << "Solving " << A << " * " << X << " == " << B << endl + << "erroneously returned " << sol << endl; + } + + return result; } static unsigned matrix_misc(void) { - unsigned result = 0; - matrix m1(2,2); - symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f"); - m1.set(0,0,a).set(0,1,b); - m1.set(1,0,c).set(1,1,d); - ex tr = trace(m1); - - // check a simple trace - if (tr.compare(a+d)) { - clog << "trace of 2x2 matrix " << m1 - << " erroneously returned " << tr << endl; - ++result; - } - - // and two simple transpositions - matrix m2 = transpose(m1); - if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) { - clog << "transpose of 2x2 matrix " << m1 - << " erroneously returned " << m2 << endl; - ++result; - } - matrix m3(3,2); - m3.set(0,0,a).set(0,1,b); - m3.set(1,0,c).set(1,1,d); - m3.set(2,0,e).set(2,1,f); - if (transpose(transpose(m3)) != m3) { - clog << "transposing 3x2 matrix " << m3 << " twice" - << " erroneously returned " << transpose(transpose(m3)) << endl; - ++result; - } - - // produce a runtime-error by inverting a singular matrix and catch it - matrix m4(2,2); - matrix m5; - bool caught = false; - try { - m5 = inverse(m4); - } catch (std::runtime_error err) { - caught = true; - } - if (!caught) { - cerr << "singular 2x2 matrix " << m4 - << " erroneously inverted to " << m5 << endl; - ++result; - } - - return result; + unsigned result = 0; + matrix m1(2,2); + symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f"); + m1.set(0,0,a).set(0,1,b); + m1.set(1,0,c).set(1,1,d); + ex tr = trace(m1); + + // check a simple trace + if (tr.compare(a+d)) { + clog << "trace of 2x2 matrix " << m1 + << " erroneously returned " << tr << endl; + ++result; + } + + // and two simple transpositions + matrix m2 = transpose(m1); + if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) { + clog << "transpose of 2x2 matrix " << m1 + << " erroneously returned " << m2 << endl; + ++result; + } + matrix m3(3,2); + m3.set(0,0,a).set(0,1,b); + m3.set(1,0,c).set(1,1,d); + m3.set(2,0,e).set(2,1,f); + if (transpose(transpose(m3)) != m3) { + clog << "transposing 3x2 matrix " << m3 << " twice" + << " erroneously returned " << transpose(transpose(m3)) << endl; + ++result; + } + + // produce a runtime-error by inverting a singular matrix and catch it + matrix m4(2,2); + matrix m5; + bool caught = false; + try { + m5 = inverse(m4); + } catch (std::runtime_error err) { + caught = true; + } + if (!caught) { + cerr << "singular 2x2 matrix " << m4 + << " erroneously inverted to " << m5 << endl; + ++result; + } + + return result; } unsigned exam_matrices(void) { - unsigned result = 0; - - cout << "examining symbolic matrix manipulations" << flush; - clog << "----------symbolic matrix manipulations:" << endl; - - result += matrix_determinants(); cout << '.' << flush; - result += matrix_invert1(); cout << '.' << flush; - result += matrix_invert2(); cout << '.' << flush; - result += matrix_invert3(); cout << '.' << flush; - result += matrix_solve2(); cout << '.' << flush; - result += matrix_misc(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + unsigned result = 0; + + cout << "examining symbolic matrix manipulations" << flush; + clog << "----------symbolic matrix manipulations:" << endl; + + result += matrix_determinants(); cout << '.' << flush; + result += matrix_invert1(); cout << '.' << flush; + result += matrix_invert2(); cout << '.' << flush; + result += matrix_invert3(); cout << '.' << flush; + result += matrix_solve2(); cout << '.' << flush; + result += matrix_misc(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/exam_misc.cpp b/check/exam_misc.cpp index 86bd12ab..8052bda9 100644 --- a/check/exam_misc.cpp +++ b/check/exam_misc.cpp @@ -1,6 +1,7 @@ /** @file exam_misc.cpp * */ + /* * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany * @@ -25,27 +26,27 @@ #define VECSIZE 30 static unsigned exam_expand_subs(void) { - unsigned result = 0; - symbol a1("a1"); - symbol a[VECSIZE]; - ex e, aux; - - a[1] = a1; - for (unsigned i=0; i 5*5^(1/2) * are carried out properly. */ static unsigned exam_numeric5(void) { - 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; + 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; } unsigned exam_numeric(void) { - 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; - - 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; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp index 2c6fe6ae..c4383c7c 100644 --- a/check/exam_paranoia.cpp +++ b/check/exam_paranoia.cpp @@ -30,30 +30,30 @@ // actually introduced the second.) static unsigned exam_paranoia1(void) { - unsigned result = 0; - symbol x("x"), y("y"), z("z"); - ex e, f, g; - - e = x * y * z; - f = y * z; - g = e / f; - - // In the first one expand did not do any job at all: - if (!g.expand().is_equal(x)) { - clog << "e = x*y*z; f = y*z; expand(e/f) erroneously returned " - << g.expand() << endl; - ++result; - } - - // This one somehow used to return 0: - e = pow(x + 1, -1); - if (!e.expand().is_equal(e)) { - clog << "expand(pow(x + 1, -1)) erroneously returned " - << e.expand() << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x("x"), y("y"), z("z"); + ex e, f, g; + + e = x * y * z; + f = y * z; + g = e / f; + + // In the first one expand did not do any job at all: + if (!g.expand().is_equal(x)) { + clog << "e = x*y*z; f = y*z; expand(e/f) erroneously returned " + << g.expand() << endl; + ++result; + } + + // This one somehow used to return 0: + e = pow(x + 1, -1); + if (!e.expand().is_equal(e)) { + clog << "expand(pow(x + 1, -1)) erroneously returned " + << e.expand() << endl; + ++result; + } + + return result; } // And here the second oops which showed up until May 17th 1999. It had to do @@ -61,35 +61,35 @@ static unsigned exam_paranoia1(void) // had the names as given here: static unsigned exam_paranoia2(void) { - unsigned result = 0; - symbol x("x"), y("y"), z("z"); - ex e, f, g; - - e = x + z*x; - f = e*y; - g = f - e*y; - - // After .expand(), g should be zero: - if (!g.expand().is_zero()) { - clog << "e = (x + z*x); f = e*y; expand(f - e*y) erroneously returned " - << g.expand() << endl; - ++result; - } - // After .eval(), g should be zero: - if (!g.eval().is_zero()) { - clog << "e = (x + z*x); f = e*y; eval(f - e*y) erroneously returned " - << g.eval() << endl; - ++result; - } - // This actually worked already back in April 1999. - // But we are *very* paranoic! - if (!g.expand().eval().is_zero()) { - clog << "e = (x + z*x); f = e*y; eval(expand(f - e*y)) erroneously returned " - << g.expand().eval() << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x("x"), y("y"), z("z"); + ex e, f, g; + + e = x + z*x; + f = e*y; + g = f - e*y; + + // After .expand(), g should be zero: + if (!g.expand().is_zero()) { + clog << "e = (x + z*x); f = e*y; expand(f - e*y) erroneously returned " + << g.expand() << endl; + ++result; + } + // After .eval(), g should be zero: + if (!g.eval().is_zero()) { + clog << "e = (x + z*x); f = e*y; eval(f - e*y) erroneously returned " + << g.eval() << endl; + ++result; + } + // This actually worked already back in April 1999. + // But we are *very* paranoic! + if (!g.expand().eval().is_zero()) { + clog << "e = (x + z*x); f = e*y; eval(expand(f - e*y)) erroneously returned " + << g.expand().eval() << endl; + ++result; + } + + return result; } // The third bug was introduced on May 18th 1999, discovered on May 19 and @@ -97,110 +97,110 @@ static unsigned exam_paranoia2(void) // other numbers: static unsigned exam_paranoia3(void) { - unsigned result = 0; - symbol x("x"), y("y"); - ex e, f; - - e = x*y - y; - f = e.subs(x == 2); - - if (!f.is_equal(y)) { - clog << "e = x*y - y; f = e.subs(x == 2) erroneously returned " - << f << endl; - ++result; - } - if (!f.eval().is_equal(y)) { - clog << "e = x*y - y; eval(e.subs(x == 2)) erroneously returned " - << f.eval() << endl; - ++result; - } - if (!f.expand().is_equal(y)) { - clog << "e = x*y - y; expand(e.subs(x == 2)) erroneously returned " - << f.expand() << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x("x"), y("y"); + ex e, f; + + e = x*y - y; + f = e.subs(x == 2); + + if (!f.is_equal(y)) { + clog << "e = x*y - y; f = e.subs(x == 2) erroneously returned " + << f << endl; + ++result; + } + if (!f.eval().is_equal(y)) { + clog << "e = x*y - y; eval(e.subs(x == 2)) erroneously returned " + << f.eval() << endl; + ++result; + } + if (!f.expand().is_equal(y)) { + clog << "e = x*y - y; expand(e.subs(x == 2)) erroneously returned " + << f.expand() << endl; + ++result; + } + + return result; } // The fourth bug was also discovered on May 19th 1999 and fixed immediately: static unsigned exam_paranoia4(void) { - unsigned result = 0; - symbol x("x"); - ex e, f, g; - - e = pow(x, 2) + x + 1; - f = pow(x, 2) + x + 1; - g = e - f; - - if (!g.is_zero()) { - clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned " - << g << endl; - ++result; - } - if (!g.is_zero()) { - clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g.eval() erroneously returned " - << g.eval() << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x("x"); + ex e, f, g; + + e = pow(x, 2) + x + 1; + f = pow(x, 2) + x + 1; + g = e - f; + + if (!g.is_zero()) { + clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned " + << g << endl; + ++result; + } + if (!g.is_zero()) { + clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g.eval() erroneously returned " + << g.eval() << endl; + ++result; + } + + return result; } // The fifth oops was discovered on May 20th 1999 and fixed a day later: static unsigned exam_paranoia5(void) { - unsigned result = 0; - symbol x("x"), y("y"); + unsigned result = 0; + symbol x("x"), y("y"); - ex e, f; - e = pow(x*y + 1, 2); - f = pow(x, 2) * pow(y, 2) + 2*x*y + 1; + ex e, f; + e = pow(x*y + 1, 2); + f = pow(x, 2) * pow(y, 2) + 2*x*y + 1; - if (!(e-f).expand().is_zero()) { - clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned " - << (e-f).expand() << endl; - ++result; - } + if (!(e-f).expand().is_zero()) { + clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned " + << (e-f).expand() << endl; + ++result; + } - return result; + return result; } // This one was discovered on Jun 1st 1999 and fixed the same day: static unsigned exam_paranoia6(void) { - unsigned result = 0; - symbol x("x"); - - ex e, f; - e = pow(x, -5); - f = e.denom(); - - if (!f.is_equal(pow(x, 5))) { - clog << "e = pow(x, -5); f = e.denom(); f was " << f << " (should be x^5)" << endl; - ++result; - } - return result; + unsigned result = 0; + symbol x("x"); + + ex e, f; + e = pow(x, -5); + f = e.denom(); + + if (!f.is_equal(pow(x, 5))) { + clog << "e = pow(x, -5); f = e.denom(); f was " << f << " (should be x^5)" << endl; + ++result; + } + return result; } // This one was introduced on June 1st 1999 by some aggressive manual // optimization. Discovered and fixed on June 2nd. static unsigned exam_paranoia7(void) { - unsigned result = 0; - symbol x("x"), y("y"); - - ex e = y + y*x + 2; - ex f = expand(pow(e, 2) - (e*y*(x + 1))); - - if (f.nops() > 3) { - clog << "e=y+y*x+2; f=expand(pow(e,2)-(e*y*(x+1))) has " - << f.nops() << " arguments instead of 3 ( f==" - << f << " )" << endl; - ++result; - } - return result; + unsigned result = 0; + symbol x("x"), y("y"); + + ex e = y + y*x + 2; + ex f = expand(pow(e, 2) - (e*y*(x + 1))); + + if (f.nops() > 3) { + clog << "e=y+y*x+2; f=expand(pow(e,2)-(e*y*(x+1))) has " + << f.nops() << " arguments instead of 3 ( f==" + << f << " )" << endl; + ++result; + } + return result; } // This one was a result of the rewrite of mul::max_coefficient when we @@ -208,23 +208,23 @@ static unsigned exam_paranoia7(void) // 1999. Fixed on Oct 4th. static unsigned exam_paranoia8(void) { - unsigned result = 0; - symbol x("x"); - - ex e = -x / (x+1); - ex f; - - try { - f = e.normal(); - if (!f.is_equal(e)) { - clog << "normal(-x/(x+1)) returns " << f << " instead of -x/(x+1)\n"; - ++result; - } - } catch (const exception &err) { - clog << "normal(-x/(x+1) throws " << err.what() << endl; - ++result; - } - return result; + unsigned result = 0; + symbol x("x"); + + ex e = -x / (x+1); + ex f; + + try { + f = e.normal(); + if (!f.is_equal(e)) { + clog << "normal(-x/(x+1)) returns " << f << " instead of -x/(x+1)\n"; + ++result; + } + } catch (const exception &err) { + clog << "normal(-x/(x+1) throws " << err.what() << endl; + ++result; + } + return result; } // This one was a result of a modification to frac_cancel() & Co. to avoid @@ -234,17 +234,17 @@ static unsigned exam_paranoia8(void) // 2000 and fixed on Jan 31th. static unsigned exam_paranoia9(void) { - unsigned result = 0; - symbol x("x"); + unsigned result = 0; + symbol x("x"); - ex e = (exp(-x)-2*x*exp(-x)+pow(x,2)/2*exp(-x))/exp(-x); - ex f = e.normal(); + ex e = (exp(-x)-2*x*exp(-x)+pow(x,2)/2*exp(-x))/exp(-x); + ex f = e.normal(); - if (!f.is_equal(1-2*x+pow(x,2)/2)) { - clog << "normal(" << e << ") returns " << f << " instead of 1-2*x+1/2*x^2\n"; - ++result; - } - return result; + if (!f.is_equal(1-2*x+pow(x,2)/2)) { + clog << "normal(" << e << ") returns " << f << " instead of 1-2*x+1/2*x^2\n"; + ++result; + } + return result; } // I have no idea when this broke. It has been working long ago, before 0.4.0 @@ -253,23 +253,23 @@ static unsigned exam_paranoia9(void) // It was fixed that same day. static unsigned exam_paranoia10(void) { - unsigned result = 0; - - ex b = numeric(2); - ex e = numeric(3,2); - ex r; - - try { - r = pow(b,e).eval(); - if (!(r-2*sqrt(ex(2))).is_zero()) { - clog << "2^(3/2) erroneously returned " << r << " instead of 2*sqrt(2)" << endl; - ++result; - } - } catch (const exception &err) { - clog << "2^(3/2) throws " << err.what() << endl; - ++result; - } - return result; + unsigned result = 0; + + ex b = numeric(2); + ex e = numeric(3,2); + ex r; + + try { + r = pow(b,e).eval(); + if (!(r-2*sqrt(ex(2))).is_zero()) { + clog << "2^(3/2) erroneously returned " << r << " instead of 2*sqrt(2)" << endl; + ++result; + } + } catch (const exception &err) { + clog << "2^(3/2) throws " << err.what() << endl; + ++result; + } + return result; } // After the rewriting of basic::normal() & Co. to return {num, den} lists, @@ -278,7 +278,7 @@ static unsigned exam_paranoia10(void) // child (did you get this? Well, never mind...). Fixed on Feb 21th 2000. static unsigned exam_paranoia11(void) { - unsigned result = 0; + unsigned result = 0; symbol x("x"); ex e = ((-5-2*x)-((2-5*x)/(-2+x))*(3+2*x))/(5-4*x); @@ -289,7 +289,7 @@ static unsigned exam_paranoia11(void) clog << "normal(" << e << ") returns " << f << " instead of " << d << endl; ++result; } - return result; + return result; } // This one returned 0 because add::normal() incorrectly assumed that if the @@ -297,16 +297,16 @@ static unsigned exam_paranoia11(void) // be +/-1). Fixed on Aug 2nd 2000. static unsigned exam_paranoia12(void) { - unsigned result = 0; + unsigned result = 0; symbol x("x"); - + ex e = 2-2*(1+x)/(-1-x); ex f = e.normal(); ex d = 4; - + if (!(f - d).expand().is_zero()) { clog << "normal(" << e << ") returns " << f - << " instead of " << d << endl; + << " instead of " << d << endl; ++result; } return result; @@ -316,53 +316,53 @@ static unsigned exam_paranoia12(void) // input polynomials against 0. Fixed on Aug 4th 2000. static unsigned exam_paranoia13(void) { - unsigned result = 0; + unsigned result = 0; symbol a("a"), b("b"), c("c"); - + ex e = (b*a-c*a)/(4-a); ex d = (c*a-b*a)/(a-4); - - try { - ex f = e.normal(); - if (!(f - d).expand().is_zero()) { - clog << "normal(" << e << ") returns " << f - << " instead of " << d << endl; - ++result; - } - } catch (const exception &err) { - clog << "normal(" << e << ") throws " << err.what() << endl; - ++result; - } + + try { + ex f = e.normal(); + if (!(f - d).expand().is_zero()) { + clog << "normal(" << e << ") returns " << f + << " instead of " << d << endl; + ++result; + } + } catch (const exception &err) { + clog << "normal(" << e << ") throws " << err.what() << endl; + ++result; + } return result; } unsigned exam_paranoia(void) { - unsigned result = 0; - - cout << "examining several historic failures just out of paranoia" << flush; - clog << "----------several historic failures:" << endl; - - result += exam_paranoia1(); cout << '.' << flush; - result += exam_paranoia2(); cout << '.' << flush; - result += exam_paranoia3(); cout << '.' << flush; - result += exam_paranoia4(); cout << '.' << flush; - result += exam_paranoia5(); cout << '.' << flush; - result += exam_paranoia6(); cout << '.' << flush; - result += exam_paranoia7(); cout << '.' << flush; - result += exam_paranoia8(); cout << '.' << flush; - result += exam_paranoia9(); cout << '.' << flush; - result += exam_paranoia10(); cout << '.' << flush; - result += exam_paranoia11(); cout << '.' << flush; - result += exam_paranoia12(); cout << '.' << flush; - result += exam_paranoia13(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + unsigned result = 0; + + cout << "examining several historic failures just out of paranoia" << flush; + clog << "----------several historic failures:" << endl; + + result += exam_paranoia1(); cout << '.' << flush; + result += exam_paranoia2(); cout << '.' << flush; + result += exam_paranoia3(); cout << '.' << flush; + result += exam_paranoia4(); cout << '.' << flush; + result += exam_paranoia5(); cout << '.' << flush; + result += exam_paranoia6(); cout << '.' << flush; + result += exam_paranoia7(); cout << '.' << flush; + result += exam_paranoia8(); cout << '.' << flush; + result += exam_paranoia9(); cout << '.' << flush; + result += exam_paranoia10(); cout << '.' << flush; + result += exam_paranoia11(); cout << '.' << flush; + result += exam_paranoia12(); cout << '.' << flush; + result += exam_paranoia13(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/exam_polygcd.cpp b/check/exam_polygcd.cpp index 34aea62a..5680bec0 100644 --- a/check/exam_polygcd.cpp +++ b/check/exam_polygcd.cpp @@ -228,27 +228,27 @@ static unsigned poly_gcd7(void) unsigned exam_polygcd(void) { - unsigned result = 0; - + unsigned result = 0; + cout << "examining polynomial GCD computation" << flush; clog << "----------polynomial GCD computation:" << endl; - + result += poly_gcd1(); cout << '.' << flush; result += poly_gcd2(); cout << '.' << flush; result += poly_gcd3(); cout << '.' << flush; - result += poly_gcd3p(); cout << '.' << flush; // takes extremely long (PRS "worst" case) + result += poly_gcd3p(); cout << '.' << flush; // takes extremely long (PRS "worst" case) result += poly_gcd4(); cout << '.' << flush; result += poly_gcd5(); cout << '.' << flush; result += poly_gcd5p(); cout << '.' << flush; result += poly_gcd6(); cout << '.' << flush; result += poly_gcd7(); cout << '.' << flush; - + if (!result) { cout << " passed " << endl; - clog << "(no output)" << endl; - } else { + clog << "(no output)" << endl; + } else { cout << " failed " << endl; - } - + } + return result; } diff --git a/check/exam_powerlaws.cpp b/check/exam_powerlaws.cpp index aba615fd..f685ffb1 100644 --- a/check/exam_powerlaws.cpp +++ b/check/exam_powerlaws.cpp @@ -25,284 +25,284 @@ static unsigned exam_powerlaws1(void) { - // (x^a)^b = x^(a*b) - - symbol x("x"); - symbol a("a"); - symbol b("b"); - - ex e1=power(power(x,a),b); - if (!(is_ex_exactly_of_type(e1,power) && - is_ex_exactly_of_type(e1.op(0),power) && - is_ex_exactly_of_type(e1.op(0).op(0),symbol) && - is_ex_exactly_of_type(e1.op(0).op(1),symbol) && - is_ex_exactly_of_type(e1.op(1),symbol) && - e1.is_equal(power(power(x,a),b)) )) { - clog << "(x^a)^b, x,a,b symbolic wrong" << endl; - clog << "returned: " << e1 << endl; - return 1; - } - - ex e2=e1.subs(a==1); - if (!(is_ex_exactly_of_type(e2,power) && - is_ex_exactly_of_type(e2.op(0),symbol) && - is_ex_exactly_of_type(e2.op(1),symbol) && - e2.is_equal(power(x,b)) )) { - clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl; - clog << "returned: " << e2 << endl; - return 1; - } - - ex e3=e1.subs(a==-1); - if (!(is_ex_exactly_of_type(e3,power) && - is_ex_exactly_of_type(e3.op(0),power) && - is_ex_exactly_of_type(e3.op(0).op(0),symbol) && - is_ex_exactly_of_type(e3.op(0).op(1),numeric) && - is_ex_exactly_of_type(e3.op(1),symbol) && - e3.is_equal(power(power(x,-1),b)) )) { - clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl; - clog << "returned: " << e3 << endl; - return 1; - } - - ex e4=e1.subs(lst(a==-1,b==2.5)); - if (!(is_ex_exactly_of_type(e4,power) && - is_ex_exactly_of_type(e4.op(0),power) && - is_ex_exactly_of_type(e4.op(0).op(0),symbol) && - is_ex_exactly_of_type(e4.op(0).op(1),numeric) && - is_ex_exactly_of_type(e4.op(1),numeric) && - e4.is_equal(power(power(x,-1),2.5)) )) { - clog << "(x^a)^b, x symbolic, a==-1, b==2.5 wrong" << endl; - clog << "returned: " << e4 << endl; - return 1; - } - - ex e5=e1.subs(lst(a==-0.9,b==2.5)); - if (!(is_ex_exactly_of_type(e5,power) && - is_ex_exactly_of_type(e5.op(0),symbol) && - is_ex_exactly_of_type(e5.op(1),numeric) && - e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) { - clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl; - clog << "returned: " << e5 << endl; - return 1; - } - - ex e6=e1.subs(lst(a==numeric(3)+numeric(5.3)*I,b==-5)); - if (!(is_ex_exactly_of_type(e6,power) && - is_ex_exactly_of_type(e6.op(0),symbol) && - is_ex_exactly_of_type(e6.op(1),numeric) && - e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) { - clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl; - clog << "returned: " << e6 << endl; - return 1; - } - - return 0; + // (x^a)^b = x^(a*b) + + symbol x("x"); + symbol a("a"); + symbol b("b"); + + ex e1=power(power(x,a),b); + if (!(is_ex_exactly_of_type(e1,power) && + is_ex_exactly_of_type(e1.op(0),power) && + is_ex_exactly_of_type(e1.op(0).op(0),symbol) && + is_ex_exactly_of_type(e1.op(0).op(1),symbol) && + is_ex_exactly_of_type(e1.op(1),symbol) && + e1.is_equal(power(power(x,a),b)) )) { + clog << "(x^a)^b, x,a,b symbolic wrong" << endl; + clog << "returned: " << e1 << endl; + return 1; + } + + ex e2=e1.subs(a==1); + if (!(is_ex_exactly_of_type(e2,power) && + is_ex_exactly_of_type(e2.op(0),symbol) && + is_ex_exactly_of_type(e2.op(1),symbol) && + e2.is_equal(power(x,b)) )) { + clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl; + clog << "returned: " << e2 << endl; + return 1; + } + + ex e3=e1.subs(a==-1); + if (!(is_ex_exactly_of_type(e3,power) && + is_ex_exactly_of_type(e3.op(0),power) && + is_ex_exactly_of_type(e3.op(0).op(0),symbol) && + is_ex_exactly_of_type(e3.op(0).op(1),numeric) && + is_ex_exactly_of_type(e3.op(1),symbol) && + e3.is_equal(power(power(x,-1),b)) )) { + clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl; + clog << "returned: " << e3 << endl; + return 1; + } + + ex e4=e1.subs(lst(a==-1,b==2.5)); + if (!(is_ex_exactly_of_type(e4,power) && + is_ex_exactly_of_type(e4.op(0),power) && + is_ex_exactly_of_type(e4.op(0).op(0),symbol) && + is_ex_exactly_of_type(e4.op(0).op(1),numeric) && + is_ex_exactly_of_type(e4.op(1),numeric) && + e4.is_equal(power(power(x,-1),2.5)) )) { + clog << "(x^a)^b, x symbolic, a==-1, b==2.5 wrong" << endl; + clog << "returned: " << e4 << endl; + return 1; + } + + ex e5=e1.subs(lst(a==-0.9,b==2.5)); + if (!(is_ex_exactly_of_type(e5,power) && + is_ex_exactly_of_type(e5.op(0),symbol) && + is_ex_exactly_of_type(e5.op(1),numeric) && + e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) { + clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl; + clog << "returned: " << e5 << endl; + return 1; + } + + ex e6=e1.subs(lst(a==numeric(3)+numeric(5.3)*I,b==-5)); + if (!(is_ex_exactly_of_type(e6,power) && + is_ex_exactly_of_type(e6.op(0),symbol) && + is_ex_exactly_of_type(e6.op(1),numeric) && + e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) { + clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl; + clog << "returned: " << e6 << endl; + return 1; + } + + return 0; } static unsigned exam_powerlaws2(void) { - // (a*x)^b = a^b * x^b - - symbol x("x"); - symbol a("a"); - symbol b("b"); - - ex e1=power(a*x,b); - if (!(is_ex_exactly_of_type(e1,power) && - is_ex_exactly_of_type(e1.op(0),mul) && - (e1.op(0).nops()==2) && - is_ex_exactly_of_type(e1.op(0).op(0),symbol) && - is_ex_exactly_of_type(e1.op(0).op(1),symbol) && - is_ex_exactly_of_type(e1.op(1),symbol) && - e1.is_equal(power(a*x,b)) )) { - clog << "(a*x)^b, x,a,b symbolic wrong" << endl; - clog << "returned: " << e1 << endl; - return 1; - } - - ex e2=e1.subs(a==3); - if (!(is_ex_exactly_of_type(e2,power) && - is_ex_exactly_of_type(e2.op(0),mul) && - (e2.op(0).nops()==2) && - is_ex_exactly_of_type(e2.op(0).op(0),symbol) && - is_ex_exactly_of_type(e2.op(0).op(1),numeric) && - is_ex_exactly_of_type(e2.op(1),symbol) && - e2.is_equal(power(3*x,b)) )) { - clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl; - clog << "returned: " << e2 << endl; - return 1; - } - - ex e3=e1.subs(b==-3); - if (!(is_ex_exactly_of_type(e3,mul) && - (e3.nops()==2) && - is_ex_exactly_of_type(e3.op(0),power) && - is_ex_exactly_of_type(e3.op(1),power) && - e3.is_equal(power(a,-3)*power(x,-3)) )) { - clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl; - clog << "returned: " << e3 << endl; - return 1; - } - - ex e4=e1.subs(b==4.5); - if (!(is_ex_exactly_of_type(e4,power) && - is_ex_exactly_of_type(e4.op(0),mul) && - (e4.op(0).nops()==2) && - is_ex_exactly_of_type(e4.op(0).op(0),symbol) && - is_ex_exactly_of_type(e4.op(0).op(1),symbol) && - is_ex_exactly_of_type(e4.op(1),numeric) && - e4.is_equal(power(a*x,4.5)) )) { - clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl; - clog << "returned: " << e4 << endl; - return 1; - } - - ex e5=e1.subs(lst(a==3.2,b==3+numeric(5)*I)); - if (!(is_ex_exactly_of_type(e5,mul) && - (e5.nops()==2) && - is_ex_exactly_of_type(e5.op(0),power) && - is_ex_exactly_of_type(e5.op(1),numeric) && - e5.is_equal(power(x,3+numeric(5)*I)* - power(numeric(3.2),3+numeric(5)*I)) )) { - clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl; - clog << "returned: " << e5 << endl; - return 1; - } - - ex e6=e1.subs(lst(a==-3.2,b==3+numeric(5)*I)); - if (!(is_ex_exactly_of_type(e6,mul) && - (e6.nops()==2) && - is_ex_exactly_of_type(e6.op(0),power) && - is_ex_exactly_of_type(e6.op(1),numeric) && - e6.is_equal(power(-x,3+numeric(5)*I)* - power(numeric(3.2),3+numeric(5)*I)) )) { - clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl; - clog << "returned: " << e6 << endl; - return 1; - } - - ex e7=e1.subs(lst(a==3+numeric(5)*I,b==3.2)); - if (!(is_ex_exactly_of_type(e7,power) && - is_ex_exactly_of_type(e7.op(0),mul) && - (e7.op(0).nops()==2) && - is_ex_exactly_of_type(e7.op(0).op(0),symbol) && - is_ex_exactly_of_type(e7.op(0).op(1),numeric) && - is_ex_exactly_of_type(e7.op(1),numeric) && - e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) { - clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl; - clog << "returned: " << e7 << endl; - return 1; - } - - return 0; + // (a*x)^b = a^b * x^b + + symbol x("x"); + symbol a("a"); + symbol b("b"); + + ex e1=power(a*x,b); + if (!(is_ex_exactly_of_type(e1,power) && + is_ex_exactly_of_type(e1.op(0),mul) && + (e1.op(0).nops()==2) && + is_ex_exactly_of_type(e1.op(0).op(0),symbol) && + is_ex_exactly_of_type(e1.op(0).op(1),symbol) && + is_ex_exactly_of_type(e1.op(1),symbol) && + e1.is_equal(power(a*x,b)) )) { + clog << "(a*x)^b, x,a,b symbolic wrong" << endl; + clog << "returned: " << e1 << endl; + return 1; + } + + ex e2=e1.subs(a==3); + if (!(is_ex_exactly_of_type(e2,power) && + is_ex_exactly_of_type(e2.op(0),mul) && + (e2.op(0).nops()==2) && + is_ex_exactly_of_type(e2.op(0).op(0),symbol) && + is_ex_exactly_of_type(e2.op(0).op(1),numeric) && + is_ex_exactly_of_type(e2.op(1),symbol) && + e2.is_equal(power(3*x,b)) )) { + clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl; + clog << "returned: " << e2 << endl; + return 1; + } + + ex e3=e1.subs(b==-3); + if (!(is_ex_exactly_of_type(e3,mul) && + (e3.nops()==2) && + is_ex_exactly_of_type(e3.op(0),power) && + is_ex_exactly_of_type(e3.op(1),power) && + e3.is_equal(power(a,-3)*power(x,-3)) )) { + clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl; + clog << "returned: " << e3 << endl; + return 1; + } + + ex e4=e1.subs(b==4.5); + if (!(is_ex_exactly_of_type(e4,power) && + is_ex_exactly_of_type(e4.op(0),mul) && + (e4.op(0).nops()==2) && + is_ex_exactly_of_type(e4.op(0).op(0),symbol) && + is_ex_exactly_of_type(e4.op(0).op(1),symbol) && + is_ex_exactly_of_type(e4.op(1),numeric) && + e4.is_equal(power(a*x,4.5)) )) { + clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl; + clog << "returned: " << e4 << endl; + return 1; + } + + ex e5=e1.subs(lst(a==3.2,b==3+numeric(5)*I)); + if (!(is_ex_exactly_of_type(e5,mul) && + (e5.nops()==2) && + is_ex_exactly_of_type(e5.op(0),power) && + is_ex_exactly_of_type(e5.op(1),numeric) && + e5.is_equal(power(x,3+numeric(5)*I)* + power(numeric(3.2),3+numeric(5)*I)) )) { + clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl; + clog << "returned: " << e5 << endl; + return 1; + } + + ex e6=e1.subs(lst(a==-3.2,b==3+numeric(5)*I)); + if (!(is_ex_exactly_of_type(e6,mul) && + (e6.nops()==2) && + is_ex_exactly_of_type(e6.op(0),power) && + is_ex_exactly_of_type(e6.op(1),numeric) && + e6.is_equal(power(-x,3+numeric(5)*I)* + power(numeric(3.2),3+numeric(5)*I)) )) { + clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl; + clog << "returned: " << e6 << endl; + return 1; + } + + ex e7=e1.subs(lst(a==3+numeric(5)*I,b==3.2)); + if (!(is_ex_exactly_of_type(e7,power) && + is_ex_exactly_of_type(e7.op(0),mul) && + (e7.op(0).nops()==2) && + is_ex_exactly_of_type(e7.op(0).op(0),symbol) && + is_ex_exactly_of_type(e7.op(0).op(1),numeric) && + is_ex_exactly_of_type(e7.op(1),numeric) && + e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) { + clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl; + clog << "returned: " << e7 << endl; + return 1; + } + + return 0; } static unsigned exam_powerlaws3(void) { - // numeric evaluation + // numeric evaluation - ex e1 = power(numeric(4),numeric(1,2)); - if (e1 != 2) { - clog << "4^(1/2) wrongly returned " << e1 << endl; - return 1; - } - - ex e2 = power(numeric(27),numeric(2,3)); - if (e2 != 9) { - clog << "27^(2/3) wrongly returned " << e2 << endl; - return 1; - } - - ex e3 = power(numeric(5),numeric(1,2)); - if (!(is_ex_exactly_of_type(e3,power) && - e3.op(0).is_equal(numeric(5)) && - e3.op(1).is_equal(numeric(1,2)))) { - clog << "5^(1/2) wrongly returned " << e3 << endl; - return 1; - } - - ex e4 = power(numeric(5),evalf(numeric(1,2))); - if (!(is_ex_exactly_of_type(e4,numeric))) { - clog << "5^(0.5) wrongly returned " << e4 << endl; - return 1; - } - - ex e5 = power(evalf(numeric(5)),numeric(1,2)); - if (!(is_ex_exactly_of_type(e5,numeric))) { - clog << "5.0^(1/2) wrongly returned " << e5 << endl; - return 1; - } - - return 0; + ex e1 = power(numeric(4),numeric(1,2)); + if (e1 != 2) { + clog << "4^(1/2) wrongly returned " << e1 << endl; + return 1; + } + + ex e2 = power(numeric(27),numeric(2,3)); + if (e2 != 9) { + clog << "27^(2/3) wrongly returned " << e2 << endl; + return 1; + } + + ex e3 = power(numeric(5),numeric(1,2)); + if (!(is_ex_exactly_of_type(e3,power) && + e3.op(0).is_equal(numeric(5)) && + e3.op(1).is_equal(numeric(1,2)))) { + clog << "5^(1/2) wrongly returned " << e3 << endl; + return 1; + } + + ex e4 = power(numeric(5),evalf(numeric(1,2))); + if (!(is_ex_exactly_of_type(e4,numeric))) { + clog << "5^(0.5) wrongly returned " << e4 << endl; + return 1; + } + + ex e5 = power(evalf(numeric(5)),numeric(1,2)); + if (!(is_ex_exactly_of_type(e5,numeric))) { + clog << "5.0^(1/2) wrongly returned " << e5 << endl; + return 1; + } + + return 0; } static unsigned exam_powerlaws4(void) { - // test for mul::eval() - - symbol a("a"); - symbol b("b"); - symbol c("c"); - - ex f1 = power(a*b,ex(1)/ex(2)); - ex f2 = power(a*b,ex(3)/ex(2)); - ex f3 = c; - - exvector v; - v.push_back(f1); - v.push_back(f2); - v.push_back(f3); - ex e1 = mul(v); - if (e1!=a*a*b*b*c) { - clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl; - return 1; - } - - return 0; + // test for mul::eval() + + symbol a("a"); + symbol b("b"); + symbol c("c"); + + ex f1 = power(a*b,ex(1)/ex(2)); + ex f2 = power(a*b,ex(3)/ex(2)); + ex f3 = c; + + exvector v; + v.push_back(f1); + v.push_back(f2); + v.push_back(f3); + ex e1 = mul(v); + if (e1!=a*a*b*b*c) { + clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl; + return 1; + } + + return 0; } static unsigned exam_powerlaws5(void) { - // cabinet of slightly pathological cases - - symbol a("a"); - - ex e1 = pow(1,a); - if (e1 != 1) { - clog << "1^a wrongly returned " << e1 << endl; - return 1; - } - - ex e2 = pow(0,a); - if (!(is_ex_exactly_of_type(e2,power))) { - clog << "0^a was evaluated to " << e2 - << " though nothing is known about a." << endl; - return 1; - } - - return 0; + // cabinet of slightly pathological cases + + symbol a("a"); + + ex e1 = pow(1,a); + if (e1 != 1) { + clog << "1^a wrongly returned " << e1 << endl; + return 1; + } + + ex e2 = pow(0,a); + if (!(is_ex_exactly_of_type(e2,power))) { + clog << "0^a was evaluated to " << e2 + << " though nothing is known about a." << endl; + return 1; + } + + return 0; } unsigned exam_powerlaws(void) { - unsigned result = 0; - - cout << "examining power laws" << flush; - clog << "----------power laws:" << endl; - - result += exam_powerlaws1(); cout << '.' << flush; - result += exam_powerlaws2(); cout << '.' << flush; - result += exam_powerlaws3(); cout << '.' << flush; - result += exam_powerlaws4(); cout << '.' << flush; - result += exam_powerlaws5(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - - return result; + unsigned result = 0; + + cout << "examining power laws" << flush; + clog << "----------power laws:" << endl; + + result += exam_powerlaws1(); cout << '.' << flush; + result += exam_powerlaws2(); cout << '.' << flush; + result += exam_powerlaws3(); cout << '.' << flush; + result += exam_powerlaws4(); cout << '.' << flush; + result += exam_powerlaws5(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + + return result; } diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp index a6ff7047..133f01b0 100644 --- a/check/exam_pseries.cpp +++ b/check/exam_pseries.cpp @@ -1,4 +1,4 @@ -/** @file exam_pseries.cpp +/** @File exam_pseries.cpp * * Series expansion test (Laurent and Taylor series). */ @@ -26,222 +26,222 @@ static symbol x("x"); static unsigned check_series(const ex &e, const ex &point, const ex &d, int order = 8) { - ex es = e.series(x==point, order); - ex ep = ex_to_pseries(es).convert_to_poly(); - if (!(ep - d).is_zero()) { - clog << "series expansion of " << e << " at " << point - << " erroneously returned " << ep << " (instead of " << d - << ")" << endl; - (ep-d).printtree(clog); - return 1; - } - return 0; + ex es = e.series(x==point, order); + ex ep = ex_to_pseries(es).convert_to_poly(); + if (!(ep - d).is_zero()) { + clog << "series expansion of " << e << " at " << point + << " erroneously returned " << ep << " (instead of " << d + << ")" << endl; + (ep-d).printtree(clog); + return 1; + } + return 0; } // Series expansion static unsigned exam_series1(void) { - unsigned result = 0; - ex e, d; - - e = sin(x); - d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, 0, d); - - e = cos(x); - d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8)); - result += check_series(e, 0, d); - - e = exp(x); - d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, 0, d); - - e = pow(1 - x, -1); - d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8)); - result += check_series(e, 0, d); - - e = x + pow(x, -1); - d = x + pow(x, -1); - result += check_series(e, 0, d); - - e = x + pow(x, -1); - d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8)); - result += check_series(e, 1, d); - - e = pow(x + pow(x, 3), -1); - d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7)); - result += check_series(e, 0, d); - - e = pow(pow(x, 2) + pow(x, 4), -1); - d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6)); - result += check_series(e, 0, d); - - e = pow(sin(x), -2); - d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5)); - result += check_series(e, 0, d); - - e = sin(x) / cos(x); - d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8)); - result += check_series(e, 0, d); - - e = cos(x) / sin(x); - d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6)); - result += check_series(e, 0, d); - - e = pow(numeric(2), x); - ex t = log(2) * x; - d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, 0, d.expand()); - - e = pow(Pi, x); - t = log(Pi) * x; - d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, 0, d.expand()); - - return result; + unsigned result = 0; + ex e, d; + + e = sin(x); + d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8)); + result += check_series(e, 0, d); + + e = cos(x); + d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8)); + result += check_series(e, 0, d); + + e = exp(x); + d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8)); + result += check_series(e, 0, d); + + e = pow(1 - x, -1); + d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8)); + result += check_series(e, 0, d); + + e = x + pow(x, -1); + d = x + pow(x, -1); + result += check_series(e, 0, d); + + e = x + pow(x, -1); + d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8)); + result += check_series(e, 1, d); + + e = pow(x + pow(x, 3), -1); + d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7)); + result += check_series(e, 0, d); + + e = pow(pow(x, 2) + pow(x, 4), -1); + d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6)); + result += check_series(e, 0, d); + + e = pow(sin(x), -2); + d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5)); + result += check_series(e, 0, d); + + e = sin(x) / cos(x); + d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8)); + result += check_series(e, 0, d); + + e = cos(x) / sin(x); + d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6)); + result += check_series(e, 0, d); + + e = pow(numeric(2), x); + ex t = log(2) * x; + d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8)); + result += check_series(e, 0, d.expand()); + + e = pow(Pi, x); + t = log(Pi) * x; + d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8)); + result += check_series(e, 0, d.expand()); + + return result; } // Series addition static unsigned exam_series2(void) { - unsigned result = 0; - ex e, d; - - e = pow(sin(x), -1).series(x==0, 8) + pow(sin(-x), -1).series(x==0, 12); - d = Order(pow(x, 6)); - result += check_series(e, 0, d); - - return result; + unsigned result = 0; + ex e, d; + + e = pow(sin(x), -1).series(x==0, 8) + pow(sin(-x), -1).series(x==0, 12); + d = Order(pow(x, 6)); + result += check_series(e, 0, d); + + return result; } // Series multiplication static unsigned exam_series3(void) { - unsigned result = 0; - ex e, d; - - e = sin(x).series(x==0, 8) * pow(sin(x), -1).series(x==0, 12); - d = 1 + Order(pow(x, 7)); - result += check_series(e, 0, d); - - return result; + unsigned result = 0; + ex e, d; + + e = sin(x).series(x==0, 8) * pow(sin(x), -1).series(x==0, 12); + d = 1 + Order(pow(x, 7)); + result += check_series(e, 0, d); + + return result; } // Order term handling static unsigned exam_series4(void) { - unsigned result = 0; - ex e, d; + unsigned result = 0; + ex e, d; - e = 1 + x + pow(x, 2) + pow(x, 3); - d = Order(1); - result += check_series(e, 0, d, 0); - d = 1 + Order(x); - result += check_series(e, 0, d, 1); - d = 1 + x + Order(pow(x, 2)); - result += check_series(e, 0, d, 2); - d = 1 + x + pow(x, 2) + Order(pow(x, 3)); - result += check_series(e, 0, d, 3); - d = 1 + x + pow(x, 2) + pow(x, 3); - result += check_series(e, 0, d, 4); - return result; + e = 1 + x + pow(x, 2) + pow(x, 3); + d = Order(1); + result += check_series(e, 0, d, 0); + d = 1 + Order(x); + result += check_series(e, 0, d, 1); + d = 1 + x + Order(pow(x, 2)); + result += check_series(e, 0, d, 2); + d = 1 + x + pow(x, 2) + Order(pow(x, 3)); + result += check_series(e, 0, d, 3); + d = 1 + x + pow(x, 2) + pow(x, 3); + result += check_series(e, 0, d, 4); + return result; } // Series expansion of tgamma(-1) static unsigned exam_series5(void) { - ex e = tgamma(2*x); - ex d = pow(x+1,-1)*numeric(1,4) + - pow(x+1,0)*(numeric(3,4) - - numeric(1,2)*Euler) + - pow(x+1,1)*(numeric(7,4) - - numeric(3,2)*Euler + - numeric(1,2)*pow(Euler,2) + - numeric(1,12)*pow(Pi,2)) + - pow(x+1,2)*(numeric(15,4) - - numeric(7,2)*Euler - - numeric(1,3)*pow(Euler,3) + - numeric(1,4)*pow(Pi,2) + - numeric(3,2)*pow(Euler,2) - - numeric(1,6)*pow(Pi,2)*Euler - - numeric(2,3)*zeta(3)) + - pow(x+1,3)*(numeric(31,4) - pow(Euler,3) - - numeric(15,2)*Euler + - numeric(1,6)*pow(Euler,4) + - numeric(7,2)*pow(Euler,2) + - numeric(7,12)*pow(Pi,2) - - numeric(1,2)*pow(Pi,2)*Euler - - numeric(2)*zeta(3) + - numeric(1,6)*pow(Euler,2)*pow(Pi,2) + - numeric(1,40)*pow(Pi,4) + - numeric(4,3)*zeta(3)*Euler) + - Order(pow(x+1,4)); - return check_series(e, -1, d, 4); + ex e = tgamma(2*x); + ex d = pow(x+1,-1)*numeric(1,4) + + pow(x+1,0)*(numeric(3,4) - + numeric(1,2)*Euler) + + pow(x+1,1)*(numeric(7,4) - + numeric(3,2)*Euler + + numeric(1,2)*pow(Euler,2) + + numeric(1,12)*pow(Pi,2)) + + pow(x+1,2)*(numeric(15,4) - + numeric(7,2)*Euler - + numeric(1,3)*pow(Euler,3) + + numeric(1,4)*pow(Pi,2) + + numeric(3,2)*pow(Euler,2) - + numeric(1,6)*pow(Pi,2)*Euler - + numeric(2,3)*zeta(3)) + + pow(x+1,3)*(numeric(31,4) - pow(Euler,3) - + numeric(15,2)*Euler + + numeric(1,6)*pow(Euler,4) + + numeric(7,2)*pow(Euler,2) + + numeric(7,12)*pow(Pi,2) - + numeric(1,2)*pow(Pi,2)*Euler - + numeric(2)*zeta(3) + + numeric(1,6)*pow(Euler,2)*pow(Pi,2) + + numeric(1,40)*pow(Pi,4) + + numeric(4,3)*zeta(3)*Euler) + + Order(pow(x+1,4)); + return check_series(e, -1, d, 4); } - + // Series expansion of tan(x==Pi/2) static unsigned exam_series6(void) { - ex e = tan(x*Pi/2); - ex d = pow(x-1,-1)/Pi*(-2) + pow(x-1,1)*Pi/6 + pow(x-1,3)*pow(Pi,3)/360 - +pow(x-1,5)*pow(Pi,5)/15120 + pow(x-1,7)*pow(Pi,7)/604800 - +Order(pow(x-1,8)); - return check_series(e,1,d,8); + ex e = tan(x*Pi/2); + ex d = pow(x-1,-1)/Pi*(-2) + pow(x-1,1)*Pi/6 + pow(x-1,3)*pow(Pi,3)/360 + +pow(x-1,5)*pow(Pi,5)/15120 + pow(x-1,7)*pow(Pi,7)/604800 + +Order(pow(x-1,8)); + return check_series(e,1,d,8); } // Series expansion of log(sin(x==0)) static unsigned exam_series7(void) { - ex e = log(sin(x)); - ex d = log(x) - pow(x,2)/6 - pow(x,4)/180 - pow(x,6)/2835 - +Order(pow(x,8)); - return check_series(e,0,d,8); + ex e = log(sin(x)); + ex d = log(x) - pow(x,2)/6 - pow(x,4)/180 - pow(x,6)/2835 + +Order(pow(x,8)); + return check_series(e,0,d,8); } // Series expansion of Li2(sin(x==0)) static unsigned exam_series8(void) { - ex e = Li2(sin(x)); - ex d = x + pow(x,2)/4 - pow(x,3)/18 - pow(x,4)/48 - - 13*pow(x,5)/1800 - pow(x,6)/360 - 23*pow(x,7)/21168 - + Order(pow(x,8)); - return check_series(e,0,d,8); + ex e = Li2(sin(x)); + ex d = x + pow(x,2)/4 - pow(x,3)/18 - pow(x,4)/48 + - 13*pow(x,5)/1800 - pow(x,6)/360 - 23*pow(x,7)/21168 + + Order(pow(x,8)); + return check_series(e,0,d,8); } // Series expansion of Li2((x==2)^2), caring about branch-cut static unsigned exam_series9(void) { - ex e = Li2(pow(x,2)); - ex d = Li2(4) + (-log(3) + I*Pi*csgn(I-I*pow(x,2))) * (x-2) - + (numeric(-2,3) + log(3)/4 - I*Pi/4*csgn(I-I*pow(x,2))) * pow(x-2,2) - + (numeric(11,27) - log(3)/12 + I*Pi/12*csgn(I-I*pow(x,2))) * pow(x-2,3) - + (numeric(-155,648) + log(3)/32 - I*Pi/32*csgn(I-I*pow(x,2))) * pow(x-2,4) - + Order(pow(x-2,5)); - return check_series(e,2,d,5); + ex e = Li2(pow(x,2)); + ex d = Li2(4) + (-log(3) + I*Pi*csgn(I-I*pow(x,2))) * (x-2) + + (numeric(-2,3) + log(3)/4 - I*Pi/4*csgn(I-I*pow(x,2))) * pow(x-2,2) + + (numeric(11,27) - log(3)/12 + I*Pi/12*csgn(I-I*pow(x,2))) * pow(x-2,3) + + (numeric(-155,648) + log(3)/32 - I*Pi/32*csgn(I-I*pow(x,2))) * pow(x-2,4) + + Order(pow(x-2,5)); + return check_series(e,2,d,5); } unsigned exam_pseries(void) { - unsigned result = 0; - - cout << "examining series expansion" << flush; - clog << "----------series expansion:" << endl; - - result += exam_series1(); cout << '.' << flush; - result += exam_series2(); cout << '.' << flush; - result += exam_series3(); cout << '.' << flush; - result += exam_series4(); cout << '.' << flush; - result += exam_series5(); cout << '.' << flush; - result += exam_series6(); cout << '.' << flush; - result += exam_series7(); cout << '.' << flush; - result += exam_series8(); cout << '.' << flush; - result += exam_series9(); cout << '.' << flush; - - if (!result) { - cout << " passed " << endl; - clog << "(no output)" << endl; - } else { - cout << " failed " << endl; - } - return result; + unsigned result = 0; + + cout << "examining series expansion" << flush; + clog << "----------series expansion:" << endl; + + result += exam_series1(); cout << '.' << flush; + result += exam_series2(); cout << '.' << flush; + result += exam_series3(); cout << '.' << flush; + result += exam_series4(); cout << '.' << flush; + result += exam_series5(); cout << '.' << flush; + result += exam_series6(); cout << '.' << flush; + result += exam_series7(); cout << '.' << flush; + result += exam_series8(); cout << '.' << flush; + result += exam_series9(); cout << '.' << flush; + + if (!result) { + cout << " passed " << endl; + clog << "(no output)" << endl; + } else { + cout << " failed " << endl; + } + return result; } diff --git a/check/exams.cpp b/check/exams.cpp index 5aa327ae..01e8f25b 100644 --- a/check/exams.cpp +++ b/check/exams.cpp @@ -27,102 +27,102 @@ int main() { - unsigned result = 0; - - try { - result += exam_paranoia(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_numeric(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_powerlaws(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_inifcns(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_differentiation(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_polygcd(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_normalization(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_pseries(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_matrices(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_lsolve(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_noncommut(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += exam_misc(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - if (result) { - cout << "Error: something went wrong. "; - if (result == 1) { - cout << "(one failure)" << endl; - } else { - cout << "(" << result << " individual failures)" << endl; - } - cout << "please check exams.out against exams.ref for more details." - << endl << "happy debugging!" << endl; - } - - return result; + unsigned result = 0; + + try { + result += exam_paranoia(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_numeric(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_powerlaws(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_inifcns(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_differentiation(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_polygcd(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_normalization(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_pseries(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_matrices(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_lsolve(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_noncommut(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += exam_misc(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + if (result) { + cout << "Error: something went wrong. "; + if (result == 1) { + cout << "(one failure)" << endl; + } else { + cout << "(" << result << " individual failures)" << endl; + } + cout << "please check exams.out against exams.ref for more details." + << endl << "happy debugging!" << endl; + } + + return result; } diff --git a/check/genex.cpp b/check/genex.cpp index 18c2ded4..b61e8e83 100644 --- a/check/genex.cpp +++ b/check/genex.cpp @@ -35,12 +35,12 @@ using namespace GiNaC; const ex dense_univariate_poly(const symbol & x, unsigned degree) { - ex unipoly; - - for (unsigned i=0; i<=degree; ++i) - unipoly += numeric((rand()-RAND_MAX/2))*pow(x,i); - - return unipoly; + ex unipoly; + + for (unsigned i=0; i<=degree; ++i) + unipoly += numeric((rand()-RAND_MAX/2))*pow(x,i); + + return unipoly; } /* Create a dense bivariate random polynomial in x1 and x2. @@ -49,110 +49,110 @@ dense_univariate_poly(const symbol & x, unsigned degree) const ex dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree) { - ex bipoly; - - for (unsigned i1=0; i1<=degree; ++i1) - for (unsigned i2=0; i2<=degree-i1; ++i2) - bipoly += numeric((rand()-RAND_MAX/2))*pow(x1,i1)*pow(x2,i2); - - return bipoly; + ex bipoly; + + for (unsigned i1=0; i1<=degree; ++i1) + for (unsigned i2=0; i2<=degree-i1; ++i2) + bipoly += numeric((rand()-RAND_MAX/2))*pow(x1,i1)*pow(x2,i2); + + return bipoly; } /* Chose a randum symbol or number from the argument list. */ const ex random_symbol(const symbol & x, - const symbol & y, - const symbol & z, - bool rational = true, - bool complex = false) + const symbol & y, + const symbol & z, + bool rational = true, + bool complex = false) { - ex e; - switch (abs(rand()) % 4) { - case 0: - e = x; - break; - case 1: - e = y; - break; - case 2: - e = z; - break; - case 3: { - int c1; - do { c1 = rand()%20 - 10; } while (!c1); - int c2; - do { c2 = rand()%20 - 10; } while (!c2); - if (!rational) - c2 = 1; - e = numeric(c1, c2); - if (complex && !(rand()%5)) - e = e*I; - break; - } - } - return e; + ex e; + switch (abs(rand()) % 4) { + case 0: + e = x; + break; + case 1: + e = y; + break; + case 2: + e = z; + break; + case 3: { + int c1; + do { c1 = rand()%20 - 10; } while (!c1); + int c2; + do { c2 = rand()%20 - 10; } while (!c2); + if (!rational) + c2 = 1; + e = numeric(c1, c2); + if (complex && !(rand()%5)) + e = e*I; + break; + } + } + return e; } /* Create a sparse random tree in three symbols. */ const ex sparse_tree(const symbol & x, - const symbol & y, - const symbol & z, - int level, - bool trig = false, // true includes trigonomatric functions - bool rational = true, // false excludes coefficients in Q - bool complex = false) // true includes complex numbers + const symbol & y, + const symbol & z, + int level, + bool trig = false, // true includes trigonomatric functions + bool rational = true, // false excludes coefficients in Q + bool complex = false) // true includes complex numbers { - if (level == 0) - return random_symbol(x,y,z,rational,complex); - switch (abs(rand()) % 10) { - case 0: - case 1: - case 2: - case 3: - return add(sparse_tree(x,y,z,level-1, trig, rational), - sparse_tree(x,y,z,level-1, trig, rational)); - case 4: - case 5: - case 6: - return mul(sparse_tree(x,y,z,level-1, trig, rational), - sparse_tree(x,y,z,level-1, trig, rational)); - case 7: - case 8: { - ex powbase; - do { - powbase = sparse_tree(x,y,z,level-1, trig, rational); - } while (powbase.is_zero()); - return pow(powbase, abs(rand() % 4)); - break; - } - case 9: - if (trig) { - switch (abs(rand()) % 4) { - case 0: - return sin(sparse_tree(x,y,z,level-1, trig, rational)); - case 1: - return cos(sparse_tree(x,y,z,level-1, trig, rational)); - case 2: - return exp(sparse_tree(x,y,z,level-1, trig, rational)); - case 3: { - ex logex; - do { - ex logarg; - do { - logarg = sparse_tree(x,y,z,level-1, trig, rational); - } while (logarg.is_zero()); - // Keep the evaluator from accidentally plugging an - // unwanted I in the tree: - if (!complex && logarg.info(info_flags::negative)) - logarg = -logarg; - logex = log(logarg); - } while (logex.is_zero()); - return logex; - break; - } - } - } else - return random_symbol(x,y,z,rational,complex); - } + if (level == 0) + return random_symbol(x,y,z,rational,complex); + switch (abs(rand()) % 10) { + case 0: + case 1: + case 2: + case 3: + return add(sparse_tree(x,y,z,level-1, trig, rational), + sparse_tree(x,y,z,level-1, trig, rational)); + case 4: + case 5: + case 6: + return mul(sparse_tree(x,y,z,level-1, trig, rational), + sparse_tree(x,y,z,level-1, trig, rational)); + case 7: + case 8: { + ex powbase; + do { + powbase = sparse_tree(x,y,z,level-1, trig, rational); + } while (powbase.is_zero()); + return pow(powbase, abs(rand() % 4)); + break; + } + case 9: + if (trig) { + switch (abs(rand()) % 4) { + case 0: + return sin(sparse_tree(x,y,z,level-1, trig, rational)); + case 1: + return cos(sparse_tree(x,y,z,level-1, trig, rational)); + case 2: + return exp(sparse_tree(x,y,z,level-1, trig, rational)); + case 3: { + ex logex; + do { + ex logarg; + do { + logarg = sparse_tree(x,y,z,level-1, trig, rational); + } while (logarg.is_zero()); + // Keep the evaluator from accidentally plugging an + // unwanted I in the tree: + if (!complex && logarg.info(info_flags::negative)) + logarg = -logarg; + logex = log(logarg); + } while (logex.is_zero()); + return logex; + break; + } + } + } else + return random_symbol(x,y,z,rational,complex); + } } diff --git a/check/time_dennyfliegner.cpp b/check/time_dennyfliegner.cpp index db9ff80c..e20fdded 100644 --- a/check/time_dennyfliegner.cpp +++ b/check/time_dennyfliegner.cpp @@ -26,73 +26,71 @@ #include "times.h" -#define VECSIZE 200 - static unsigned expand_subs(unsigned size) { - unsigned result = 0; - // create a vector of size symbols named "a0", "a1", ... - vector a; - for (unsigned i=0; i a; + for (unsigned i=0; i sizes; - vector times; - timer breitling; - - sizes.push_back(25); - sizes.push_back(50); - sizes.push_back(100); - sizes.push_back(200); - - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { - breitling.start(); - result += expand_subs(*i); - times.push_back(breitling.read()); - cout << '.' << flush; - } - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - // print the report: - cout << endl << " size: "; - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) - cout << '\t' << (*i); - cout << endl << " time/s:"; - for (vector::iterator i=times.begin(); i!=times.end(); ++i) - cout << '\t' << int(1000*(*i))*0.001; - cout << endl; - - return result; + unsigned result = 0; + + cout << "timing commutative expansion and substitution" << flush; + clog << "-------commutative expansion and substitution:" << endl; + + vector sizes; + vector times; + timer breitling; + + sizes.push_back(25); + sizes.push_back(50); + sizes.push_back(100); + sizes.push_back(200); + + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { + breitling.start(); + result += expand_subs(*i); + times.push_back(breitling.read()); + cout << '.' << flush; + } + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + // print the report: + cout << endl << " size: "; + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) + cout << '\t' << (*i); + cout << endl << " time/s:"; + for (vector::iterator i=times.begin(); i!=times.end(); ++i) + cout << '\t' << int(1000*(*i))*0.001; + cout << endl; + + return result; } diff --git a/check/time_gammaseries.cpp b/check/time_gammaseries.cpp index 0bfbe302..59a08a9f 100644 --- a/check/time_gammaseries.cpp +++ b/check/time_gammaseries.cpp @@ -24,62 +24,62 @@ unsigned tgammaseries(unsigned order) { - unsigned result = 0; - symbol x; - - ex myseries = series(tgamma(x),x==0,order); - // compute the last coefficient numerically: - ex last_coeff = myseries.coeff(x,order-1).evalf(); - // compute a bound for that coefficient using a variation of the leading - // term in Stirling's formula: - ex bound = evalf(exp(ex(-.57721566490153286*(order-1)))/(order-1)); - if (evalf(abs((last_coeff-pow(-1,order))/bound)) > numeric(1)) { - clog << "The " << order-1 - << "th order coefficient in the power series expansion of tgamma(0) was erroneously found to be " - << last_coeff << ", violating a simple estimate." << endl; - ++result; - } - - return result; + unsigned result = 0; + symbol x; + + ex myseries = series(tgamma(x),x==0,order); + // compute the last coefficient numerically: + ex last_coeff = myseries.coeff(x,order-1).evalf(); + // compute a bound for that coefficient using a variation of the leading + // term in Stirling's formula: + ex bound = evalf(exp(ex(-.57721566490153286*(order-1)))/(order-1)); + if (evalf(abs((last_coeff-pow(-1,order))/bound)) > numeric(1)) { + clog << "The " << order-1 + << "th order coefficient in the power series expansion of tgamma(0) was erroneously found to be " + << last_coeff << ", violating a simple estimate." << endl; + ++result; + } + + return result; } unsigned time_gammaseries(void) { - unsigned result = 0; - - cout << "timing Laurent series expansion of Gamma function" << flush; - clog << "-------Laurent series expansion of Gamma function:" << endl; - - vector sizes; - vector times; - timer omega; - - sizes.push_back(10); - sizes.push_back(15); - sizes.push_back(20); - sizes.push_back(25); - - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { - omega.start(); - result += tgammaseries(*i); - times.push_back(omega.read()); - cout << '.' << flush; - } - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - // print the report: - cout << endl << " order: "; - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) - cout << '\t' << (*i); - cout << endl << " time/s:"; - for (vector::iterator i=times.begin(); i!=times.end(); ++i) - cout << '\t' << int(1000*(*i))*0.001; - cout << endl; - - return result; + unsigned result = 0; + + cout << "timing Laurent series expansion of Gamma function" << flush; + clog << "-------Laurent series expansion of Gamma function:" << endl; + + vector sizes; + vector times; + timer omega; + + sizes.push_back(10); + sizes.push_back(15); + sizes.push_back(20); + sizes.push_back(25); + + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { + omega.start(); + result += tgammaseries(*i); + times.push_back(omega.read()); + cout << '.' << flush; + } + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + // print the report: + cout << endl << " order: "; + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) + cout << '\t' << (*i); + cout << endl << " time/s:"; + for (vector::iterator i=times.begin(); i!=times.end(); ++i) + cout << '\t' << int(1000*(*i))*0.001; + cout << endl; + + return result; } diff --git a/check/time_lw_A.cpp b/check/time_lw_A.cpp index 2fa40c92..ac56eed4 100644 --- a/check/time_lw_A.cpp +++ b/check/time_lw_A.cpp @@ -25,42 +25,42 @@ static unsigned test(void) { - for (int i=1; i<=99; ++i) - factorial(1000+i)/factorial(900+i); - ex rat(factorial(1100)/factorial(1000)); - - if (abs(evalf(rat)-numeric(".13280014101512E303"))>numeric("1.0E289")) { - clog << "1100!/1000! erroneously returned " << rat << endl; - return 1; - } - return 0; + for (int i=1; i<=99; ++i) + factorial(1000+i)/factorial(900+i); + ex rat(factorial(1100)/factorial(1000)); + + if (abs(evalf(rat)-numeric(".13280014101512E303"))>numeric("1.0E289")) { + clog << "1100!/1000! erroneously returned " << rat << endl; + return 1; + } + return 0; } unsigned time_lw_A(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test A (divide factorials)" << flush; - clog << "-------Lewis-Wester test A (divide factorials)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(1000*(time/count))*0.001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test A (divide factorials)" << flush; + clog << "-------Lewis-Wester test A (divide factorials)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(1000*(time/count))*0.001 << 's' << endl; + + return result; } diff --git a/check/time_lw_B.cpp b/check/time_lw_B.cpp index 3431203d..6fd60dae 100644 --- a/check/time_lw_B.cpp +++ b/check/time_lw_B.cpp @@ -25,43 +25,43 @@ static unsigned test(void) { - numeric s; - - for (int i=1; i<=1000; ++i) - s += numeric(i).inverse(); - - if (abs(s.evalf()-numeric("7.4854708605503449"))>numeric("2.0E-16")) { - clog << "sum(1/i,i=1..1000) erroneously returned " << s << endl; - return 1; - } - return 0; + numeric s; + + for (int i=1; i<=1000; ++i) + s += numeric(i).inverse(); + + if (abs(s.evalf()-numeric("7.4854708605503449"))>numeric("2.0E-16")) { + clog << "sum(1/i,i=1..1000) erroneously returned " << s << endl; + return 1; + } + return 0; } unsigned time_lw_B(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test B (sum of rational numbers)" << flush; - clog << "-------Lewis-Wester test B (sum of rational numbers)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(1000*(time/count))*0.001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test B (sum of rational numbers)" << flush; + clog << "-------Lewis-Wester test B (sum of rational numbers)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(1000*(time/count))*0.001 << 's' << endl; + + return result; } diff --git a/check/time_lw_C.cpp b/check/time_lw_C.cpp index 8885d870..ce2921c8 100644 --- a/check/time_lw_C.cpp +++ b/check/time_lw_C.cpp @@ -25,47 +25,47 @@ static unsigned test(void) { - numeric x(13*17*31); - numeric y(13*19*29); - - for (int i=1; i<200; ++i) - gcd(pow(x,300+(i%181)),pow(y,200+(i%183))); - - ex lastgcd = gcd(pow(x,300+(200%181)),pow(y,200+(200%183))); - if (lastgcd != numeric("53174994123961114423610399251974962981084780166115806651505844915220196792416194060680805428433601792982500430324916963290494659936522782673704312949880308677990050199363768068005367578752699785180694630122629259539608472261461289805919741933")) { - clog << "gcd(" << x << "^" << 300+(200%181) << "," - << y << "^" << 200+(200%183) << ") erroneously returned " - << lastgcd << endl; - return 1; - } - return 0; + numeric x(13*17*31); + numeric y(13*19*29); + + for (int i=1; i<200; ++i) + gcd(pow(x,300+(i%181)),pow(y,200+(i%183))); + + ex lastgcd = gcd(pow(x,300+(200%181)),pow(y,200+(200%183))); + if (lastgcd != numeric("53174994123961114423610399251974962981084780166115806651505844915220196792416194060680805428433601792982500430324916963290494659936522782673704312949880308677990050199363768068005367578752699785180694630122629259539608472261461289805919741933")) { + clog << "gcd(" << x << "^" << 300+(200%181) << "," + << y << "^" << 200+(200%183) << ") erroneously returned " + << lastgcd << endl; + return 1; + } + return 0; } unsigned time_lw_C(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test C (gcd of big integers)" << flush; - clog << "-------Lewis-Wester test C (gcd of big integers)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(1000*(time/count))*0.001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test C (gcd of big integers)" << flush; + clog << "-------Lewis-Wester test C (gcd of big integers)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(1000*(time/count))*0.001 << 's' << endl; + + return result; } diff --git a/check/time_lw_D.cpp b/check/time_lw_D.cpp index 1f63e58f..f226f51a 100644 --- a/check/time_lw_D.cpp +++ b/check/time_lw_D.cpp @@ -25,45 +25,45 @@ static unsigned test(void) { - ex s; - symbol y("y"); - symbol t("t"); - - for (int i=1; i<=10; ++i) - s += i*y*pow(t,i)/pow(y + i*t,i); - - if (s.nops()!=10) { - clog << "something very strange happened" << endl; - return 1; - } - return 0; + ex s; + symbol y("y"); + symbol t("t"); + + for (int i=1; i<=10; ++i) + s += i*y*pow(t,i)/pow(y + i*t,i); + + if (s.nops()!=10) { + clog << "something very strange happened" << endl; + return 1; + } + return 0; } unsigned time_lw_D(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test D (sum of rational fcns)" << flush; - clog << "-------Lewis-Wester test D (sum of rational fcns)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(100000*(time/count))*0.00001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test D (sum of rational fcns)" << flush; + clog << "-------Lewis-Wester test D (sum of rational fcns)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(100000*(time/count))*0.00001 << 's' << endl; + + return result; } diff --git a/check/time_lw_E.cpp b/check/time_lw_E.cpp index 5dc004d2..95a4fd4a 100644 --- a/check/time_lw_E.cpp +++ b/check/time_lw_E.cpp @@ -25,45 +25,45 @@ static unsigned test(void) { - ex s; - symbol y("y"); - symbol t("t"); - - for (int i=1; i<=10; ++i) - s += i*y*pow(t,i)/pow(y + abs(5-i)*t,i); - - if (s.nops()!=10) { - clog << "something very strange happened" << endl; - return 1; - } - return 0; + ex s; + symbol y("y"); + symbol t("t"); + + for (int i=1; i<=10; ++i) + s += i*y*pow(t,i)/pow(y + abs(5-i)*t,i); + + if (s.nops()!=10) { + clog << "something very strange happened" << endl; + return 1; + } + return 0; } unsigned time_lw_E(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test E (sum of rational fcns)" << flush; - clog << "-------Lewis-Wester test E (sum of rational fcns)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(100000*(time/count))*0.00001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test E (sum of rational fcns)" << flush; + clog << "-------Lewis-Wester test E (sum of rational fcns)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(100000*(time/count))*0.00001 << 's' << endl; + + return result; } diff --git a/check/time_lw_F.cpp b/check/time_lw_F.cpp index ce8920e5..de8aa92b 100644 --- a/check/time_lw_F.cpp +++ b/check/time_lw_F.cpp @@ -25,44 +25,44 @@ static unsigned test(void) { - symbol x("x"); - symbol y("y"); + symbol x("x"); + symbol y("y"); - ex p = expand(pow(pow(x,2)-3*x*y+pow(y,2),4)*pow(3*x-7*y+2,5)); - ex q = expand(pow(pow(x,2)-3*x*y+pow(y,2),3)*pow(3*x-7*y-2,6)); - ex result = gcd(p,q); - if (result!=expand(pow(pow(x,2)-3*x*y+pow(y,2),3))) { - clog << "gcd((x^2-3*x*y+y^2)^4*(3*x-7*y+2)^5),(x^2-3*x*y+y^2)^3*(3*x-7*y-2)^6)) erroneously returned " << result << endl; - return 1; - } - return 0; + ex p = expand(pow(pow(x,2)-3*x*y+pow(y,2),4)*pow(3*x-7*y+2,5)); + ex q = expand(pow(pow(x,2)-3*x*y+pow(y,2),3)*pow(3*x-7*y-2,6)); + ex result = gcd(p,q); + if (result!=expand(pow(pow(x,2)-3*x*y+pow(y,2),3))) { + clog << "gcd(expand((x^2-3*x*y+y^2)^4*(3*x-7*y+2)^5),expand((x^2-3*x*y+y^2)^3*(3*x-7*y-2)^6)) erroneously returned " << result << endl; + return 1; + } + return 0; } unsigned time_lw_F(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test F (gcd of 2-var polys)" << flush; - clog << "-------Lewis-Wester test F (gcd of 2-var polys)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(1000*(time/count))*0.001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test F (gcd of 2-var polys)" << flush; + clog << "-------Lewis-Wester test F (gcd of 2-var polys)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(1000*(time/count))*0.001 << 's' << endl; + + return result; } diff --git a/check/time_lw_G.cpp b/check/time_lw_G.cpp index 75dd4261..1bae9c88 100644 --- a/check/time_lw_G.cpp +++ b/check/time_lw_G.cpp @@ -25,47 +25,47 @@ static unsigned test(void) { - symbol x("x"); - symbol y("y"); - symbol z("z"); + symbol x("x"); + symbol y("y"); + symbol z("z"); - ex p = expand(pow(7*y*pow(x*z,2)-3*x*y*z+11*(x+1)*pow(y,2)+5*z+1,4) - *pow(3*x-7*y+2*z-3,5)); - ex q = expand(pow(7*y*pow(x*z,2)-3*x*y*z+11*(x+1)*pow(y,2)+5*z+1,3) - *pow(3*x-7*y+2*z+3,6)); - ex result = gcd(p,q); - if (result.expand()!=expand(pow(7*y*pow(x*z,2)-3*x*y*z+11*(x+1)*pow(y,2)+5*z+1,3))) { - clog << "gcd(expand((7*y*x^2*z^2-3*x*y*z+11*(x+1)*y^2+5*z+1)^4*(3*x-7*y+2*z-3)^5),expand((7*y*x^2*z^2-3*x*y*z+11*(x+1)*y^2+5*z+1)^3*(3*x-7*y+2*z+3)^6)) erroneously returned " << result << endl; - return 1; - } - return 0; + ex p = expand(pow(7*y*pow(x*z,2)-3*x*y*z+11*(x+1)*pow(y,2)+5*z+1,4) + *pow(3*x-7*y+2*z-3,5)); + ex q = expand(pow(7*y*pow(x*z,2)-3*x*y*z+11*(x+1)*pow(y,2)+5*z+1,3) + *pow(3*x-7*y+2*z+3,6)); + ex result = gcd(p,q); + if (result.expand()!=expand(pow(7*y*pow(x*z,2)-3*x*y*z+11*(x+1)*pow(y,2)+5*z+1,3))) { + clog << "gcd(expand((7*y*x^2*z^2-3*x*y*z+11*(x+1)*y^2+5*z+1)^4*(3*x-7*y+2*z-3)^5),expand((7*y*x^2*z^2-3*x*y*z+11*(x+1)*y^2+5*z+1)^3*(3*x-7*y+2*z+3)^6)) erroneously returned " << result << endl; + return 1; + } + return 0; } unsigned time_lw_G(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test G (gcd of 3-var polys)" << flush; - clog << "-------Lewis-Wester test G (gcd of 3-var polys)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(1000*(time/count))*0.001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test G (gcd of 3-var polys)" << flush; + clog << "-------Lewis-Wester test G (gcd of 3-var polys)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(1000*(time/count))*0.001 << 's' << endl; + + return result; } diff --git a/check/time_lw_H.cpp b/check/time_lw_H.cpp index a7c5c852..8ea4b47f 100644 --- a/check/time_lw_H.cpp +++ b/check/time_lw_H.cpp @@ -25,45 +25,45 @@ static unsigned test(void) { - matrix h80(80,80); - - for (unsigned r=0; r<80; ++r) - for (unsigned c=0; c<80; ++c) - h80.set(r,c,numeric(1,r+c+1)); - ex det = h80.determinant(); - - if (abs(det.evalf()-numeric(".10097939769690107E-3789"))>numeric("1.E-3800")) { - clog << "determinant of 80x80 erroneously returned " << det << endl; - return 1; - } - return 0; + matrix h80(80,80); + + for (unsigned r=0; r<80; ++r) + for (unsigned c=0; c<80; ++c) + h80.set(r,c,numeric(1,r+c+1)); + ex det = h80.determinant(); + + if (abs(det.evalf()-numeric(".10097939769690107E-3789"))>numeric("1.E-3800")) { + clog << "determinant of 80x80 erroneously returned " << det << endl; + return 1; + } + return 0; } unsigned time_lw_H(void) { - unsigned result = 0; - unsigned count = 0; - timer rolex; - double time = .0; - - cout << "timing Lewis-Wester test H (det of 80x80 Hilbert)" << flush; - clog << "-------Lewis-Wester test H (det of 80x80 Hilbert)" << endl; - - rolex.start(); - // correct for very small times: - do { - result = test(); - ++count; - } while ((time=rolex.read())<0.1 && !result); - cout << '.' << flush; - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - cout << int(1000*(time/count))*0.001 << 's' << endl; - - return result; + unsigned result = 0; + unsigned count = 0; + timer rolex; + double time = .0; + + cout << "timing Lewis-Wester test H (det of 80x80 Hilbert)" << flush; + clog << "-------Lewis-Wester test H (det of 80x80 Hilbert)" << endl; + + rolex.start(); + // correct for very small times: + do { + result = test(); + ++count; + } while ((time=rolex.read())<0.1 && !result); + cout << '.' << flush; + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + cout << int(1000*(time/count))*0.001 << 's' << endl; + + return result; } diff --git a/check/time_lw_IJKL.cpp b/check/time_lw_IJKL.cpp index 62f69b66..41c07702 100644 --- a/check/time_lw_IJKL.cpp +++ b/check/time_lw_IJKL.cpp @@ -25,68 +25,68 @@ static unsigned test(unsigned n) { - unsigned result = 0; - timer cartier; - char name = (n==40?'I':(n==70?'K':'?')); - - cout << "timing Lewis-Wester test " << name - << " (invert rank " << n << " Hilbert)" << flush; - clog << "-------Lewis-Wester test " << name - << " (invert rank " << n << " Hilbert)" << endl; - - // Create a rank n Hilbert matrix: - matrix H(n,n); - for (unsigned r=0; r sizes; - vector times; - timer longines; - - sizes.push_back(5); - sizes.push_back(6); - sizes.push_back(7); - sizes.push_back(8); - - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { - int count = 1; - longines.start(); - result += toeplitz_det(*i); - // correct for very small times: - while (longines.read()<0.1) { - toeplitz_det(*i); - ++count; - } - times.push_back(longines.read()/count); - cout << '.' << flush; - } - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - // print the report: - cout << endl << " dim: "; - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) - cout << '\t' << *i << 'x' << *i; - cout << endl << " time/s:"; - for (vector::iterator i=times.begin(); i!=times.end(); ++i) - cout << '\t' << int(1000*(*i))*0.001; - cout << endl; - - return result; + unsigned result = 0; + + cout << "timing determinant of polyvariate symbolic Toeplitz matrices" << flush; + clog << "-------determinant of polyvariate symbolic Toeplitz matrices:" << endl; + + vector sizes; + vector times; + timer longines; + + sizes.push_back(5); + sizes.push_back(6); + sizes.push_back(7); + sizes.push_back(8); + + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { + int count = 1; + longines.start(); + result += toeplitz_det(*i); + // correct for very small times: + while (longines.read()<0.1) { + toeplitz_det(*i); + ++count; + } + times.push_back(longines.read()/count); + cout << '.' << flush; + } + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + // print the report: + cout << endl << " dim: "; + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) + cout << '\t' << *i << 'x' << *i; + cout << endl << " time/s:"; + for (vector::iterator i=times.begin(); i!=times.end(); ++i) + cout << '\t' << int(1000*(*i))*0.001; + cout << endl; + + return result; } diff --git a/check/time_vandermonde.cpp b/check/time_vandermonde.cpp index cfef6ea2..69f8934a 100644 --- a/check/time_vandermonde.cpp +++ b/check/time_vandermonde.cpp @@ -28,77 +28,77 @@ static unsigned vandermonde_det(unsigned size) { - unsigned result = 0; - symbol a("a"); - - // construct Vandermonde matrix: - matrix M(size,size); - for (unsigned ro=0; ro sizes; - vector times; - timer swatch; - - sizes.push_back(4); - sizes.push_back(6); - sizes.push_back(8); - sizes.push_back(10); - - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { - int count = 1; - swatch.start(); - result += vandermonde_det(*i); - // correct for very small times: - while (swatch.read()<0.02) { - vandermonde_det(*i); - ++count; - } - times.push_back(swatch.read()/count); - cout << '.' << flush; - } - - if (!result) { - cout << " passed "; - clog << "(no output)" << endl; - } else { - cout << " failed "; - } - // print the report: - cout << endl << " dim: "; - for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) - cout << '\t' << *i << 'x' << *i; - cout << endl << " time/s:"; - for (vector::iterator i=times.begin(); i!=times.end(); ++i) - cout << '\t' << int(1000*(*i))*0.001; - cout << endl; - - return result; + unsigned result = 0; + + cout << "timing determinant of univariate symbolic Vandermonde matrices" << flush; + clog << "-------determinant of univariate symbolic Vandermonde matrices:" << endl; + + vector sizes; + vector times; + timer swatch; + + sizes.push_back(4); + sizes.push_back(6); + sizes.push_back(8); + sizes.push_back(10); + + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { + int count = 1; + swatch.start(); + result += vandermonde_det(*i); + // correct for very small times: + while (swatch.read()<0.02) { + vandermonde_det(*i); + ++count; + } + times.push_back(swatch.read()/count); + cout << '.' << flush; + } + + if (!result) { + cout << " passed "; + clog << "(no output)" << endl; + } else { + cout << " failed "; + } + // print the report: + cout << endl << " dim: "; + for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) + cout << '\t' << *i << 'x' << *i; + cout << endl << " time/s:"; + for (vector::iterator i=times.begin(); i!=times.end(); ++i) + cout << '\t' << int(1000*(*i))*0.001; + cout << endl; + + return result; } diff --git a/check/timer.cpp b/check/timer.cpp index c4c21a30..ec2e557c 100644 --- a/check/timer.cpp +++ b/check/timer.cpp @@ -24,43 +24,43 @@ timer::timer(void) : on(false) { - getrusage(RUSAGE_SELF, &used1); - getrusage(RUSAGE_SELF, &used2); + getrusage(RUSAGE_SELF, &used1); + getrusage(RUSAGE_SELF, &used2); } void timer::start(void) { - on = true; - getrusage(RUSAGE_SELF, &used1); - getrusage(RUSAGE_SELF, &used2); + on = true; + getrusage(RUSAGE_SELF, &used1); + getrusage(RUSAGE_SELF, &used2); } void timer::stop(void) { - on = false; - getrusage(RUSAGE_SELF, &used2); + on = false; + getrusage(RUSAGE_SELF, &used2); } void timer::reset(void) { - getrusage(RUSAGE_SELF, &used1); - getrusage(RUSAGE_SELF, &used2); + getrusage(RUSAGE_SELF, &used1); + getrusage(RUSAGE_SELF, &used2); } double timer::read(void) { - double elapsed; - if (this->running()) - getrusage(RUSAGE_SELF, &used2); - elapsed = ((used2.ru_utime.tv_sec - used1.ru_utime.tv_sec) + - (used2.ru_stime.tv_sec - used1.ru_stime.tv_sec) + - (used2.ru_utime.tv_usec - used1.ru_utime.tv_usec) / 1e6 + - (used2.ru_stime.tv_usec - used1.ru_stime.tv_usec) / 1e6); - // round to 10ms for safety: - return 0.01*int(elapsed*100+0.5); + double elapsed; + if (this->running()) + getrusage(RUSAGE_SELF, &used2); + elapsed = ((used2.ru_utime.tv_sec - used1.ru_utime.tv_sec) + + (used2.ru_stime.tv_sec - used1.ru_stime.tv_sec) + + (used2.ru_utime.tv_usec - used1.ru_utime.tv_usec) / 1e6 + + (used2.ru_stime.tv_usec - used1.ru_stime.tv_usec) / 1e6); + // round to 10ms for safety: + return 0.01*int(elapsed*100+0.5); } bool timer::running(void) { - return on; + return on; } diff --git a/check/times.cpp b/check/times.cpp index 9542c84e..19031dd9 100644 --- a/check/times.cpp +++ b/check/times.cpp @@ -27,151 +27,151 @@ int main() { - unsigned result = 0; - - try { - result += time_dennyfliegner(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_gammaseries(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_vandermonde(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_toeplitz(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_A(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_B(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_C(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_D(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_E(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_F(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_G(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_H(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_IJKL(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_M1(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_O(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_P(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_Pprime(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_Q(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - try { - result += time_lw_Qprime(); - } catch (const exception &e) { - cout << "Error: caught exception " << e.what() << endl; - ++result; - } - - if (result) { - cout << "Error: something went wrong. "; - if (result == 1) { - cout << "(one failure)" << endl; - } else { - cout << "(" << result << " individual failures)" << endl; - } - cout << "please check times.out against times.ref for more details." - << endl << "happy debugging!" << endl; - } - - return result; + unsigned result = 0; + + try { + result += time_dennyfliegner(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_gammaseries(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_vandermonde(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_toeplitz(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_A(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_B(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_C(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_D(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_E(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_F(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_G(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_H(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_IJKL(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_M1(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_O(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_P(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_Pprime(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_Q(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + try { + result += time_lw_Qprime(); + } catch (const exception &e) { + cout << "Error: caught exception " << e.what() << endl; + ++result; + } + + if (result) { + cout << "Error: something went wrong. "; + if (result == 1) { + cout << "(one failure)" << endl; + } else { + cout << "(" << result << " individual failures)" << endl; + } + cout << "please check times.out against times.ref for more details." + << endl << "happy debugging!" << endl; + } + + return result; } diff --git a/check/times.h b/check/times.h index 90aac753..cc792bde 100644 --- a/check/times.h +++ b/check/times.h @@ -37,15 +37,15 @@ using namespace GiNaC; class timer { public: - timer(); - void start(void); - void stop(void); - void reset(void); - double read(void); - bool running(void); + timer(); + void start(void); + void stop(void); + void reset(void); + double read(void); + bool running(void); private: - bool on; - struct rusage used1, used2; + bool on; + struct rusage used1, used2; }; // prototypes for all individual timings should be unsigned fcn():