- Same shit as Christian did yesterday in ginac/.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Sat, 12 Aug 2000 21:59:10 +0000 (21:59 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Sat, 12 Aug 2000 21:59:10 +0000 (21:59 +0000)
42 files changed:
check/check_inifcns.cpp
check/check_lsolve.cpp
check/check_matrices.cpp
check/check_numeric.cpp
check/checks.cpp
check/checks.h
check/exam_differentiation.cpp
check/exam_inifcns.cpp
check/exam_lsolve.cpp
check/exam_matrices.cpp
check/exam_misc.cpp
check/exam_noncommut.cpp
check/exam_normalization.cpp
check/exam_numeric.cpp
check/exam_paranoia.cpp
check/exam_polygcd.cpp
check/exam_powerlaws.cpp
check/exam_pseries.cpp
check/exams.cpp
check/genex.cpp
check/time_dennyfliegner.cpp
check/time_gammaseries.cpp
check/time_lw_A.cpp
check/time_lw_B.cpp
check/time_lw_C.cpp
check/time_lw_D.cpp
check/time_lw_E.cpp
check/time_lw_F.cpp
check/time_lw_G.cpp
check/time_lw_H.cpp
check/time_lw_IJKL.cpp
check/time_lw_M1.cpp
check/time_lw_O.cpp
check/time_lw_P.cpp
check/time_lw_Pprime.cpp
check/time_lw_Q.cpp
check/time_lw_Qprime.cpp
check/time_toeplitz.cpp
check/time_vandermonde.cpp
check/timer.cpp
check/times.cpp
check/times.h

index 9cebcee..bec2ed9 100644 (file)
 /* 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;
 }
index 1873bb2..4a84a13 100644 (file)
 #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;
 }
index 62b74ba..4a768b9 100644 (file)
  * 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;
 }
index 811ba52..36a0d1c 100644 (file)
@@ -31,90 +31,90 @@ using namespace GiNaC;
  * 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;
 }
index b782738..87782fc 100644 (file)
 
 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;
 }
index 5827f22..e585eab 100644 (file)
@@ -37,8 +37,8 @@ using namespace GiNaC;
 const ex dense_univariate_poly(const symbol & x, unsigned degree);
 const ex dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree);
 const ex sparse_tree(const symbol & x, const symbol & y, const symbol & z,
-                     int level,
-                     bool trig = false, bool rational = true, bool complex = false);
+                                        int level,
+                                        bool trig = false, bool rational = true, bool complex = false);
 
 // prototypes for all individual checks should be unsigned fcn():
 unsigned check_numeric();
index 1a5aabb..13e5471 100644 (file)
 #include "exams.h"
 
 static unsigned check_diff(const ex &e, const symbol &x,
-                           const ex &d, unsigned nth=1)
+                                                  const ex &d, unsigned nth=1)
 {
-    ex ed = e.diff(x, nth);
-    if ((ed - d).compare(ex(0)) != 0) {
-        switch (nth) {
-        case 0:
-            clog << "zeroth ";
-            break;
-        case 1:
-            break;
-        case 2:
-            clog << "second ";
-            break;
-        case 3:
-            clog << "third ";
-            break;
-        default:
-            clog << nth << "th ";
-        }
-        clog << "derivative of " << e << " by " << x << " returned "
-             << ed << " instead of " << d << endl;
-        clog << "returned:" << endl;
-        ed.printtree(clog);
-        clog << endl << "instead of" << endl;
-        d.printtree(clog);
+       ex ed = e.diff(x, nth);
+       if ((ed - d).compare(ex(0)) != 0) {
+               switch (nth) {
+               case 0:
+                       clog << "zeroth ";
+                       break;
+               case 1:
+                       break;
+               case 2:
+                       clog << "second ";
+                       break;
+               case 3:
+                       clog << "third ";
+                       break;
+               default:
+                       clog << nth << "th ";
+               }
+               clog << "derivative of " << e << " by " << x << " returned "
+                        << ed << " instead of " << d << endl;
+               clog << "returned:" << endl;
+               ed.printtree(clog);
+               clog << endl << "instead of" << endl;
+               d.printtree(clog);
 
-        return 1;
-    }
-    return 0;
+               return 1;
+       }
+       return 0;
 }
 
 // Simple (expanded) polynomials
 static unsigned exam_differentiation1(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y");
-    ex e1, e2, e, d;
-    
-    // construct bivariate polynomial e to be diff'ed:
-    e1 = pow(x, -2) * 3 + pow(x, -1) * 5 + 7 + x * 11 + pow(x, 2) * 13;
-    e2 = pow(y, -2) * 5 + pow(y, -1) * 7 + 11 + y * 13 + pow(y, 2) * 17;
-    e = (e1 * e2).expand();
-    
-    // d e / dx:
-    d = ex("121-55/x^2-66/x^3-30/x^3/y^2-42/x^3/y-78/x^3*y-102/x^3*y^2-25/x^2/y^2-35/x^2/y-65/x^2*y-85/x^2*y^2+77/y+143*y+187*y^2+130*x/y^2+182/y*x+338*x*y+442*x*y^2+55/y^2+286*x",lst(x,y));
-    result += check_diff(e, x, d);
-    
-    // d e / dy:
-    d = ex("91-30/x^2/y^3-21/x^2/y^2+39/x^2+102/x^2*y-50/x/y^3-35/x/y^2+65/x+170/x*y-77*x/y^2+143*x+374*x*y-130/y^3*x^2-91/y^2*x^2+169*x^2+442*x^2*y-110/y^3*x-70/y^3+238*y-49/y^2",lst(x,y));
-    result += check_diff(e, y, d);
-    
-    // d^2 e / dx^2:
-    d = ex("286+90/x^4/y^2+126/x^4/y+234/x^4*y+306/x^4*y^2+50/x^3/y^2+70/x^3/y+130/x^3*y+170/x^3*y^2+130/y^2+182/y+338*y+442*y^2+198/x^4+110/x^3",lst(x,y));
-    result += check_diff(e, x, d, 2);
-    
-    // d^2 e / dy^2:
-    d = ex("238+90/x^2/y^4+42/x^2/y^3+102/x^2+150/x/y^4+70/x/y^3+170/x+330*x/y^4+154*x/y^3+374*x+390*x^2/y^4+182*x^2/y^3+442*x^2+210/y^4+98/y^3",lst(x,y));
-    result += check_diff(e, y, d, 2);
-    
-    return result;
+       unsigned result = 0;
+       symbol x("x"), y("y");
+       ex e1, e2, e, d;
+       
+       // construct bivariate polynomial e to be diff'ed:
+       e1 = pow(x, -2) * 3 + pow(x, -1) * 5 + 7 + x * 11 + pow(x, 2) * 13;
+       e2 = pow(y, -2) * 5 + pow(y, -1) * 7 + 11 + y * 13 + pow(y, 2) * 17;
+       e = (e1 * e2).expand();
+       
+       // d e / dx:
+       d = ex("121-55/x^2-66/x^3-30/x^3/y^2-42/x^3/y-78/x^3*y-102/x^3*y^2-25/x^2/y^2-35/x^2/y-65/x^2*y-85/x^2*y^2+77/y+143*y+187*y^2+130*x/y^2+182/y*x+338*x*y+442*x*y^2+55/y^2+286*x",lst(x,y));
+       result += check_diff(e, x, d);
+       
+       // d e / dy:
+       d = ex("91-30/x^2/y^3-21/x^2/y^2+39/x^2+102/x^2*y-50/x/y^3-35/x/y^2+65/x+170/x*y-77*x/y^2+143*x+374*x*y-130/y^3*x^2-91/y^2*x^2+169*x^2+442*x^2*y-110/y^3*x-70/y^3+238*y-49/y^2",lst(x,y));
+       result += check_diff(e, y, d);
+       
+       // d^2 e / dx^2:
+       d = ex("286+90/x^4/y^2+126/x^4/y+234/x^4*y+306/x^4*y^2+50/x^3/y^2+70/x^3/y+130/x^3*y+170/x^3*y^2+130/y^2+182/y+338*y+442*y^2+198/x^4+110/x^3",lst(x,y));
+       result += check_diff(e, x, d, 2);
+       
+       // d^2 e / dy^2:
+       d = ex("238+90/x^2/y^4+42/x^2/y^3+102/x^2+150/x/y^4+70/x/y^3+170/x+330*x/y^4+154*x/y^3+374*x+390*x^2/y^4+182*x^2/y^3+442*x^2+210/y^4+98/y^3",lst(x,y));
+       result += check_diff(e, y, d, 2);
+       
+       return result;
 }
 
 // Trigonometric functions
 static unsigned exam_differentiation2(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y"), a("a"), b("b");
-    ex e1, e2, e, d;
-    
-    // construct expression e to be diff'ed:
-    e1 = y*pow(x, 2) + a*x + b;
-    e2 = sin(e1);
-    e = b*pow(e2, 2) + y*e2 + a;
-    
-    d = 2*b*e2*cos(e1)*(2*x*y + a) + y*cos(e1)*(2*x*y + a);
-    result += check_diff(e, x, d);
-    
-    d = 2*b*pow(cos(e1),2)*pow(2*x*y + a, 2) + 4*b*y*e2*cos(e1)
-        - 2*b*pow(e2,2)*pow(2*x*y + a, 2) - y*e2*pow(2*x*y + a, 2)
-        + 2*pow(y,2)*cos(e1);
-    result += check_diff(e, x, d, 2);
-    
-    d = 2*b*e2*cos(e1)*pow(x, 2) + e2 + y*cos(e1)*pow(x, 2);
-    result += check_diff(e, y, d);
+       unsigned result = 0;
+       symbol x("x"), y("y"), a("a"), b("b");
+       ex e1, e2, e, d;
+       
+       // construct expression e to be diff'ed:
+       e1 = y*pow(x, 2) + a*x + b;
+       e2 = sin(e1);
+       e = b*pow(e2, 2) + y*e2 + a;
+       
+       d = 2*b*e2*cos(e1)*(2*x*y + a) + y*cos(e1)*(2*x*y + a);
+       result += check_diff(e, x, d);
+       
+       d = 2*b*pow(cos(e1),2)*pow(2*x*y + a, 2) + 4*b*y*e2*cos(e1)
+               - 2*b*pow(e2,2)*pow(2*x*y + a, 2) - y*e2*pow(2*x*y + a, 2)
+               + 2*pow(y,2)*cos(e1);
+       result += check_diff(e, x, d, 2);
+       
+       d = 2*b*e2*cos(e1)*pow(x, 2) + e2 + y*cos(e1)*pow(x, 2);
+       result += check_diff(e, y, d);
 
-    d = 2*b*pow(cos(e1),2)*pow(x,4) - 2*b*pow(e2,2)*pow(x,4)
-        + 2*cos(e1)*pow(x,2) - y*e2*pow(x,4);
-    result += check_diff(e, y, d, 2);
-    
-    // construct expression e to be diff'ed:
-    e2 = cos(e1);
-    e = b*pow(e2, 2) + y*e2 + a;
-    
-    d = -2*b*e2*sin(e1)*(2*x*y + a) - y*sin(e1)*(2*x*y + a);
-    result += check_diff(e, x, d);
-    
-    d = 2*b*pow(sin(e1),2)*pow(2*y*x + a,2) - 4*b*e2*sin(e1)*y 
-        - 2*b*pow(e2,2)*pow(2*y*x + a,2) - y*e2*pow(2*y*x + a,2)
-        - 2*pow(y,2)*sin(e1);
-    result += check_diff(e, x, d, 2);
-    
-    d = -2*b*e2*sin(e1)*pow(x,2) + e2 - y*sin(e1)*pow(x, 2);
-    result += check_diff(e, y, d);
-    
-    d = -2*b*pow(e2,2)*pow(x,4) + 2*b*pow(sin(e1),2)*pow(x,4)
-        - 2*sin(e1)*pow(x,2) - y*e2*pow(x,4);
-    result += check_diff(e, y, d, 2);
+       d = 2*b*pow(cos(e1),2)*pow(x,4) - 2*b*pow(e2,2)*pow(x,4)
+               + 2*cos(e1)*pow(x,2) - y*e2*pow(x,4);
+       result += check_diff(e, y, d, 2);
+       
+       // construct expression e to be diff'ed:
+       e2 = cos(e1);
+       e = b*pow(e2, 2) + y*e2 + a;
+       
+       d = -2*b*e2*sin(e1)*(2*x*y + a) - y*sin(e1)*(2*x*y + a);
+       result += check_diff(e, x, d);
+       
+       d = 2*b*pow(sin(e1),2)*pow(2*y*x + a,2) - 4*b*e2*sin(e1)*y 
+               - 2*b*pow(e2,2)*pow(2*y*x + a,2) - y*e2*pow(2*y*x + a,2)
+               - 2*pow(y,2)*sin(e1);
+       result += check_diff(e, x, d, 2);
+       
+       d = -2*b*e2*sin(e1)*pow(x,2) + e2 - y*sin(e1)*pow(x, 2);
+       result += check_diff(e, y, d);
+       
+       d = -2*b*pow(e2,2)*pow(x,4) + 2*b*pow(sin(e1),2)*pow(x,4)
+               - 2*sin(e1)*pow(x,2) - y*e2*pow(x,4);
+       result += check_diff(e, y, d, 2);
 
        return result;
 }
-    
+       
 // exp function
 static unsigned exam_differentiation3(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y"), a("a"), b("b");
-    ex e1, e2, e, d;
+       unsigned result = 0;
+       symbol x("x"), y("y"), a("a"), b("b");
+       ex e1, e2, e, d;
 
-    // construct expression e to be diff'ed:
-    e1 = y*pow(x, 2) + a*x + b;
-    e2 = exp(e1);
-    e = b*pow(e2, 2) + y*e2 + a;
-    
-    d = 2*b*pow(e2, 2)*(2*x*y + a) + y*e2*(2*x*y + a);
-    result += check_diff(e, x, d);
-    
-    d = 4*b*pow(e2,2)*pow(2*y*x + a,2) + 4*b*pow(e2,2)*y
-        + 2*pow(y,2)*e2 + y*e2*pow(2*y*x + a,2);
-    result += check_diff(e, x, d, 2);
-    
-    d = 2*b*pow(e2,2)*pow(x,2) + e2 + y*e2*pow(x,2);
-    result += check_diff(e, y, d);
-    
-    d = 4*b*pow(e2,2)*pow(x,4) + 2*e2*pow(x,2) + y*e2*pow(x,4);
-    result += check_diff(e, y, d, 2);
+       // construct expression e to be diff'ed:
+       e1 = y*pow(x, 2) + a*x + b;
+       e2 = exp(e1);
+       e = b*pow(e2, 2) + y*e2 + a;
+       
+       d = 2*b*pow(e2, 2)*(2*x*y + a) + y*e2*(2*x*y + a);
+       result += check_diff(e, x, d);
+       
+       d = 4*b*pow(e2,2)*pow(2*y*x + a,2) + 4*b*pow(e2,2)*y
+               + 2*pow(y,2)*e2 + y*e2*pow(2*y*x + a,2);
+       result += check_diff(e, x, d, 2);
+       
+       d = 2*b*pow(e2,2)*pow(x,2) + e2 + y*e2*pow(x,2);
+       result += check_diff(e, y, d);
+       
+       d = 4*b*pow(e2,2)*pow(x,4) + 2*e2*pow(x,2) + y*e2*pow(x,4);
+       result += check_diff(e, y, d, 2);
 
        return result;
 }
@@ -165,29 +165,29 @@ static unsigned exam_differentiation3(void)
 // log functions
 static unsigned exam_differentiation4(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y"), a("a"), b("b");
-    ex e1, e2, e, d;
-    
-    // construct expression e to be diff'ed:
-    e1 = y*pow(x, 2) + a*x + b;
-    e2 = log(e1);
-    e = b*pow(e2, 2) + y*e2 + a;
-    
-    d = 2*b*e2*(2*x*y + a)/e1 + y*(2*x*y + a)/e1;
-    result += check_diff(e, x, d);
-    
-    d = 2*b*pow((2*x*y + a),2)*pow(e1,-2) + 4*b*y*e2/e1
-        - 2*b*e2*pow(2*x*y + a,2)*pow(e1,-2) + 2*pow(y,2)/e1
-        - y*pow(2*x*y + a,2)*pow(e1,-2);
-    result += check_diff(e, x, d, 2);
-    
-    d = 2*b*e2*pow(x,2)/e1 + e2 + y*pow(x,2)/e1;
-    result += check_diff(e, y, d);
-    
-    d = 2*b*pow(x,4)*pow(e1,-2) - 2*b*e2*pow(e1,-2)*pow(x,4)
-        + 2*pow(x,2)/e1 - y*pow(x,4)*pow(e1,-2);
-    result += check_diff(e, y, d, 2);
+       unsigned result = 0;
+       symbol x("x"), y("y"), a("a"), b("b");
+       ex e1, e2, e, d;
+       
+       // construct expression e to be diff'ed:
+       e1 = y*pow(x, 2) + a*x + b;
+       e2 = log(e1);
+       e = b*pow(e2, 2) + y*e2 + a;
+       
+       d = 2*b*e2*(2*x*y + a)/e1 + y*(2*x*y + a)/e1;
+       result += check_diff(e, x, d);
+       
+       d = 2*b*pow((2*x*y + a),2)*pow(e1,-2) + 4*b*y*e2/e1
+               - 2*b*e2*pow(2*x*y + a,2)*pow(e1,-2) + 2*pow(y,2)/e1
+               - y*pow(2*x*y + a,2)*pow(e1,-2);
+       result += check_diff(e, x, d, 2);
+       
+       d = 2*b*e2*pow(x,2)/e1 + e2 + y*pow(x,2)/e1;
+       result += check_diff(e, y, d);
+       
+       d = 2*b*pow(x,4)*pow(e1,-2) - 2*b*e2*pow(e1,-2)*pow(x,4)
+               + 2*pow(x,2)/e1 - y*pow(x,4)*pow(e1,-2);
+       result += check_diff(e, y, d, 2);
 
        return result;
 }
@@ -195,102 +195,102 @@ static unsigned exam_differentiation4(void)
 // Functions with two variables
 static unsigned exam_differentiation5(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y"), a("a"), b("b");
-    ex e1, e2, e, d;
-    
-    // test atan2
-    e1 = y*pow(x, 2) + a*x + b;
-    e2 = x*pow(y, 2) + b*y + a;
-    e = atan2(e1,e2);
-    /*
-    d = pow(y,2)*(-b-y*pow(x,2)-a*x)/(pow(b+y*pow(x,2)+a*x,2)+pow(x*pow(y,2)+b*y+a,2))
-        +(2*y*x+a)/((x*pow(y,2)+b*y+a)*(1+pow(b*y*pow(x,2)+a*x,2)/pow(x*pow(y,2)+b*y+a,2)));
-    */
-    /*
-    d = ((a+2*y*x)*pow(y*b+pow(y,2)*x+a,-1)-(a*x+b+y*pow(x,2))*
-         pow(y*b+pow(y,2)*x+a,-2)*pow(y,2))*
-        pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1);
-    */
-    /*
-    d = pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1)
-        *pow(y*b+pow(y,2)*x+a,-1)*(a+2*y*x)
-        +pow(y,2)*(-a*x-b-y*pow(x,2))*
-        pow(pow(y*b+pow(y,2)*x+a,2)+pow(a*x+b+y*pow(x,2),2),-1);
-    */
-    d = pow(y,2)*pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)*
-        (-b-y*pow(x,2)-x*a)+
-        pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)*
-        (y*b+pow(y,2)*x+a)*(2*y*x+a);
-    result += check_diff(e, x, d);
-    
-    return result;
+       unsigned result = 0;
+       symbol x("x"), y("y"), a("a"), b("b");
+       ex e1, e2, e, d;
+       
+       // test atan2
+       e1 = y*pow(x, 2) + a*x + b;
+       e2 = x*pow(y, 2) + b*y + a;
+       e = atan2(e1,e2);
+       /*
+       d = pow(y,2)*(-b-y*pow(x,2)-a*x)/(pow(b+y*pow(x,2)+a*x,2)+pow(x*pow(y,2)+b*y+a,2))
+               +(2*y*x+a)/((x*pow(y,2)+b*y+a)*(1+pow(b*y*pow(x,2)+a*x,2)/pow(x*pow(y,2)+b*y+a,2)));
+       */
+       /*
+       d = ((a+2*y*x)*pow(y*b+pow(y,2)*x+a,-1)-(a*x+b+y*pow(x,2))*
+                pow(y*b+pow(y,2)*x+a,-2)*pow(y,2))*
+               pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1);
+       */
+       /*
+       d = pow(1+pow(a*x+b+y*pow(x,2),2)*pow(y*b+pow(y,2)*x+a,-2),-1)
+               *pow(y*b+pow(y,2)*x+a,-1)*(a+2*y*x)
+               +pow(y,2)*(-a*x-b-y*pow(x,2))*
+               pow(pow(y*b+pow(y,2)*x+a,2)+pow(a*x+b+y*pow(x,2),2),-1);
+       */
+       d = pow(y,2)*pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)*
+               (-b-y*pow(x,2)-x*a)+
+               pow(pow(b+y*pow(x,2)+x*a,2)+pow(y*b+pow(y,2)*x+a,2),-1)*
+               (y*b+pow(y,2)*x+a)*(2*y*x+a);
+       result += check_diff(e, x, d);
+       
+       return result;
 }
 
 // Series
 static unsigned exam_differentiation6(void)
 {
-    symbol x("x");
-    ex e, d, ed;
-    
-    e = sin(x).series(x==0, 8);
-    d = cos(x).series(x==0, 7);
-    ed = e.diff(x);
-    ed = series_to_poly(ed);
-    d = series_to_poly(d);
-    
-    if ((ed - d).compare(ex(0)) != 0) {
-        clog << "derivative of " << e << " by " << x << " returned "
-             << ed << " instead of " << d << ")" << endl;
-        return 1;
-    }
-    return 0;
+       symbol x("x");
+       ex e, d, ed;
+       
+       e = sin(x).series(x==0, 8);
+       d = cos(x).series(x==0, 7);
+       ed = e.diff(x);
+       ed = series_to_poly(ed);
+       d = series_to_poly(d);
+       
+       if ((ed - d).compare(ex(0)) != 0) {
+               clog << "derivative of " << e << " by " << x << " returned "
+                        << ed << " instead of " << d << ")" << endl;
+               return 1;
+       }
+       return 0;
 }
 
 // Hashing can help a lot, if differentiation is done cleverly
 static unsigned exam_differentiation7(void)
 {
-    symbol x("x");
-    ex P = x + pow(x,3);
-    ex e = (P.diff(x) / P).diff(x, 2);
-    ex d = 6/P - 18*x/pow(P,2) - 54*pow(x,3)/pow(P,2) + 2/pow(P,3)
-        +18*pow(x,2)/pow(P,3) + 54*pow(x,4)/pow(P,3) + 54*pow(x,6)/pow(P,3);
-    
-    if (!(e-d).expand().is_zero()) {
-        clog << "expanded second derivative of " << (P.diff(x) / P) << " by " << x
-             << " returned " << e.expand() << " instead of " << d << endl;
-        return 1;
-    }
-    if (e.nops() > 3) {
-        clog << "second derivative of " << (P.diff(x) / P) << " by " << x
-             << " has " << e.nops() << " operands.  "
-             << "The result is still correct but not optimal: 3 are enough!  "
-             << "(Hint: maybe the product rule for objects of class mul should be more careful about assembling the result?)" << endl;
-        return 1;
-    }
-    return 0;
+       symbol x("x");
+       ex P = x + pow(x,3);
+       ex e = (P.diff(x) / P).diff(x, 2);
+       ex d = 6/P - 18*x/pow(P,2) - 54*pow(x,3)/pow(P,2) + 2/pow(P,3)
+               +18*pow(x,2)/pow(P,3) + 54*pow(x,4)/pow(P,3) + 54*pow(x,6)/pow(P,3);
+       
+       if (!(e-d).expand().is_zero()) {
+               clog << "expanded second derivative of " << (P.diff(x) / P) << " by " << x
+                        << " returned " << e.expand() << " instead of " << d << endl;
+               return 1;
+       }
+       if (e.nops() > 3) {
+               clog << "second derivative of " << (P.diff(x) / P) << " by " << x
+                        << " has " << e.nops() << " operands.  "
+                        << "The result is still correct but not optimal: 3 are enough!  "
+                        << "(Hint: maybe the product rule for objects of class mul should be more careful about assembling the result?)" << endl;
+               return 1;
+       }
+       return 0;
 }
 
 unsigned exam_differentiation(void)
 {
-    unsigned result = 0;
-    
-    cout << "examining symbolic differentiation" << flush;
-    clog << "----------symbolic differentiation:" << endl;
-    
-    result += exam_differentiation1();  cout << '.' << flush;
-    result += exam_differentiation2();  cout << '.' << flush;
-    result += exam_differentiation3();  cout << '.' << flush;
-    result += exam_differentiation4();  cout << '.' << flush;
-    result += exam_differentiation5();  cout << '.' << flush;
-    result += exam_differentiation6();  cout << '.' << flush;
-    result += exam_differentiation7();  cout << '.' << flush;
-    
-    if (!result) {
-        cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
-        cout << " failed " << endl;
-    }
-    return result;
+       unsigned result = 0;
+       
+       cout << "examining symbolic differentiation" << flush;
+       clog << "----------symbolic differentiation:" << endl;
+       
+       result += exam_differentiation1();  cout << '.' << flush;
+       result += exam_differentiation2();  cout << '.' << flush;
+       result += exam_differentiation3();  cout << '.' << flush;
+       result += exam_differentiation4();  cout << '.' << flush;
+       result += exam_differentiation5();  cout << '.' << flush;
+       result += exam_differentiation6();  cout << '.' << flush;
+       result += exam_differentiation7();  cout << '.' << flush;
+       
+       if (!result) {
+               cout << " passed " << endl;
+               clog << "(no output)" << endl;
+       } else {
+               cout << " failed " << endl;
+       }
+       return result;
 }
