Wishlist (relational, pseries, matrix)

[Richard B. Kreckel]

Hi,

On Tue, 24 Apr 2001, Pearu Peterson wrote:
> 1) Having relational::get_operator() would be useful (read: necessary) for
>    extensions that can perform symbolic operations like
>    (a < b) && (b > a)    ->    a < b

There are currently no intentions to target such operations within
GiNaC and I don't know anything about the algorithms involved once you get
to the nontrivial cases.  Be warned: RJF has recently informed  somebody
on sci.math.symbolic that this is still an ongoing field of research.  But
for cases of interfacing like yours, yes, we should add something.  Patch
suggestion?  Maybe relational should even be sub-classed???

> 2) If a,b are compatible pseries then
>       ex(a) * ex(b)
>    would automatically produce
>       ex(a.mul_series(b))
>    The same for operations +,-,scalar_mul,scalar_div.
>    If a,b are incompatible, then apply current
>    behaviour (that is, do nothing).
> 
>    I have implemented this behaviour in the GiNaC->Python interface but
>    it does not cover cases where multiplication is carried out inside
>    GiNaC functions. I am thinking of applications where, for example, 
>    a matrix consists of compatible pseries and the task is to find its
>    inverse.

To do this in an automated way would probably become quiet involved from
the complexity point of view.  What is so wrong about reexpanding the
pseries?  What other cases are there where this happens internally?  In
the matrix case: can't you just iterate through it and reexpand the
elements and assemble a new matrix in case there really _are_ series
objects involved?

> 3) If a,b are compatible matrices (for a given operation) then
>      ex(a) * ex(b)   ->  ex(a.mul(b))
>    The same for operations +,-,scalar_mul,scalar_div.

Matrices as derived from basic obey the equivalence relation established
by the .compare().  This has some funny consequences.  Let A be any n x m
matrix.  Let B be the same n x m matrix.  Compute A-B.  You'll get 0, yet
not the n x m zero matrix but a plain scalar ex(numeric(0)) by the rules
of add::eval().  Incidentally, this happens in Maple too more or less,
where for any matrix A you have A - A   ->  0.

Since multiplication of matrices is O(n^x) where 2<x<3 we had decided a
while ago that the rule you are suggesting should not be applied but
matrices should be left as unevaluated until somebody triggers full matrix
evaluation, maybe by calling a (not yet present) ex::evalm(). This may be
somewhat confusing for the novice but should really be considered a
feature since it allows us to produce an n x n matrix A, later produce B =
pow(A,-1) and then A*B will evaluate to 1 without expensive matrix
inversion going on under the hood -- unless of course one says B =
pow(A,-1).evalm() and then A*B.

Regards
     -richy.
-- 
Richard Kreckel
<Richard.Kreckel at Uni-Mainz.DE>
<http://wwwthep.physik.uni-mainz.de/~kreckel/>