- introduced info_flag::algebraic.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Sat, 29 Apr 2000 02:15:55 +0000 (02:15 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Sat, 29 Apr 2000 02:15:55 +0000 (02:15 +0000)
- matrix::fraction_free_elimination(bool) finished and used by matrix::determinant().
- introduced some comments.

ginac/add.cpp
ginac/basic.cpp
ginac/flags.h
ginac/matrix.cpp
ginac/matrix.h
ginac/mul.cpp
ginac/normal.cpp
ginac/numeric.cpp
ginac/power.cpp

index 717279b..8a7e904 100644 (file)
@@ -290,21 +290,28 @@ void add::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons
 
 bool add::info(unsigned inf) const
 {
-    // TODO: optimize
-    if (inf==info_flags::polynomial ||
-        inf==info_flags::integer_polynomial ||
-        inf==info_flags::cinteger_polynomial ||
-        inf==info_flags::rational_polynomial ||
-        inf==info_flags::crational_polynomial ||
-        inf==info_flags::rational_function) {
-        for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
-            if (!(recombine_pair_to_ex(*it).info(inf)))
-                return false;
+    switch (inf) {
+        case info_flags::polynomial:
+        case info_flags::integer_polynomial:
+        case info_flags::cinteger_polynomial:
+        case info_flags::rational_polynomial:
+        case info_flags::crational_polynomial:
+        case info_flags::rational_function: {
+            for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+                if (!(recombine_pair_to_ex(*i).info(inf)))
+                    return false;
+            }
+            return overall_coeff.info(inf);
+        }
+        case info_flags::algebraic: {
+            for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+                if ((recombine_pair_to_ex(*i).info(inf)))
+                    return true;
+            }
+            return false;
         }
-        return overall_coeff.info(inf);
-    } else {
-        return inherited::info(inf);
     }
+    return inherited::info(inf);
 }
 
 int add::degree(const symbol & s) const
index 2a4a3e3..83b777b 100644 (file)
@@ -284,12 +284,14 @@ ex basic::collect(const symbol & s) const
 /* Perform automatic symbolic evaluations on expression. */
 ex basic::eval(int level) const
 {
+    // There is nothing to do for basic objects:
     return this->hold();
 }
 
 /** Evaluate object numerically. */
 ex basic::evalf(int level) const
 {
+    // There is nothing to do for basic objects:
     return *this;
 }
 
index 16228a2..8b9caed 100644 (file)
@@ -45,7 +45,7 @@ public:
 /** Possible attributes an object can have. */
 class info_flags {
 public:
-    enum { 
+    enum {
            // answered by class numeric
            numeric,
            real,
@@ -62,7 +62,7 @@ public:
            even,
            odd,
            prime,
-
+           
            // answered by class relation
            relation,
            relation_equal,
@@ -71,16 +71,16 @@ public:
            relation_less_or_equal,
            relation_greater,
            relation_greater_or_equal,
-
+           
            // answered by class symbol
            symbol,
-
+           
            // answered by class lst
            list,
-
+           
            // answered by class exprseq
            exprseq,
-
+           
            // answered by classes numeric, symbol, add, mul, power
            polynomial,
            integer_polynomial,
@@ -88,17 +88,18 @@ public:
            rational_polynomial,
            crational_polynomial,
            rational_function,
-
+           algebraic,
+           
            // answered by class indexed
            indexed,      // class can carry indices
            has_indices,  // object has at least one index
-
+           
            // answered by class idx
            idx,
-
+           
            // answered by class coloridx
            coloridx,
-
+           
            // answered by class lorentzidx
            lorentzidx
     };
index f76aa9d..c5d056c 100644 (file)
@@ -460,24 +460,28 @@ ex matrix::determinant(void) const
     bool normal_flag = false;
     unsigned sparse_count = 0;  // count non-zero elements
     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
-        if (!(*r).is_zero()) {
+        if (!(*r).is_zero())
             ++sparse_count;
-        }
-        if (!(*r).info(info_flags::numeric)) {
+        if (!(*r).info(info_flags::numeric))
             numeric_flag = false;
-        }
         if ((*r).info(info_flags::rational_function) &&
-            !(*r).info(info_flags::crational_polynomial)) {
+            !(*r).info(info_flags::crational_polynomial))
             normal_flag = true;
-        }
     }
     
     if (numeric_flag)  // purely numeric matrix
         return determinant_numeric();
