/* 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;
}
#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<m)&&(ro<n); ++ro) {
- for (unsigned co=0; co<n; ++co)
- A.set(ro,co,dense_univariate_poly(a,degree));
- for (unsigned co=0; co<p; ++co)
- B.set(ro,co,dense_univariate_poly(a,degree));
- }
- // repeat excessive rows of A and B to avoid excessive construction of
- // overdetermined linear systems
- for (unsigned ro=n; ro<m; ++ro) {
- for (unsigned co=0; co<n; ++co)
- A.set(ro,co,A(ro-1,co));
- for (unsigned co=0; co<p; ++co)
- B.set(ro,co,B(ro-1,co));
- }
- // create a vector of n*p symbols all named "xrc" where r and c are ints
- vector<symbol> x;
- matrix X(n,p);
- for (unsigned i=0; i<n; ++i) {
- for (unsigned j=0; j<p; ++j) {
- char buf[4];
- ostrstream(buf,sizeof(buf)) << i << j << ends;
- x.push_back(symbol(string("x")+buf));
- X.set(i,j,x[p*i+j]);
- }
- }
- matrix sol(n,p);
- // Solve the system A*X==B:
- try {
- sol = A.solve(X, B);
- } catch (const exception & err) { // catch runtime_error
- // Presumably, the coefficient matrix A was degenerate
- string errwhat = err.what();
- if (errwhat == "matrix::solve(): inconsistent linear system")
- return 0;
- else
- clog << "caught exception: " << errwhat << endl;
- throw;
- }
-
- // check the result with our original matrix:
- bool errorflag = false;
- for (unsigned ro=0; ro<m; ++ro) {
- for (unsigned pco=0; pco<p; ++pco) {
- ex e = 0;
- for (unsigned co=0; co<n; ++co)
- e += A(ro,co)*sol(co,pco);
- if (!(e-B(ro,pco)).normal().is_zero())
- errorflag = true;
- }
- }
- if (errorflag) {
- clog << "Our solve method claims that A*X==B, with matrices" << endl
- << "A == " << A << endl
- << "X == " << sol << endl
- << "B == " << B << endl;
- return 1;
- }
-
- return 0;
+ 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<m)&&(ro<n); ++ro) {
+ for (unsigned co=0; co<n; ++co)
+ A.set(ro,co,dense_univariate_poly(a,degree));
+
for (unsigned co=0; co<p; ++co)
+ B.set(ro,co,dense_univariate_poly(a,degree));
+ }
+ // repeat excessive rows of A and B to avoid excessive construction of
+ // overdetermined linear systems
+ for (unsigned ro=n; ro<m; ++ro) {
+ for (unsigned co=0; co<n; ++co)
+ A.set(ro,co,A(ro-1,co));
+
for (unsigned co=0; co<p; ++co)
+ B.set(ro,co,B(ro-1,co));
+ }
+ // create a vector of n*p symbols all named "xrc" where r and c are ints
+ vector<symbol> x;
+ matrix X(n,p);
+ for (unsigned i=0; i<n; ++i) {
+
for (unsigned j=0; j<p; ++j) {
+ char buf[4];
+ ostrstream(buf,sizeof(buf)) << i << j << ends;
+ x.push_back(symbol(string("x")+buf));
+ X.set(i,j,x[p*i+j]);
+ }
+ }
+ matrix sol(n,p);
+ // Solve the system A*X==B:
+ try {
+ sol = A.solve(X, B);
+ } catch (const exception & err) { // catch runtime_error
+ // Presumably, the coefficient matrix A was degenerate
+ string errwhat = err.what();
+ if (errwhat == "matrix::solve(): inconsistent linear system")
+ return 0;
+ else
+ clog << "caught exception: " << errwhat << endl;
+ throw;
+ }
+
+ // check the result with our original matrix:
+ bool errorflag = false;
+ for (unsigned ro=0; ro<m; ++ro) {
+
for (unsigned pco=0; pco<p; ++pco) {
+ ex e = 0;
+ for (unsigned co=0; co<n; ++co)
+ e += A(ro,co)*sol(co,pco);
+ if (!(e-B(ro,pco)).normal().is_zero())
+ errorflag = true;
+ }
+ }
+ if (errorflag) {
+ clog << "Our solve method claims that A*X==B, with matrices" << endl
+ << "A == " << A << endl
+ << "X == " << sol << endl
+ << "B == " << B << endl;
+ return 1;
+ }
+
+ return 0;
}
static unsigned check_inifcns_lsolve(unsigned n)
{
- unsigned result = 0;
-
- for (int repetition=0; repetition<100; ++repetition) {
- // create two size n vectors of symbols, one for the coefficients
- // a[0],..,a[n], one for indeterminates x[0]..x[n]:
- vector<symbol> a;
- vector<symbol> x;
- for (unsigned i=0; i<n; ++i) {
- char buf[3];
- ostrstream(buf,sizeof(buf)) << i << ends;
- a.push_back(symbol(string("a")+buf));
- x.push_back(symbol(string("x")+buf));
- }
- lst eqns; // equation list
- lst vars; // variable list
- ex sol; // solution
- // Create a random linear system...
- for (unsigned i=0; i<n; ++i) {
- ex lhs = rand()%201-100;
- ex rhs = rand()%201-100;
- for (unsigned j=0; j<n; ++j) {
- // ...with small coefficients to give degeneracy a chance...
- lhs += a[j]*(rand()%21-10);
- rhs += x[j]*(rand()%21-10);
- }
- eqns.append(lhs==rhs);
- vars.append(x[i]);
- }
- // ...solve it...
- sol = lsolve(eqns, vars);
-
- // ...and check the solution:
- if (sol.nops() == 0) {
- // no solution was found
- // is the coefficient matrix really, really, really degenerate?
- matrix coeffmat(n,n);
- for (unsigned ro=0; ro<n; ++ro)
- for (unsigned co=0; co<n; ++co)
- coeffmat.set(ro,co,eqns.op(co).rhs().coeff(a[co],1));
- if (!coeffmat.determinant().is_zero()) {
- ++result;
- clog << "solution of the system " << eqns << " for " << vars
- << " was not found" << endl;
- }
- } else {
- // insert the solution into rhs of out equations
- bool errorflag = false;
- for (unsigned i=0; i<n; ++i)
- if (eqns.op(i).rhs().subs(sol) != eqns.op(i).lhs())
- errorflag = true;
- if (errorflag) {
- ++result;
- clog << "solution of the system " << eqns << " for " << vars
- << " erroneously returned " << sol << endl;
- }
- }
- }
-
- return result;
+ unsigned result = 0;
+
+ for (int repetition=0; repetition<100; ++repetition) {
+ // create two size n vectors of symbols, one for the coefficients
+ // a[0],..,a[n], one for indeterminates x[0]..x[n]:
+ vector<symbol> a;
+ vector<symbol> x;
+ for (unsigned i=0; i<n; ++i) {
+ char buf[3];
+ ostrstream(buf,sizeof(buf)) << i << ends;
+ a.push_back(symbol(string("a")+buf));
+ x.push_back(symbol(string("x")+buf));
+ }
+ lst eqns; // equation list
+ lst vars; // variable list
+ ex sol; // solution
+ // Create a random linear system...
+ for (unsigned i=0; i<n; ++i) {
+ ex lhs = rand()%201-100;
+ ex rhs = rand()%201-100;
+ for (unsigned j=0; j<n; ++j) {
+ // ...with small coefficients to give degeneracy a chance...
+ lhs += a[j]*(rand()%21-10);
+ rhs += x[j]*(rand()%21-10);
+ }
+ eqns.append(lhs==rhs);
+ vars.append(x[i]);
+ }
+ // ...solve it...
+ sol = lsolve(eqns, vars);
+
+ // ...and check the solution:
+ if (sol.nops() == 0) {
+ // no solution was found
+ // is the coefficient matrix really, really, really degenerate?
+ matrix coeffmat(n,n);
+ for (unsigned ro=0; ro<n; ++ro)
+ for (unsigned co=0; co<n; ++co)
+ coeffmat.set(ro,co,eqns.op(co).rhs().coeff(a[co],1));
+ if (!coeffmat.determinant().is_zero()) {
+ ++result;
+ clog << "solution of the system " << eqns << " for " << vars
+ << " was not found" << endl;
+ }
+ } else {
+ // insert the solution into rhs of out equations
+ bool errorflag = false;
+ for (unsigned i=0; i<n; ++i)
+ if (eqns.op(i).rhs().subs(sol) != eqns.op(i).lhs())
+ errorflag = true;
+ if (errorflag) {
+ ++result;
+ clog << "solution of the system " << eqns << " for " << vars
+ << " erroneously returned " << sol << endl;
+ }
+ }
+ }
+
+ return result;
}
unsigned check_lsolve(void)
{
- unsigned result = 0;
-
- cout << "checking linear solve" << flush;
- clog << "---------linear solve:" << endl;
-
- // solve some numeric linear systems
- for (unsigned n=1; n<12; ++n)
- result += check_matrix_solve(n, n, 1, 0);
- cout << '.' << flush;
- // solve some underdetermined numeric systems
- for (unsigned n=1; n<12; ++n)
- result += check_matrix_solve(n+1, n, 1, 0);
- cout << '.' << flush;
- // solve some overdetermined numeric systems
- for (unsigned n=1; n<12; ++n)
- result += check_matrix_solve(n, n+1, 1, 0);
- cout << '.' << flush;
- // solve some multiple numeric systems
- for (unsigned n=1; n<12; ++n)
- result += check_matrix_solve(n, n, n/3+1, 0);
- cout << '.' << flush;
- // solve some symbolic linear systems
- for (unsigned n=1; n<7; ++n)
- result += check_matrix_solve(n, n, 1, 2);
- cout << '.' << flush;
-
- // check lsolve, the wrapper function around matrix::solve()
- result += check_inifcns_lsolve(2); cout << '.' << flush;
- result += check_inifcns_lsolve(3); cout << '.' << flush;
- result += check_inifcns_lsolve(4); cout << '.' << flush;
- result += check_inifcns_lsolve(5); cout << '.' << flush;
- result += check_inifcns_lsolve(6); cout << '.' << flush;
-
- if (!result) {
- cout << " passed " << endl;
- clog << "(no output)" << endl;
- } else {
- cout << " failed " << endl;
- }
-
- return result;
+ unsigned result = 0;
+
+ cout << "checking linear solve" << flush;
+ clog << "---------linear solve:" << endl;
+
+ // solve some numeric linear systems
+ for (unsigned n=1; n<12; ++n)
+ result += check_matrix_solve(n, n, 1, 0);
+ cout << '.' << flush;
+ // solve some underdetermined numeric systems
+ for (unsigned n=1; n<12; ++n)
+ result += check_matrix_solve(n+1, n, 1, 0);
+ cout << '.' << flush;
+ // solve some overdetermined numeric systems
+ for (unsigned n=1; n<12; ++n)
+ result += check_matrix_solve(n, n+1, 1, 0);
+ cout << '.' << flush;
+ // solve some multiple numeric systems
+ for (unsigned n=1; n<12; ++n)
+ result += check_matrix_solve(n, n, n/3+1, 0);
+ cout << '.' << flush;
+ // solve some symbolic linear systems
+ for (unsigned n=1; n<7; ++n)
+ result += check_matrix_solve(n, n, 1, 2);
+ cout << '.' << flush;
+
+ // check lsolve, the wrapper function around matrix::solve()
+ result += check_inifcns_lsolve(2); cout << '.' << flush;
+ result += check_inifcns_lsolve(3); cout << '.' << flush;
+ result += check_inifcns_lsolve(4); cout << '.' << flush;
+ result += check_inifcns_lsolve(5); cout << '.' << flush;
+ result += check_inifcns_lsolve(6); cout << '.' << flush;
+
+ if (!result) {
+ cout << " passed " << endl;
+ clog << "(no output)" << endl;
+ } else {
+ cout << " failed " << endl;
+ }
+
+ return result;
}
* an integral domain. */
static unsigned integdom_matrix_determinants(void)
{
- unsigned result = 0;
- symbol a("a");
-
- for (unsigned size=3; size<20; ++size) {
- matrix A(size,size);
- // populate one element in each row:
- for (unsigned r=0; r<size-1; ++r)
- A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
- // set the last row to a linear combination of two other lines
- // to guarantee that the determinant is zero:
- for (unsigned c=0; c<size; ++c)
- A.set(size-1,c,A(0,c)-A(size-2,c));
- if (!A.determinant().is_zero()) {
- clog << "Determinant of " << size << "x" << size << " matrix "
- << endl << A << endl
- << "was not found to vanish!" << endl;
- ++result;
- }
- }
-
- return result;
+ unsigned result = 0;
+ symbol a("a");
+
+ for (unsigned size=3; size<20; ++size) {
+ matrix A(size,size);
+ // populate one element in each row:
+ for (unsigned r=0; r<size-1; ++r)
+ A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
+ // set the last row to a linear combination of two other lines
+ // to guarantee that the determinant is zero:
+ for (unsigned c=0; c<size; ++c)
+ A.set(size-1,c,A(0,c)-A(size-2,c));
+ if (!A.determinant().is_zero()) {
+ clog << "Determinant of " << size << "x" << size << " matrix "
+ << endl << A << endl
+ << "was not found to vanish!" << endl;
+ ++result;
+ }
+ }
+
+ return result;
}
/* determinants of some symbolic matrices with multivariate rational function
* coefficients. */
static unsigned rational_matrix_determinants(void)
{
- unsigned result = 0;
- symbol a("a"), b("b"), c("c");
-
- for (unsigned size=3; size<8; ++size) {
- matrix A(size,size);
- for (unsigned r=0; r<size-1; ++r) {
- // populate one or two elements in each row:
- for (unsigned ec=0; ec<2; ++ec) {
- ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
- ex denom;
- do {
- denom = sparse_tree(a, b, c, rand()%2, false, false, false);
- } while (denom.is_zero());
- A.set(r,unsigned(rand()%size),numer/denom);
- }
- }
- // set the last row to a linear combination of two other lines
- // to guarantee that the determinant is zero:
- for (unsigned co=0; co<size; ++co)
- A.set(size-1,co,A(0,co)-A(size-2,co));
- if (!A.determinant().is_zero()) {
- clog << "Determinant of " << size << "x" << size << " matrix "
- << endl << A << endl
- << "was not found to vanish!" << endl;
- ++result;
- }
- }
-
- return result;
+ unsigned result = 0;
+ symbol a("a"), b("b"), c("c");
+
+ for (unsigned size=3; size<8; ++size) {
+ matrix A(size,size);
+ for (unsigned r=0; r<size-1; ++r) {
+ // populate one or two elements in each row:
+ for (unsigned ec=0; ec<2; ++ec) {
+ ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
+ ex denom;
+ do {
+ denom = sparse_tree(a, b, c, rand()%2, false, false, false);
+ } while (denom.is_zero());
+ A.set(r,unsigned(rand()%size),numer/denom);
+ }
+ }
+ // set the last row to a linear combination of two other lines
+ // to guarantee that the determinant is zero:
+ for (unsigned co=0; co<size; ++co)
+ A.set(size-1,co,A(0,co)-A(size-2,co));
+ if (!A.determinant().is_zero()) {
+ clog << "Determinant of " << size << "x" << size << " matrix "
+ << endl << A << endl
+ << "was not found to vanish!" << endl;
+ ++result;
+ }
+ }
+
+ return result;
}
/* Some quite funny determinants with functions and stuff like that inside. */
static unsigned funny_matrix_determinants(void)
{
- unsigned result = 0;
- symbol a("a"), b("b"), c("c");
-
- for (unsigned size=3; size<7; ++size) {
- matrix A(size,size);
- for (unsigned co=0; co<size-1; ++co) {
- // populate one or two elements in each row:
- for (unsigned ec=0; ec<2; ++ec) {
- ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
- ex denom;
- do {
- denom = sparse_tree(a, b, c, rand()%2, false, true, false);
- } while (denom.is_zero());
- A.set(unsigned(rand()%size),co,numer/denom);
- }
- }
- // set the last column to a linear combination of two other columns
- // to guarantee that the determinant is zero:
- for (unsigned ro=0; ro<size; ++ro)
- A.set(ro,size-1,A(ro,0)-A(ro,size-2));
- if (!A.determinant().is_zero()) {
- clog << "Determinant of " << size << "x" << size << " matrix "
- << endl << A << endl
- << "was not found to vanish!" << endl;
- ++result;
- }
- }
-
- return result;
+ unsigned result = 0;
+ symbol a("a"), b("b"), c("c");
+
+ for (unsigned size=3; size<7; ++size) {
+ matrix A(size,size);
+ for (unsigned co=0; co<size-1; ++co) {
+ // populate one or two elements in each row:
+ for (unsigned ec=0; ec<2; ++ec) {
+ ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
+ ex denom;
+ do {
+ denom = sparse_tree(a, b, c, rand()%2, false, true, false);
+ } while (denom.is_zero());
+ A.set(unsigned(rand()%size),co,numer/denom);
+ }
+ }
+ // set the last column to a linear combination of two other columns
+ // to guarantee that the determinant is zero:
+ for (unsigned ro=0; ro<size; ++ro)
+ A.set(ro,size-1,A(ro,0)-A(ro,size-2));
+ if (!A.determinant().is_zero()) {
+ clog << "Determinant of " << size << "x" << size << " matrix "
+ << endl << A << endl
+ << "was not found to vanish!" << endl;
+ ++result;
+ }
+ }
+
+ return result;
}
/* compare results from different determinant algorithms.*/
static unsigned compare_matrix_determinants(void)
{
- unsigned result = 0;
- symbol a("a");
-
- for (unsigned size=2; size<7; ++size) {
- matrix A(size,size);
- for (unsigned co=0; co<size; ++co) {
- for (unsigned ro=0; ro<size; ++ro) {
- // populate some elements
- ex elem = 0;
- if (rand()%(size/2) == 0)
- elem = sparse_tree(a, a, a, rand()%3, false, true, false);
- A.set(ro,co,elem);
- }
- }
- ex det_gauss = A.determinant(determinant_algo::gauss);
- ex det_laplace = A.determinant(determinant_algo::laplace);
- ex det_divfree = A.determinant(determinant_algo::divfree);
- ex det_bareiss = A.determinant(determinant_algo::bareiss);
- if ((det_gauss-det_laplace).normal() != 0 ||
- (det_bareiss-det_laplace).normal() != 0 ||
- (det_divfree-det_laplace).normal() != 0) {
- clog << "Determinant of " << size << "x" << size << " matrix "
- << endl << A << endl
- << "is inconsistent between different algorithms:" << endl
- << "Gauss elimination: " << det_gauss << endl
- << "Minor elimination: " << det_laplace << endl
- << "Division-free elim.: " << det_divfree << endl
- << "Fraction-free elim.: " << det_bareiss << endl;
- ++result;
- }
- }
-
- return result;
+ unsigned result = 0;
+ symbol a("a");
+
+ for (unsigned size=2; size<7; ++size) {
+ matrix A(size,size);
+ for (unsigned co=0; co<size; ++co) {
+ for (unsigned ro=0; ro<size; ++ro) {
+ // populate some elements
+ ex elem = 0;
+ if (rand()%(size/2) == 0)
+ elem = sparse_tree(a, a, a, rand()%3, false, true, false);
+ A.set(ro,co,elem);
+ }
+ }
+ ex det_gauss = A.determinant(determinant_algo::gauss);
+ ex det_laplace = A.determinant(determinant_algo::laplace);
+ ex det_divfree = A.determinant(determinant_algo::divfree);
+ ex det_bareiss = A.determinant(determinant_algo::bareiss);
+ if ((det_gauss-det_laplace).normal() != 0 ||
+ (det_bareiss-det_laplace).normal() != 0 ||
+ (det_divfree-det_laplace).normal() != 0) {
+ clog << "Determinant of " << size << "x" << size << " matrix "
+ << endl << A << endl
+ << "is inconsistent between different algorithms:" << endl
+ << "Gauss elimination: " << det_gauss << endl
+ << "Minor elimination: " << det_laplace << endl
+ << "Division-free elim.: " << det_divfree << endl
+ << "Fraction-free elim.: " << det_bareiss << endl;
+ ++result;
+ }
+ }
+
+ return result;
}
static unsigned symbolic_matrix_inverse(void)
{
- unsigned result = 0;
- symbol a("a"), b("b"), c("c");
-
- for (unsigned size=2; size<5; ++size) {
- matrix A(size,size);
- do {
- for (unsigned co=0; co<size; ++co) {
- for (unsigned ro=0; ro<size; ++ro) {
- // populate some elements
- ex elem = 0;
- if (rand()%(size/2) == 0)
- elem = sparse_tree(a, b, c, rand()%2, false, true, false);
- A.set(ro,co,elem);
- }
- }
- } while (A.determinant() == 0);
- matrix B = A.inverse();
- matrix C = A.mul(B);
- bool ok = true;
- for (unsigned ro=0; ro<size; ++ro)
- for (unsigned co=0; co<size; ++co)
- if (C(ro,co).normal() != (ro==co?1:0))
- ok = false;
- if (!ok) {
- clog << "Inverse of " << size << "x" << size << " matrix "
- << endl << A << endl
- << "erroneously returned: "
- << endl << B << endl;
- ++result;
- }
- }
-
- return result;
+ unsigned result = 0;
+ symbol a("a"), b("b"), c("c");
+
+ for (unsigned size=2; size<5; ++size) {
+ matrix A(size,size);
+ do {
+ for (unsigned co=0; co<size; ++co) {
+ for (unsigned ro=0; ro<size; ++ro) {
+ // populate some elements
+ ex elem = 0;
+ if (rand()%(size/2) == 0)
+ elem = sparse_tree(a, b, c, rand()%2, false, true, false);
+ A.set(ro,co,elem);
+ }
+ }
+ } while (A.determinant() == 0);
+ matrix B = A.inverse();
+ matrix C = A.mul(B);
+ bool ok = true;
+ for (unsigned ro=0; ro<size; ++ro)
+ for (unsigned co=0; co<size; ++co)
+ if (C(ro,co).normal() != (ro==co?1:0))
+ ok = false;
+ if (!ok) {
+ clog << "Inverse of " << size << "x" << size << " matrix "
+ << endl << A << endl
+ << "erroneously returned: "
+ << endl << B << endl;
+ ++result;
+ }
+ }
+
+ return result;
}
unsigned check_matrices(void)
{
- unsigned result = 0;
-
- cout << "checking symbolic matrix manipulations" << flush;
- clog << "---------symbolic matrix manipulations:" << endl;
-
- result += integdom_matrix_determinants(); cout << '.' << flush;
- result += rational_matrix_determinants(); cout << '.' << flush;
- result += funny_matrix_determinants(); cout << '.' << flush;
- result += compare_matrix_determinants(); cout << '.' << flush;
- result += symbolic_matrix_inverse(); cout << '.' << flush;
-
- if (!result) {
- cout << " passed " << endl;
- clog << "(no output)" << endl;
- } else {
- cout << " failed " << endl;
- }
-
- return result;
+ unsigned result = 0;
+
+ cout << "checking symbolic matrix manipulations" << flush;
+ clog << "---------symbolic matrix manipulations:" << endl;
+
+ result += integdom_matrix_determinants(); cout << '.' << flush;
+ result += rational_matrix_determinants(); cout << '.' << flush;
+ result += funny_matrix_determinants(); cout << '.' << flush;
+ result += compare_matrix_determinants(); cout << '.' << flush;
+ result += symbolic_matrix_inverse(); cout << '.' << flush;
+
+ if (!result) {
+ cout << " passed " << endl;
+ clog << "(no output)" << endl;
+ } else {
+ cout << " failed " << endl;
+ }
+
+ return result;
}
* conversions. */
static unsigned check_numeric1(void)
{
- unsigned result = 0;
- bool errorflag = false;
- int re_q, im_q;
-
- // Check some numerator and denominator calculations:
- for (int i=0; i<200; ++i) {
- do { re_q = rand(); } while (re_q == 0);
- do { im_q = rand(); } while (im_q == 0);
- numeric r(rand()-RAND_MAX/2, re_q);
- numeric i(rand()-RAND_MAX/2, im_q);
- numeric z = r + I*i;
- numeric p = numer(z);
- numeric q = denom(z);
- numeric res = p/q;
- if (res != z) {
- clog << z << " erroneously transformed into "
- << p << "/" << q << " by numer() and denom()" << endl;
- errorflag = true;
- }
- }
- if (errorflag)
- ++result;
-
- return result;
+ unsigned result = 0;
+ bool errorflag = false;
+ int re_q, im_q;
+
+ // Check some numerator and denominator calculations:
+ for (int rep=0; rep<200; ++rep) {
+ do { re_q = rand(); } while (re_q == 0);
+ do { im_q = rand(); } while (im_q == 0);
+ numeric r(rand()-RAND_MAX/2, re_q);
+ numeric i(rand()-RAND_MAX/2, im_q);
+ numeric z = r + I*i;
+ numeric p = numer(z);
+ numeric q = denom(z);
+ numeric res = p/q;
+ if (res != z) {
+ clog << z << " erroneously transformed into "
+ << p << "/" << q << " by numer() and denom()" << endl;
+ errorflag = true;
+ }
+ }
+ if (errorflag)
+ ++result;
+
+ return result;
}
static unsigned check_numeric2(void)
{
- 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 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;
}
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;
}
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();
#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;
}
// 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;
}
// 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;
}
/* 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
* 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;
}
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;
}
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;
}
/** @file exam_misc.cpp
*
*/
+
/*
* GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
*
#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<VECSIZE; ++i) {
- e = e + a[i];
- }
-
- // prepare aux so it will swallow anything but a1^2:
- aux = -e + a[0] + a[1];
- e = expand(subs(expand(pow(e, 2)), a[0] == aux));
-
- if (e != pow(a1,2)) {
- clog << "Denny Fliegner's quick consistency check erroneously returned "
- << e << "." << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+ symbol a1("a1");
+ symbol a[VECSIZE];
+ ex e, aux;
+
+ a[1] = a1;
+ for (unsigned i=0; i<VECSIZE; ++i) {
+ e = e + a[i];
+ }
+
+ // prepare aux so it will swallow anything but a1^2:
+ aux = -e + a[0] + a[1];
+ e = expand(subs(expand(pow(e, 2)), a[0] == aux));
+
+ if (e != pow(a1,2)) {
+ clog << "Denny Fliegner's quick consistency check erroneously returned "
+ << e << "." << endl;
+ ++result;
+ }
+
+ return result;
}
/* A simple modification of Denny Fliegner's three step consistency test:
* after which e should return 0 (without expanding). */
static unsigned exam_expand_subs2(void)
{
- unsigned result = 0;
- symbol a("a"), b("b");
- ex e, f;
-
- e = pow(a+b,200).expand();
- f = e.subs(a == -b);
-
- if (f != 0) {
- clog << "e = pow(a+b,200).expand(); f = e.subs(a == -b); erroneously returned "
- << f << " instead of simplifying to 0." << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+ symbol a("a"), b("b");
+ ex e, f;
+
+ e = pow(a+b,200).expand();
+ f = e.subs(a == -b);
+
+ if (f != 0) {
+ clog << "e = pow(a+b,200).expand(); f = e.subs(a == -b); erroneously returned "
+ << f << " instead of simplifying to 0." << endl;
+ ++result;
+ }
+
+ return result;
}
unsigned exam_misc(void)
{
- unsigned result = 0;
-
- cout << "examining miscellaneous other things" << flush;
- clog << "----------miscellaneous other things:" << endl;
-
- result += exam_expand_subs(); cout << '.' << flush;
- result += exam_expand_subs2(); cout << '.' << flush;
-
- if (!result) {
- cout << " passed " << endl;
- clog << "(no output)" << endl;
- } else {
- cout << " failed " << endl;
- }
-
- return result;
+ unsigned result = 0;
+
+ cout << "examining miscellaneous other things" << flush;
+ clog << "----------miscellaneous other things:" << endl;
+
+ result += exam_expand_subs(); cout << '.' << flush;
+ result += exam_expand_subs2(); cout << '.' << flush;
+
+ if (!result) {
+ cout << " passed " << endl;
+ clog << "(no output)" << endl;
+ } else {
+ cout << " failed " << endl;
+ }
+
+ return result;
}
static unsigned lortensor_check1(void)
{
- // checks simple identities of the metric tensor!
-
- unsigned result = 0;
- lorentzidx mu("mu"), nu("nu");
- ex e1, e2, e3, e4, e5, e6;
- e1 = lortensor_g(mu,nu);
- e2 = lortensor_g(nu,mu);
- e3 = e1 - e2; // g(~mu,~nu) - g(~nu,~mu) = 0 !
- e4 = lortensor_g(mu,mu.toggle_covariant());
- e5 = lortensor_g(mu.toggle_covariant(),mu);
- e6 = e4 - e5; // g(~mu,_mu) - g(_mu,~mu) = 0!
- if (!e3.is_zero()) {
- clog << e1 << "-" << e2 << " erroneously returned "
- << e3 << " instead of 0" << endl;
- ++result;
- }
- if (!e6.is_zero()) {
- clog << e4 << "-" << e5 << " erroneously returned "
- << e6 << " instead of 0" << endl;
- ++result;
- }
+ // checks simple identities of the metric tensor!
+
+ unsigned result = 0;
+ lorentzidx mu("mu"), nu("nu");
+ ex e1, e2, e3, e4, e5, e6;
+ e1 = lortensor_g(mu,nu);
+ e2 = lortensor_g(nu,mu);
+ e3 = e1 - e2; // g(~mu,~nu) - g(~nu,~mu) = 0 !
+ e4 = lortensor_g(mu,mu.toggle_covariant());
+ e5 = lortensor_g(mu.toggle_covariant(),mu);
+ e6 = e4 - e5; // g(~mu,_mu) - g(_mu,~mu) = 0!
+ if (!e3.is_zero()) {
+ clog << e1 << "-" << e2 << " erroneously returned "
+ << e3 << " instead of 0" << endl;
+ ++result;
+ }
+ if (!e6.is_zero()) {
+ clog << e4 << "-" << e5 << " erroneously returned "
+ << e6 << " instead of 0" << endl;
+ ++result;
+ }
- return result;
+ return result;
}
static unsigned lortensor_check2(void)
{
- // checks simple contraction properties of an arbitrary (symmetric!) rankn lortensor!
-
- unsigned result = 0;
- lorentzidx mu("mu"), nu("nu"), rho("rho");
- ex e1, e2, e3, e4, e5, e6, e7, e8, e9, e10;
- e1 = lortensor_g(mu,nu);
- e2 = lortensor_g(nu,mu);
- e3 = lortensor_rank1("p",mu.toggle_covariant());
- e4 = lortensor_rank1("p",nu);
- e5 = e3 * e1 - e3 * e2; // p_mu g(~mu,~nu) - p_mu g(~nu,~mu) = 0!
- e6 = simplify_lortensor(e3 * e1) - e4; // p~nu - p~nu = 0!
- e7 = lortensor_g(nu,rho);
- e8 = lortensor_rank2("F",mu.toggle_covariant(),nu.toggle_covariant());
- e9 = lortensor_rank2("F",mu.toggle_covariant(),rho);
- e10 = simplify_lortensor(e8 * e7) - e9; // F(_mu,_nu) g(~nu,~rho) - F(_mu,~rho) = 0!
- if (!e5.is_zero()) {
- clog << e3 << "*" << e1 << "-" << e3 << "*" << e2 << " erroneously returned "
- << e5 << " instead of 0" << endl;
- ++result;
- }
- if (!e6.is_zero()) {
- clog << " simplify_lortensor(e3 * e1)" << "-" << e4 << " erroneously returned"
- << e6 << " instead of 0" << endl;
- ++result;
- }
- if (!e10.is_zero()) {
- clog << " simplify_lortensor(e8 * e7)" << "-" << e9 << " erroneously returned"
- << e10 << " instead of 0" << endl;
- ++result;
- }
+ // checks simple contraction properties of an arbitrary (symmetric!) rankn lortensor!
+
+ unsigned result = 0;
+ lorentzidx mu("mu"), nu("nu"), rho("rho");
+ ex e1, e2, e3, e4, e5, e6, e7, e8, e9, e10;
+ e1 = lortensor_g(mu,nu);
+ e2 = lortensor_g(nu,mu);
+ e3 = lortensor_rank1("p",mu.toggle_covariant());
+ e4 = lortensor_rank1("p",nu);
+ e5 = e3 * e1 - e3 * e2; // p_mu g(~mu,~nu) - p_mu g(~nu,~mu) = 0!
+ e6 = simplify_lortensor(e3 * e1) - e4; // p~nu - p~nu = 0!
+ e7 = lortensor_g(nu,rho);
+ e8 = lortensor_rank2("F",mu.toggle_covariant(),nu.toggle_covariant());
+ e9 = lortensor_rank2("F",mu.toggle_covariant(),rho);
+ e10 = simplify_lortensor(e8 * e7) - e9; // F(_mu,_nu) g(~nu,~rho) - F(_mu,~rho) = 0!
+ if (!e5.is_zero()) {
+ clog << e3 << "*" << e1 << "-" << e3 << "*" << e2 << " erroneously returned "
+ << e5 << " instead of 0" << endl;
+ ++result;
+ }
+ if (!e6.is_zero()) {
+ clog << " simplify_lortensor(e3 * e1)" << "-" << e4 << " erroneously returned"
+ << e6 << " instead of 0" << endl;
+ ++result;
+ }
+ if (!e10.is_zero()) {
+ clog << " simplify_lortensor(e8 * e7)" << "-" << e9 << " erroneously returned"
+ << e10 << " instead of 0" << endl;
+ ++result;
+ }
- return result;
+ return result;
}
unsigned exam_noncommut(void)
{
- unsigned result = 0;
-
- cout << "examining behaviour of noncommutative objects" << flush;
- clog << "----------behaviour of noncommutative objects:" << endl;
-
- result += lortensor_check1(); cout << '.' << flush;
- result += lortensor_check2(); cout << '.' << flush;
-
- if (!result) {
- cout << " passed " << endl;
- clog << "(no output)" << endl;
- } else {
- cout << " failed " << endl;
- }
-
- return result;
+ unsigned result = 0;
+
+ cout << "examining behaviour of noncommutative objects" << flush;
+ clog << "----------behaviour of noncommutative objects:" << endl;
+
+ result += lortensor_check1(); cout << '.' << flush;
+ result += lortensor_check2(); cout << '.' << flush;
+
+ if (!result) {
+ cout << " passed " << endl;
+ clog << "(no output)" << endl;
+ } else {
+ cout << " failed " << endl;
+ }
+
+ return result;
}
static unsigned check_normal(const ex &e, const ex &d)
{
- ex en = e.normal();
- if (en.compare(d) != 0) {
- clog << "normal form of " << e << " erroneously returned "
- << en << " (should be " << d << ")" << endl;
- return 1;
- }
- return 0;
+ ex en = e.normal();
+ if (en.compare(d) != 0) {
+ clog << "normal form of " << e << " erroneously returned "
+ << en << " (should be " << d << ")" << endl;
+ return 1;
+ }
+ return 0;
}
static unsigned exam_normal1(void)
{
- unsigned result = 0;
- ex e, d;
-
- // Expansion
- e = pow(x, 2) - (x+1)*(x-1) - 1;
- d = ex(0);
- result += check_normal(e, d);
-
- // Expansion inside functions
- e = sin(x*(x+1)-x) + 1;
- d = sin(pow(x, 2)) + 1;
- result += check_normal(e, d);
-
- // Fraction addition
- e = 2/x + y/3;
- d = (x*y + 6) / (x*3);
- result += check_normal(e, d);
-
- e = pow(x, -1) + x/(x+1);
- d = (pow(x, 2)+x+1)/(x*(x+1));
- result += check_normal(e, d);
+ unsigned result = 0;
+ ex e, d;
+
+ // Expansion
+ e = pow(x, 2) - (x+1)*(x-1) - 1;
+ d = 0;
+ result += check_normal(e, d);
+
+ // Expansion inside functions
+ e = sin(x*(x+1)-x) + 1;
+ d = sin(pow(x, 2)) + 1;
+ result += check_normal(e, d);
+
+ // Fraction addition
+ e = 2/x + y/3;
+ d = (x*y + 6) / (x*3);
+ result += check_normal(e, d);
+
+ e = pow(x, -1) + x/(x+1);
+ d = (pow(x, 2)+x+1)/(x*(x+1));
+ result += check_normal(e, d);
- return result;
+ return result;
}
static unsigned exam_normal2(void)
{
- unsigned result = 0;
- ex e, d;
-
- // Fraction cancellation
- e = numeric(1)/2 * z * (2*x + 2*y);
- d = z * (x + y);
- result += check_normal(e, d);
-
- e = numeric(1)/6 * z * (3*x + 3*y) * (2*x + 2*w);
- d = z * (x + y) * (x + w);
- result += check_normal(e, d);
-
- e = (3*x + 3*y) * (w/3 + z/3);
- d = (x + y) * (w + z);
- result += check_normal(e, d);
-
- e = (pow(x, 2) - pow(y, 2)) / pow(x-y, 3);
- d = (x + y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
- result += check_normal(e, d);
-
- e = (pow(x, -1) + x) / (pow(x , 2) * 2 + 2);
- d = pow(x * 2, -1);
- result += check_normal(e, d);
-
- // Fraction cancellation with rational coefficients
- e = (pow(x, 2) - pow(y, 2)) / pow(x/2 - y/2, 3);
- d = (8 * x + 8 * y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
- result += check_normal(e, d);
-
- // Fraction cancellation with rational coefficients
- e = z/5 * (x/7 + y/10) / (x/14 + y/20);
- d = 2*z/5;
- result += check_normal(e, d);
-
- return result;
+ unsigned result = 0;
+ ex e, d;
+
+ // Fraction cancellation
+ e = numeric(1)/2 * z * (2*x + 2*y);
+ d = z * (x + y);
+ result += check_normal(e, d);
+
+ e = numeric(1)/6 * z * (3*x + 3*y) * (2*x + 2*w);
+ d = z * (x + y) * (x + w);
+ result += check_normal(e, d);
+
+ e = (3*x + 3*y) * (w/3 + z/3);
+ d = (x + y) * (w + z);
+ result += check_normal(e, d);
+
+ e = (pow(x, 2) - pow(y, 2)) / pow(x-y, 3);
+ d = (x + y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
+ result += check_normal(e, d);
+
+ e = (pow(x, -1) + x) / (pow(x , 2) * 2 + 2);
+ d = pow(x * 2, -1);
+ result += check_normal(e, d);
+
+ // Fraction cancellation with rational coefficients
+ e = (pow(x, 2) - pow(y, 2)) / pow(x/2 - y/2, 3);
+ d = (8 * x + 8 * y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
+ result += check_normal(e, d);
+
+ // Fraction cancellation with rational coefficients
+ e = z/5 * (x/7 + y/10) / (x/14 + y/20);
+ d = 2*z/5;
+ result += check_normal(e, d);
+
+ return result;
}
static unsigned exam_normal3(void)
{
- unsigned result = 0;
- ex e, d;
-
- // Distribution of powers
- e = pow(x/y, 2);
- d = pow(x, 2) / pow(y, 2);
- result += check_normal(e, d);
-
- // Distribution of powers (integer, distribute) and fraction addition
- e = pow(pow(x, -1) + x, 2);
- d = pow(pow(x, 2) + 1, 2) / pow(x, 2);
- result += check_normal(e, d);
-
- // Distribution of powers (non-integer, don't distribute) and fraction addition
- e = pow(pow(x, -1) + x, numeric(1)/2);
- d = pow((pow(x, 2) + 1) / x, numeric(1)/2);
- result += check_normal(e, d);
-
- return result;
+ unsigned result = 0;
+ ex e, d;
+
+ // Distribution of powers
+ e = pow(x/y, 2);
+ d = pow(x, 2) / pow(y, 2);
+ result += check_normal(e, d);
+
+ // Distribution of powers (integer, distribute) and fraction addition
+ e = pow(pow(x, -1) + x, 2);
+ d = pow(pow(x, 2) + 1, 2) / pow(x, 2);
+ result += check_normal(e, d);
+
+ // Distribution of powers (non-integer, don't distribute) and fraction addition
+ e = pow(pow(x, -1) + x, numeric(1)/2);
+ d = pow((pow(x, 2) + 1) / x, numeric(1)/2);
+ result += check_normal(e, d);
+
+ return result;
}
static unsigned exam_normal4(void)
{
- unsigned result = 0;
- ex e, d;
-
- // Replacement of functions with temporary symbols and fraction cancellation
- e = pow(sin(x), 2) - pow(cos(x), 2);
- e /= sin(x) + cos(x);
- d = sin(x) - cos(x);
- result += check_normal(e, d);
-
- // Replacement of non-integer powers with temporary symbols
- e = (pow(numeric(2), numeric(1)/2) * x + x) / x;
- d = pow(numeric(2), numeric(1)/2) + 1;
- result += check_normal(e, d);
-
- // Replacement of complex numbers with temporary symbols
- e = (x + y + x*I + y*I) / (x + y);
- d = 1 + I;
- result += check_normal(e, d);
-
- e = (pow(x, 2) + pow(y, 2)) / (x + y*I);
- d = e;
- result += check_normal(e, d);
-
- // More complex rational function
- e = (pow(x-y*2,4)/pow(pow(x,2)-pow(y,2)*4,2)+1)*(x+y*2)*(y+z)/(pow(x,2)+pow(y,2)*4);
- d = (y*2 + z*2) / (x + y*2);
- result += check_normal(e, d);
-
- return result;
+ unsigned result = 0;
+ ex e, d;
+
+ // Replacement of functions with temporary symbols and fraction cancellation
+ e = pow(sin(x), 2) - pow(cos(x), 2);
+ e /= sin(x) + cos(x);
+ d = sin(x) - cos(x);
+ result += check_normal(e, d);
+
+ // Replacement of non-integer powers with temporary symbols
+ e = (pow(numeric(2), numeric(1)/2) * x + x) / x;
+ d = pow(numeric(2), numeric(1)/2) + 1;
+ result += check_normal(e, d);
+
+ // Replacement of complex numbers with temporary symbols
+ e = (x + y + x*I + y*I) / (x + y);
+ d = 1 + I;
+ result += check_normal(e, d);
+
+ e = (pow(x, 2) + pow(y, 2)) / (x + y*I);
+ d = e;
+ result += check_normal(e, d);
+
+ // More complex rational function
+ e = (pow(x-y*2,4)/pow(pow(x,2)-pow(y,2)*4,2)+1)*(x+y*2)*(y+z)/(pow(x,2)+pow(y,2)*4);
+ d = (y*2 + z*2) / (x + y*2);
+ result += check_normal(e, d);
+
+ return result;
}
unsigned exam_normalization(void)
{
- unsigned result = 0;
-
- cout << "examining rational function normalization" << flush;
- clog << "----------rational function normalization:" << endl;
-
- result += exam_normal1(); cout << '.' << flush;
- result += exam_normal2(); cout << '.' << flush;
- result += exam_normal3(); cout << '.' << flush;
- result += exam_normal4(); cout << '.' << flush;
-
- if (!result) {
- cout << " passed " << endl;
- clog << "(no output)" << endl;
- } else {
- cout << " failed " << endl;
- }
-
- return result;
+ unsigned result = 0;
+
+ cout << "examining rational function normalization" << flush;
+ clog << "----------rational function normalization:" << endl;
+
+ result += exam_normal1(); cout << '.' << flush;
+ result += exam_normal2(); cout << '.' << flush;
+ result += exam_normal3(); cout << '.' << flush;
+ result += exam_normal4(); cout << '.' << flush;
+
+ if (!result) {
+ cout << " passed " << endl;
+ clog << "(no output)" << endl;
+ } else {
+ cout << " failed " << endl;
+ }
+
+ return result;
}
* conversions. */
static unsigned exam_numeric1(void)
{
- unsigned result = 0;
- numeric test_int1(42);
- numeric test_int2(5);
- numeric test_rat1 = test_int1; test_rat1 /= test_int2;
- test_rat1 = -test_rat1; // -42/5
- numeric test_crat = test_rat1+I*test_int2; // 5*I-42/5
- symbol a("a");
- ex e1, e2;
-
- if (!test_int1.is_integer()) {
- clog << test_int1
- << " erroneously not recognized as integer" << endl;
- ++result;
- }
- if (!test_int1.is_rational()) {
- clog << test_int1
- << " erroneously not recognized as rational" << endl;
- ++result;
- }
-
- if (!test_rat1.is_rational()) {
- clog << test_rat1
- << " erroneously not recognized as rational" << endl;
- ++result;
- }
- if (test_rat1.is_integer()) {
- clog << test_rat1
- << " erroneously recognized as integer" << endl;
- ++result;
- }
-
- if (!test_crat.is_crational()) {
- clog << test_crat
- << " erroneously not recognized as complex rational" << endl;
- ++result;
- }
-
- int i = numeric(1984).to_int();
- if (i-1984) {
- clog << "conversion of " << i
- << " from numeric to int failed" << endl;
- ++result;
- }
-
- e1 = test_int1;
- if (!e1.info(info_flags::posint)) {
- clog << "expression " << e1
- << " erroneously not recognized as positive integer" << endl;
- ++result;
- }
-
- e2 = test_int1 + a;
- if (ex_to_numeric(e2).is_integer()) {
- clog << "expression " << e2
- << " erroneously recognized as integer" << endl;
- ++result;
- }
-
- // The next two were two actual bugs in CLN till June, 12, 1999:
- test_rat1 = numeric(3)/numeric(2);
- test_rat1 += test_rat1;
- if (!test_rat1.is_integer()) {
- clog << "3/2 + 3/2 erroneously not integer 3 but instead "
- << test_rat1 << endl;
- ++result;
- }
- test_rat1 = numeric(3)/numeric(2);
- numeric test_rat2 = test_rat1 + numeric(1); // 5/2
- test_rat2 -= test_rat1; // 1
- if (!test_rat2.is_integer()) {
- clog << "5/2 - 3/2 erroneously not integer 1 but instead "
- << test_rat2 << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+ numeric test_int1(42);
+ numeric test_int2(5);
+ numeric test_rat1 = test_int1; test_rat1 /= test_int2;
+ test_rat1 = -test_rat1; // -42/5
+ numeric test_crat = test_rat1+I*test_int2; // 5*I-42/5
+ symbol a("a");
+ ex e1, e2;
+
+ if (!test_int1.is_integer()) {
+ clog << test_int1
+ << " erroneously not recognized as integer" << endl;
+ ++result;
+ }
+ if (!test_int1.is_rational()) {
+ clog << test_int1
+ << " erroneously not recognized as rational" << endl;
+ ++result;
+ }
+
+ if (!test_rat1.is_rational()) {
+ clog << test_rat1
+ << " erroneously not recognized as rational" << endl;
+ ++result;
+ }
+ if (test_rat1.is_integer()) {
+ clog << test_rat1
+ << " erroneously recognized as integer" << endl;
+ ++result;
+ }
+
+ if (!test_crat.is_crational()) {
+ clog << test_crat
+ << " erroneously not recognized as complex rational" << endl;
+ ++result;
+ }
+
+ int i = numeric(1984).to_int();
+ if (i-1984) {
+ clog << "conversion of " << i
+ << " from numeric to int failed" << endl;
+ ++result;
+ }
+
+ e1 = test_int1;
+ if (!e1.info(info_flags::posint)) {
+ clog << "expression " << e1
+ << " erroneously not recognized as positive integer" << endl;
+ ++result;
+ }
+
+ e2 = test_int1 + a;
+ if (ex_to_numeric(e2).is_integer()) {
+ clog << "expression " << e2
+ << " erroneously recognized as integer" << endl;
+ ++result;
+ }
+
+ // The next two were two actual bugs in CLN till June, 12, 1999:
+ test_rat1 = numeric(3)/numeric(2);
+ test_rat1 += test_rat1;
+ if (!test_rat1.is_integer()) {
+ clog << "3/2 + 3/2 erroneously not integer 3 but instead "
+ << test_rat1 << endl;
+ ++result;
+ }
+ test_rat1 = numeric(3)/numeric(2);
+ numeric test_rat2 = test_rat1 + numeric(1); // 5/2
+ test_rat2 -= test_rat1; // 1
+ if (!test_rat2.is_integer()) {
+ clog << "5/2 - 3/2 erroneously not integer 1 but instead "
+ << test_rat2 << endl;
+ ++result;
+ }
+
+ return result;
}
/* We had some fun with a bug in CLN that caused it to loop forever when
* the original bug in CLN was finally killed on September 2nd. */
static unsigned exam_numeric2(void)
{
- unsigned result = 0;
-
- ex zero = numeric(0);
- ex two = numeric(2);
- ex three = numeric(3);
-
- // The hang in this code was the reason for the original workaround
- if (pow(two,two/three)==42) {
- clog << "pow(2,2/3) erroneously returned 42" << endl;
- ++result; // cannot happen
- }
-
- // Actually, this used to raise a FPE after introducing the workaround
- if (two*zero!=zero) {
- clog << "2*0 erroneously returned " << two*zero << endl;
- ++result;
- }
-
- // And this returned a cl_F due to the implicit call of numeric::power()
- ex six = two*three;
- if (!six.info(info_flags::integer)) {
- clog << "2*3 erroneously returned the non-integer " << six << endl;
- ++result;
- }
-
- // The fix in the workaround left a whole which was fixed hours later...
- ex another_zero = pow(zero,numeric(1)/numeric(2));
- if (!another_zero.is_zero()) {
- clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+
+ ex zero = numeric(0);
+ ex two = numeric(2);
+ ex three = numeric(3);
+
+ // The hang in this code was the reason for the original workaround
+ if (pow(two,two/three)==42) {
+ clog << "pow(2,2/3) erroneously returned 42" << endl;
+ ++result; // cannot happen
+ }
+
+ // Actually, this used to raise a FPE after introducing the workaround
+ if (two*zero!=zero) {
+ clog << "2*0 erroneously returned " << two*zero << endl;
+ ++result;
+ }
+
+ // And this returned a cl_F due to the implicit call of numeric::power()
+ ex six = two*three;
+ if (!six.info(info_flags::integer)) {
+ clog << "2*3 erroneously returned the non-integer " << six << endl;
+ ++result;
+ }
+
+ // The fix in the workaround left a whole which was fixed hours later...
+ ex another_zero = pow(zero,numeric(1)/numeric(2));
+ if (!another_zero.is_zero()) {
+ clog << "pow(0,1/2) erroneously returned" << another_zero << endl;
+ ++result;
+ }
+
+ return result;
}
/* Assorted tests to ensure some crucial functions behave exactly as specified
* in the documentation. */
static unsigned exam_numeric3(void)
{
- unsigned result = 0;
- numeric calc_rem, calc_quo;
- numeric a, b;
-
- // check if irem(a, b) and irem(a, b, q) really behave like Maple's
- // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
- // These overloaded routines indeed need to be checked separately since
- // internally they might be doing something completely different:
- a = 23; b = 4; calc_rem = irem(a, b);
- if (calc_rem != 3) {
- clog << "irem(" << a << "," << b << ") erroneously returned "
- << calc_rem << endl;
- ++result;
- }
- a = 23; b = -4; calc_rem = irem(a, b);
- if (calc_rem != 3) {
- clog << "irem(" << a << "," << b << ") erroneously returned "
- << calc_rem << endl;
- ++result;
- }
- a = -23; b = 4; calc_rem = irem(a, b);
- if (calc_rem != -3) {
- clog << "irem(" << a << "," << b << ") erroneously returned "
- << calc_rem << endl;
- ++result;
- }
- a = -23; b = -4; calc_rem = irem(a, b);
- if (calc_rem != -3) {
- clog << "irem(" << a << "," << b << ") erroneously returned "
- << calc_rem << endl;
- ++result;
- }
- // and now the overloaded irem(a,b,q):
- a = 23; b = 4; calc_rem = irem(a, b, calc_quo);
- if (calc_rem != 3 || calc_quo != 5) {
- clog << "irem(" << a << "," << b << ",q) erroneously returned "
- << calc_rem << " with q=" << calc_quo << endl;
- ++result;
- }
- a = 23; b = -4; calc_rem = irem(a, b, calc_quo);
- if (calc_rem != 3 || calc_quo != -5) {
- clog << "irem(" << a << "," << b << ",q) erroneously returned "
- << calc_rem << " with q=" << calc_quo << endl;
- ++result;
- }
- a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
- if (calc_rem != -3 || calc_quo != -5) {
- clog << "irem(" << a << "," << b << ",q) erroneously returned "
- << calc_rem << " with q=" << calc_quo << endl;
- ++result;
- }
- a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
- if (calc_rem != -3 || calc_quo != 5) {
- clog << "irem(" << a << "," << b << ",q) erroneously returned "
- << calc_rem << " with q=" << calc_quo << endl;
- ++result;
- }
- // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's
- // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
- // These overloaded routines indeed need to be checked separately since
- // internally they might be doing something completely different:
- a = 23; b = 4; calc_quo = iquo(a, b);
- if (calc_quo != 5) {
- clog << "iquo(" << a << "," << b << ") erroneously returned "
- << calc_quo << endl;
- ++result;
- }
- a = 23; b = -4; calc_quo = iquo(a, b);
- if (calc_quo != -5) {
- clog << "iquo(" << a << "," << b << ") erroneously returned "
- << calc_quo << endl;
- ++result;
- }
- a = -23; b = 4; calc_quo = iquo(a, b);
- if (calc_quo != -5) {
- clog << "iquo(" << a << "," << b << ") erroneously returned "
- << calc_quo << endl;
- ++result;
- }
- a = -23; b = -4; calc_quo = iquo(a, b);
- if (calc_quo != 5) {
- clog << "iquo(" << a << "," << b << ") erroneously returned "
- << calc_quo << endl;
- ++result;
- }
- // and now the overloaded iquo(a,b,r):
- a = 23; b = 4; calc_quo = iquo(a, b, calc_rem);
- if (calc_quo != 5 || calc_rem != 3) {
- clog << "iquo(" << a << "," << b << ",r) erroneously returned "
- << calc_quo << " with r=" << calc_rem << endl;
- ++result;
- }
- a = 23; b = -4; calc_quo = iquo(a, b, calc_rem);
- if (calc_quo != -5 || calc_rem != 3) {
- clog << "iquo(" << a << "," << b << ",r) erroneously returned "
- << calc_quo << " with r=" << calc_rem << endl;
- ++result;
- }
- a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
- if (calc_quo != -5 || calc_rem != -3) {
- clog << "iquo(" << a << "," << b << ",r) erroneously returned "
- << calc_quo << " with r=" << calc_rem << endl;
- ++result;
- }
- a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
- if (calc_quo != 5 || calc_rem != -3) {
- clog << "iquo(" << a << "," << b << ",r) erroneously returned "
- << calc_quo << " with r=" << calc_rem << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+ numeric calc_rem, calc_quo;
+ numeric a, b;
+
+ // check if irem(a, b) and irem(a, b, q) really behave like Maple's
+ // irem(a, b) and irem(a, b, 'q') as advertised in our documentation.
+ // These overloaded routines indeed need to be checked separately since
+ // internally they might be doing something completely different:
+ a = 23; b = 4; calc_rem = irem(a, b);
+ if (calc_rem != 3) {
+ clog << "irem(" << a << "," << b << ") erroneously returned "
+ << calc_rem << endl;
+ ++result;
+ }
+ a = 23; b = -4; calc_rem = irem(a, b);
+ if (calc_rem != 3) {
+ clog << "irem(" << a << "," << b << ") erroneously returned "
+ << calc_rem << endl;
+ ++result;
+ }
+ a = -23; b = 4; calc_rem = irem(a, b);
+ if (calc_rem != -3) {
+ clog << "irem(" << a << "," << b << ") erroneously returned "
+ << calc_rem << endl;
+ ++result;
+ }
+ a = -23; b = -4; calc_rem = irem(a, b);
+ if (calc_rem != -3) {
+ clog << "irem(" << a << "," << b << ") erroneously returned "
+ << calc_rem << endl;
+ ++result;
+ }
+ // and now the overloaded irem(a,b,q):
+ a = 23; b = 4; calc_rem = irem(a, b, calc_quo);
+ if (calc_rem != 3 || calc_quo != 5) {
+ clog << "irem(" << a << "," << b << ",q) erroneously returned "
+ << calc_rem << " with q=" << calc_quo << endl;
+ ++result;
+ }
+ a = 23; b = -4; calc_rem = irem(a, b, calc_quo);
+ if (calc_rem != 3 || calc_quo != -5) {
+ clog << "irem(" << a << "," << b << ",q) erroneously returned "
+ << calc_rem << " with q=" << calc_quo << endl;
+ ++result;
+ }
+ a = -23; b = 4; calc_rem = irem(a, b, calc_quo);
+ if (calc_rem != -3 || calc_quo != -5) {
+ clog << "irem(" << a << "," << b << ",q) erroneously returned "
+ << calc_rem << " with q=" << calc_quo << endl;
+ ++result;
+ }
+ a = -23; b = -4; calc_rem = irem(a, b, calc_quo);
+ if (calc_rem != -3 || calc_quo != 5) {
+ clog << "irem(" << a << "," << b << ",q) erroneously returned "
+ << calc_rem << " with q=" << calc_quo << endl;
+ ++result;
+ }
+ // check if iquo(a, b) and iquo(a, b, r) really behave like Maple's
+ // iquo(a, b) and iquo(a, b, 'r') as advertised in our documentation.
+ // These overloaded routines indeed need to be checked separately since
+ // internally they might be doing something completely different:
+ a = 23; b = 4; calc_quo = iquo(a, b);
+ if (calc_quo != 5) {
+ clog << "iquo(" << a << "," << b << ") erroneously returned "
+ << calc_quo << endl;
+ ++result;
+ }
+ a = 23; b = -4; calc_quo = iquo(a, b);
+ if (calc_quo != -5) {
+ clog << "iquo(" << a << "," << b << ") erroneously returned "
+ << calc_quo << endl;
+ ++result;
+ }
+ a = -23; b = 4; calc_quo = iquo(a, b);
+ if (calc_quo != -5) {
+ clog << "iquo(" << a << "," << b << ") erroneously returned "
+ << calc_quo << endl;
+ ++result;
+ }
+ a = -23; b = -4; calc_quo = iquo(a, b);
+ if (calc_quo != 5) {
+ clog << "iquo(" << a << "," << b << ") erroneously returned "
+ << calc_quo << endl;
+ ++result;
+ }
+ // and now the overloaded iquo(a,b,r):
+ a = 23; b = 4; calc_quo = iquo(a, b, calc_rem);
+ if (calc_quo != 5 || calc_rem != 3) {
+ clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+ << calc_quo << " with r=" << calc_rem << endl;
+ ++result;
+ }
+ a = 23; b = -4; calc_quo = iquo(a, b, calc_rem);
+ if (calc_quo != -5 || calc_rem != 3) {
+ clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+ << calc_quo << " with r=" << calc_rem << endl;
+ ++result;
+ }
+ a = -23; b = 4; calc_quo = iquo(a, b, calc_rem);
+ if (calc_quo != -5 || calc_rem != -3) {
+ clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+ << calc_quo << " with r=" << calc_rem << endl;
+ ++result;
+ }
+ a = -23; b = -4; calc_quo = iquo(a, b, calc_rem);
+ if (calc_quo != 5 || calc_rem != -3) {
+ clog << "iquo(" << a << "," << b << ",r) erroneously returned "
+ << calc_quo << " with r=" << calc_rem << endl;
+ ++result;
+ }
+
+ return result;
}
/* Now we perform some less trivial checks about several functions which should
* return exact numbers if possible. */
static unsigned exam_numeric4(void)
{
- unsigned result = 0;
- bool passed;
-
- // square roots of squares of integers:
- passed = true;
- for (int i=0; i<42; ++i) {
- if (!sqrt(numeric(i*i)).is_integer()) {
- passed = false;
- }
- }
- if (!passed) {
- clog << "One or more square roots of squares of integers did not return exact integers" << endl;
- ++result;
- }
-
- // square roots of squares of rationals:
- passed = true;
- for (int num=0; num<41; ++num) {
- for (int den=1; den<42; ++den) {
- if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational()) {
- passed = false;
- }
- }
- }
- if (!passed) {
- clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+ bool passed;
+
+ // square roots of squares of integers:
+ passed = true;
+ for (int i=0; i<42; ++i) {
+ if (!sqrt(numeric(i*i)).is_integer()) {
+ passed = false;
+ }
+ }
+ if (!passed) {
+ clog << "One or more square roots of squares of integers did not return exact integers" << endl;
+ ++result;
+ }
+
+ // square roots of squares of rationals:
+ passed = true;
+ for (int num=0; num<41; ++num) {
+ for (int den=1; den<42; ++den) {
+ if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational()) {
+ passed = false;
+ }
+ }
+ }
+ if (!passed) {
+ clog << "One or more square roots of squares of rationals did not return exact integers" << endl;
+ ++result;
+ }
+
+ return result;
}
/* This test examines that simplifications of the form 5^(3/2) -> 5*5^(1/2)
* are carried out properly. */
static unsigned exam_numeric5(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;
}
// 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
// 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
// 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
// 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
// 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
// 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,
// 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);
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
// 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;
// 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;
}
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;
}
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;
}
-/** @file exam_pseries.cpp
+/** @File exam_pseries.cpp
*
* Series expansion test (Laurent and Taylor series). */
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;
}
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;
}
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.
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);
+ }
}
#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<symbol> a;
- for (unsigned i=0; i<size; ++i) {
- char buf[4];
- ostrstream(buf,sizeof(buf)) << i << ends;
- a.push_back(symbol(string("a")+buf));
- }
- ex e, aux;
-
- for (unsigned i=0; i<size; ++i)
- e = e + a[i];
-
- // prepare aux so it will swallow anything but a1^2:
- aux = -e + a[0] + a[1];
- e = pow(e,2).expand().subs(a[0]==aux).expand();
-
- if (e != pow(a[1],2)) {
- clog << "Denny Fliegner's quick consistency check erroneously returned "
- << e << "." << endl;
- ++result;
- }
-
- return result;
+ unsigned result = 0;
+ // create a vector of size symbols named "a0", "a1", ...
+ vector<symbol> a;
+ for (unsigned i=0; i<size; ++i) {
+ char buf[4];
+ ostrstream(buf,sizeof(buf)) << i << ends;
+ a.push_back(symbol(string("a")+buf));
+ }
+ ex e, aux;
+
+ for (unsigned i=0; i<size; ++i)
+ e = e + a[i];
+
+ // prepare aux so it will swallow anything but a1^2:
+ aux = -e + a[0] + a[1];
+ e = pow(e,2).expand().subs(a[0]==aux).expand();
+
+ if (e != pow(a[1],2)) {
+ clog << "Denny Fliegner's quick consistency check erroneously returned "
+ << e << "." << endl;
+ ++result;
+ }
+
+ return result;
}
unsigned time_dennyfliegner(void)
{
- unsigned result = 0;
-
- cout << "timing commutative expansion and substitution" << flush;
- clog << "-------commutative expansion and substitution:" << endl;
-
- vector<unsigned> sizes;
- vector<double> times;
- timer breitling;
-
- sizes.push_back(25);
- sizes.push_back(50);
- sizes.push_back(100);
- sizes.push_back(200);
-
- for (vector<unsigned>::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<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i)
- cout << '\t' << (*i);
- cout << endl << " time/s:";
- for (vector<double>::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<unsigned> sizes;
+ vector<double> times;
+ timer breitling;
+
+ sizes.push_back(25);
+ sizes.push_back(50);
+ sizes.push_back(100);
+ sizes.push_back(200);
+
+ for (vector<unsigned>::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<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i)
+ cout << '\t' << (*i);
+ cout << endl << " time/s:";
+ for (vector<double>::iterator i=times.begin(); i!=times.end(); ++i)
+ cout << '\t' << int(1000*(*i))*0.001;
+ cout << endl;
+
+ return result;
}
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<unsigned> sizes;
- vector<double> times;
- timer omega;
-
- sizes.push_back(10);
- sizes.push_back(15);
- sizes.push_back(20);
- sizes.push_back(25);
-
- for (vector<unsigned>::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<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i)
- cout << '\t' << (*i);
- cout << endl << " time/s:";
- for (vector<double>::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<unsigned> sizes;
+ vector<double> times;
+ timer omega;
+
+ sizes.push_back(10);
+ sizes.push_back(15);
+ sizes.push_back(20);
+ sizes.push_back(25);
+
+ for (vector<unsigned>::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<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i)
+ cout << '\t' << (*i);
+ cout << endl << " time/s:";
+ for (vector<double>::iterator i=times.begin(); i!=times.end(); ++i)
+ cout << '\t' << int(1000*(*i))*0.001;
+ cout << endl;
+
+ return result;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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<n; ++r)
- for (unsigned c=0; c<n; ++c)
- H.set(r,c,numeric(1,r+c+1));
- // invert it:
- cartier.start();
- matrix Hinv(n,n);
- Hinv = H.inverse();
- cout << ". passed ";
- clog << "(no output)" << endl;
- cout << int(1000*cartier.read())*0.001 << 's' << endl;
-
- // check result:
- name = (n==40?'J':(n==70?'L':'?'));
-
- cout << "timing Lewis-Wester test " << name
- << " (check rank " << n << " Hilbert)" << flush;
- clog << "-------Lewis-Wester test " << name
- << " (check rank " << n << " Hilbert)" << endl;
-
- cartier.reset();
- matrix identity = H.mul(Hinv);
- bool correct = true;
- for (unsigned r=0; r<n; ++r)
- for (unsigned c=0; c<n; ++c) {
- if (r==c) {
- if (identity(r,c)!=1)
- correct = false;
- } else {
- if (identity(r,c)!=0)
- correct = false;
- }
- }
- if (correct) {
- cout << ". passed ";
- clog << "(no output)" << endl;
- } else {
- cout << ". failed ";
- ++result;
- }
- cout << int(1000*cartier.read())*0.001 << 's' << endl;
- return result;
+ 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<n; ++r)
+ for (unsigned c=0; c<n; ++c)
+ H.set(r,c,numeric(1,r+c+1));
+ // invert it:
+ cartier.start();
+ matrix Hinv(n,n);
+ Hinv = H.inverse();
+ cout << ". passed ";
+ clog << "(no output)" << endl;
+ cout << int(1000*cartier.read())*0.001 << 's' << endl;
+
+ // check result:
+ name = (n==40?'J':(n==70?'L':'?'));
+
+ cout << "timing Lewis-Wester test " << name
+ << " (check rank " << n << " Hilbert)" << flush;
git://www.ginac.de / ginac.git / commitdiff