What is PyGiNaC, anyway?¶
Vladimir V Kisil¶
Introduction¶
PyGiNaC is a Python package that provides an interface to the C++ library GiNaC, which is an open framework for symbolic computation within C++. PyGiNaC is implemented with the help of the Boost.Python library. At the moment, the package is more or less in an development state, i.e.
the GiNaC classes are only partially exposed (yet most common methods are covered);
parts of regression test suite are unconverted; and
no actual documentation exists (but who is reading the full documentation nowadays if you can quickly ask at StackOverflow?)
In short: many things are already usable and further improvements are possible.
A glimpse into usage¶
Despite of being not-so-complete, PyGiNaC can do some fancy stuff for you. For example, solving a linear system of equations in the Python intepreter can be as simple as
from ginac import *
x = symbol('x')
y = symbol('y')
lsolve([3*x + 5*y == 2, 5*x+y == -3], [x,y])
The result is returned as a dictionary. To see it in a meaningful form we convert GiNaC symbol and numeric objects to strings:
soln = lsolve([3*x + 5*y == 2, 5*x+y == -3], [x,y])
[f'{str(x)} : {str(soln[x])}' for x in soln]
Power series of functions are also handled:
x = symbol('x')
str(sin(x).series(x==0, 8))
The same result is much more readable if we upgrade the output to $\LaTeX$ pretty-printing:
from IPython.display import Latex
latex_on()
Latex(f'${sin(x).series(x==0, 8)}$')
Here is a simple example of algebraic expansion:
x=realsymbol("x")
e=pow(x+1,2)
Latex(f'${e.expand()}$')
A less obvious example of an algebraic expansion and simplification with exact arithmetic is a modified version of one of Ramanujan's identities (this example is ripped off from GiNaC's regression test suite):
e1 = pow(1 + pow(3, numeric(1,5)) - pow(3, numeric(2,5)),3)
e2 = (e1 - 10 + 5*pow(3, numeric(3,5)))
display(Latex(f'e2 is ${e2}$'))
f'e2 expands to {e2.expand()}'
Above numeric(3,5)
is the fraction $\frac{3}{5}$.
Further examples can be found below in the Appendix.
Similarity and differences between GiNaC and PyGiNaC¶
Currently, (at least partially) exposed GiNaC classes are:
- Basic (the GiNaC root class)
- Constant
- Numeric
- Symbol, realsymbol, possymbol
- Expression, which are sums, products, non-commutative products, or powers of other expressions or basic objects
- Relational (see, however, the caution below)
- Matrix
- Function
- Indices
- Tensor
- Clifford unit
- Series
- Integral
- Lists are converted bidirectionally to Python3 lists
- Wildcard
- Symmetry
- Parser
- More technical: flags, exmap, expairs, expairseq, symtab
Names of classes and methods in PyGiNaC are as close as possible to their prototypes in GiNaC. In most cases you can use GiNaC Tutorial as your user manual of PyGiNaC. Some inevitable differences come from the gap between C++ and Python, for example:
- GiNaC variables (like
DIGITS
) cannot be directly assigned asDIGITS=50
. Instead you need to call special helper functions for these, e.g.set_digits(150)
as demonstrated below. - Namespace syntax in C++ like
subs_options::no_pattern
shall be replaced by Pythonishsubs_options.no_pattern
. - List initialisation in C++ like
lst l = {a, b, c}
shall be replaced byl=[a, b, c]
. - There are no semicolons at the end of statements and you need to use the proper Python identification in your scripts.
- Due to the way Python handle the operator
==
(and alike) a relational created from the operator may have the left and right parts swapped:
x = realsymbol("x")
p = possymbol("p")
print(x == p, " vs ", relational(x, p, relop.eq))
print(p == x, " vs ", relational(p, x, relop.eq))
For this reason if the order of parts is important, e.g. in
series()
orsubs()
methods, the proper constructorrelational()
needs to be used. Forsubs()
method some other alternatives are present, see the next item.
- Due to the previously discussed issues with
==
-operator and cumbersomeness of the constructorrelational()
some additional alternatives are provided. The most convenient is Python dictionaries. For example, C++ statement
e.subs(lst{x==2, y==a+3, z=2*b-1});
with the Python dictionary:
e.subs({x : 2, y : a+3, z : 2*b-1})
See an example below. Another opportunity is to use two parallel Python lists for substitutions:
e.subs([x, y, z], [2, a+3, 2*b-1])
If a GiNaC method modifies its argument, then the respective PyGiNaC wrapper returns a tuple: the first element of the tuple is the return of GiNaC method, the following element(s) are values of the modified arguments.
For example, the GiNaC method
ex ex::to_polynomial(exmap & m);
returns an expression converted to a polynomial using the substitutions from
exmap m
andm
becomes updated after the function call. Therefore the respective PyGiNaC call is
x = realsymbol("x")
e = (x+1/x)**3 +sin(x)
m = {}
[a, m] = e.to_polynomial(m)
Latex(f'Polynomial form: ${a}$ with substitutions '+', '.join(map(lambda k: f'${k}$ : ${m[k]}$', m.keys())))
(Non-)Downloads¶
You can try PyGiNaC without any local installation from two cloud services. Since Python wrapper for MoebInv libraries are built on top of PyGiNaC, the full access to PyGiNaC is provided in the following cloud Jupyter notebooks:
Effective, on those service (or any supported Ubuntu/Debian box indeed) it is enough to execute the following cell or its content as a shell script.
%%bash
DISTRO=` grep -e "^ID=" /etc/os-release | cut -d= -f2 `
# Execute this cell on Ubuntu only
if echo "ubuntu debian" | grep -qw "$DISTRO" ; then \
echo "Continue on $DISTRO " ; \
else \
echo 'Does not look like an Ubuntu/Debian box, exiting' ; \
exit ; \
fi
# Check if the software is already installed
if dpkg -l python3-pyginac > /dev/null ; then \
echo 'The package is already installed, exiting' ; \
exit ; \
fi
# Install signature key of the Ubuntu/Debian repository
## apt-key is now obsolete
curl -L https://sourceforge.net/projects/moebinv/files/binary/$DISTRO/moebinv.gpg.key | \
gpg --dearmor > /etc/apt/trusted.gpg.d/moebinv-archive-keyring.gpg
## An alternative installation with a keyring package.
#cd /tmp;\
#wget --backups=1 https://sourceforge.net/projects/moebinv/files/binary/$DISTRO/moebinv-archive-keyring_0.2_all.deb && \
#dpkg -i moebinv-archive-keyring_0.2_all.deb
# Add Ubuntu/Debian repository to known sources
CODENAME=`grep -e "VERSION_CODENAME" /etc/os-release | cut -d= -f2`
echo "deb https://sourceforge.net/projects/moebinv/files/binary/$DISTRO $CODENAME main" > \
/etc/apt/sources.list.d/moebinv.list
# Update the contents of the repository
apt update
# Install required packages and their dependencies
apt-get -q -y install python3-pyginac
If you want to have PyGiNaC installed locally you can use:
If the pre-compiled Debian/Ubuntu packages is not working for your system you can compile the binary file as described in the next section.
Compiling from source¶
Prerequisites¶
To run PyGiNaC you need to have the following software installed:
Python 3 (tested for 3.6 or higher)
The boost libraries (tested for 1.65.0 or higher)
GiNaC (tested for 1.7.2 or higher)
In addition to the above, to compile PyGiNaC the following are needed
GNU g++ (tested for 8 or higher).
Compiling PyGiNaC takes considerable memory (but not so high for modern computers), although generally not as much is needed to run it.
Compilation and test¶
Once you have all the dependencies listed above installed, issue the command
$ make
from the source directory. This will build the module in-place. The script run
can be used to start an interactive Python session that will be able to use PyGiNaC by
$ ./run python3
or, to run a Python script
$ ./run python3 some_script_file.py
For example you can run the collection of self-test by
$ ./run python3 bin/checkall.py
That's it. Have fun. :)
Installation¶
To install PyGiNaC globally run
$ make install
as the root from the source directory. Optional installation prefix, e.g. /opt
can be specified as follows:
$ make install DESTDIR=/opt
Building Debian package¶
If you are on Debian/Ubuntu system and with to re-build the binary file for some reasons you can do this with the standard command
$ debuild -us -uc
The relevant Debian packaging infrastructure need to be installed for this, of course.
Mailing list¶
There is no currently a PyGiNaC related mailing list. Feel free to write an email to the curent maintainer at kisilv@maths.leeds.ac.uk.
TODO¶
The following were identified by founding fathers as further targets:
Wrap more of GiNaC classes and objects.
Pythonize the GiNaC regression tests.
Create a map_function system, probably based on Python function objects. [See a partial solution through Python
map()
function and iterables below in App.C.]Add extra member access functions for higher-level access of containerish types, like power.basis() and power.exponent(), for example.
Make every function capable of taking more than one argument also take arguments in keyword form.
Prepare some documention with examples written in Python instead of ginsh or C++.
Feel free to contribute to this or other worthy developments.
History and contributors¶
The current implementation of PyGiNaC has, up to our knowledge, two predecessors: a version written by Pearu Peterson many years ago, and another one by Ondrej Certik.
The present version of PyGiNaC is originally written by Jonathan Brandmeyer and later co-authored by Matti Peltomäki. Patches have been submitted by Ondrej Certik. Here is the historic site.
Vladimir V. Kisil is the current maintainer of the PyGiNaC code as a subproject of MoebInv -- C++ libraries for symbolic, numeric and graphical manipulations in non-Euclidean geometry.
Appendix A: from Ginsh to PyGiNaC¶
This section is taken from the GiNaC Tutorial "2.2 What it can do for you" and shows how you can migrate your interactive usage from Ginsh/C++
to PyGiNaC
.
After invoking Python3/IPython/Jupyter shell one can test and experiment with PyGiNaC's features much like in other Computer Algebra Systems and mix it with arbitrary Python3 programming constructs like loops or conditionals. In IPython/Jupyter you can additionally benefit from extra features and magics, e.g. pretty-printed mathematics output.
Arbitrary precision arithmetic¶
(Py)GiNaC can manipulate arbitrary precision integers in a very fast way. Rational numbers are automatically converted to fractions of coprime integers:
x=pow(3,150)
str(x)
Note that a statement x=3**150
would produce a Python long integer instance. We need to use the dedicated function pow()
that the PyGiNaC will take the precedence from the interpreter.
A slightly different techniques to create a (Py)GiNaC numeric is:
y=numeric(3)**149
str(y)
The next two results are exact numbers:
str(x/y)
To pretty-print the next output we use $\LaTeX$ facilities:
Latex(f'${y/x}$')
These may be compared to the ordinary Python arithmetic:
(3**149)/(3**150)
Exact numbers are always retained as exact numbers and only evaluated as floating point numbers if requested. For instance, with numeric radicals is dealt pretty much as with symbols. Products of sums of them can be expanded:
a=symbol("a")
Latex(f'${expand((1+a**numeric(1,5)-a**numeric(2,5))**3)}$')
Latex(f'${expand((1+3**numeric(1,5)-3**numeric(2,5))**3)}$')
Latex(f'${evalf((1+3**numeric(1,5)-3**numeric(2,5))**3)}$')
The function evalf()
that was used above converts any number in (Py)GiNaC's expressions into floating point numbers. This can be done to arbitrary predefined accuracy:
Latex(f'${evalf(numeric(1,7))}$')
Now we change the required number of evaluated digits:
set_digits(150)
Latex(f'${evalf(numeric(1,7))}$')
Exact numbers other than rationals that can be manipulated in (Py)GiNaC include predefined constants like Archimedes' $\pi$, called Pi
in (Py)GiNaC. They can both be used in symbolic manipulations (as an exact number) as well as in numeric expressions (as an inexact number):
set_digits(15)
x=symbol("x")
a=Pi**2+x
Latex(f'${a}$')
Latex(f'${evalf(a)}$')
Latex(f'${evalf(a.subs({x : 2}))}$')
Latex(f'${evalf(a.subs({x : 2}))}$')
(Py)GiNaC does not provide the Euler constant $e$, because it is primary needed as a base of the exponent function, see the next subsection.
Mathematical functions¶
(Py)GiNaC is aware of main mathematical functions and can manipulate them either in the exact or an approximate manner. For example, for the above mentioned exponential function $\exp(x)=e^x$, (Py)GiNaC knows the Euler identity:
X=exp(I*Pi)+1
str(X)
Of course, the Euler constant $e$ can be created as exp(1)
:
E=exp(1)
Latex(f'${E}$')
It has the expected value:
str(evalf(E))
But (Py)GiNaC is reluctant to make a reduction based on the Euler identity for power function with complex exponent:
Latex(f'${eval(pow(E,I*Pi)+1)}$')
To see the reason think about the identity $e^0 = e^{2\pi i}$ which would imply $(e^0)^i =( e^{2\pi i})^i$.
Built-in functions can be evaluated to exact numbers if this is possible.
Latex(f'${cos(42*Pi)}$')
Latex(f'${cos(42*Pi).eval()}$')
Conversions that can be safely performed are done immediately; conversions that are not generally valid are not done:
Latex(f'${cos(acos(x))}$')
Latex(f'${cos(acos(x)).eval()}$')
However we have:
Latex(f'${acos(cos(x)).eval()}$')
Note that converting the last input to $x$ would allow one to conclude that $42 \pi$ is equal to 0.
Linear Algebra¶
Linear equation systems can be solved along with basic linear algebra manipulations over symbolic expressions. In (Py)GiNaC offers a matrix class
. We start from a single equation:
a=symbol("a")
x=symbol("x")
y=symbol("y")
z=symbol("z")
soln=lsolve([a+x*y==z], [x])
Latex(f'${soln[x]}$')
A pair of linear equations with two variables:
soln=lsolve([3*x+5*y == 7, -2*x+10*y == -5], [x, y])
for t in soln:
display(Latex(f'${t}$ : ${soln[t]}$'))
A matrix can be created from a list of lists of row elements:
M = matrix([ [1, 3], [-3, 2] ])
Latex(f'${M.determinant()}$')
The characteristic polynomial (note that lambda
is a reserved keyword in Python3):
lam=symbol("lambda")
Latex(f'${M.charpoly(lam)}$')
Matrix operations can be called in a usual way:
A = matrix([ [1, 1], [2, -1] ])
Latex(f'${A+2*M}$')
However their evaluation is postponed until an explicit request by the dedicated function evalm()
for matrix evaluation:
Latex(f'${evalm(A+2*M)}$')
Matrix arithmetic is also performed in the exact manner:
a=symbol("a")
b=symbol("b")
B = matrix([ [0, 0, a], [b, 1, -b], [-1/a, 0, 0] ])
Latex(f'${evalm(B**(2**12345))}$')
Polynomials and rational functions¶
Multivariate polynomials and rational functions may be expanded, collected and normalized (i.e. converted to a ratio of two coprime polynomials):
a = x**4 + 2*x**2*y**2 + 4*x**3*y + 12*x*y**3 - 3*y**4
Latex(f'${a}$')
b = x**2 + 4*x*y - y**2
Latex(f'${b}$')
Latex(f'${expand(a*b)}$')
Latex(f'${collect(a+b,x)}$')
Latex(f'${collect(a+b,y)}$')
Latex(f'${normal(a/b)}$')
Latex(f'${diff(tan(x),x)}$')
Any expression can be expanded as Taylor or Laurent series in a very natural syntax (the second argument of series is a relation defining the evaluation point, the third specifies the order):
x_is_0 = relational(x, 0, relop.eq)
Latex(f'${series(sin(x), x_is_0, 4)}$')
Another one:
Latex(f'${series(1/tan(x), x_is_0, 4)}$')
Or a bit more involved
Latex(f'${series(tgamma(x), x_is_0, 3)}$')
In necessary a floating point evaluation can be called as well:
Latex(f'${series(tgamma(x), x_is_0, 3).evalf()}$')
(If the last output is identical to the previous one, this shall be due to a GiNaC bug up to the version $\leq$ 1.7.7).
Latex(f'${series(tgamma(2*sin(x)-2), relational(x, Pi/2, relop.eq), 6)}$')
What about integration? Let us try:
t=realsymbol("t")
Int=integral(x, 0, t, x*x+sin(x))
display(Latex(f'${Int}$'))
Latex(f'${Int.eval_integ()}$')
That is: (Py)GiNaC is not very useful at symbolically evaluating integrals, it can do it for polynomials only. However, (Py)GiNaC is aware of the Fundamental Theorem of Calculus:
Latex(f'${Int.diff(t)}$')
Furthermore, (Py)GiNaC is not so bad at the numeric evaluation of definite integrals, see the next subsection.
Latex(f'${integral(x, 0, 1, x*x+sin(x)).evalf()}$')
See GiNaC Tutorial for fine-tuning of numerical integration.
Often, functions don't have roots in closed form. Nevertheless, it's quite easy to compute a solution numerically, to arbitrary precision (note, that for an equation we are not worried if lhs and rhs of the relational will be swapped, cf. relational caution):
set_digits(50)
str(fsolve(cos(x)==x,x,0,2))
f=exp(sin(x))-x
X=fsolve(f,x,-10,10)
str(X)
str(f.subs({x : X}))
Notice how the final result above differs slightly from zero by about $6\cdot 10^{-58}$. This is because with 50 decimal digits precision the root cannot be represented more accurately than X. Such inaccuracies are to be expected when computing with finite floating point values.
If you ever wanted to convert units in C or C++ and found this is cumbersome, here is the solution. Symbolic types can always be used as tags for different types of objects. Converting from wrong units to the metric system is now easy:
set_digits(10)
m=symbol("m")
kg=symbol("kg")
inch=.0254*m
lb=.45359237*kg
Latex(f'${(200*lb/inch**2).evalf()}$')
Parser and integration with other packages¶
(Py)GiNaC has some attractive or unique features (e.g. the universal Clifford algebras support), but it is not a one-for-all package. You may want to integrate it with other mathematical software. This can be done in a way usually used by human being: by exchange of strings of characters. (Py)GiNaC is literate enough to parse such strings possibly produced by either a user or other programme. This can be used as simple as this:
reader=parser()
e=reader("Pi^2+sin(2*t)")
Latex(f'${e}$')
Note C++-style notation for powers, the (Py)GiNaC expects these rather than pythonish expressions. The output looks nice but you may be disappointed by the next line:
str(e.diff(t))
In fact, our reader
was not aware of the previously defined symbol $t$ and created a new symbol still represented by the same letter. To avoid such a confusion the reader
need to be aware of desirable substitutions from strings to existing symbols (or even expressions):
reader=parser({"t" : t})
e=reader("Pi^2+sin(2*t)")
display(Latex(f'e=${e}$'))
Latex(r'$\frac{de}{dt}='+f'{e.diff(t)}$')
We can disable parsing of unknown symbols (for example, you want treat an unexpected string in the input as an error):
import sys
reader.set_strict(True)
try:
str(reader("1+2*z"))
except ValueError:
print("There is an undefined symbol:", sys.exc_info()[1])
and enable it back:
reader.set_strict(False)
str(reader("1+2*z"))
You can obtain all parser-generated symbols with 'get_syms()' method:
D=reader.get_syms()
[f'"{x}" : {str(D[x])}' for x in D]
The dictionary of known symbols can be updated as needed:
D.update({"w" : log(y)})
reader.set_syms(D)
Latex(f'${reader("t+s+w")}$')
Parser provides a convenient communication tool with other mathematical packages or user inputs. Fine-tuning can be achieved by additional string manipulations.
class vector3D:
x = 0
y = 0
z = 0
def __init__(self, x_, y_, z_):
self.x = x_
self.y = y_
self.z = z_
def __str__(self):
return "(" + str(self.x) + "," + str(self.y) + "," + str(self.z) + ")"
def __copy(self):
return vector3D(self.x, self.y, self.z)
def __neg__(self):
return vector3D(-self.x, -self.y, -self.z)
def __add__(self, other):
return vector3D(self.x + other.x, self.y + other.y, self.z + other.z)
def __sub__(self, other):
return self.__copy() + -other
def __mul__(self, number):
return vector3D(number * self.x, number * self.y, number * self.z)
def __rmul__(self, number):
return self.__mul__(number)
def dot(self, other):
return self.x*other.x + self.y*other.y + self.z*other.z
def cross(self, other):
return vector3D(self.y*other.z - self.z*other.y, -self.x*other.z + self.z*other.x, self.x*other.y - self.y*other.x)
Within the pure Python the class vector3D
can be used in numeric computations for vector algebra. If PyGiNaC module is imported the same class shall handle vectors with symbolic components as well. However, to make it even more convenient we can add some further methods, e.g. normal()
, expand()
, is_zero()
and may be even diff()
:
class vector3D:
x = 0
y = 0
z = 0
def __init__(self, x_, y_, z_):
self.x = x_
self.y = y_
self.z = z_
def __str__(self):
return "(" + str(self.x) + "," + str(self.y) + "," + str(self.z) + ")"
def __copy(self):
return vector3D(self.x, self.y, self.z)
def __neg__(self):
return vector3D(-self.x, -self.y, -self.z)
def __add__(self, other):
return vector3D(self.x + other.x, self.y + other.y, self.z + other.z)
def __sub__(self, other):
return self.__copy() + -other
def __mul__(self, number):
return vector3D(number * self.x, number * self.y, number * self.z)
def __rmul__(self, number):
return self.__mul__(number)
def dot(self, other):
return self.x*other.x + self.y*other.y + self.z*other.z
def cross(self, other):
return vector3D(self.y*other.z - self.z*other.y, -self.x*other.z + self.z*other.x, self.x*other.y - self.y*other.x)
#
# Some PyGiNac specific addtional merhods for symbolic computations
#
def normal(self):
return vector3D(self.x.normal(), self.y.normal(), self.z.normal())
def expand(self):
return vector3D(self.x.expand(), self.y.expand(), self.z.expand())
def is_zero(self):
return self.x.is_zero() and self.y.is_zero() and self.z.is_zero()
def diff(self, var):
return vector3D(self.x.diff(var), self.y.diff(var), self.z.diff(var))
Note, that there is no need to modify any of the already existing method to work with PyGiNaC.
Now we can check some properties of vector algebra by symbolic calculations in PyGiNaC.
v1 = vector3D( realsymbol("x1"), realsymbol("y1"), realsymbol("z1"))
v2 = vector3D( realsymbol("x2"), realsymbol("y2"), realsymbol("z2"))
v3 = vector3D( realsymbol("x3"), realsymbol("y3"), realsymbol("z3"))
"Cross product is associative: %s" % \
(v1.cross(v2).cross(v3) - v1.cross(v2.cross(v3))).normal().is_zero()
"Mixed product is given by the determinant: %s" % \
(v1.cross(v2).dot(v3) \
- matrix([[v1.x, v1.y, v1.z], \
[v2.x, v2.y, v2.z], \
[v3.x, v3.y, v3.z]]).determinant()).normal().is_zero()
"Cross product is orthogonal to the first factor: %s" % \
v1.cross(v2).dot(v1).expand().is_zero()
t = realsymbol("t")
vt = vector3D(cos(t), sin(t), numeric(0))
"Velocity of rotation is orthogonal to the radius-vector: %s" % \
vt.diff(t).dot(vt).is_zero()
Hopefully, the above examples are sufficient to give an idea how PyGiNaC can be used. Further advice can be found in GiNaC tutorial, PyGINaC test suit and MoebInv notebooks, see references below.
e = cos(x).series(x==0, 5)
Latex(f'In series ${e}:\\\\$' + '$\\\\$'.join(map(lambda t: f'Term ${t}$ has degree {t.degree(x)}', e)))
e = integral(t, 0, x, t**2+ t/2)
Latex(f'Find the value of $${e}$$ for $x={numeric(1,4)}$')
In this case you may prefer to use flatlatex package if your console supports Unicode. An example of usage is:
import flatlatex
FLC = flatlatex.converter()
def Latex(e):
e = e.replace('$$', '$') # replace displayed match swith by in-line math switch
l = e.split('$') # split into text and math chunks
is_text = True
s = ''
for t in l: # Iterate through alternating text/math chunks
if is_text:
s += t
else:
s += FLC.convert(t)
is_text = not is_text
return s
Latex(f'Find the value of $${e}$$ for $x={numeric(1,4)}$')
Note that with such redefinition of Latex()
function no alteration to other notebook's cells is required. Thus your notebook can be quickly adopted to both worlds just by enabling/disabling the above code.
import pandas
pandas.to_datetime('today')