-    
     // Does anybody really know when a matrix is sparse?
-    if (4*sparse_count<row*col) {  // < row/2 nonzero elements average in a row
-        return determinant_bareiss(normal_flag);
+    // Maybe <~row/2.2 nonzero elements average in a row?
+    if (5*sparse_count<=row*col) {
+        // copy *this:
+        matrix tmp(*this);
+        int sign;
+        sign = tmp.fraction_free_elimination(true);
+        if (normal_flag)
+            return (sign*tmp.m[row*col-1]).normal();
+        else
+            return (sign*tmp.m[row*col-1]).expand();
     }
     
     // Now come the minor expansion schemes.  We always develop such that the
@@ -564,7 +568,7 @@ ex matrix::charpoly(const symbol & lambda) const
     
     // The pure numeric case is traditionally rather common.  Hence, it is
     // trapped and we use Leverrier's algorithm which goes as row^3 for
-    // every coefficient.  The expensive section is the matrix multiplication.
+    // every coefficient.  The expensive part is the matrix multiplication.
     if (numeric_flag) {
         matrix B(*this);
         ex c = B.trace();
@@ -669,7 +673,7 @@ ex matrix::ffe_get(unsigned r, unsigned c) const
 matrix matrix::fraction_free_elim(const matrix & vars,
                                   const matrix & rhs) const
 {
-    // FIXME: use implementation of matrix::fraction_free_elim
+    // FIXME: use implementation of matrix::fraction_free_elimination
     if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
         throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
     
@@ -712,18 +716,13 @@ matrix matrix::fraction_free_elim(const matrix & vars,
             divisor = a.ffe_get(r,k);
             r++;
         }
-    }
-    // optionally compute the determinant for square or augmented matrices
-    // if (r==m+1) { det = sign*divisor; } else { det = 0; }
-    
-    /*
-    for (unsigned r=1; r<=m; ++r) {
-        for (unsigned c=1; c<=n; ++c) {
-            cout << a.ffe_get(r,c) << "\t";
-        }
-        cout << " | " <<  b.ffe_get(r,1) << endl;
-    }
-    */
+    }    
+//     for (unsigned r=1; r<=m; ++r) {
+//         for (unsigned c=1; c<=n; ++c) {
+//             cout << a.ffe_get(r,c) << "\t";
+//         }
+//         cout << " | " <<  b.ffe_get(r,1) << endl;
+//     }
     
 #ifdef DO_GINAC_ASSERT
     // test if we really have an upper echelon matrix
@@ -731,7 +730,7 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned r=1; r<=m; ++r) {
         int zero_in_this_row=0;
         for (unsigned c=1; c<=n; ++c) {
-            if (a.ffe_get(r,c).is_equal(_ex0()))
+            if (a.ffe_get(r,c).is_zero())
                zero_in_this_row++;
             else
                 break;
@@ -741,12 +740,6 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     }
 #endif // def DO_GINAC_ASSERT
     
-    /*
-    cout << "after" << endl;
-    cout << "a=" << a << endl;
-    cout << "b=" << b << endl;
-    */
-    
     // assemble solution
     matrix sol(n,1);
     unsigned last_assigned_sol = n+1;
@@ -1015,112 +1008,6 @@ ex matrix::determinant_minor(void) const
     return det;
 }
 
