Re: [Axiom-developer] Axiom domains and Aldor return types

[William Sit]

It appears that the previous file may have some display alignment problems (in
the code indentations). So here is hopefully a better version. My apologies for
duplicates.

William
-----


Hi Steve:

Your example is very interesting and it took me quite some time to
understand it (I am slow in learning and I tossed and turned several
times between endorsing it to rejecting it).  My tentative short
answers are that (1) Yes, this can be implemented in current Axiom and
(2) There are big hidden problems (see (II, III) below) in
implementing the Aldor version!

But: there is no such thing as a simple explanation.

So please bear with me.  To avoid getting into residue class rings,
let me simplify and abstract the hypothesis under which the example
makes sense.  If you do want the analysis for residue class ring (and
on how I arrive at the above conclusions), just keep reading.

Given a certain domain R of category A, and another domain P of
category B, let's assume that there is at least one uniform way of
constructing, for each p:P, a new domain S (belonging to some category
C that depends on R and p).  For a given R, there may be several such
uniform methods to construct the same mathematical object S; and
certainly for different R's, there would be other ways of constructing
S.  For example, if R' (of category A', which is a subcategory of A)
has some special properties, the construction may be different.  In
Axiom, each construction of S requires a domain constructor:

   Method1(R:A, p:P)==S:C(R,p)
   Method2(R:A, p:P)==S:C(R,p)
   Method3(R':A, p:P)==S:C(R',p)
   Method4(R':A', p:P)==S:C(R',p)

In Axiom, because the signatures of the first three methods would be
the same, the domain constructors need different names.

Now suppose we want to implement an algorithm for S and this is to be
done categorically, independent of the actual construction of S, but
depend on a function foo which is defined for all domains in the
category C(R,p).  Currently, in Axiom, this would be a package:

   Algorithm(R:A, p:P, S:C(R,p)) == add
        bar(...) == ... foo(...)$S ...

and a typical calling sequence would be:

   bar(...)$Algorithm(R,p, Method1(R,p))

Now let's see how we may do the same thing a la Steve's example in
Aldor.  We would create a new category to encapsulate the various
methods for constructing S.

ComputationMethod: Category == Join(...) with
  method:(p:P) -> C(%,p)

and the algorithm for S would be in a package:

  Algorithm(R:Join(A,ComputationMethod)) ==
      add { bar(p:P, ...)== ... foo(...)$method(p) ...}

and a typical calling sequence would be:

  bar(p:P, ...)$Algorithm(R)

There are several problems:

(I) The algorithm, which is really for S, is now an algorithm for R.
It hides the existence of S completely.  It also hides the relevance
of p:P.

(II) For each R:A for which a method to construct S for p:P exist, we
must *retroactively* make R:ComputationMethod.  This requires
modifying the original construction of R.  If the construction of R is
from a domain constructor T, say R := T(...) where

       T(...):A == ...

and if Method1 is used to construct S, we must now modify this to

      T'(...):Join(A, ComputationMethod) == T(...) add
         method(p:P) == Method1(%,p)

We must create a new constructor T' and not modify T because not every
T(...) is to become a domain in ComputationMethod.  So even though we
eliminated one domain constructor (Method1, assuming this is inlined),
we need one new domain constructor to "wrap" it.  If T is actually the
construction of R from scratch (that is, R is T, and (...) is empty;
such as Integer), then inlining retroactively Method1 would require
recompiling the whole world.

(III) Now if R has several methods, then we must *rename* R using
three different names in (II), even though R is still the same object
in category A.  This would be very undesirable if R is Integer.

[Remark:  It seems to me the original algebra developers in Axiom
avoided at all cost to modify existing domains because recompiling the
world took a long time in those days.  They typically added packages
and extended domains to cover any shortcoming in the original
implementations.]

(IV) The "advantage" at first glance, seems to be operator
overloading:  using one single name "method" to refer to all the
Method[i] in constructions of S.  But as (II) shows, this "powerful"
and categorical construction is nothing but a wrapper, which is the
the reverse of what I proposed for Martin's example:  instead of
lifting a function g(n,k) to the level of a package, where the
dependence on parameters in the signature belongs, the construction of
ComputationMethod pushed the level of a domain constructor (which is
what each Method[i] is) to the level of a function.  I don't think
that is a convincing example of the need for allowing dependence on
parameters for functions.

*************************************************************
*  It is clear that the two set ups are equivalent and
*  the translation is bidirectional. The Axiom implementation
*  is more natural and does not have the disadvantages.
*************************************************************
The power of that construction in Aldor, as I pointed out in another
email, is allowing such dependence in functions on the SOURCE side,
not on the TARGET side of the map. In short notation:

       F(a:A, b:B(a), c:C(a,b)):D(a,b,c)

is a powerful signature, not because the a, b, c appears on the Target
D, but because they appear in other parameters on the Source side. If
A, B, C represent three axes in 3-space, then F is analogous to an
iterated triple integral over a non-rectangular region in 3-space,
whereas

       F(a:A, b:B, c:C):D(a,b,c)

is like the same over a rectangular region (a cuboid) in 3-space.

If we may borrow the way Matlab does these computation (Matlab, for
Matrix Laboratory, is constrained in these numerical triple
integration computations because it must define grid points in 3-space
in a 3-dimensional matrix), we have to "zero" out the complement of
the conditionals b:B(a) and c:C(a,b).  This would be difficult without
an over category B' and C' of which B(a) and C(a,b) are subcategories.
But we can always create such over categories and it would not be
difficult then to restrict the inputs by using "if b has B(a)" and "if
c has C(a,b)".  In fact, Axiom does this all the time with
non-parametrized conditionals like "if R has CharacteristicZero".
Even though I can't recall any attributes being parametrized, there is
no reason why Axiom cannot support such.  If I understand correctly,
attributes are defined as Categories.  Let me leave that as an
exercise for someone to find out.  :-).  I am more convinced now that
the signature limitation in Axiom does not actually limit its power,
but does restrict its freedom of expression in some situations.