index 823e25e..baeacfd 100644 (file)
 /* Assorted tests on other transcendental functions. */
 static unsigned inifcns_consist_trans(void)
 {
-    unsigned result = 0;
-    symbol x("x");
-    ex chk;
-    
-    chk = asin(1)-acos(0);
-    if (!chk.is_zero()) {
-        clog << "asin(1)-acos(0) erroneously returned " << chk
-             << " instead of 0" << endl;
-        ++result;
-    }
-    
-    // arbitrary check of type sin(f(x)):
-    chk = pow(sin(acos(x)),2) + pow(sin(asin(x)),2)
-        - (1+pow(x,2))*pow(sin(atan(x)),2);
-    if (chk != 1-pow(x,2)) {
-        clog << "sin(acos(x))^2 + sin(asin(x))^2 - (1+x^2)*sin(atan(x))^2 "
-             << "erroneously returned " << chk << " instead of 1-x^2" << endl;
-        ++result;
-    }
-    
-    // arbitrary check of type cos(f(x)):
-    chk = pow(cos(acos(x)),2) + pow(cos(asin(x)),2)
-        - (1+pow(x,2))*pow(cos(atan(x)),2);
-    if (!chk.is_zero()) {
-        clog << "cos(acos(x))^2 + cos(asin(x))^2 - (1+x^2)*cos(atan(x))^2 "
-             << "erroneously returned " << chk << " instead of 0" << endl;
-        ++result;
-    }
-    
-    // arbitrary check of type tan(f(x)):
-    chk = tan(acos(x))*tan(asin(x)) - tan(atan(x));
-    if (chk != 1-x) {
-        clog << "tan(acos(x))*tan(asin(x)) - tan(atan(x)) "
-             << "erroneously returned " << chk << " instead of -x+1" << endl;
-        ++result;
-    }
-    
-    // arbitrary check of type sinh(f(x)):
-    chk = -pow(sinh(acosh(x)),2).expand()*pow(sinh(atanh(x)),2)
-        - pow(sinh(asinh(x)),2);
-    if (!chk.is_zero()) {
-        clog << "expand(-(sinh(acosh(x)))^2)*(sinh(atanh(x))^2) - sinh(asinh(x))^2 "
-             << "erroneously returned " << chk << " instead of 0" << endl;
-        ++result;
-    }
-    
-    // arbitrary check of type cosh(f(x)):
-    chk = (pow(cosh(asinh(x)),2) - 2*pow(cosh(acosh(x)),2))
-        * pow(cosh(atanh(x)),2);
-    if (chk != 1) {
-        clog << "(cosh(asinh(x))^2 - 2*cosh(acosh(x))^2) * cosh(atanh(x))^2 "
-             << "erroneously returned " << chk << " instead of 1" << endl;
-        ++result;
-    }
-    
-    // arbitrary check of type tanh(f(x)):
-    chk = (pow(tanh(asinh(x)),-2) - pow(tanh(acosh(x)),2)).expand()
-        * pow(tanh(atanh(x)),2);
-    if (chk != 2) {
-        clog << "expand(tanh(acosh(x))^2 - tanh(asinh(x))^(-2)) * tanh(atanh(x))^2 "
-             << "erroneously returned " << chk << " instead of 2" << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       symbol x("x");
+       ex chk;
+       
+       chk = asin(1)-acos(0);
+       if (!chk.is_zero()) {
+               clog << "asin(1)-acos(0) erroneously returned " << chk
+                        << " instead of 0" << endl;
+               ++result;
+       }
+       
+       // arbitrary check of type sin(f(x)):
+       chk = pow(sin(acos(x)),2) + pow(sin(asin(x)),2)
+               - (1+pow(x,2))*pow(sin(atan(x)),2);
+       if (chk != 1-pow(x,2)) {
+               clog << "sin(acos(x))^2 + sin(asin(x))^2 - (1+x^2)*sin(atan(x))^2 "
+                        << "erroneously returned " << chk << " instead of 1-x^2" << endl;
+               ++result;
+       }
+       
+       // arbitrary check of type cos(f(x)):
+       chk = pow(cos(acos(x)),2) + pow(cos(asin(x)),2)
+               - (1+pow(x,2))*pow(cos(atan(x)),2);
+       if (!chk.is_zero()) {
+               clog << "cos(acos(x))^2 + cos(asin(x))^2 - (1+x^2)*cos(atan(x))^2 "
+                        << "erroneously returned " << chk << " instead of 0" << endl;
+               ++result;
+       }
+       
+       // arbitrary check of type tan(f(x)):
+       chk = tan(acos(x))*tan(asin(x)) - tan(atan(x));
+       if (chk != 1-x) {
+               clog << "tan(acos(x))*tan(asin(x)) - tan(atan(x)) "
+                        << "erroneously returned " << chk << " instead of -x+1" << endl;
+               ++result;
+       }
+       
+       // arbitrary check of type sinh(f(x)):
+       chk = -pow(sinh(acosh(x)),2).expand()*pow(sinh(atanh(x)),2)
+               - pow(sinh(asinh(x)),2);
+       if (!chk.is_zero()) {
+               clog << "expand(-(sinh(acosh(x)))^2)*(sinh(atanh(x))^2) - sinh(asinh(x))^2 "
+                        << "erroneously returned " << chk << " instead of 0" << endl;
+               ++result;
+       }
+       
+       // arbitrary check of type cosh(f(x)):
+       chk = (pow(cosh(asinh(x)),2) - 2*pow(cosh(acosh(x)),2))
+               * pow(cosh(atanh(x)),2);
+       if (chk != 1) {
+               clog << "(cosh(asinh(x))^2 - 2*cosh(acosh(x))^2) * cosh(atanh(x))^2 "
+                        << "erroneously returned " << chk << " instead of 1" << endl;
+               ++result;
+       }
+       
+       // arbitrary check of type tanh(f(x)):
+       chk = (pow(tanh(asinh(x)),-2) - pow(tanh(acosh(x)),2)).expand()
+               * pow(tanh(atanh(x)),2);
+       if (chk != 2) {
+               clog << "expand(tanh(acosh(x))^2 - tanh(asinh(x))^(-2)) * tanh(atanh(x))^2 "
+                        << "erroneously returned " << chk << " instead of 2" << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 /* Simple tests on the tgamma function.  We stuff in arguments where the results
  * exists in closed form and check if it's ok. */
 static unsigned inifcns_consist_gamma(void)
 {
-    unsigned result = 0;
-    ex e;
-    
-    e = tgamma(ex(1));
-    for (int i=2; i<8; ++i)
-        e += tgamma(ex(i));
-    if (e != numeric(874)) {
-        clog << "tgamma(1)+...+tgamma(7) erroneously returned "
-             << e << " instead of 874" << endl;
-        ++result;
-    }
-    
-    e = tgamma(ex(1));
-    for (int i=2; i<8; ++i)
-        e *= tgamma(ex(i));    
-    if (e != numeric(24883200)) {
-        clog << "tgamma(1)*...*tgamma(7) erroneously returned "
-             << e << " instead of 24883200" << endl;
-        ++result;
-    }
-    
-    e = tgamma(ex(numeric(5, 2)))*tgamma(ex(numeric(9, 2)))*64;
-    if (e != 315*Pi) {
-        clog << "64*tgamma(5/2)*tgamma(9/2) erroneously returned "
-             << e << " instead of 315*Pi" << endl;
-        ++result;
-    }
-    
-    e = tgamma(ex(numeric(-13, 2)));
-    for (int i=-13; i<7; i=i+2)
-        e += tgamma(ex(numeric(i, 2)));
-    e = (e*tgamma(ex(numeric(15, 2)))*numeric(512));
-    if (e != numeric(633935)*Pi) {
-        clog << "512*(tgamma(-13/2)+...+tgamma(5/2))*tgamma(15/2) erroneously returned "
-             << e << " instead of 633935*Pi" << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       ex e;
+       
+       e = tgamma(ex(1));
+       for (int i=2; i<8; ++i)
+               e += tgamma(ex(i));
+       if (e != numeric(874)) {
+               clog << "tgamma(1)+...+tgamma(7) erroneously returned "
+                        << e << " instead of 874" << endl;
+               ++result;
+       }
+       
+       e = tgamma(ex(1));
+       for (int i=2; i<8; ++i)
+               e *= tgamma(ex(i));     
+       if (e != numeric(24883200)) {
+               clog << "tgamma(1)*...*tgamma(7) erroneously returned "
+                        << e << " instead of 24883200" << endl;
+               ++result;
+       }
+       
+       e = tgamma(ex(numeric(5, 2)))*tgamma(ex(numeric(9, 2)))*64;
+       if (e != 315*Pi) {
+               clog << "64*tgamma(5/2)*tgamma(9/2) erroneously returned "
+                        << e << " instead of 315*Pi" << endl;
+               ++result;
+       }
+       
+       e = tgamma(ex(numeric(-13, 2)));
+       for (int i=-13; i<7; i=i+2)
+               e += tgamma(ex(numeric(i, 2)));
+       e = (e*tgamma(ex(numeric(15, 2)))*numeric(512));
+       if (e != numeric(633935)*Pi) {
+               clog << "512*(tgamma(-13/2)+...+tgamma(5/2))*tgamma(15/2) erroneously returned "
+                        << e << " instead of 633935*Pi" << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 /* Simple tests on the Psi-function (aka polygamma-function).  We stuff in
    arguments where the result exists in closed form and check if it's ok. */
 static unsigned inifcns_consist_psi(void)
 {
-    unsigned result = 0;
-    symbol x;
-    ex e, f;
-    
-    // We check psi(1) and psi(1/2) implicitly by calculating the curious
-    // little identity tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) == 2*log(2).
-    e += (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1));
-    e -= (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1,2));
-    if (e!=2*log(2)) {
-        clog << "tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) erroneously returned "
-             << e << " instead of 2*log(2)" << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       symbol x;
+       ex e, f;
+       
+       // We check psi(1) and psi(1/2) implicitly by calculating the curious
+       // little identity tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) == 2*log(2).
+       e += (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1));
+       e -= (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1,2));
+       if (e!=2*log(2)) {
+               clog << "tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) erroneously returned "
+                        << e << " instead of 2*log(2)" << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 /* Simple tests on the Riemann Zeta function.  We stuff in arguments where the
@@ -164,47 +164,47 @@ static unsigned inifcns_consist_psi(void)
  * the Bernoulli numbers as a side effect. */
 static unsigned inifcns_consist_zeta(void)
 {
-    unsigned result = 0;
-    ex e;
-    
-    for (int i=0; i<13; i+=2)
-        e += zeta(i)/pow(Pi,i);
-    if (e!=numeric(-204992279,638512875)) {
-        clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned "
-             << e << " instead of -204992279/638512875" << endl;
-        ++result;
-    }
-    
-    e = 0;
-    for (int i=-1; i>-16; i--)
-        e += zeta(i);
-    if (e!=numeric(487871,1633632)) {
-        clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned "
-             << e << " instead of 487871/1633632" << endl;
-        ++result;
-    }
-    
-    return result;
+       unsigned result = 0;
+       ex e;
+       
+       for (int i=0; i<13; i+=2)
+               e += zeta(i)/pow(Pi,i);
+       if (e!=numeric(-204992279,638512875)) {
+               clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned "
+                        << e << " instead of -204992279/638512875" << endl;
+               ++result;
+       }
+       
+       e = 0;
+       for (int i=-1; i>-16; i--)
+               e += zeta(i);
+       if (e!=numeric(487871,1633632)) {
+               clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned "
+                        << e << " instead of 487871/1633632" << endl;
+               ++result;
+       }
+       
+       return result;
 }
 
 unsigned exam_inifcns(void)
 {
-    unsigned result = 0;
-    
-    cout << "examining consistency of symbolic functions" << flush;
-    clog << "----------consistency of symbolic functions:" << endl;
-    
-    result += inifcns_consist_trans();  cout << '.' << flush;
-    result += inifcns_consist_gamma();  cout << '.' << flush;
-    result += inifcns_consist_psi();  cout << '.' << flush;
-    result += inifcns_consist_zeta();  cout << '.' << flush;
-    
-    if (!result) {
-        cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
-        cout << " failed " << endl;
-    }
-    
-    return result;
+       unsigned result = 0;
+       
+       cout << "examining consistency of symbolic functions" << flush;
+       clog << "----------consistency of symbolic functions:" << endl;
+       
+       result += inifcns_consist_trans();  cout << '.' << flush;
+       result += inifcns_consist_gamma();  cout << '.' << flush;
+       result += inifcns_consist_psi();  cout << '.' << flush;
+       result += inifcns_consist_zeta();  cout << '.' << flush;
+       
+       if (!result) {
+               cout << " passed " << endl;
+               clog << "(no output)" << endl;
+       } else {
+               cout << " failed " << endl;
+       }
+       
+       return result;
 }
index 04eb438..cafda13 100644 (file)
 
 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;
 }
index 9f17bb3..a032bfe 100644 (file)
 
 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;
 }
index 86bd12a..8052bda 100644 (file)
@@ -1,6 +1,7 @@
 /** @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:
@@ -55,38 +56,38 @@ static unsigned exam_expand_subs(void)
  *  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;
 }
index 634fad7..3cf251b 100644 (file)
 
 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;
 }
index 98e1d7c..50d717d 100644 (file)
@@ -26,155 +26,155 @@ static symbol w("w"), x("x"), y("y"), z("z");
 
 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;
 }
index e37c565..04cf354 100644 (file)
  * 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
@@ -112,233 +112,233 @@ static unsigned exam_numeric1(void)
  * 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;
 }
index 2c6fe6a..c4383c7 100644 (file)
 // actually introduced the second.)
 static unsigned exam_paranoia1(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y"), z("z");
-    ex e, f, g;
-
-    e = x * y * z;
-    f = y * z;
-    g = e / f;
-
-    // In the first one expand did not do any job at all:
-    if (!g.expand().is_equal(x)) {
-        clog << "e = x*y*z; f = y*z; expand(e/f) erroneously returned "
-             << g.expand() << endl;
-        ++result;
-    }
-
-    // This one somehow used to return 0:
-    e = pow(x + 1, -1);
-    if (!e.expand().is_equal(e)) {
-        clog << "expand(pow(x + 1, -1)) erroneously returned "
-             << e.expand() << endl;
-        ++result;
-    }
-
-    return result;
+       unsigned result = 0;
+       symbol x("x"), y("y"), z("z");
+       ex e, f, g;
+
+       e = x * y * z;
+       f = y * z;
+       g = e / f;
+
+       // In the first one expand did not do any job at all:
+       if (!g.expand().is_equal(x)) {
+               clog << "e = x*y*z; f = y*z; expand(e/f) erroneously returned "
+                        << g.expand() << endl;
+               ++result;
+       }
+
+       // This one somehow used to return 0:
+       e = pow(x + 1, -1);
+       if (!e.expand().is_equal(e)) {
+               clog << "expand(pow(x + 1, -1)) erroneously returned "
+                        << e.expand() << endl;
+               ++result;
+       }
+
+       return result;
 }
 
 // And here the second oops which showed up until May 17th 1999.  It had to do
@@ -61,35 +61,35 @@ static unsigned exam_paranoia1(void)
 // had the names as given here:
 static unsigned exam_paranoia2(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y"), z("z");
-    ex e, f, g;
-
-    e = x + z*x;
-    f = e*y;
-    g = f - e*y;
-
-    // After .expand(), g should be zero:
-    if (!g.expand().is_zero()) {
-        clog << "e = (x + z*x); f = e*y; expand(f - e*y) erroneously returned "
-             << g.expand() << endl;
-        ++result;
-    }
-    // After .eval(), g should be zero:
-    if (!g.eval().is_zero()) {
-        clog << "e = (x + z*x); f = e*y; eval(f - e*y) erroneously returned "
-             << g.eval() << endl;
-        ++result;
-    }
-    // This actually worked already back in April 1999.
-    // But we are *very* paranoic!
-    if (!g.expand().eval().is_zero()) {
-        clog << "e = (x + z*x); f = e*y; eval(expand(f - e*y)) erroneously returned "
-             << g.expand().eval() << endl;
-        ++result;
-    }
-
-    return result;
+       unsigned result = 0;
+       symbol x("x"), y("y"), z("z");
+       ex e, f, g;
+
+       e = x + z*x;
+       f = e*y;
+       g = f - e*y;
+
+       // After .expand(), g should be zero:
+       if (!g.expand().is_zero()) {
+               clog << "e = (x + z*x); f = e*y; expand(f - e*y) erroneously returned "
+                        << g.expand() << endl;
+               ++result;
+       }
+       // After .eval(), g should be zero:
+       if (!g.eval().is_zero()) {
+               clog << "e = (x + z*x); f = e*y; eval(f - e*y) erroneously returned "
+                        << g.eval() << endl;
+               ++result;
+       }
+       // This actually worked already back in April 1999.
+       // But we are *very* paranoic!
+       if (!g.expand().eval().is_zero()) {
+               clog << "e = (x + z*x); f = e*y; eval(expand(f - e*y)) erroneously returned "
+                        << g.expand().eval() << endl;
+               ++result;
+       }
+
+       return result;
 }
 
 // The third bug was introduced on May 18th 1999, discovered on May 19 and
@@ -97,110 +97,110 @@ static unsigned exam_paranoia2(void)
 // other numbers:
 static unsigned exam_paranoia3(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y");
-    ex e, f;
-
-    e = x*y - y;
-    f = e.subs(x == 2);
-
-    if (!f.is_equal(y)) {
-        clog << "e = x*y - y; f = e.subs(x == 2) erroneously returned "
-             << f << endl;
-        ++result;
-    }
-    if (!f.eval().is_equal(y)) {
-        clog << "e = x*y - y; eval(e.subs(x == 2)) erroneously returned "
-             << f.eval() << endl;
-        ++result;
-    }
-    if (!f.expand().is_equal(y)) {
-        clog << "e = x*y - y; expand(e.subs(x == 2)) erroneously returned "
-             << f.expand() << endl;
-        ++result;
-    }
-
-    return result;
+       unsigned result = 0;
+       symbol x("x"), y("y");
+       ex e, f;
+
+       e = x*y - y;
+       f = e.subs(x == 2);
+
+       if (!f.is_equal(y)) {
+               clog << "e = x*y - y; f = e.subs(x == 2) erroneously returned "
+                        << f << endl;
+               ++result;
+       }
+       if (!f.eval().is_equal(y)) {
+               clog << "e = x*y - y; eval(e.subs(x == 2)) erroneously returned "
+                        << f.eval() << endl;
+               ++result;
+       }
+       if (!f.expand().is_equal(y)) {
+               clog << "e = x*y - y; expand(e.subs(x == 2)) erroneously returned "
+                        << f.expand() << endl;
+               ++result;
+       }
+
+       return result;
 }
 
 // The fourth bug was also discovered on May 19th 1999 and fixed immediately:
 static unsigned exam_paranoia4(void)
 {
-    unsigned result = 0;
-    symbol x("x");
-    ex e, f, g;
-
-    e = pow(x, 2) + x + 1;
-    f = pow(x, 2) + x + 1;
-    g = e - f;
-
-    if (!g.is_zero()) {
-        clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned "
-             << g << endl;
-        ++result;
-    }
-    if (!g.is_zero()) {
-        clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g.eval() erroneously returned "
-             << g.eval() << endl;
-        ++result;
-    }
-
-    return result;
+       unsigned result = 0;
+       symbol x("x");
+       ex e, f, g;
+
+       e = pow(x, 2) + x + 1;
+       f = pow(x, 2) + x + 1;
+       g = e - f;
+
+       if (!g.is_zero()) {
+               clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned "
+                        << g << endl;
+               ++result;
+       }
+       if (!g.is_zero()) {
+               clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g.eval() erroneously returned "
+                        << g.eval() << endl;
+               ++result;
+       }
+
+       return result;
 }
 
 // The fifth oops was discovered on May 20th 1999 and fixed a day later:
 static unsigned exam_paranoia5(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y");
+       unsigned result = 0;
+       symbol x("x"), y("y");
 
-    ex e, f;
-    e = pow(x*y + 1, 2);
-    f = pow(x, 2) * pow(y, 2) + 2*x*y + 1;
+       ex e, f;
+       e = pow(x*y + 1, 2);
+       f = pow(x, 2) * pow(y, 2) + 2*x*y + 1;
 
-    if (!(e-f).expand().is_zero()) {
-        clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned "
-             << (e-f).expand() << endl;
-        ++result;
-    }
+       if (!(e-f).expand().is_zero()) {
+               clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned "
+                        << (e-f).expand() << endl;
+               ++result;
+       }
 
-    return result;
+       return result;
 }
 
 // This one was discovered on Jun 1st 1999 and fixed the same day:
 static unsigned exam_paranoia6(void)
 {
-    unsigned result = 0;
-    symbol x("x");
-
-    ex e, f;
-    e = pow(x, -5);
-    f = e.denom();
-
-    if (!f.is_equal(pow(x, 5))) {
-        clog << "e = pow(x, -5); f = e.denom(); f was " << f << " (should be x^5)" << endl;
-        ++result;
-    }
-    return result;
+       unsigned result = 0;
+       symbol x("x");
+
+       ex e, f;
+       e = pow(x, -5);
+       f = e.denom();
+
+       if (!f.is_equal(pow(x, 5))) {
+               clog << "e = pow(x, -5); f = e.denom(); f was " << f << " (should be x^5)" << endl;
+               ++result;
+       }
+       return result;
 }
 
 // This one was introduced on June 1st 1999 by some aggressive manual
 // optimization. Discovered and fixed on June 2nd.
 static unsigned exam_paranoia7(void)
 {
-    unsigned result = 0;
-    symbol x("x"), y("y");
-
-    ex e = y + y*x + 2;
-    ex f = expand(pow(e, 2) - (e*y*(x + 1)));
-
-    if (f.nops() > 3) {
-        clog << "e=y+y*x+2; f=expand(pow(e,2)-(e*y*(x+1))) has "
-             << f.nops() << " arguments instead of 3 ( f=="
-             << f << " )" << endl;
-        ++result;
-    }
-    return result;
+       unsigned result = 0;
+       symbol x("x"), y("y");
+
+       ex e = y + y*x + 2;
+       ex f = expand(pow(e, 2) - (e*y*(x + 1)));
+
+       if (f.nops() > 3) {
+               clog << "e=y+y*x+2; f=expand(pow(e,2)-(e*y*(x+1))) has "
+                        << f.nops() << " arguments instead of 3 ( f=="
+                        << f << " )" << endl;
+               ++result;
+       }
+       return result;
 }
 
 // This one was a result of the rewrite of mul::max_coefficient when we
@@ -208,23 +208,23 @@ static unsigned exam_paranoia7(void)
 // 1999. Fixed on Oct 4th.
 static unsigned exam_paranoia8(void)
 {
-    unsigned result = 0;
-    symbol x("x");
-
-    ex e = -x / (x+1);
-    ex f;
-    
-    try {
-        f = e.normal();
-        if (!f.is_equal(e)) {
-            clog << "normal(-x/(x+1)) returns " << f << " instead of -x/(x+1)\n";
-            ++result;
-        }
-    } catch (const exception &err) {
-        clog << "normal(-x/(x+1) throws " << err.what() << endl;
-        ++result;
-    }
-    return result;
+       unsigned result = 0;
+       symbol x("x");
+
+       ex e = -x / (x+1);
+       ex f;
+       
+       try {
+               f = e.normal();
+               if (!f.is_equal(e)) {
+                       clog << "normal(-x/(x+1)) returns " << f << " instead of -x/(x+1)\n";
+                       ++result;
+               }
+       } catch (const exception &err) {
+               clog << "normal(-x/(x+1) throws " << err.what() << endl;
+               ++result;
+       }
+       return result;
 }
 
 // This one was a result of a modification to frac_cancel() & Co. to avoid
@@ -234,17 +234,17 @@ static unsigned exam_paranoia8(void)
 // 2000 and fixed on Jan 31th.
 static unsigned exam_paranoia9(void)
 {
-    unsigned result = 0;
-    symbol x("x");
+       unsigned result = 0;
+       symbol x("x");
 
-    ex e = (exp(-x)-2*x*exp(-x)+pow(x,2)/2*exp(-x))/exp(-x);
-    ex f = e.normal();
+       ex e = (exp(-x)-2*x*exp(-x)+pow(x,2)/2*exp(-x))/exp(-x);
+       ex f = e.normal();
 
-    if (!f.is_equal(1-2*x+pow(x,2)/2)) {
-        clog << "normal(" << e << ") returns " << f << " instead of 1-2*x+1/2*x^2\n";
-        ++result;
-    }
-    return result;
+       if (!f.is_equal(1-2*x+pow(x,2)/2)) {
+               clog << "normal(" << e << ") returns " << f << " instead of 1-2*x+1/2*x^2\n";
+               ++result;
+       }
+       return result;
 }
 
 // I have no idea when this broke.  It has been working long ago, before 0.4.0
@@ -253,23 +253,23 @@ static unsigned exam_paranoia9(void)
 // It was fixed that same day.
 static unsigned exam_paranoia10(void)
 {
-    unsigned result = 0;
-    
-    ex b = numeric(2);
-    ex e = numeric(3,2);
-    ex r;
-    
-    try {
-        r = pow(b,e).eval();
-        if (!(r-2*sqrt(ex(2))).is_zero()) {
-            clog << "2^(3/2) erroneously returned " << r << " instead of 2*sqrt(2)" << endl;
-            ++result;
-        }
-    } catch (const exception &err) {
-        clog << "2^(3/2) throws " << err.what() << endl;
-        ++result;
-    }
-    return result;
+       unsigned result = 0;
+       
+       ex b = numeric(2);
+       ex e = numeric(3,2);
+       ex r;
+       
+       try {
+               r = pow(b,e).eval();
+               if (!(r-2*sqrt(ex(2))).is_zero()) {
+                       clog << "2^(3/2) erroneously returned " << r << " instead of 2*sqrt(2)" << endl;
+                       ++result;
+               }
+       } catch (const exception &err) {
+               clog << "2^(3/2) throws " << err.what() << endl;
+               ++result;
+       }
+       return result;
 }
 
 // After the rewriting of basic::normal() & Co. to return {num, den} lists,
@@ -278,7 +278,7 @@ static unsigned exam_paranoia10(void)
 // child (did you get this? Well, never mind...). Fixed on Feb 21th 2000.
 static unsigned exam_paranoia11(void)
 {
-    unsigned result = 0;
+       unsigned result = 0;
        symbol x("x");
 
        ex e = ((-5-2*x)-((2-5*x)/(-2+x))*(3+2*x))/(5-4*x);
@@ -289,7 +289,7 @@ static unsigned exam_paranoia11(void)
                clog << "normal(" << e << ") returns " << f << " instead of " << d << endl;
                ++result;
        }
-    return result;
+       return result;
 }
 
 // This one returned 0 because add::normal() incorrectly assumed that if the
@@ -297,16 +297,16 @@ static unsigned exam_paranoia11(void)
 // be +/-1). Fixed on Aug 2nd 2000.
 static unsigned exam_paranoia12(void)
 {
-    unsigned result = 0;
+       unsigned result = 0;
        symbol x("x");
-
+       
        ex e = 2-2*(1+x)/(-1-x);
        ex f = e.normal();
        ex d = 4;
-
+       
        if (!(f - d).expand().is_zero()) {
                clog << "normal(" << e << ") returns " << f
-             << " instead of " << d << endl;
+                        << " instead of " << d << endl;
                ++result;
        }
        return result;
@@ -316,53 +316,53 @@ static unsigned exam_paranoia12(void)
 // input polynomials against 0. Fixed on Aug 4th 2000.
 static unsigned exam_paranoia13(void)
 {
-    unsigned result = 0;
+       unsigned result = 0;
        symbol a("a"), b("b"), c("c");
-    
+       
        ex e = (b*a-c*a)/(4-a);
        ex d = (c*a-b*a)/(a-4);
-    
-    try {
-        ex f = e.normal();    
-        if (!(f - d).expand().is_zero()) {
-            clog << "normal(" << e << ") returns " << f
-                 << " instead of " << d << endl;
-            ++result;
-        }
-    } catch (const exception &err) {
-        clog << "normal(" << e << ") throws " << err.what() << endl;
-        ++result;
-    }
+       
+       try {
+               ex f = e.normal();      
+               if (!(f - d).expand().is_zero()) {
+                       clog << "normal(" << e << ") returns " << f
+                                << " instead of " << d << endl;
+                       ++result;
+               }
+       } catch (const exception &err) {
+               clog << "normal(" << e << ") throws " << err.what() << endl;
+               ++result;
+       }
        return result;
 }
 
 unsigned exam_paranoia(void)
 {
-    unsigned result = 0;
-    
-    cout << "examining several historic failures just out of paranoia" << flush;
-    clog << "----------several historic failures:" << endl;
-    
-    result += exam_paranoia1();  cout << '.' << flush;
-    result += exam_paranoia2();  cout << '.' << flush;
-    result += exam_paranoia3();  cout << '.' << flush;
-    result += exam_paranoia4();  cout << '.' << flush;
-    result += exam_paranoia5();  cout << '.' << flush;
-    result += exam_paranoia6();  cout << '.' << flush;
-    result += exam_paranoia7();  cout << '.' << flush;
-    result += exam_paranoia8();  cout << '.' << flush;
-    result += exam_paranoia9();  cout << '.' << flush;
-    result += exam_paranoia10();  cout << '.' << flush;
-    result += exam_paranoia11();  cout << '.' << flush;
-    result += exam_paranoia12();  cout << '.' << flush;
-    result += exam_paranoia13();  cout << '.' << flush;
-    
-    if (!result) {
-        cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
-        cout << " failed " << endl;
-    }
-    
-    return result;
+       unsigned result = 0;
+       
+       cout << "examining several historic failures just out of paranoia" << flush;
+       clog << "----------several historic failures:" << endl;
+       
+       result += exam_paranoia1();  cout << '.' << flush;
+       result += exam_paranoia2();  cout << '.' << flush;
+       result += exam_paranoia3();  cout << '.' << flush;
+       result += exam_paranoia4();  cout << '.' << flush;
+       result += exam_paranoia5();  cout << '.' << flush;
+       result += exam_paranoia6();  cout << '.' << flush;
+       result += exam_paranoia7();  cout << '.' << flush;
+       result += exam_paranoia8();  cout << '.' << flush;
+       result += exam_paranoia9();  cout << '.' << flush;
+       result += exam_paranoia10();  cout << '.' << flush;
+       result += exam_paranoia11();  cout << '.' << flush;
+       result += exam_paranoia12();  cout << '.' << flush;
+       result += exam_paranoia13();  cout << '.' << flush;
+       
+       if (!result) {
+               cout << " passed " << endl;
+               clog << "(no output)" << endl;
+       } else {
+               cout << " failed " << endl;
+       }
+       
+       return result;
 }
index 34aea62..5680bec 100644 (file)
@@ -228,27 +228,27 @@ static unsigned poly_gcd7(void)
 
 unsigned exam_polygcd(void)
 {
-    unsigned result = 0;
-    
+       unsigned result = 0;
+       
        cout << "examining polynomial GCD computation" << flush;
        clog << "----------polynomial GCD computation:" << endl;
-    
+       
        result += poly_gcd1();  cout << '.' << flush;
        result += poly_gcd2();  cout << '.' << flush;
        result += poly_gcd3();  cout << '.' << flush;
-    result += poly_gcd3p();     cout << '.' << flush; // takes extremely long (PRS "worst" case)
+       result += poly_gcd3p();  cout << '.' << flush; // takes extremely long (PRS "worst" case)
        result += poly_gcd4();  cout << '.' << flush;
        result += poly_gcd5();  cout << '.' << flush;
        result += poly_gcd5p();  cout << '.' << flush;
        result += poly_gcd6();  cout << '.' << flush;
        result += poly_gcd7();  cout << '.' << flush;
-    
+       
        if (!result) {
                cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
+               clog << "(no output)" << endl;
+       } else {
                cout << " failed " << endl;
-    }
-    
+       }
+       
        return result;
 }
index aba615f..f685ffb 100644 (file)
 
 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;
 }
index a6ff704..133f01b 100644 (file)
@@ -1,4 +1,4 @@
-/** @file exam_pseries.cpp
+/** @File exam_pseries.cpp
  *
  *  Series expansion test (Laurent and Taylor series). */
 
@@ -26,222 +26,222 @@ static symbol x("x");
 
 static unsigned check_series(const ex &e, const ex &point, const ex &d, int order = 8)
 {
-    ex es = e.series(x==point, order);
-    ex ep = ex_to_pseries(es).convert_to_poly();
-    if (!(ep - d).is_zero()) {
-        clog << "series expansion of " << e << " at " << point
-             << " erroneously returned " << ep << " (instead of " << d
-             << ")" << endl;
-        (ep-d).printtree(clog);
-        return 1;
-    }
-    return 0;
+       ex es = e.series(x==point, order);
+       ex ep = ex_to_pseries(es).convert_to_poly();
+       if (!(ep - d).is_zero()) {
+               clog << "series expansion of " << e << " at " << point
+                        << " erroneously returned " << ep << " (instead of " << d
+                        << ")" << endl;
+               (ep-d).printtree(clog);
+               return 1;
+       }
+       return 0;
 }
 
 // Series expansion
 static unsigned exam_series1(void)
 {
-    unsigned result = 0;
-    ex e, d;
-    
-    e = sin(x);
-    d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
-    result += check_series(e, 0, d);
-    
-    e = cos(x);
-    d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
-    result += check_series(e, 0, d);
-    
-    e = exp(x);
-    d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8));
-    result += check_series(e, 0, d);
-    
-    e = pow(1 - x, -1);
-    d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8));
-    result += check_series(e, 0, d);
-    
-    e = x + pow(x, -1);
-    d = x + pow(x, -1);
-    result += check_series(e, 0, d);
-    
-    e = x + pow(x, -1);
-    d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8));
-    result += check_series(e, 1, d);
-    
-    e = pow(x + pow(x, 3), -1);
-    d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
-    result += check_series(e, 0, d);
-    
-    e = pow(pow(x, 2) + pow(x, 4), -1);
-    d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
-    result += check_series(e, 0, d);
-    
-    e = pow(sin(x), -2);
-    d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
-    result += check_series(e, 0, d);
-    
-    e = sin(x) / cos(x);
-    d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
-    result += check_series(e, 0, d);
-    
-    e = cos(x) / sin(x);
-    d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
-    result += check_series(e, 0, d);
-    
-    e = pow(numeric(2), x);
-    ex t = log(2) * x;
-    d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
-    result += check_series(e, 0, d.expand());
-    
-    e = pow(Pi, x);
-    t = log(Pi) * x;
-    d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
-    result += check_series(e, 0, d.expand());
-    
-    return result;
+       unsigned result = 0;
+       ex e, d;
+       
+       e = sin(x);
+       d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
+       result += check_series(e, 0, d);
+       
+       e = cos(x);
+       d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
+       result += check_series(e, 0, d);
+       
+       e = exp(x);
+       d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8));
+       result += check_series(e, 0, d);
+       
+       e = pow(1 - x, -1);
+       d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8));
+       result += check_series(e, 0, d);
+       
+       e = x + pow(x, -1);
+       d = x + pow(x, -1);
+       result += check_series(e, 0, d);
+       
+       e = x + pow(x, -1);
+       d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8));
+       result += check_series(e, 1, d);
+       
+       e = pow(x + pow(x, 3), -1);
+       d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
+       result += check_series(e, 0, d);
+       
+       e = pow(pow(x, 2) + pow(x, 4), -1);
+       d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
+       result += check_series(e, 0, d);
+       
+       e = pow(sin(x), -2);
+       d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
+       result += check_series(e, 0, d);
+       
+       e = sin(x) / cos(x);
+       d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
+       result += check_series(e, 0, d);
+       
+       e = cos(x) / sin(x);
+       d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
+       result += check_series(e, 0, d);
+       
+       e = pow(numeric(2), x);
+       ex t = log(2) * x;
+       d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
+       result += check_series(e, 0, d.expand());
+       
+       e = pow(Pi, x);
+       t = log(Pi) * x;
+       d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
+       result += check_series(e, 0, d.expand());
+       
+       return result;
 }
 
 // Series addition
 static unsigned exam_series2(void)
 {
-    unsigned result = 0;
-    ex e, d;
-    
-    e = pow(sin(x), -1).series(x==0, 8) + pow(sin(-x), -1).series(x==0, 12);
-    d = Order(pow(x, 6));
-    result += check_series(e, 0, d);
-    
-    return result;
+       unsigned result = 0;
+       ex e, d;
+       
+       e = pow(sin(x), -1).series(x==0, 8) + pow(sin(-x), -1).series(x==0, 12);
+       d = Order(pow(x, 6));
+       result += check_series(e, 0, d);
+       
+       return result;
 }
 
 // Series multiplication
 static unsigned exam_series3(void)
 {
-    unsigned result = 0;
-    ex e, d;
-    
-    e = sin(x).series(x==0, 8) * pow(sin(x), -1).series(x==0, 12);
-    d = 1 + Order(pow(x, 7));
-    result += check_series(e, 0, d);
-    
-    return result;
+       unsigned result = 0;
+       ex e, d;
+       
+       e = sin(x).series(x==0, 8) * pow(sin(x), -1).series(x==0, 12);
+       d = 1 + Order(pow(x, 7));
+       result += check_series(e, 0, d);
+       
+       return result;
 }
 
 // Order term handling
 static unsigned exam_series4(void)
 {
-    unsigned result = 0;
-    ex e, d;
+       unsigned result = 0;
+       ex e, d;
 
-    e = 1 + x + pow(x, 2) + pow(x, 3);
-    d = Order(1);
-    result += check_series(e, 0, d, 0);
-    d = 1 + Order(x);
-    result += check_series(e, 0, d, 1);
-    d = 1 + x + Order(pow(x, 2));
-    result += check_series(e, 0, d, 2);
-    d = 1 + x + pow(x, 2) + Order(pow(x, 3));
-    result += check_series(e, 0, d, 3);
-    d = 1 + x + pow(x, 2) + pow(x, 3);
-    result += check_series(e, 0, d, 4);
-    return result;
+       e = 1 + x + pow(x, 2) + pow(x, 3);
+       d = Order(1);
+       result += check_series(e, 0, d, 0);
+       d = 1 + Order(x);
+       result += check_series(e, 0, d, 1);
+       d = 1 + x + Order(pow(x, 2));
+       result += check_series(e, 0, d, 2);
+       d = 1 + x + pow(x, 2) + Order(pow(x, 3));
+       result += check_series(e, 0, d, 3);
+       d = 1 + x + pow(x, 2) + pow(x, 3);
+       result += check_series(e, 0, d, 4);
+       return result;
 }
 
 // Series expansion of tgamma(-1)
 static unsigned exam_series5(void)
 {
-    ex e = tgamma(2*x);
-    ex d = pow(x+1,-1)*numeric(1,4) +
-           pow(x+1,0)*(numeric(3,4) -
-                       numeric(1,2)*Euler) +
-           pow(x+1,1)*(numeric(7,4) -
-                       numeric(3,2)*Euler +
-                       numeric(1,2)*pow(Euler,2) +
-                       numeric(1,12)*pow(Pi,2)) +
-           pow(x+1,2)*(numeric(15,4) -
-                       numeric(7,2)*Euler -
-                       numeric(1,3)*pow(Euler,3) +
-                       numeric(1,4)*pow(Pi,2) +
-                       numeric(3,2)*pow(Euler,2) -
-                       numeric(1,6)*pow(Pi,2)*Euler -
-                       numeric(2,3)*zeta(3)) +
-           pow(x+1,3)*(numeric(31,4) - pow(Euler,3) -
-                       numeric(15,2)*Euler +
-                       numeric(1,6)*pow(Euler,4) +
-                       numeric(7,2)*pow(Euler,2) +
-                       numeric(7,12)*pow(Pi,2) -
-                       numeric(1,2)*pow(Pi,2)*Euler -
-                       numeric(2)*zeta(3) +
-                       numeric(1,6)*pow(Euler,2)*pow(Pi,2) +
-                       numeric(1,40)*pow(Pi,4) +
-                       numeric(4,3)*zeta(3)*Euler) +
-           Order(pow(x+1,4));
-    return check_series(e, -1, d, 4);
+       ex e = tgamma(2*x);
+       ex d = pow(x+1,-1)*numeric(1,4) +
+                  pow(x+1,0)*(numeric(3,4) -
+                                          numeric(1,2)*Euler) +
+                  pow(x+1,1)*(numeric(7,4) -
+                                          numeric(3,2)*Euler +
+                                          numeric(1,2)*pow(Euler,2) +
+                                          numeric(1,12)*pow(Pi,2)) +
+                  pow(x+1,2)*(numeric(15,4) -
+                                          numeric(7,2)*Euler -
+                                          numeric(1,3)*pow(Euler,3) +
+                                          numeric(1,4)*pow(Pi,2) +
+                                          numeric(3,2)*pow(Euler,2) -
+                                          numeric(1,6)*pow(Pi,2)*Euler -
+                                          numeric(2,3)*zeta(3)) +
+                  pow(x+1,3)*(numeric(31,4) - pow(Euler,3) -
+                                          numeric(15,2)*Euler +
+                                          numeric(1,6)*pow(Euler,4) +
+                                          numeric(7,2)*pow(Euler,2) +
+                                          numeric(7,12)*pow(Pi,2) -
+                                          numeric(1,2)*pow(Pi,2)*Euler -
+                                          numeric(2)*zeta(3) +
+                                          numeric(1,6)*pow(Euler,2)*pow(Pi,2) +
+                                          numeric(1,40)*pow(Pi,4) +
+                                          numeric(4,3)*zeta(3)*Euler) +
+                  Order(pow(x+1,4));
+       return check_series(e, -1, d, 4);
 }
