Mathemagix

[Joris van der Hoeven <address@concealed>]

Hi all,

There are a few design issues with the Mathemagix libraries
which continue to bother me. Maybe discussing them with you,
would help me to make the appropriate decisions.

First of all, there is the question whether the standard bracket
notation [1,2,3] should rather be used for lists, arrays or vectors.

Currently, we use the notation for lists. However, this is not very
nice if we want to use the compact notation [1,2;3,4] for matrices.
Also, it has turned out that at the typical places where one might
expect lists to be used (e.g. the Mathemagix compiler), the map
primitive has proven to be more useful than recursive programming.
Finally, vectors are better than lists when raw performance is required,
and also better in the case of parallelism. For these reasons, I am
currently inclined to use [1,2,3] to denote vectors rather than lists.

Now there is a second subtlety with arrays and vectors. Mathemagix
currently implements both arrays and vectors. The difference is that
vectors can be customized using variants and that they support
a constructor C -> Vector C, which also accounts for some overhead.
Indeed, using Vector C as a ring requires at least the existence
of a homomorphism Integer -> Vector C, for which the size of
the vector is not yet known (the same holds for matrices by the way).

Currently, this is implemented by regarding Vector C as some kind
of union of C and Vector C, using a special flag "bool scalar_flag"
which indicates whether the vector is in fact in C or a real vector.
Of course, this accounts for some overhead for most operations,
which does not exist for arrays. Do we want to keep separate types
array and vector, or try to modify vector in such a way that array
becomes redundant and the overhead disappears. For instance,
one might decide that a vector of size 1 behaves like a scalar
when added to a vector of another size and only undertake special
action if the sizes of operands don't coincide. We might also use
negative sizes for scalar vectors so as to save some space.

A second main and recurrent concern is the way to customize
implementations, such as the use of different implementations
of Polynomial C (dense, sparse, Bernstein). In synaps this is done
by templating containers over namespaces. In mathemagix, we have
so far used the concept of a special template argument called the variant,
using which a class can be customized a posteriori.

In the Mathemagix language, the ideal solution would be to
declare the correctness assumptions for every routine.
For instance, if we have a routine foo which does something
interesting on polynomials, we might write

        forall (C: Ring, P: Polynomial_ring C) foo (p: P): P == ...

We see here that we *explicitly* have to abstract the routine
so as to work for all kinds of polynomial rings instead of
a particular one. This is more or less equivalent to specifying
a variant template parameter, except that the customization has
not to be foreseen in any of the implementations of polynomials
themselves. Rather, the signatures of the different implementations
have to match Polynomial_ring (and wrappers might be added to force this).
Of course, in C++, we might have written

        template<class C,class P> P foo (const P& p) { ... }

or

        template<class P> P foo (const P& p) { ... }

where C equals P::C_type for some typedef inside P.
This leads me to the opinion that we probably should do something
similar in C++: never try to "prepare the ground" for future
customization of some class, but rather allow for the customization
of the implementation itself.

In fact, it seems that the variant mechanism still remains useful,
but only for the customization of underlying implementation issues,
much in the same way as synaps specifies namespaces for the underlying
arithmetic. Notice that having classes as template arguments is nicer
than having namespaces, because classes can be templated themselves.
For instance, consider

        template<class C, class V= dense_polynomial_variant<C> >
        class dense_polynomial { ... };

In this example, the default underlying arithmetic for dense_polynomial<C>
naturally depends on C.

An advantage of namespaces is that they may be extended a posteriori.
However, in view what has been said above, we do not really want to do this:
strict application of the philosophy would rather mean that the variant
(class or namespace) should only be used in order to customize
the implementation in consideration (dense_polynomial in the above example).
Indeed, V corresponds to mathematical or implementation hypothesis
under which the implementation of a new function or class is valid.
Such hypothesis should never be extended: instead a new hypothesis class
should be created.

Of course, carrying out this strategy is somewhat verbose (thank you C++).
Mathemagix style factorizations:

        forall (C: Ring, Pol: Polynomial_ring C) {
           foo1 ...
           foo2 ...
        }

can be mimicked using the usual #define #undef cruft:

        #define TMPL template<class Pol>
        #define C (typename Pol::scalar_type)

        TMPL foo1 ...
        TMPL foo2 ...

        #undef C
        #undef TMPL

or by regrouping things in a class with static member functions.

Summarizing, what about the following general practice:

  1) For types whose implementations depend on lower layers,
     use a variant class to select particular implementations
     of the lower layers.

  2) Never use variant classes to factor different implementations
     of the same kind of object; rather give the different implementations
     different names and template over the implementation when using
     them in an even higher level.

Another benefit of this practice might be that it is actually quite
straightforward to upgrade an existing code for, say, polynomials
to work with other implementations: it suffices to systematically
use templates and replace polynomial<C> by a template class Pol.

Best wishes, Joris