> How could one write the bar()$Foo above with the current axiom
> language? All you can do is write a package/domain parameterize
> by a commutative ring R and a representative of R, and write the
> exports dependently on R's type via `if R has ...'  constructs.
> This is what I meant in the previous email about having a
> RESCLASS  package/domain which needed to know too much about the
> algebra.

In Axiom, as well as in Aldor, the conditional "if R has ..." is never
meant to be algorithmic.  It is a declarative that a domain has a
certain attribute.  Such a declarative is never verified by code.  It
is one of the "trust me" situations.  We do NOT need to know about
every ring when we use these conditionals, because for whatever actual
domain R is inputted, the *programmer* will have to declare his/her
knowledge about these conditionals in the domain constructor for R, as
in, for example, "Join(CharacteristicZero, ..." There is no algorithm
to test if an arbitrary ring has characteristic zero.  One just knows
from the mathematics of particular rings.  You used these conditionals
also in defining the ResidueClassRing(R,p) category and it does not
mean you know which commutative ring R has SourceOfPrimes or has
implemented prime?(p).  But for the rings R for which you do know, you
declare them to have SourceOfPrimes, and you implement prime?(p) in
the domain constructor that gives R.

Below, I'll give more explanation how I arrive at my conclusion above.
Stanzas marked between ==== and **** are mine (and those lines between
**** and ==== are Steve's).  In each stanza, I put down a paraphrase
of Steve's code (in a combination of code, math and English, with
added background information), and an analog taken from Axiom (which
helped me notice what the map residueClassRing really is --- it became
obvious only after such analysis).  The way to follow my analysis
would be to read the paraphrase from top down, then read the analogy
from bottom up (if you don't, you may find some notations not yet
defined).  Then read the above abstract discussion again.

I skipped the Axiom constructions for all the constructors in Steve's
example, except for Foo.  I think it is clear from the more general
discussion above how to complete the conversion.

By the way, the analogy would fit the abstract situation above too.  R
is a ring, P is List Symbol, S is the polynomial ring R[p], and the
Algorithm is GroebnerBasis.  The methods are DMP, HDMP, GDMP.

William
------
Steve Wilson wrote:

The following is an example with a view towards generic modular
computations. Aldor has a category (approximately):

---------
 ModularComputation: Category == CommutativeRing with {
         residueClassRing: (p: %) -> ResidueClassRing(%,p);
                 ....
          }
---------
============

Paraphrase:

A ModularComputation domain R is a commutative ring with a map

   residueClassRing: R -> ResidueClassRing category

In such a domain R, there is an algorithm to construct the residue
class ring for any prime ideal p in R.  Mathematically, the residue
class ring of R with respect to a prime ideal p is (R_p)/(p R_p),
where R_p is the localization of R at the prime ideal p.  When lifted
to R_p, the prime ideal p becomes p R_p, the unique maximal ideal of
R_p (which is called a local ring).  The residue class ring is then
formed by the "modding out" the maximal ideal.