-    
+       
 // Series expansion of tan(x==Pi/2)
 static unsigned exam_series6(void)
 {
-    ex e = tan(x*Pi/2);
-    ex d = pow(x-1,-1)/Pi*(-2) + pow(x-1,1)*Pi/6 + pow(x-1,3)*pow(Pi,3)/360
-          +pow(x-1,5)*pow(Pi,5)/15120 + pow(x-1,7)*pow(Pi,7)/604800
-          +Order(pow(x-1,8));
-    return check_series(e,1,d,8);
+       ex e = tan(x*Pi/2);
+       ex d = pow(x-1,-1)/Pi*(-2) + pow(x-1,1)*Pi/6 + pow(x-1,3)*pow(Pi,3)/360
+                 +pow(x-1,5)*pow(Pi,5)/15120 + pow(x-1,7)*pow(Pi,7)/604800
+                 +Order(pow(x-1,8));
+       return check_series(e,1,d,8);
 }
 
 // Series expansion of log(sin(x==0))
 static unsigned exam_series7(void)
 {
-    ex e = log(sin(x));
-    ex d = log(x) - pow(x,2)/6 - pow(x,4)/180 - pow(x,6)/2835
-          +Order(pow(x,8));
-    return check_series(e,0,d,8);
+       ex e = log(sin(x));
+       ex d = log(x) - pow(x,2)/6 - pow(x,4)/180 - pow(x,6)/2835
+                 +Order(pow(x,8));
+       return check_series(e,0,d,8);
 }
 
 // Series expansion of Li2(sin(x==0))
 static unsigned exam_series8(void)
 {
-    ex e = Li2(sin(x));
-    ex d = x + pow(x,2)/4 - pow(x,3)/18 - pow(x,4)/48
-           - 13*pow(x,5)/1800 - pow(x,6)/360 - 23*pow(x,7)/21168
-           + Order(pow(x,8));
-    return check_series(e,0,d,8);
+       ex e = Li2(sin(x));
+       ex d = x + pow(x,2)/4 - pow(x,3)/18 - pow(x,4)/48
+                  - 13*pow(x,5)/1800 - pow(x,6)/360 - 23*pow(x,7)/21168
+                  + Order(pow(x,8));
+       return check_series(e,0,d,8);
 }
 
 // Series expansion of Li2((x==2)^2), caring about branch-cut
 static unsigned exam_series9(void)
 {
-    ex e = Li2(pow(x,2));
-    ex d = Li2(4) + (-log(3) + I*Pi*csgn(I-I*pow(x,2))) * (x-2)
-           + (numeric(-2,3) + log(3)/4 - I*Pi/4*csgn(I-I*pow(x,2))) * pow(x-2,2)
-           + (numeric(11,27) - log(3)/12 + I*Pi/12*csgn(I-I*pow(x,2))) * pow(x-2,3)
-           + (numeric(-155,648) + log(3)/32 - I*Pi/32*csgn(I-I*pow(x,2))) * pow(x-2,4)
-           + Order(pow(x-2,5));
-    return check_series(e,2,d,5);
+       ex e = Li2(pow(x,2));
+       ex d = Li2(4) + (-log(3) + I*Pi*csgn(I-I*pow(x,2))) * (x-2)
+                  + (numeric(-2,3) + log(3)/4 - I*Pi/4*csgn(I-I*pow(x,2))) * pow(x-2,2)
+                  + (numeric(11,27) - log(3)/12 + I*Pi/12*csgn(I-I*pow(x,2))) * pow(x-2,3)
+                  + (numeric(-155,648) + log(3)/32 - I*Pi/32*csgn(I-I*pow(x,2))) * pow(x-2,4)
+                  + Order(pow(x-2,5));
+       return check_series(e,2,d,5);
 }
 
 unsigned exam_pseries(void)
 {
-    unsigned result = 0;
-    
-    cout << "examining series expansion" << flush;
-    clog << "----------series expansion:" << endl;
-    
-    result += exam_series1();  cout << '.' << flush;
-    result += exam_series2();  cout << '.' << flush;
-    result += exam_series3();  cout << '.' << flush;
-    result += exam_series4();  cout << '.' << flush;
-    result += exam_series5();  cout << '.' << flush;
-    result += exam_series6();  cout << '.' << flush;
-    result += exam_series7();  cout << '.' << flush;
-    result += exam_series8();  cout << '.' << flush;
-    result += exam_series9();  cout << '.' << flush;
-    
-    if (!result) {
-        cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
-        cout << " failed " << endl;
-    }
-    return result;
+       unsigned result = 0;
+       
+       cout << "examining series expansion" << flush;
+       clog << "----------series expansion:" << endl;
+       
+       result += exam_series1();  cout << '.' << flush;
+       result += exam_series2();  cout << '.' << flush;
+       result += exam_series3();  cout << '.' << flush;
+       result += exam_series4();  cout << '.' << flush;
+       result += exam_series5();  cout << '.' << flush;
+       result += exam_series6();  cout << '.' << flush;
+       result += exam_series7();  cout << '.' << flush;
+       result += exam_series8();  cout << '.' << flush;
+       result += exam_series9();  cout << '.' << flush;
+       
+       if (!result) {
+               cout << " passed " << endl;
+               clog << "(no output)" << endl;
+       } else {
+               cout << " failed " << endl;
+       }
+       return result;
 }
