Axiom Notes - sympy/sympy GitHub Wiki

The question is if we ("we" meaning maybe a GSoC student) can implement the rest of the Risch algorithm by translating it from Axiom, where Bronstein implemented the full thing. Here are my (Aaron's) notes on playing with Axiom and trying to understand the Risch code there.

To use Axiom, use the Fricas fork. There is a GitHub mirror of the sources at https://github.com/fricas/fricas.

Compiling Fricas requires sbcl.

I personally couldn't get it to compile on macOS, but the binary at http://fricas.sourceforge.net/download.html worked.

The best resource to understand the Axiom source code is the Axiom book. I highly recommend reading sections 5 and 6 of the book to understand the Aldor language that Axiom is written it.

For whatever reason, the Axiom book uses ** for exponentiation but it does not work in Fricas. You have to use ^ instead.

The Risch algorithm is defined in several files (mainly the ones that start with int, rde, and efstruc.

Axiom's entrypoint to the Risch algorihtm is the integrate function. It works very similar to SymPy:

(4) -> f := (exp(x) + 1)/(2*x^2 - 1)

          x
        %e  + 1
   (4)  -------
          2
        2x  - 1
                                                    Type: Expression(Integer)
(5) -> integrate(f, x)

   (5)
      +-+                                             +-+ 2
     \|2                                             \|2
     ----       2      +-+               +-+         ----         +-+
       2     (2x  + 1)\|2  - 4x       - \|2  + 2x      2         \|2  + 2x
   %e    log(------------------) + Ei(-----------)(%e    )  - Ei(---------)
                     2                     2                         2
                   2x  - 1
   ------------------------------------------------------------------------
                                          +-+
                                         \|2
                                         ----
                                    +-+    2
                                  2\|2 %e
                                         Type: Union(Expression(Integer),...)

Every time you evaluate an expression interactively in Axiom, it shows you the type of the expression. Polynomials expand automatically:

(6) -> (x + 1)^2

         2
   (6)  x  + 2x + 1
                                                    Type: Polynomial(Integer)

Polynomial(Integer) means a polynomial with integer coefficients. Integer is the ring in question. The field of rational numbers is Fraction(Integer).

Commands to the Axiom system start with a ), for instance, )help.

Axiom uses a global namespace of functions. Functions are defined by their name and arity. For example, f(x), and f(x, y) can both exist at the same time. The same function can exist in many forms for different types of the inputs. Functions are either defined for a specific type (like f(x: Integer): Integer = x + 1, or defined generically (like f(x) == x + 1). If it's defined generically a different version is compiled for each type the first time it is used with that type.

The operators ::, @, and $ after an expression manipulate the type. From the book (section 0.4.5):

:: means explicitly convert X to type T if possible. $ means use the available operators for type T to compute X. @ means choose operators to compute X so that the result is of type T.

In particular, $ is used to call a function from a specific package. This is needed to access functions that aren't "exposed". You can also expose all functions from a module using )set expose add con (con is short for constructor)

(8) -> )set expose add con RDEaux
   RDEaux is now explicitly exposed in frame frame1

Note that a single file can contain multiple packages. The RDEaux package comes from src/algebra/intpar.spad, and is defined in the middle of the file (search for )abbrev package RDEAUX RDEaux).

Packages take parameters, which template the types used in the package. For instance, RDEaux takes a parameter F, which is the field used in the package (like RDEaux(Fraction(Integer))). You can use multi_SPDE(...)$RDEaux(Fraction(Integer)) for instance (or just expose RDEaux as above).

Axiom does not seem to have tests for these internal functions, only for the publicly exposed functions (integrate). I did not yet figure out how to use a debugger to see what internal functions are being called from integrate.

The functions in Axiom roughly correspond to the layout of Bronstein's book, but they don't follow the pseudocode exactly. I haven't yet found an exact duplicate function Axiom <-> SymPy that can be used to

Types convert automatically in most cases.

Aldor compiles to lisp. I'm not sure how easy it is to manipulate/debug at the lisp level.

Read sections 5 and 6 of the book. They explain the syntax really well and concisely.

Fricas requires using ^ instead of **. It looks like you can easily fix this by running

x**y == x^y

Functions of one argument can be called with spaces. f x is the same as f(x). This does not work for functions of more than one argument.

: defined type information. f(x: Integer): Fraction(Integer) is a function that takes an integer and returns a rational number.

See above for information on ::, $, and @.

+-> is an anonymous function (lambda). For instance, x +-> x + 1 is a function that adds one to the argument.

:= defines variables. == defines functions. = is a symbolic equality (like Eq in SymPy).

=> exits the current block. If the block is the function it returns the value. It can also be used as a "continue" when in a loop. For instance, from the book (5.4.4):

i := 1
repeat
  i := i + 1
  i > 3 => i
  output(i)

=> breaks out of the "block", which is just the lines

i := i + 1
i > 3 => i
output(i)

but not the repeat. So it will loop forever (but only call output once).

This is clearer if you write it in one line:

i := 1
repeat (i := i + 1; i>3 => i; output i)

The (... ; ...; ...) bit is the "block".

=> is often used to do if x: return ... constructs. For example, from intpar.spad:

multi_SPDE(a, b, lc, d, der) ==
    d <$Z 0 => [[0, c]$RSOL for c in lc]
    every?(zero?, lc) => [[0, 0]$RSOL for c in lc]
    ...

[1, 2], and (1, 2) are called "lists" and "tuples", resp., just like in Python (this is important when looking at the type information, for instance, List(UP) means a list of elements of type UP).

Undefined names become variables automatically.

++ creates a comment. This is where you should look for hints as to what the internal functions do, mathematically.

• How are the der derivation functions defined in the Risch code? I tried for instance,

  )set expose add con RDEaux
  a := t^2 + x*t*2 + x^2
  b := t^2/x^2 + (2/x - 1)*t
  c := t^2/x^2 + (2/x - 1)*t
  multi_SPDE(a, b, [c], 0, t +-> t)
  

  Trying to emulate the test from SymPy, but it crashed Fricas. (I'm not sure if this is a correct correspondence either, as c must be a list in multi_SPDE)

  The arguments a, etc. are expected to be of type UP. Those are univariate polynomials of an anonymous variable (like PurePoly in SymPy). They should probably be created by something like univariate(a, t). (t stands here for the 'kernel' exp(x).) The derivative should be of type UP -> UP. There is an example in intef.spad (an expintegrate call in expint):

           (x1 : UP) : UP+->differentiate(x1,
             (x2 : F) : F+->differentiate(x2, x),
             monomial(differentiate(eta, x), 1))
  

  Here eta is the argument of t = exp(x), so differentiate(eta, x) is equal to 1, and monomial(differentiate(eta, x), 1) is the derivative of exp(x) (i.e. exp(x)) as a polynomial. (x2 : F) : F+->differentiate(x2, x) says that the coefficients of the polynomial are differentiated with respect to x. (The outer differentiate is a three-argument function of polynomials.)

  (I have not found an explanation for the repeated appearance of UP and F. Maybe the latter is the expected result type.)

• Is it possible to implement the rest of the Risch algorithm in SymPy by copying the code from Fricas?