The notation p:% above is technically incorrect and should be
something like p:  Ideals(%), and for R in the IntegerCategory
(below), there is a coercion from a prime integer p to the prime ideal
(p), at least in the envisioned situation.  There are other rings,
typically polynomial rings, where there is an algorithm to detect
prime ideals (using a primary decomposition algorithm and for these
rings, the residue class ring can also be constructed).

Analogy:

PolynomialComputation: Category == CommutativeRing with {

   polyomialRing: (v:List Symbol) -> POLYCAT(v, %)

There is really no need for the map polynomialRing because it
encapsulates a domain constructor, which must be implemented for each
instance of v and R.  Examples of this map polynomialRing in action
are:

   DMP(v,R): POLYCAT(v,R)
   HDMP(v,R): POLYCAT(v,R)

Each of these domain constructor is equivalent to an instance of the
map polynomialRing and therefore, a domain of the category
PolynomialComputation.

************


So any domain satisfying Modular computation is a CommutativeRing R,
which exports a function which takes a representative p of R and
returns something which satisfies ResidueClassRing(R,p).

---------
ResidueClassRing(R: CommutativeRing, p: R): Category ==
   CommutativeRing with {
      modularRep: R -> %;
      canonicalPreImage: % -> R;
      if R has EuclideanDomain then {
         symmetricPreImage: % -> R;
         if R has SourceOfPrimes and prime?(p)
         then Field;
     } }
---------
=========

Paraphrase:

A ResideClassRing S is constructed from a base commutative ring R and
a prime ideal p of R. It is a commutative ring with these operations:

  modularRep : R -> S
    canonicalPreImage: S -> R
      if R is a Euclidean domain, then there are more operations:
          symmetricPreImage: S -> R
          if we know the prime ideals of R, and p is a prime ideal,
          then S is a field

ResidueClassRing is a category constructor because there may be
several representations of S, given one R and one prime ideal p of R.
In order to perform this construction, we would need to have
constructed already Ideal(R), the domain of ideals of R, and a
function that can decide whether a given ideal is prime or not.  The
localization construction is already in Axiom (the FRAC constructor is
a special case where the prime ideal is (0)).  The modulo operation is
available for polynomial rings using Groebner basis method (in
PolynomialIdeal), and of course also for Integer using plain old
division.  These two are the most important examples.  For this
discussion, the if-clauses are not relevant.

Analogy:

POLYCAT(v: List Symbol, R: CommutativeRing):Category ==
  CommutativeRing with {
      coerce: R -> %;
      retract: % -> R;
        ...
                 }

**********


Here we use the notion of SourceOfPrimes until someone figures out a
meaningful way to represent a MaximalIdeal generally:).

Aldor has an IntegerCategory, roughly:

---------
IntegerCategory: Category ==
    Join(IntegerType, CharacteristicZero, EuclideanDomain,
                 ModularComputation, SourceOfPrimes,
                              GeneralExponentCategory, Specializable,
Parsable) with {
       ...
       default {
          residueClassRing(p:%):ResidueClassRing(%, p) ==
            IntegerMod(%,p);
       ...
          } }
---------
=========

Paraphrase

A domain R of category IntegerCategory is a Euclidean domain of
characteristic zero, etc., where we know the prime ideals, and we know
how to construct a ResidueClassRing S given any prime ideal p in R.

A default way to construct S is via IntegerMod(R, p) when R is an
IntegerCategory domain.  The construction IntegerMod(R, p) is assumed
known and efficient.  This default construction does not really apply
always, for example, when R is a polynomial ring.  But this is outside
the scope of the current discussion.

Analogy:

PolynomialConstructable: Category ==
  Join( ...) with {
    default {
        polynomialRing(v:List Symbol): POLYCAT(v,%) == DMP(v, %);
            ...
              }  }


**********

And IntegerMod is an efficient implementation:

---------
   IntegerMod(Z:IntegerCategory, p:Z):ResidueClassRing(Z, p) == add {
... }
---------

=========
Paraphrase:

IntegerMod is a domain constructor that is actually implementable
because we supposedly know how to construct the residue class ring for
a domain R (or Z) of the IntegerCategory and a prime ideal p of R.
IntegerMod represents ONE way of construction for S:= (R_p)/(p R_p).

Analogy:

DMP (DistributedMultivariatePolynomial) is a domain constructor in
Axiom that implements a polynomial ring with coefficient ring R and a
list of indeterminates v.  It uses the pure lexicographic ordering on
monomials.

   DMP(v:List Symbol, R: Ring): POLYCAT(v,R) == Join(...) add ...