index 5aa327a..01e8f25 100644 (file)
 
 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;
 }
index 18c2ded..b61e8e8 100644 (file)
@@ -35,12 +35,12 @@ using namespace GiNaC;
 const ex
 dense_univariate_poly(const symbol & x, unsigned degree)
 {
-    ex unipoly;
-    
-    for (unsigned i=0; i<=degree; ++i)
-        unipoly += numeric((rand()-RAND_MAX/2))*pow(x,i);
-    
-    return unipoly;
+       ex unipoly;
+       
+       for (unsigned i=0; i<=degree; ++i)
+               unipoly += numeric((rand()-RAND_MAX/2))*pow(x,i);
+       
+       return unipoly;
 }
 
 /* Create a dense bivariate random polynomial in x1 and x2.
@@ -49,110 +49,110 @@ dense_univariate_poly(const symbol & x, unsigned degree)
 const ex
 dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree)
 {
-    ex bipoly;
-    
-    for (unsigned i1=0; i1<=degree; ++i1)
-        for (unsigned i2=0; i2<=degree-i1; ++i2)
-            bipoly += numeric((rand()-RAND_MAX/2))*pow(x1,i1)*pow(x2,i2);
-    
-    return bipoly;
+       ex bipoly;
+       
+       for (unsigned i1=0; i1<=degree; ++i1)
+               for (unsigned i2=0; i2<=degree-i1; ++i2)
+                       bipoly += numeric((rand()-RAND_MAX/2))*pow(x1,i1)*pow(x2,i2);
+       
+       return bipoly;
 }
 
 /* Chose a randum symbol or number from the argument list. */
 const ex
 random_symbol(const symbol & x,
-              const symbol & y,
-              const symbol & z,
-              bool rational = true,
-              bool complex = false)
+                         const symbol & y,
+                         const symbol & z,
+                         bool rational = true,
+                         bool complex = false)
 {
-    ex e;
-    switch (abs(rand()) % 4) {
-        case 0:
-            e = x;
-            break;
-        case 1:
-            e = y;
-            break;
-        case 2:
-            e = z;
-            break;
-        case 3: {
-            int c1;
-            do { c1 = rand()%20 - 10; } while (!c1);
-            int c2;
-            do { c2 = rand()%20 - 10; } while (!c2);
-            if (!rational)
-                c2 = 1;
-            e = numeric(c1, c2);
-            if (complex && !(rand()%5))
-                e = e*I;
-            break;
-        }
-    }
-    return e;
+       ex e;
+       switch (abs(rand()) % 4) {
+               case 0:
+                       e = x;
+                       break;
+               case 1:
+                       e = y;
+                       break;
+               case 2:
+                       e = z;
+                       break;
+               case 3: {
+                       int c1;
+                       do { c1 = rand()%20 - 10; } while (!c1);
+                       int c2;
+                       do { c2 = rand()%20 - 10; } while (!c2);
+                       if (!rational)
+                               c2 = 1;
+                       e = numeric(c1, c2);
+                       if (complex && !(rand()%5))
+                               e = e*I;
+                       break;
+               }
+       }
+       return e;
 }
 
 /* Create a sparse random tree in three symbols. */
 const ex
 sparse_tree(const symbol & x,
-            const symbol & y,
-            const symbol & z,
-            int level,
-            bool trig = false,    // true includes trigonomatric functions
-            bool rational = true, // false excludes coefficients in Q
-            bool complex = false) // true includes complex numbers
+                       const symbol & y,
+                       const symbol & z,
+                       int level,
+                       bool trig = false,      // true includes trigonomatric functions
+                       bool rational = true, // false excludes coefficients in Q
+                       bool complex = false) // true includes complex numbers
 {
-    if (level == 0)
-        return random_symbol(x,y,z,rational,complex);
-    switch (abs(rand()) % 10) {
-        case 0:
-        case 1:
-        case 2:
-        case 3:
-            return add(sparse_tree(x,y,z,level-1, trig, rational),
-                       sparse_tree(x,y,z,level-1, trig, rational));
-        case 4:
-        case 5:
-        case 6:
-            return mul(sparse_tree(x,y,z,level-1, trig, rational),
-                       sparse_tree(x,y,z,level-1, trig, rational));
-        case 7:
-        case 8: {
-            ex powbase;
-            do {
-                powbase = sparse_tree(x,y,z,level-1, trig, rational);
-            } while (powbase.is_zero());
-            return pow(powbase, abs(rand() % 4));
-            break;
-        }
-        case 9:
-            if (trig) {
-                switch (abs(rand()) % 4) {
-                    case 0:
-                        return sin(sparse_tree(x,y,z,level-1, trig, rational));
-                    case 1:
-                        return cos(sparse_tree(x,y,z,level-1, trig, rational));
-                    case 2:
-                        return exp(sparse_tree(x,y,z,level-1, trig, rational));
-                    case 3: {
-                        ex logex;
-                        do {
-                            ex logarg;
-                            do {
-                                logarg = sparse_tree(x,y,z,level-1, trig, rational);
-                            } while (logarg.is_zero());
-                            // Keep the evaluator from accidentally plugging an
-                            // unwanted I in the tree:
-                            if (!complex && logarg.info(info_flags::negative))
-                                logarg = -logarg;
-                            logex = log(logarg);
-                        } while (logex.is_zero());
-                        return logex;
-                        break;
-                    }
-                }
-            } else
-                return random_symbol(x,y,z,rational,complex);
-    }
+       if (level == 0)
+               return random_symbol(x,y,z,rational,complex);
+       switch (abs(rand()) % 10) {
+               case 0:
+               case 1:
+               case 2:
+               case 3:
+                       return add(sparse_tree(x,y,z,level-1, trig, rational),
+                                          sparse_tree(x,y,z,level-1, trig, rational));
+               case 4:
+               case 5:
+               case 6:
+                       return mul(sparse_tree(x,y,z,level-1, trig, rational),
+                                          sparse_tree(x,y,z,level-1, trig, rational));
+               case 7:
+               case 8: {
+                       ex powbase;
+                       do {
+                               powbase = sparse_tree(x,y,z,level-1, trig, rational);
+                       } while (powbase.is_zero());
+                       return pow(powbase, abs(rand() % 4));
+                       break;
+               }
+               case 9:
+                       if (trig) {
+                               switch (abs(rand()) % 4) {
+                                       case 0:
+                                               return sin(sparse_tree(x,y,z,level-1, trig, rational));
+                                       case 1:
+                                               return cos(sparse_tree(x,y,z,level-1, trig, rational));
+                                       case 2:
+                                               return exp(sparse_tree(x,y,z,level-1, trig, rational));
+                                       case 3: {
+                                               ex logex;
+                                               do {
+                                                       ex logarg;
+                                                       do {
+                                                               logarg = sparse_tree(x,y,z,level-1, trig, rational);
+                                                       } while (logarg.is_zero());
+                                                       // Keep the evaluator from accidentally plugging an
+                                                       // unwanted I in the tree:
+                                                       if (!complex && logarg.info(info_flags::negative))
+                                                               logarg = -logarg;
+                                                       logex = log(logarg);
+                                               } while (logex.is_zero());
+                                               return logex;
+                                               break;
+                                       }
+                               }
+                       } else
+                               return random_symbol(x,y,z,rational,complex);
+       }
 }
index db9ff80..e20fdde 100644 (file)
 
 #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;
 }
index 0bfbe30..59a08a9 100644 (file)
 
 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;
 }
index 2fa40c9..ac56eed 100644 (file)
 
 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;
 }
index 3431203..6fd60da 100644 (file)
 
 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;
 }
index 8885d87..ce2921c 100644 (file)
 
 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;
 }
index 1f63e58..f226f51 100644 (file)
 
 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;
 }
index 5dc004d..95a4fd4 100644 (file)
 
 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;
 }
index ce8920e..de8aa92 100644 (file)
 
 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;
 }
index 75dd426..1bae9c8 100644 (file)
 
 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;
 }
index a7c5c85..8ea4b47 100644 (file)
 
 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;
 }
index 62f69b6..41c0770 100644 (file)
 
 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