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

[William Sit]

See the attached file. (sorry)

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


*******