-/** Helper function to divide rational functions, as needed in any Bareiss
- *  elimination scheme over a quotient field.
- *  
- *  @see divide() */
-bool rat_divide(const ex & a, const ex & b, ex & q, bool check_args = true)
-{
-    q = _ex0();
-    if (b.is_zero())
-        throw(std::overflow_error("rat_divide(): division by zero"));
-    if (a.is_zero())
-        return true;
-    if (is_ex_exactly_of_type(b, numeric)) {
-        q = a / b;
-        return true;
-    } else if (is_ex_exactly_of_type(a, numeric))
-        return false;
-    ex a_n = a.numer();
-    ex a_d = a.denom();
-    ex b_n = b.numer();
-    ex b_d = b.denom();
-    ex n;  // new numerator
-    ex d;  // new denominator
-    bool check = true;
-    check &= divide(a_n, b_n, n, check_args);
-    check &= divide(a_d, b_d, d, check_args);
-    q = n/d;
-    return check;
-}
-
-
-/** Determinant computed by using fraction free elimination.  This
- *  routine is only called internally by matrix::determinant().
- *
- *  @param normalize may be set to false only in integral domains. */
-ex matrix::determinant_bareiss(bool normalize) const
-{
-    if (rows()==1)
-        return m[0];
-    
-    int sign = 1;
-    ex divisor = 1;
-    ex dividend;
-    
-    // we populate a tmp matrix to subsequently operate on, it should
-    // be normalized even though this algorithm doesn't need GCDs since
-    // the elements of *this might be unnormalized, which complicates
-    // things:
-    matrix tmp(*this);
-    exvector::const_iterator i = m.begin();
-    exvector::iterator ti = tmp.m.begin();
-    for (; i!= m.end(); ++i, ++ti) {
-        if (normalize)
-            (*ti) = (*i).normal();
-        else
-            (*ti) = (*i);
-    }
-    
-    for (unsigned r1=0; r1<row-1; ++r1) {
-        int indx = tmp.pivot(r1);
-        if (indx==-1)
-            return _ex0();
-        if (indx>0)
-            sign = -sign;
-        if (r1>0) {
-            divisor = tmp.m[(r1-1)*col+(r1-1)].expand();
-            // delete the elements we don't need anymore:
-            for (unsigned c=0; c<col; ++c)
-                tmp.m[(r1-1)*col+c] = _ex0();
-        }
-        for (unsigned r2=r1+1; r2<row; ++r2) {
-            for (unsigned c=r1+1; c<col; ++c) {
-                lst srl;  // symbol replacement list for .to_rational()
-                dividend = (tmp.m[r1*tmp.col+r1]*tmp.m[r2*tmp.col+c]
-                           -tmp.m[r2*tmp.col+r1]*tmp.m[r1*tmp.col+c]).expand();
-                if (normalize) {
-#ifdef DO_GINAC_ASSERT
-                    GINAC_ASSERT(rat_divide(dividend.to_rational(srl),
-                                            divisor.to_rational(srl),
-                                            tmp.m[r2*tmp.col+c],true));
-#else
-                    rat_divide(dividend.to_rational(srl),
-                               divisor.to_rational(srl),
-                               tmp.m[r2*tmp.col+c],false);
-#endif                    
-                }
-                else {
-#ifdef DO_GINAC_ASSERT
-                    GINAC_ASSERT(divide(dividend.to_rational(srl),
-                                        divisor.to_rational(srl),
-                                        tmp.m[r2*tmp.col+c],true));
-#else
-                    divide(dividend.to_rational(srl),
-                           divisor.to_rational(srl),
-                           tmp.m[r2*tmp.col+c],false);
-#endif                    
-                }
-                tmp.m[r2*tmp.col+c] = tmp.m[r2*tmp.col+c].subs(srl);
-            }
-            for (unsigned c=0; c<=r1; ++c)
-                tmp.m[r2*tmp.col+c] = _ex0();
-        }
-    }
-    
-    return sign*tmp.m[tmp.row*tmp.col-1];
-}
-
 
 /** Perform the steps of an ordinary Gaussian elimination to bring the matrix
  *  into an upper echelon form.
@@ -1177,52 +1064,128 @@ int matrix::division_free_elimination(void)
 
 
 /** Perform the steps of Bareiss' one-step fraction free elimination to bring
- *  the matrix into an upper echelon form.
- *
+ *  the matrix into an upper echelon form.  Fraction free elimination means
+ *  that divide is used straightforwardly, without computing GCDs first.  This
+ *  is possible, since we know the divisor at each step.
+ *  
+ *  @param det may be set to true to save a lot of space if one is only
+ *  interested in the last element (i.e. for calculating determinants), the
+ *  others are set to zero in this case.
  *  @return sign is 1 if an even number of rows was swapped, -1 if an odd
  *  number of rows was swapped and 0 if the matrix is singular. */
