- ++version_major.

- added matrix::pow() to handle integer exponents with the least amount of
  multiplications possible and...
- ...added power::evalm() to actually use it.
- classhierarchy.fig: add class wildcard.
- some cleanups.

diff --git a/COPYING b/COPYING
index 12065f35cb96440fa294b2c0d202a2254b76d322..d406b2738a5dfb2e9dba704c242a0ee739fd7b88 100644 (file)
--- a/COPYING
+++ b/COPYING
@@ -292,7 +292,7 @@ convey the exclusion of warranty; and each file should have at least
 the "copyright" line and a pointer to where the full notice is found.
 
     <one line to give the program's name and a brief idea of what it does.>
-    Copyright (C) 19yy  <name of author>
+    Copyright (C) <year>  <name of author>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -313,7 +313,7 @@ Also add information on how to contact you by electronic and paper mail.
 If the program is interactive, make it output a short notice like this
 when it starts in an interactive mode:
 
-    Gnomovision version 69, Copyright (C) 19yy name of author
+    Gnomovision version 69, Copyright (C) year  name of author
     Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
     This is free software, and you are welcome to redistribute it
     under certain conditions; type `show c' for details.

diff --git a/NEWS b/NEWS
index d5d29512a9ccdf74aaebe4d2f1eca3860781982a..af7ec2202cb0764b49703ca437aea3fa6b13e996 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
 This file records noteworthy changes.
 
-0.8.4 (<date>)
+0.9.0 (<date>)
 * In the output and in ginsh, lists are now delimited by { } braces, and
   matrices are delimited by single [ ] brackets.
 * simplify_indexed() renames dummy indices so, e.g., "a.i*a.i+a.j*a.j" gets
@@ -21,6 +21,7 @@ This file records noteworthy changes.
   products.
 * Positive integer powers of non-commutative expressions are automatically
   expanded.
+* Several little bugfixes.
 
 0.8.3 (11 May 2001)
 * color and clifford classes are functional and documented.

diff --git a/configure.in b/configure.in
index c4f45821d1a9892350da892a6d0e24f794e9e359..0bf0cc21f1c0f5704dee6e20e3525a85a1ea7d16 100644 (file)
--- a/configure.in
+++ b/configure.in
@@ -18,8 +18,8 @@ dnl autoconf sees "AC_MAJOR_VERSION" and complains about an undefined macro
 dnl (don't we all *love* M4?)...
 
 GINACLIB_MAJOR_VERSION=0
-GINACLIB_MINOR_VERSION=8
-GINACLIB_MICRO_VERSION=3
+GINACLIB_MINOR_VERSION=9
+GINACLIB_MICRO_VERSION=0
 GINACLIB_INTERFACE_AGE=0
 GINACLIB_BINARY_AGE=0
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION

diff --git a/doc/tutorial/classhierarchy.fig b/doc/tutorial/classhierarchy.fig
index 0a3ab0e0de1e4069573c57ef94687346eaf46838..936dc1589de56a20978b513ec1ba57fb145d2d63 100644 (file)
--- a/doc/tutorial/classhierarchy.fig
+++ b/doc/tutorial/classhierarchy.fig
@@ -9,13 +9,6 @@ Single
 1200 2
 5 1 1 1 0 7 50 0 -1 4.000 0 0 1 0 1912.500 4102.500 945 720 1980 585 2880 720
        1 0 1.00 60.00 120.00
-6 5580 1080 5895 2025
-2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
-        5850 1080 5850 1980 5580 1980 5580 1080 5850 1080
-2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
-        5895 1125 5895 2025 5625 2025 5625 1125 5895 1125
-4 1 0 50 0 14 10 4.7124 4 135 540 5670 1530 symbol\001
--6
 6 675 1125 945 2115
 2 2 0 1 0 7 51 0 20 0.000 0 0 -1 0 0 5
         675 2115 675 1125 945 1125 945 2115 675 2115
@@ -49,20 +42,6 @@ Single
         4365 2250 4635 2250 4635 1485 4365 1485 4365 2250
 4 0 0 50 0 14 10 4.7124 4 105 540 4410 1575 matrix\001
 -6
-6 4770 1260 5085 2205
-2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
-        5085 1305 5085 2205 4815 2205 4815 1305 5085 1305
-2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
-        5040 1260 5040 2160 4770 2160 4770 1260 5040 1260
-4 1 0 50 0 14 10 4.7124 4 105 630 4860 1710 numeric\001
--6
-6 5175 1170 5490 2115
-2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
-        5445 1170 5445 2070 5175 2070 5175 1170 5445 1170
-2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
-        5490 1215 5490 2115 5220 2115 5220 1215 5490 1215
-4 1 0 50 0 14 10 4.7124 4 105 720 5265 1620 constant\001
--6
 6 2430 1215 2700 2025
 2 2 4 1 0 7 51 0 20 4.000 0 0 -1 0 0 5
         2430 2025 2700 2025 2700 1215 2430 1215 2430 2025
@@ -226,6 +205,34 @@ Single
         3870 3825 3870 4725 3600 4725 3600 3825 3870 3825
 4 1 0 50 0 14 10 4.7124 4 105 810 3690 4275 tensdelta\001
 -6
+6 5580 1170 5895 2115
+2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
+        5850 1170 5850 2070 5580 2070 5580 1170 5850 1170
+2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
+        5895 1215 5895 2115 5625 2115 5625 1215 5895 1215
+4 1 0 50 0 14 10 4.7124 4 105 720 5670 1620 constant\001
+-6
+6 5985 1080 6300 2025
+2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
+        6255 1080 6255 1980 5985 1980 5985 1080 6255 1080
+2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
+        6300 1125 6300 2025 6030 2025 6030 1125 6300 1125
+4 1 0 50 0 14 10 4.7124 4 135 540 6075 1530 symbol\001
+-6
+6 5175 1260 5490 2205
+2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
+        5490 1305 5490 2205 5220 2205 5220 1305 5490 1305
+2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
+        5445 1260 5445 2160 5175 2160 5175 1260 5445 1260
+4 1 0 50 0 14 10 4.7124 4 105 630 5265 1710 numeric\001
+-6
+6 4770 1350 5085 2295
+2 4 0 0 0 7 52 0 10 0.000 0 0 7 0 0 5
+        5085 1395 5085 2295 4815 2295 4815 1395 5085 1395
+2 4 0 1 0 7 51 0 20 0.000 0 0 7 0 0 5
+        5040 1350 5040 2250 4770 2250 4770 1350 5040 1350
+4 1 0 50 0 14 10 4.7124 4 105 720 4860 1800 wildcard\001
+-6
 2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
        1 1 1.00 45.00 90.00
         810 1125 2835 855
@@ -237,13 +244,13 @@ Single
         1800 1260 3015 900
 2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
        1 1 1.00 45.00 90.00
-        5625 1080 3645 810
+        5580 1215 3645 810
 2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
        1 1 1.00 45.00 90.00
-        5220 1215 3645 855
+        5220 1260 3645 855
 2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
        1 1 1.00 45.00 90.00
-        4815 1260 3600 900
+        4815 1350 3600 900
 2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
        1 1 1.00 45.00 90.00
         4005 1485 3465 900
@@ -291,5 +298,8 @@ Single
 2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
        1 1 1.00 45.00 90.00
         3375 3015 3330 2745
+2 1 0 1 0 7 50 0 10 0.000 0 0 -1 1 0 2
+       1 1 1.00 45.00 90.00
+        6030 1125 3650 740
 4 0 0 50 0 0 16 0.0000 4 30 180 3150 4275 ...\001
 4 0 0 50 0 0 16 0.0000 4 30 180 3960 4365 ...\001

diff --git a/ginac/ex.cpp b/ginac/ex.cpp
index 1849789871585a3d96318454d8f49e86769bae2c..b8a35d3a67c78b98843906ae17b14cda9627c0a5 100644 (file)
--- a/ginac/ex.cpp
+++ b/ginac/ex.cpp
@@ -207,14 +207,12 @@ void ex::makewriteable()
  *  @see ex::ex(const basic &) */
 void ex::construct_from_basic(const basic & other)
 {
-       if ((other.flags & status_flags::evaluated)==0) {
-               // cf. copy ctor
-               const ex & tmpex = other.eval(1); // evaluate only one (top) level
+       if (!(other.flags & status_flags::evaluated)) {
+               const ex & tmpex(other.eval(1)); // evaluate only one (top) level
                bp = tmpex.bp;
-               GINAC_ASSERT(bp!=0);
                GINAC_ASSERT(bp->flags & status_flags::dynallocated);
                ++bp->refcount;
-               if ((other.flags & status_flags::dynallocated)&&(other.refcount==0))
+               if ((other.refcount==0) && (other.flags & status_flags::dynallocated))
                        delete &const_cast<basic &>(other);
        } else {
                if (other.flags & status_flags::dynallocated) {

diff --git a/ginac/function.pl b/ginac/function.pl
index b0564e6e45f105422029acccdf6c108a95323979..e2287091db4391ccac0daa8ceb1206251d3cc652 100755 (executable)
--- a/ginac/function.pl
+++ b/ginac/function.pl
@@ -678,7 +678,7 @@ void function::print(const print_context & c, unsigned level) const
                }
                c.s << ")";
 
-       } else if is_of_type(c, print_latex) {
+       } else if (is_of_type(c, print_latex)) {
                c.s << registered_functions()[serial].TeX_name;
                printseq(c, '(', ',', ')', exprseq::precedence(), function::precedence());
        } else {

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 5767d4cfb924f91e4ec5adc7c43a59cecc409f05..28c2f1d25a853faf7750220d1a82867d9f425775 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -534,7 +534,7 @@ bool matrix::contract_with(exvector::iterator self, exvector::iterator other, ex
 matrix matrix::add(const matrix & other) const
 {
        if (col != other.col || row != other.row)
-               throw (std::logic_error("matrix::add(): incompatible matrices"));
+               throw std::logic_error("matrix::add(): incompatible matrices");
        
        exvector sum(this->m);
        exvector::iterator i;
@@ -552,7 +552,7 @@ matrix matrix::add(const matrix & other) const
 matrix matrix::sub(const matrix & other) const
 {
        if (col != other.col || row != other.row)
-               throw (std::logic_error("matrix::sub(): incompatible matrices"));
+               throw std::logic_error("matrix::sub(): incompatible matrices");
        
        exvector dif(this->m);
        exvector::iterator i;
@@ -570,7 +570,7 @@ matrix matrix::sub(const matrix & other) const
 matrix matrix::mul(const matrix & other) const
 {
        if (this->cols() != other.rows())
-               throw (std::logic_error("matrix::mul(): incompatible matrices"));
+               throw std::logic_error("matrix::mul(): incompatible matrices");
        
        exvector prod(this->rows()*other.cols());
        
@@ -602,6 +602,9 @@ matrix matrix::mul(const numeric & other) const
 /** Product of matrix and scalar expression. */
 matrix matrix::mul_scalar(const ex & other) const
 {
+       if (other.return_type() != return_types::commutative)
+               throw std::runtime_error("matrix::mul_scalar(): non-commutative scalar");
+
        exvector prod(row * col);
 
        for (unsigned r=0; r<row; ++r)
@@ -612,6 +615,47 @@ matrix matrix::mul_scalar(const ex & other) const
 }
 
 
+/** Power of a matrix.  Currently handles integer exponents only. */
+matrix matrix::pow(const ex & expn) const
+{
+       if (col!=row)
+               throw (std::logic_error("matrix::pow(): matrix not square"));
+       
+       if (is_ex_exactly_of_type(expn, numeric)) {
+               // Integer cases are computed by successive multiplication, using the
+               // obvious shortcut of storing temporaries, like A^4 == (A*A)*(A*A).
+               if (expn.info(info_flags::integer)) {
+                       numeric k;
+                       matrix prod(row,col);
+                       if (expn.info(info_flags::negative)) {
+                               k = -ex_to_numeric(expn);
+                               prod = this->inverse();
+                       } else {
+                               k = ex_to_numeric(expn);
+                               prod = *this;
+                       }
+                       matrix result(row,col);
+                       for (unsigned r=0; r<row; ++r)
+                               result.set(r,r,_ex1());
+                       numeric b(1);
+                       // this loop computes the representation of k in base 2 and multiplies
+                       // the factors whenever needed:
+                       while (b.compare(k)<=0) {
+                               b *= numeric(2);
+                               numeric r(mod(k,b));
+                               if (!r.is_zero()) {
+                                       k -= r;
+                                       result = result.mul(prod);
+                               }
+                               prod = prod.mul(prod);
+                       }
+                       return result;
+               }
+       }
+       throw (std::runtime_error("matrix::pow(): don't know how to handle exponent"));
+}
+
+
 /** operator() to access elements.
  *
  *  @param ro row of element
@@ -633,6 +677,8 @@ matrix & matrix::set(unsigned ro, unsigned co, ex value)
 {
        if (ro>=row || co>=col)
                throw (std::range_error("matrix::set(): index out of range"));
+       if (value.return_type() != return_types::commutative)
+               throw std::runtime_error("matrix::set(): non-commutative argument");
     
        ensure_if_modifiable();
        m[ro*col+co] = value;

diff --git a/ginac/matrix.h b/ginac/matrix.h
index 08f99edea92015177dd14beaed3182ffb3ef9b3e..4a267885a8ef0bc85721661328e2ace1deb257ee 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -71,6 +71,7 @@ public:
        matrix mul(const matrix & other) const;
        matrix mul(const numeric & other) const;
        matrix mul_scalar(const ex & other) const;
+       matrix pow(const ex & expn) const;
        const ex & operator() (unsigned ro, unsigned co) const;
        matrix & set(unsigned ro, unsigned co, ex value);
        matrix transpose(void) const;

diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index a4bb802ec9475c6b8f87ac4becbeaaf22b13cea5..cb6794323b8cb96911b40e3dc14d2b125d3edd6e 100644 (file)
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -144,7 +144,7 @@ void mul::print(const print_context & c, unsigned level) const
                        overall_coeff.bp->print(c, precedence());
                        c.s << "*";
                }
-       
+
                // Print arguments, separated by "*" or "/"
                epvector::const_iterator it = seq.begin(), itend = seq.end();
                while (it != itend) {

diff --git a/ginac/numeric.h b/ginac/numeric.h
index 1cd7d976dfaefd1eca9f3b6d029a1c764065813e..4c7fd0b33162ab7572f8de42a75e7ee0d3747e83 100644 (file)
--- a/ginac/numeric.h
+++ b/ginac/numeric.h
@@ -200,7 +200,6 @@ const numeric doublefactorial(const numeric &n);
 const numeric binomial(const numeric &n, const numeric &k);
 const numeric bernoulli(const numeric &n);
 const numeric fibonacci(const numeric &n);
-const numeric abs(const numeric &x);
 const numeric isqrt(const numeric &x);
 const numeric sqrt(const numeric &x);
 const numeric abs(const numeric &x);

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 11fcbffeebbfc3a5ae2180ef7156c67625e365f7..eefd438f9bb334415694e313bf0018e0bae48df7 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -31,6 +31,7 @@
 #include "ncmul.h"
 #include "numeric.h"
 #include "inifcns.h"
+#include "matrix.h"
 #include "symbol.h"
 #include "print.h"
 #include "archive.h"
@@ -450,8 +451,10 @@ ex power::eval(int level) const
                        }
                }
 
-               // ^(nc,c1) -> ncmul(nc,nc,...) (c1 positive integer)
-               if (num_exponent->is_pos_integer() && ebasis.return_type() != return_types::commutative) {
+               // ^(nc,c1) -> ncmul(nc,nc,...) (c1 positive integer, unless nc is a matrix)
+               if (num_exponent->is_pos_integer() &&
+                   ebasis.return_type() != return_types::commutative &&
+                   !is_ex_of_type(ebasis,matrix)) {
                        return ncmul(exvector(num_exponent->to_int(), ebasis), true);
                }
        }
@@ -487,6 +490,18 @@ ex power::evalf(int level) const
        return power(ebasis,eexponent);
 }
 
+ex power::evalm(void) const
+{
+       ex ebasis = basis.evalm();
+       ex eexponent = exponent.evalm();
+       if (is_ex_of_type(ebasis,matrix)) {
+               if (is_ex_of_type(eexponent,numeric)) {
+                       return (new matrix(ex_to_matrix(ebasis).pow(eexponent)))->setflag(status_flags::dynallocated);
+               }
+       }
+       return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
+}
+
 ex power::subs(const lst & ls, const lst & lr, bool no_pattern) const
 {
        const ex &subsed_basis = basis.subs(ls, lr, no_pattern);

diff --git a/ginac/power.h b/ginac/power.h
index 00f02b1f09f1dd7b3b706c839dc7ad1f3aa2f60d..0a156545d1367802cd393e493b685e9d54106804 100644 (file)
--- a/ginac/power.h
+++ b/ginac/power.h
@@ -59,6 +59,7 @@ public:
        ex coeff(const ex & s, int n = 1) const;
        ex eval(int level=0) const;
        ex evalf(int level=0) const;
+       ex evalm(void) const;
        ex series(const relational & s, int order, unsigned options = 0) const;
        ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const;
        ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;