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
The Risch algorithm is defined in several files (mainly the ones that start
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.
is the ring in question. The field of rational numbers is
Commands to the Axiom system start with a
), for instance,
Axiom uses a global namespace of functions. Functions are defined by their
name and arity. For example,
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.
$ 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.
$ 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
(8) -> )set expose add con RDEaux RDEaux is now explicitly exposed in frame frame1
Note that a single file can contain multiple packages. The
src/algebra/intpar.spad, and is defined in the middle of the file
)abbrev package RDEAUX RDEaux).
Packages take parameters, which template the types used in the package. For
RDEaux takes a parameter
F, which is the field used in the
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
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
+-> 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
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
This is clearer if you write it in one line:
i := 1 repeat (i := i + 1; i>3 => i; output i)
(... ; ...; ...) bit is the "block".
=> is often used to do
if x: return ... constructs. For example, from
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
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
derderivation 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
cmust be a list in
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). (
tstands here for the 'kernel'
exp(x).) The derivative should be of type
UP -> UP. There is an example in
(x1 : UP) : UP+->differentiate(x1, (x2 : F) : F+->differentiate(x2, x), monomial(differentiate(eta, x), 1))
etais 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)) 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
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?