It may serve as the default constructor.  Other constructors use other
term orderings, for example HDMP or GDMP.

Here POLYCAT(v,R) is a specialization of an actual category
constructor from Axiom, except I have abbreviated the parameter set in
this analogy.  The full macro expansion is

    POLYCAT(v,R) ==> PolynomialCategory(R,_
        Direct Product(#(v),NonNegativeInteger),_
           OrderedVariableList(v))


Axiom Version:

   IntegerMod(R, p): T == C where
     R: IntegerCategory
     p: Ideal(R)
     T == ResidueClassRing(R, p)
     C == add { ... }

*********

Assuming this type of stuff is implemented in the library where it is
needed we can write very generic functions:

---------
   Foo(R: ModularComputation): with { ... } == add {
          bar(r: R, p:R): R == {

         elem : ResidueClassRing(R, p) :=
                  modularRep(r)$residueClassRing(p)

            -- hairy computation...

         canonicalPreImage(elem)
                } }
---------
========

Paraphrase:

Foo is a package for a ModularComputation domain R, with a function
bar: (R, R) -> R and the function bar is implemented as:

   bar(r,p) ==
        compute the elem:= modularRep(r) in S
         -- S is the ResidueClassRing for R and p.
        compute some other things
         -- (which may or may not change elem, but
         -- presumably elem remains in S
        return canonicalPreImage(elem) in R

or the function bar may be bar: (R, p:R) -> S(p) if it returns elem.


Analogy:

The function bar is just a supped-up version of Martin's example:

   g(n,k) == (k mod n):PrimeField(n)  -- assuming n is prime

so it can be implemented in Axiom. The mod function is the coerce
function in PrimeField(n)

   coerce: PositiveInteger -> PrimeField(n)

So similarly, the modularRep(r) function is a function in S

   modularRep: R -> S

is similar to a coercion from R to S and the canonicalPreImage is a
function from S to R similar to a retract:S -> R.


Axiom version:

Foo(R,p,S): T == C where
  R: ModularCategory
  p: Ideal(R)
  S: ResidueClassRing(R,p)
  Q ==> any domain constructed from R, p, S
  T == with
    bar: R -> Q
  C == add
    bar(r)==
       elem:S:=modularRep(r)$S
        -- hairy computations
       q:Q:= ...

Calling sequence in some other package, assuming R, p, S are already
defined:

   bar(r)$Foo(R,p,S)

If you want to use the default implementation when R is a domain in
IntegerCategory, you can use:

  if R has IntegerCategory then
      S:=IntegerMod(R,p)
  bar(r)$Foo(R,p,S)

Below are just some random notes (my tosses and turns):

Note here S is typically defined using one of the constructors.  If
the map residueClassRing(p) exists, then there will be a corresponding
domain constructor in Axiom.  The advantage of using residueClassRing
is that we can use ONE name for all the constructors for ALL rings R,
even if these constructions depend on R (but uniform on p).
OVERLOADING and S does not have to appear as a parameter because we
don't care how S is constructed.  It is cleaner and corresponds to the
mathematics by ignoring the data representation or implementation.  On
the other hand, by tagging S along in the package, we can use special
features in the construction of S in the computation (not really, S is
given categorically:  Example, in GroebnerBasis, we tag along S, the
Dpol, but since Dpol is just PolynomialCategory based on the other
parameters, we don't know more.  However, in actual computation,
calling groebner for example, the implementation of Dpol will come
into action.  It makes no difference to let R be a ModularComputation
because then the implementation is residueClassRing(R,p), so by
specifying R:  ModularComputation, we are singling out the
implementation, just like when we input Dpol.  So ModularComputation
is only a wrapper.  In Steve's example, each time you need functions
from S, you have to make a function call to residueClassRing(R,p),
unless, you assign S to it.  So it is only a wrapper!

********

All of this depends on the fact that we can express a dependently
typed function residueClassRing(p:R), which can be implemented by any
given domain as appropriate.  The Foo package knows all it needs to,
the Ring, and an element of the ring to get at the the quotient ring.
Of course, the bar function above could be more complex and return an
element of ResidueClassRing(R,p), etc.

How could one write the bar()$Foo above with the current axiom
language?  All you can do is write a package/domain parameterized by a
commutative ring R and a representative of R, and write the exports
dependently on R's type via `if R has ...' constructs.  This is what I
meant in the previous email about having a RESCLASS package/domain
which needed to know too much about the algebra.

=======

see discussion on top of email


*******