-int matrix::fraction_free_elimination(void)
+int matrix::fraction_free_elimination(bool det)
 {
-    ensure_if_modifiable();
+    // Method:
+    // (single-step fraction free elimination scheme, already known to Jordan)
+    //
+    // Usual division-free elimination sets m[0](r,c) = m(r,c) and then sets
+    //     m[k+1](r,c) = m[k](k,k) * m[k](r,c) - m[k](r,k) * m[k](k,c).
+    //
+    // Bareiss (fraction-free) elimination in addition divides that element
+    // by m[k-1](k-1,k-1) for k>1, where it can be shown by means of the
+    // Sylvester determinant that this really divides m[k+1](r,c).
+    //
+    // We also allow rational functions where the original prove still holds.
+    // However, we must care for numerator and denominator separately and
+    // "manually" work in the integral domains because of subtle cancellations
+    // (see below).  This blows up the bookkeeping a bit and the formula has
+    // to be modified to expand like this (N{x} stands for numerator of x,
+    // D{x} for denominator of x):
+    //     N{m[k+1](r,c)} = N{m[k](k,k)}*N{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)}
+    //                     -N{m[k](r,k)}*N{m[k](k,c)}*D{m[k](k,k)}*D{m[k](r,c)}
+    //     D{m[k+1](r,c)} = D{m[k](k,k)}*D{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)}
+    // where for k>1 we now divide N{m[k+1](r,c)} by
+    //     N{m[k-1](k-1,k-1)}
+    // and D{m[k+1](r,c)} by
+    //     D{m[k-1](k-1,k-1)}.
     
-    // first normal all elements:
-    for (exvector::iterator i=m.begin(); i!=m.end(); ++i)
-        (*i) = (*i).normal();
+    GINAC_ASSERT(det || row==col);
+    ensure_if_modifiable();
+    if (rows()==1)
+        return 1;
     
-    // FIXME: this is unfinished, once matrix::determinant_bareiss is
-    // bulletproof, some code ought to be copy from there to here.
     int sign = 1;
-    ex divisor = 1;
-    ex dividend;
-    lst srl;      // symbol replacement list for .to_rational()
+    ex divisor_n = 1;
+    ex divisor_d = 1;
+    ex dividend_n;
+    ex dividend_d;
+    
+    // We populate temporary matrices to subsequently operate on.  There is
+    // one holding numerators and another holding denominators of entries.
+    // This is a must since the evaluator (or even earlier mul's constructor)
+    // might cancel some trivial element which causes divide() to fail.  The
+    // elements are normalized first (yes, even though this algorithm doesn't
+    // need GCDs) since the elements of *this might be unnormalized, which
+    // makes things more complicated than they need to be.
+    matrix tmp_n(*this);
+    matrix tmp_d(row,col);  // for denominators, if needed
+    lst srl;  // symbol replacement list
+    exvector::iterator it = m.begin();
+    exvector::iterator tmp_n_it = tmp_n.m.begin();
+    exvector::iterator tmp_d_it = tmp_d.m.begin();
+    for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it) {
+        (*tmp_n_it) = (*it).normal().to_rational(srl);
+        (*tmp_d_it) = (*tmp_n_it).denom();
+        (*tmp_n_it) = (*tmp_n_it).numer();
+    }
     
     for (unsigned r1=0; r1<row-1; ++r1) {
-        int indx = pivot(r1);
-        if (indx==-1)
-            return 0;  // Note: leaves *this in a messy state.
-        if (indx>0)
+//         cout << "==<" << r1 << ">" << string(60,'=') << endl;
+        int indx = tmp_n.pivot(r1);
+        if (det && indx==-1)
+            return 0;  // FIXME: what to do if det is false?
+        if (indx>0) {
             sign = -sign;
-        if (r1>0)
-            divisor = this->m[(r1-1)*col+(r1-1)].expand();
+            // rows r1 and indx were swapped, so pivot matrix tmp_d:
+            for (unsigned c=0; c<col; ++c)
+                tmp_d.m[row*indx+c].swap(tmp_d.m[row*r1+c]);
+        }
+//         cout << tmp_n << endl;
+//         cout << tmp_d << endl;
+        if (r1>0) {
+            divisor_n = tmp_n.m[(r1-1)*col+(r1-1)].expand();
+            divisor_d = tmp_d.m[(r1-1)*col+(r1-1)].expand();
+            // save space by deleting no longer needed elements:
+            if (det) {
+                for (unsigned c=0; c<col; ++c) {
+                    tmp_n.m[(r1-1)*col+c] = 0;
+                    tmp_d.m[(r1-1)*col+c] = 1;
+                }
+            }
+        }
         for (unsigned r2=r1+1; r2<row; ++r2) {
             for (unsigned c=r1+1; c<col; ++c) {
-                dividend = (this->m[r1*col+r1]*this->m[r2*col+c]
-                           -this->m[r2*col+r1]*this->m[r1*col+c]).expand();
-#ifdef DO_GINAC_ASSERT
-                GINAC_ASSERT(divide(dividend.to_rational(srl),
-                                    divisor.to_rational(srl),
-                                    this->m[r2*col+c]));
-#else
-                divide(dividend.to_rational(srl),
-                       divisor.to_rational(srl),
-                       this->m[r2*col+c]);
-#endif // DO_GINAC_ASSERT
-                this->m[r2*col+c] = this->m[r2*col+c].subs(srl);
+                dividend_n = (tmp_n.m[r1*col+r1]*tmp_n.m[r2*col+c]*
+                              tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]
+                             -tmp_n.m[r2*col+r1]*tmp_n.m[r1*col+c]*
+                              tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand();
+                dividend_d = (tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]*
+                              tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand();
+//                 cout << "Element " << r2 << ',' << c << endl;
+//                 cout << "dividend_n==" << dividend_n << endl;
+//                 cout << "dividend_d==" << dividend_d << endl;
+//                 cout << " divisor_n==" << divisor_n << endl;
+//                 cout << " divisor_d==" << divisor_d << endl;
+//                 cout << string(20,'-') << endl;
+                bool check = divide(dividend_n, divisor_n,
+                                    tmp_n.m[r2*col+c],true);
+                check &= divide(dividend_d, divisor_d,
+                                tmp_d.m[r2*col+c],true);
+                GINAC_ASSERT(check);
             }
+            // fill up left hand side.
             for (unsigned c=0; c<=r1; ++c)
-                this->m[r2*col+c] = _ex0();
+                tmp_n.m[r2*col+c] = _ex0();
         }
+//         cout << tmp_n << endl;
+//         cout << tmp_d << endl;        
     }
+    // repopulate *this matrix:
+    it = m.begin();
+    tmp_n_it = tmp_n.m.begin();
+    tmp_d_it = tmp_d.m.begin();
+    for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
+        (*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl);
     
     return sign;
 }
index a36276f..9742f39 100644 (file)
@@ -97,10 +97,9 @@ public:
 protected:
     ex determinant_numeric(void) const;
     ex determinant_minor(void) const;
-    ex determinant_bareiss(bool normalize=true) const;
     int gauss_elimination(void);
-    int fraction_free_elimination(void);
     int division_free_elimination(void);
+    int fraction_free_elimination(bool det = false);
     int pivot(unsigned ro, bool symbolic=true);
 private:  // FIXME: these should be obsoleted
     void ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2);
index 860fb10..e6d4a61 100644 (file)
@@ -280,21 +280,28 @@ void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons
 
 bool mul::info(unsigned inf) const
 {
-    // TODO: optimize
-    if (inf==info_flags::polynomial ||
-        inf==info_flags::integer_polynomial ||
-        inf==info_flags::cinteger_polynomial ||
-        inf==info_flags::rational_polynomial ||
-        inf==info_flags::crational_polynomial ||
-        inf==info_flags::rational_function) {
-        for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
-            if (!(recombine_pair_to_ex(*it).info(inf)))
-                return false;
+    switch (inf) {
+        case info_flags::polynomial:
+        case info_flags::integer_polynomial:
+        case info_flags::cinteger_polynomial:
+        case info_flags::rational_polynomial:
+        case info_flags::crational_polynomial:
+        case info_flags::rational_function: {
+            for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+                if (!(recombine_pair_to_ex(*i).info(inf)))
+                    return false;
+            }
+            return overall_coeff.info(inf);
+        }
+        case info_flags::algebraic: {
+            for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+                if ((recombine_pair_to_ex(*i).info(inf)))
+                    return true;
+            }
+            return false;
         }
-        return overall_coeff.info(inf);
-    } else {
-        return inherited::info(inf);
     }
+    return inherited::info(inf);
 }
 
 typedef vector<int> intvector;
@@ -419,7 +426,7 @@ ex mul::eval(int level) const
                         ex_to_numeric(addref.overall_coeff).
                         mul_dyn(ex_to_numeric(overall_coeff))))
             ->setflag(status_flags::dynallocated |
-                      status_flags::evaluated );
+                      status_flags::evaluated);
     }
     return this->hold();
 }
index 4e81436..f93ecc8 100644 (file)
@@ -516,7 +516,8 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
         return true;
     }
 #endif
-    if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
+    if (check_args && (!a.info(info_flags::rational_polynomial) ||
+                       !b.info(info_flags::rational_polynomial)))
         throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
 
     // Find first symbol
index 49e65be..3652f4d 100644 (file)
@@ -221,12 +221,12 @@ numeric::numeric(double d) : basic(TINFO_numeric)
 
 
 numeric::numeric(const char *s) : basic(TINFO_numeric)
-{   // MISSING: treatment of complex and ints and rationals.
+{   // MISSING: treatment of complex numbers
     debugmsg("numeric constructor from string",LOGLEVEL_CONSTRUCT);
     if (strchr(s, '.'))
         value = new cl_LF(s);
     else
-        value = new cl_I(s);
+        value = new cl_R(s);
     calchash();
     setflag(status_flags::evaluated|
             status_flags::hash_calculated);
@@ -489,42 +489,44 @@ void numeric::printcsrc(ostream & os, unsigned type, unsigned upper_precedence)
 bool numeric::info(unsigned inf) const
 {
     switch (inf) {
-    case info_flags::numeric:
-    case info_flags::polynomial:
-    case info_flags::rational_function:
-        return true;
-    case info_flags::real:
-        return is_real();
-    case info_flags::rational:
-    case info_flags::rational_polynomial:
-        return is_rational();
-    case info_flags::crational:
-    case info_flags::crational_polynomial:
-        return is_crational();
-    case info_flags::integer:
-    case info_flags::integer_polynomial:
-        return is_integer();
-    case info_flags::cinteger:
-    case info_flags::cinteger_polynomial:
-        return is_cinteger();
-    case info_flags::positive:
-        return is_positive();
-    case info_flags::negative:
-        return is_negative();
-    case info_flags::nonnegative:
-        return !is_negative();
-    case info_flags::posint:
-        return is_pos_integer();
-    case info_flags::negint:
-        return is_integer() && is_negative();
-    case info_flags::nonnegint:
-        return is_nonneg_integer();
-    case info_flags::even:
-        return is_even();
-    case info_flags::odd:
-        return is_odd();
-    case info_flags::prime:
-        return is_prime();
+        case info_flags::numeric:
+        case info_flags::polynomial:
+        case info_flags::rational_function:
+            return true;
+        case info_flags::real:
+            return is_real();
+        case info_flags::rational:
+        case info_flags::rational_polynomial:
+            return is_rational();
+        case info_flags::crational:
+        case info_flags::crational_polynomial:
+            return is_crational();
+        case info_flags::integer:
+        case info_flags::integer_polynomial:
+            return is_integer();
+        case info_flags::cinteger:
+        case info_flags::cinteger_polynomial:
+            return is_cinteger();
+        case info_flags::positive:
+            return is_positive();
+        case info_flags::negative:
+            return is_negative();
+        case info_flags::nonnegative:
+            return !is_negative();
+        case info_flags::posint:
+            return is_pos_integer();
+        case info_flags::negint:
+            return is_integer() && is_negative();
+        case info_flags::nonnegint:
+            return is_nonneg_integer();
+        case info_flags::even:
+            return is_even();
+        case info_flags::odd:
+            return is_odd();
+        case info_flags::prime:
+            return is_prime();
+        case info_flags::algebraic:
+            return !is_real();
     }
     return false;
 }
index fb4edd5..8ec1945 100644 (file)
@@ -250,17 +250,20 @@ void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) co
 
 bool power::info(unsigned inf) const
 {
-    if (inf==info_flags::polynomial ||
-        inf==info_flags::integer_polynomial ||
-        inf==info_flags::cinteger_polynomial ||
-        inf==info_flags::rational_polynomial ||
-        inf==info_flags::crational_polynomial) {
-        return exponent.info(info_flags::nonnegint);
-    } else if (inf==info_flags::rational_function) {
-        return exponent.info(info_flags::integer);
-    } else {
-        return inherited::info(inf);
+    switch (inf) {
+        case info_flags::polynomial:
+        case info_flags::integer_polynomial:
+        case info_flags::cinteger_polynomial:
+        case info_flags::rational_polynomial:
+        case info_flags::crational_polynomial:
+            return exponent.info(info_flags::nonnegint);
+        case info_flags::rational_function:
+            return exponent.info(info_flags::integer);
+        case info_flags::algebraic:
+            return (!exponent.info(info_flags::integer) ||
+                    basis.info(inf));
     }
+    return inherited::info(inf);
 }
 
 unsigned power::nops() const