1 SymPy: Symbolic Computing in Python

2 Supplementary material

3 As in the paper, all examples in the supplement assume that the following has been run:
4 >>> from sympy import *
5 >>> x, y, z = symbols('x y z')

6 1 LIMITS: THE GRUNTZ ALGORITHM

7 SymPy calculates limits using the Gruntz algorithm, as described in [7]. The basic idea is as
8 follows: any limit can be converted to a limit lim f (x) by substitutions like x → x1 . Then
x→∞
9 the subexpression ω (that converges to zero as x → ∞ faster than all other subexpressions) is
10 identified in f (x), and f (x) is expanded into a series with respect to ω. Any positive powers
11 of ω converge to zero (while negative powers indicate an infinite limit) and any constant term
12 independent of ω determines the limit. When a constant term still dependends on x the Gruntz
13 algorithm is applied again until a final numerical value is obtained as the limit.
To determine the most rapidly varying subexpression, the comparability classes must first be
defined, by calculating L:
log |f (x)|
L ≡ lim (1)
x→∞ log |g(x)|

The relations <, >, and ∼ are defined as follows: f > g when L = ±∞ (it is said that f is more
rapidly varying than g, i.e., f goes to ∞ or 0 faster than g), f < g when L = 0 (f is less rapidly
varying than g) and f ∼ g when L 6= 0, ±∞ (both f and g are bounded from above and below by
suitable integral powers of the other). Note that if f > g, then f > g n for any n. Here are some
examples of comparability classes:
2 x
2 < x < ex < ex < ee

2 ∼ 3 ∼ −5
1
x ∼ x2 ∼ x3 ∼ ∼ xm ∼ −x
x
−x
ex ∼ e−x ∼ e2x ∼ ex+e
1
f (x) ∼
f (x)
The Gruntz algorithm is now illustrated with the following example:
−x 1
f (x) = ex+2e − ex + . (2)
x
14 First, the set of most rapidly varying subexpressions is determined — the so-called mrv set.
−x
15 For (2), the mrv set {ex , e−x , ex+2e } is obtained. These are all subexpressions of (2) and they
16 all belong to the same comparability class. This calculation can be done using SymPy as follows:

17 >>> from sympy.series.gruntz import mrv

18 >>> mrv(exp(x+2*exp(-x))-exp(x) + 1/x, x)[0].keys()
19 dict_keys([exp(x + 2*exp(-x)), exp(x), exp(-x)])

20 Next, an arbitrary item ω is taken from mrv set that converges to zero for x → ∞ and doesn’t
21 have subexpressions in the given mrv set. If such a term is not present in the mrv set (i.e.,
22 all terms converge to infinity instead of zero), the relation f (x) ∼ f (x)
1
can be used. In the
23 considered case, only item ω = e −x can be accepted.
 The next step is to rewrite the mrv set in terms of ω = g(x). Every element f (x) of the mrv
set is rewritten as Aω c , where

log f (x)
c = lim , A = elog f −c log g (3)
x→∞ log g(x)

Note that this step includes calculation of more simple limits, for instance
−x
log ex+2e x + 2e−x
lim = lim = −1 (4)
x→∞ log e−x x→∞ −x

24 In this example we obtain the rewritten mrv set: { ω1 , ω, ω1 e2ω }. This can be done in SymPy with
25 >>> from sympy.series.gruntz import mrv, rewrite
26 >>> m = mrv(exp(x+2*exp(-x))-exp(x) + 1/x, x)
27 >>> w = Symbol('w')
28 >>> rewrite(m[1], m[0], x, w)[0]
29 1/x + exp(2*w)/w - 1/w

Then the rewritten subexpressions are substituted back into f (x) in (2) and the result is expanded
with respect to ω:
1 1 1 2ω 1
f (x) = − + e = 2 + + 2ω + O(ω 2 ) (5)
x ω ω x
Since ω is from the mrv set, then in the limit as x → ∞, ω → 0, and so 2ω + O(ω 2 ) → 0 in (5):

1 1 1 2ω 1 1
f (x) = − + e = 2 + + 2ω + O(ω 2 ) → 2 + (6)
x ω ω x x
30 In this example the result (2 + x1 ) still depends on x, so the above procedure is repeated until
31 just a value independent of x is obtained. This is the final limit. In the above case the limit is 2,
32 as can be verified by SymPy:

33 >>> limit(exp(x+2*exp(-x))-exp(x) + 1/x, x, oo)

34 2

In general, when f (x) is expanded in terms of ω, the following is obtained:

1 C−2 (x) C−1 (x)

 
f (x) = O 3
+ 2
+ +C0 (x) + C1 (x)ω + O(ω 2 ) (7)
ω ω ω | {z } | {z }
| {z } | {z } | {z } 0 0
∞ ∞ ∞

35 The positive powers of ω are zero. If there are any negative powers of ω, then the result of the
36 limit is infinity, otherwise the limit is equal to lim C0 (x). The expression C0 (x) is always simpler
x→∞
37 than original f (x), same is true for limits, arising in the rewrite stage (3), so the algorithm
38 converges. A proof of this and further details on the algorithm are given in Gruntz’s PhD
39 thesis [7].

40 2 SERIES
41 2.1 Series Expansion
42 SymPy is able to calculate the symbolic series expansion of an arbitrary series or expression
43 involving elementary and special functions and multiple variables. For this it has two different
44 implementations: the series method and Ring Series.
45 The first approach stores a series as an instance of the Expr class. Each function has its
46 specific implementation of its expansion, which is able to evaluate the Puiseux series expansion
47 about a specified point. For example, consider a Taylor expansion about 0:

2/16
48 >>> series(sin(x+y) + cos(x*y), x, 0, 2)
49 1 + sin(y) + x*cos(y) + O(x**2)

50 The newer and much faster approach called Ring Series makes use of the fact that a truncated
51 Taylor series is simply a polynomial. Correspondingly, they may be represented by sparse
52 polynomials which perform well in a under a wide range of cases. Ring Series also gives the user
53 the freedom to choose the type of coefficients to use, resulting in faster operations on certain
54 types.
55 For this, several low-level methods for expansion of trigonometric, hyperbolic and other
56 elementary operations (like series inversion, calculating the nth root, etc.) are implemented
57 using variants of the Newton Method [Brent and Zimmermann]. All these support Puiseux series
58 expansion. The following example demonstrates the use of an elementary function that calculates
59 the Taylor expansion of the sine of a series.

60 >>> from sympy.polys.ring_series import rs_sin

61 >>> R, t = ring('t', QQ)
62 >>> rs_sin(t**2 + t, t, 5)
63 -1/2*t**4 - 1/6*t**3 + t**2 + t

64 The function sympy.polys.rs_series makes use of these elementary functions to expand

65 an arbitrary SymPy expression. It does so by following a recursive strategy of expanding the
66 lowermost functions first and then composing them recursively to calculate the desired expansion.
67 Currently, it only supports expansion about 0 and is under active development. Ring Series is
68 several times faster than the default implementation with the speed difference increasing with
69 the size of the series. The sympy.polys.rs_series takes as input any SymPy expression and
70 hence there is no need to explicitly create a polynomial ring. An example demonstrating its use:

71 >>> from sympy.polys.ring_series import rs_series

72 >>> from sympy.abc import a, b
73 >>> rs_series(sin(a + b), a, 4)
74 -1/2*(sin(b))*a**2 + (sin(b)) - 1/6*a**3*(cos(b)) + a*(cos(b))

75 2.2 Formal Power Series

76 SymPy can be used for computing the formal power series of a function. The implementation
77 is based on the algorithm described in the paper on formal power series [8]. The advantage of
78 this approach is that an explicit formula for the coefficients of the series expansion is generated
79 rather than just computing a few terms.
80 The following example shows how to use fps:

81 >>> f = fps(sin(x), x, x0=0)

82 >>> f.truncate(6)
83 x - x**3/6 + x**5/120 + O(x**6)
84 >>> f[15]
85 -x**15/1307674368000

86 2.3 Fourier Series

87 SymPy provides functionality to compute Fourier series of a function using the fourier_series
88 function:

89 >>> L = symbols('L')
90 >>> expr = 2 * (Heaviside(x/L) - Heaviside(x/L - 1)) - 1
91 >>> f = fourier_series(expr, (x, 0, 2*L))
92 >>> f.truncate(3)
93 4*sin(pi*x/L)/pi + 4*sin(3*pi*x/L)/(3*pi) + 4*sin(5*pi*x/L)/(5*pi)

3/16
94 3 LOGIC
95 SymPy supports construction and manipulation of boolean expressions through the sympy.logic
96 module. SymPy symbols can be used as propositional variables and subsequently be replaced
97 with True or False values. Many functions for manipulating boolean expressions have been
98 implemented in the logic module.

99 3.1 Constructing boolean expressions

100 A boolean variable can be declared as a SymPy Symbol. Python operators &, | and ~ are
101 overridden when using SymPy objects to use the SymPy functionality for logical And, Or, and
102 Not. Other logic functions are also integrated into SymPy, including Xor and Implies, which are
103 constructed with ˆ and >>, respectively. Expressions can therefore be constructed either by using
104 the shortcut operator notation or by directly creating the relevant objects: And(), Or(), Not(),
105 Xor(), Implies(), Nand(), Nor(), etc.:

106 >>> e = (x & y) | z

107 >>> e.subs({x: True, y: True, z: False})
108 True

109 3.2 CNF and DNF

110 Any boolean expression can be converted to conjunctive normal form, disjunctive normal form,
111 or negation normal form. The API also exposes methods to check if a boolean expression is in
112 any of the aforementioned forms.
113 >>> from sympy.logic.boolalg import is_dnf, is_cnf
114 >>> to_cnf((x & y) | z)
115 And(Or(x, z), Or(y, z))
116 >>> to_dnf(x & (y | z))
117 Or(And(x, y), And(x, z))
118 >>> is_cnf((x | y) & z)
119 True
120 >>> is_dnf((x & y) | z)
121 True

122 3.3 Simplification and Equivalence

123 The sympy.logic module supports simplification of given boolean expression by making deductions
124 from the expression. Equivalence of two logical expressions can also be checked. In the case of
125 equivalence, the function bool_map can be used to show which variables of the first expression
126 correspond to which variables of the second one.
127 >>> a, b, c = symbols('a b c')
128 >>> e = a & (~a | ~b) & (a | c)
129 >>> simplify(e)
130 And(Not(b), a)
131 >>> e1 = a & (b | c)
132 >>> e2 = (x & y) | (x & z)
133 >>> bool_map(e1, e2)
134 (And(Or(b, c), a), {a: x, b: y, c: z})

135 3.4 SAT solving

136 The module also supports satisfiability (SAT) checking of a given boolean expression. If an
137 expression is satisfiable, it is possible to return a variable assignment which satisfies it. The
138 API also supports listing all possible assignments. The SAT solver has a clause learning DPLL
139 algorithm implemented with a watch literal scheme and VSIDS heuristic [11].
140 >>> satisfiable(a & (~a | b) & (~b | c) & ~c)
141 False
142 >>> satisfiable(a & (~a | b) & (~b | c) & c)
143 {a: True, b: True, c: True}

4/16
144 4 DIOPHANTINE EQUATIONS
145 Diophantine equations play a central role in number theory. A Diophantine equation has the
146 form, f (x1 , x2 , . . . , xn ) = 0 where n ≥ 2 and x1 , x2 , . . . , xn are integer variables. If there are n
147 integers a1 , a2 , . . . , an such that x1 = a1 , x2 = a2 , . . . , xn = an satisfies the above equation, the
148 equation is said to be solvable.
149 Currently, the following five types of Diophantine equations can be solved using SymPy’s
150 Diophantine module (a1 , . . . , an+1 , a, b, c, d, e, f , and k are explicitly given rational constants):
151 • Linear Diophantine equations: a1 x1 + a2 x2 + · · · + an xn = b
152 • General binary quadratic equation: ax2 + bxy + cy 2 + dx + ey + f = 0
153 • Homogeneous ternary quadratic equation: ax2 + by 2 + cz 2 + dxy + eyz + f zx = 0
154 • Extended Pythagorean equation: a1 x21 + a2 x22 + · · · + an x2n = an+1 x2n+1
155 • General sum of squares: x21 + x22 + · · · + x2n = k
156 The diophantine function factors the equation it is given (if possible), solves each factor
157 separately, and combines the results to give a final solution set. The following examples illustrate
158 some of the basic functionalities of the Diophantine module.
159 >>> from sympy.solvers.diophantine import *
160 >>> diophantine(2*x + 3*y - 5)
161 set([(3*t_0 - 5, -2*t_0 + 5)])
162

163 >>> diophantine(2*x + 4*y - 3)

164 set()
165

166 >>> diophantine(x**2 - 4*x*y + 8*y**2 - 3*x + 7*y - 5)

167 set([(2, 1), (5, 1)])
168

169 >>> diophantine(x**2 - 4*x*y + 4*y**2 - 3*x + 7*y - 5)

170 set([(-2*t**2 - 7*t + 10, -t**2 - 3*t + 5)])
171

172 >>> diophantine(3*x**2 + 4*y**2 - 5*z**2 + 4*x*y - 7*y*z + 7*z*x)

173 set([(-16*p**2 + 28*p*q + 20*q**2,
174 3*p**2 + 38*p*q - 25*q**2,
175 4*p**2 - 24*p*q + 68*q**2)])
176

177 >>> x1, x2, x3, x4, x5, x6 = symbols('x1, x2, x3, x4, x5, x6')
178 >>> diophantine(9*x1**2 + 16*x2**2 + x3**2 + 49*x4**2 + 4*x5**2 - 25*x6**2)
179 set([(70*t1**2 + 70*t2**2 + 70*t3**2 + 70*t4**2 - 70*t5**2, 105*t1*t5,
180 420*t2*t5, 60*t3*t5, 210*t4*t5,
181 42*t1**2 + 42*t2**2 + 42*t3**2 + 42*t4**2 + 42*t5**2)])
182

183 >>> a, b, c, d = symbols('a:d')

184 >>> diophantine(a**2 + b**2 + c**2 + d**2 - 23)
185 set([(2, 3, 3, 1)])

186 5 SETS
187 SymPy supports representation of a wide variety of mathematical sets. This is achieved by first
188 defining abstract representations of atomic set classes and then combining and transforming
189 them using various set operations.
190 Each of the set classes inherits from the base class Set and defines methods to check
191 membership and calculate unions, intersections, and set differences. When these methods are
192 not able to evaluate to atomic set classes, they are represented as abstract unevaluated objects.
193 SymPy has the following atomic set classes:

5/16
194 • EmptySet represents the empty set ∅.
195 • UniversalSet is an abstract “universal set” of which everything is a member. The union of
196 the universal set with any set gives the universal set and the intersection gives the other
197 set itself.
198 • FiniteSet is functionally equivalent to Python’s built in set object. Its members can be
199 any SymPy object including other sets.
200 • Integers represents the set of integers Z.
201 • Naturals represents the set of natural numbers N, i.e., the set of positive integers.
202 • Naturals0 represents the set of whole numbers N0 , which are all the non-negative integers.
203 • Range represents a range of integers. A range is defined by specifying a start value, an end
204 value, and a step size. The enumeration of a Range object is functionally equivalent to
205 Python’s range except it supports infinite endpoints, allowing the representation of infinite
206 ranges.
207 • Interval represents an interval of real numbers. It is defined by giving the start and the
208 end points and by specifying if the interval is open or closed on the respective ends.

209 Other than unevaluated classes of Union, Intersection, and Complement operations, SymPy
210 has the following set classes.

211 • ProductSet defines the Cartesian product of two or more sets. The product set is useful
212 when representing higher dimensional spaces. For example, to represent a three-dimensional
213 space, SymPy uses the Cartesian product of three real sets.
214 • ImageSet represents the image of a function when applied to a particular set. The image
215 set of a function F with respect to a set S is {F (x) | x ∈ S}. SymPy uses image sets to
216 represent sets of infinite solutions of equations such as sin(x) = 0.
217 • ConditionSet represents a subset of a set whose members satisfy a particular condition.
218 The subset of set S given by the condition H is {x | H(x), x ∈ S}. SymPy uses condition
219 sets to represent the set of solutions of equations and inequalities, where the equation or
220 the inequality is the condition and the set is the domain over which it is being solved.

221 A few other classes are implemented as special cases of the classes described above. The set of
222 real numbers, Reals, is implemented as a special case of Interval. ComplexRegion is implemented
223 as a special case of ImageSet. ComplexRegion supports both polar and rectangular representation
224 of regions on the complex plane.

225 6 STATISTICS
226 The sympy.stats module provides random variable types and methods for computing of statistical
227 properties of expressions involving random variables, which can be either continuous or discrete,
228 the latter ones being further divided into finite and infinite. The variables are associated with
229 probability densities on corresponding domains and internally defined in terms of probability
230 spaces. Apart from the possibility of defining the random variables from user supplied density
231 distribution, SymPy provides definitions of most common distributions, including Uniform,
232 Poisson, Normal, Binomial, Bernoulli, and many others.
233 Properties of random expressions can be calculated using, e.g., expectation (abbreviated E)
234 and variance to calculate expectation and variance. Internally, these functions generate integrals
235 and summations, which are automatically evaluated. The evaluation can be suppressed using
236 evaluate=False keyword argument.
237 Conditions on random variables can be defined with inequalities, equalities, and logical
238 operators and their overall probabilities are obtained using P. The features can be illustrated on
239 a model of two dice throws:

6/16
240 >>> from sympy.stats import Die, P, E
241 >>> X, Y = Die("X"), Die("Y")
242 >>> P(Eq(X, 6) & Eq(Y, 6))
243 1/36
244 >>> P(X>Y)
245 5/12

246 The conditions can also be supplied as a second parameter to E, P, and other methods to calculate
247 the property given the condition:

248 >>> E(X, X+Y<5)

249 5/3

250 Using the facilities of the sympy.stats module, one can, for example, calculate the well known
251 properties of maxwellian velocity distribution

252 >>> from sympy.stats import Maxwell, density

253 >>> kT, m, x = symbols("kT m x", positive=True)
254 >>> v = Maxwell("v", sqrt(kT/m))
255 >>> E(v) # mean velocity
256 2*sqrt(2)*sqrt(kT)/(sqrt(pi)*sqrt(m))
257 >>> E(v, evaluate=False) # unevaluated mean velocity
258 Integral(sqrt(2)*m**(3/2)*v**3*exp(-m*v**2/(2*kT))/(sqrt(pi)*kT**(3/2)),
259 (v, 0, oo))
260 >>> E(m*v**2/2) # mean energy
261 3*kT/2
262 >>> solve(density(v)(x).diff(x), x)[0] # most probable velocity
263 sqrt(2)*sqrt(kT)/sqrt(m)

264 More information on the sympy.stats module can be found in [13].

265 7 CATEGORY THEORY

266 SymPy includes a module for dealing with categories—abstract mathematical objects representing
267 classes of structures as classes of objects (points) and morphisms (arrows) between the objects.
268 The module was designed with the following two goals in mind:

269 1. automatic typesetting of diagrams given by a collection of objects and of morphisms

270 between them, and

271 2. specification and semi-automatic derivation of properties using commutative diagrams.

272 As of version 1.0, SymPy only implements the first goal, while a partially working draft
273 of implementation of the second goal is available at https://github.com/scolobb/sympy/tree/
274 ct4-commutativity.
275 In order to achieve the two goals, the module sympy.categories defines several classes
276 representing some of the essential concepts: objects, morphisms, categories, and diagrams. In
277 category theory, the inner structure of objects is often discarded in the favor of studying the
278 properties of morphisms, so the class Object is essentially a synonym of the class Symbol. There
279 are several morphism classes which do not have a particular internal structure either, though an
280 exception is CompositeMorphism, which essentially stores a list of morphisms.
281 The class Diagram captures the properties of morphisms. This class stores a family of
282 morphisms, the corresponding source and target objects, and, possibly, some properties of the
283 morphisms. Generally, no restrictions are imposed on what the properties may be—for example,
284 one might use strings of the form “forall”, “exists”, “unique”, etc. Furthermore, the morphisms
285 of a diagram are grouped into premises and conclusions in order to be able to represent logical
286 implications of the form “for a collection of morphisms P with properties p : P → Ω (the premises),
287 there exists a collection of morphisms C with properties c : C → Ω (the conclusions)”, where Ω is

7/16
288 the universal collection of properties. Finally, the class Category includes a collection of diagrams
289 which are deemed commutative and which therefore define the properties of this category.
290 Automatic typesetting of diagrams takes a Diagram and produces LATEX code using the Xy-pic
291 package. Typesetting is done in two stages: layout and generation of Xy-pic code. The layout
292 stage is taken care of by the class DiagramGrid, which takes a Diagram and lays out the objects
293 in a grid, trying to reduce the average length of the arrows in the final picture. By default,
294 DiagramGrid uses a series of triangle-based heuristics to produce a rectangular grid. A linear
295 layout can also be imposed. Furthermore, groups of objects can be given; in this case, the groups
296 will be treated as atomic cells, and the member objects will be typeset independently of the
297 other objects.
298 The second phase of diagram typesetting consists in actually drawing the picture and is
299 carried out by the class XypicDiagramDrawer. An example of a diagram automatically typeset by
DiagramgGrid and XypicDiagramDrawer in given in Figure 1.

h1
h
h2
lA } )  lDl
,A /Bo D
n M
A f k nD
g

CM l
n
lC C

Figure 1. An automatically typeset commutative diagram

300

301 As far as the second main goal of the module is concerned, the principal idea consists in
302 automatically deciding whether a diagram is commutative or not, given a collection of “axioms”:
303 diagrams known to be commutative. The implementation is based on graph embeddings (injective
304 maps): whenever an embedding of a commutative diagram into a given diagram is found, one
305 concludes that the subdiagram is commutative. Deciding commutativity of the whole diagram is
306 therefore based (theoretically) on finding a “cover” of the target diagram by embeddings of the
307 axioms. The naïve implementation proved to be prohibitively slow; a better optimized version is
308 therefore in order, as well as application of heuristics.

309 8 SYMPY GAMMA

310 SymPy Gamma is a simple web application that runs on Google App Engine. It executes and
311 displays the results of SymPy expressions as well as additional related computations, in a fashion
312 similar to that of Wolfram|Alpha. For instance, entering an integer will display its prime factors,
313 digits in the base-10 expansion, and a factorization diagram. Entering a function will display its
314 docstring; in general, entering an arbitrary expression will display its derivative, integral, series
315 expansion, plot, and roots.
316 SymPy Gamma also has several features beyond just computing the results using SymPy.

317 • SymPy Gamma displays integration and differentiation steps in detail, which can be viewed
318 in Figure 2:

8/16
319

320 Figure 2. Integral steps of tan(x)

321 • SymPy Gamma displays the factor tree diagrams for different numbers.

322 • SymPy Gamma saves user search queries, and offers many such similar features for free,
323 which Wolfram|Alpha only offers to its paid users.

324 Every input query from the user on SymPy Gamma is first parsed by its own parser capable of
325 handling several different forms of function names which SymPy as a library does not support.
326 For instance, SymPy Gamma supports queries like sin x, whereas SymPy will only recognise
327 sin(x).
328 This parser converts the input query to the equivalent SymPy readable code, which is then
329 processed by SymPy, and the result is finally printed with the built-in LATEX output and rendered
330 by the SymPy Gamma web application.

331 9 SYMPY LIVE

332 SymPy Live is an online Python shell, which uses the Google App Engine to executes SymPy
333 code. It is integrated in the SymPy documentation examples at http://docs.sympy.org.
334 This is accomplished by providing a HTML/JavaScript GUI for entering source code and
335 visualization of output, and a server that evaluates the requested source code. It is an interactive
336 AJAX shell that runs SymPy code using Python on the server.

9/16
337 10 COMPARISON WITH MATHEMATICA
338 Wolfram Mathematica is a popular proprietary CAS that features highly advanced algorithms,
339 has a core written in C++ [16], and interprets its own programming language, Wolfram Language.
340 Analogous to Lisp S-expressions, Mathematica uses its own style of M-expressions, which
341 are arrays of either atoms or other M-expressions. The first element of the expression identifies
342 the type of the expression and is indexed by zero, and the first argument is indexed starting
343 with one. In SymPy, expression arguments are stored in a Python tuple (that is, an immutable
344 array), while the expression type is identified by the type of the object storing the expression.
345 Mathematica can associate attributes to its atoms. Attributes may define mathematical
346 properties and behavior of the nodes associated to the atom. In SymPy, the usage of static class
347 fields is roughly similar to Mathematica’s attributes, though other programming patterns may
348 also be used to achieve an equivalent behavior such as class inheritance.
349 Unlike SymPy, Mathematica’s expressions are mutable: one can change parts of the expression
350 tree without the need of creating a new object. The mutability of Mathematica expressions
351 allows for a lazy updating of any references to a given data structure.
352 Products in Mathematica are determined by some built in node types, such as Times, Dot,
353 and others. Times is a representation of the * operator, and is always meant to represent a
354 commutative product operator. The other notable product is Dot, which represents the . operator.
355 This product represents matrix multiplication. It is not commutative. Unlike Mathematica,
356 SymPy determines commutativity with respect to multiplication from the expression type of the
357 factors. Mathematica puts the Orderless attribute on the expression type.
358 Regarding associative expressions, SymPy handles associativity of sums and products by
359 automatically flattening them, Mathematica specifies the Flat attribute on the expression type.
360 Mathematica relies heavily on pattern matching—even the so-called equivalent of function
361 declaration is in reality the definition of a pattern generating an expression tree transformation
362 on input expressions. Mathematica’s pattern matching is sensitive to associative, commutative,
363 and one-identity properties of its expression tree nodes. SymPy has various ways to perform
364 pattern matching. All of them play a lesser role in the CAS than in Mathematica and are
365 basically available as a tool to rewrite expressions. The differential equation solver in SymPy
366 somewhat relies on pattern matching to identify differential equation types, but it is envisaged to
367 replace that strategy with analysis of Lie symmetries in the future. Mathematica’s real advantage
368 is the ability to add (at runtime) new overloading to the expression builder or specific subnodes.
369 Consider for example:

370 In[1]:= Unprotect[Plus]

371 Out[1]= {Plus}
372

373 In[2]:= Sin[x_]^2 + Cos[y_]^2 := 1

374

375 In[3]:= x + Sin[t]^2 + y + Cos[t]^2

376 Out[3]= 1 + x + y

377 This expression in Mathematica defines a substitution rule that overloads the functionality of
378 the Plus node (the node for additions in Mathematica). A symbol with a trailing underscore is
379 treated as a wildcard. Although one may wish to keep this identity unevaluated, this example
380 clearly illustrates the potential to define one’s own immediate transformation rules. In SymPy,
381 the operations constructing the addition node in the expression tree are Python class constructors
382 and cannot be modified at runtime.1 The way SymPy deals with extending the missing runtime
383 overloadability functionality is by subclassing the node types: subclasses may redefine the class
384 constructor to yield the proper extended functionality.
385 Unlike SymPy, Mathematica does not support type inheritance or polymorphism [4]. SymPy
386 relies heavily on class inheritance, but for the most part, class inheritance is used to make sure
387 that SymPy objects inherit the proper methods and implement the basic hashing system.
1 Nonetheless, Python supports monkey patching but it is a discouraged programming pattern.

10/16
388 While Mathematica interprets nested lists as matrices whenever the sublists have the same
389 length, matrices in SymPy are a type in their own right, allowing ordinary operators and functions
390 (like multiplication and exponentiation) to be used as they traditionally are in mathematics.

391 >>> exp(Matrix([[1, 1],[0, 2]])) * Matrix([a, b])

392 Matrix([
393 [E*a + b*(-E + exp(2))],
394 [ b*exp(2)]])

395 Using the standard multiplication in Mathematica performs an element-wise product and
396 calling the exponential function Exp on a matrix returns an element-wise exponentiation of its
397 elements.
398 Unevaluated expressions in Mathematica can be achieved in various ways, most commonly
399 with the HoldForm or Hold nodes, that block the evaluation of subnodes by the parser. Such a
400 node cannot be expressed in Python because of greedy evaluation. Whenever needed in SymPy,
401 it is necessary to add the parameter evaluate=False to all subnodes.
402 In Mathematica, the operator == returns a boolean whenever it is able to immediately evaluate
403 the truth of the equality, otherwise it returns an Equal expression. In SymPy, == means structural
404 equality and is always guaranteed to return a boolean expression. To express a mathematical
405 equality in SymPy it is necessary to explicitly construct an instance of the Equality class.
406 SymPy, in accordance with Python (and unlike the usual programming convention), uses **
407 to express the power operator, while Mathematica uses the more common ^.
408 SymPy’s use of floating-point numbers is similar to that of most other CASs, including
409 Maple and Maxima. By contrast, Mathematica uses a form of significance arithmetic [14]
410 for approximate numbers. This offers further protection against numerical errors, although it
411 comes with its own set of problems (for a critique of significance arithmetic, see Fateman [4]).
412 Internally, SymPy’s evalf method works similarly to Mathematica’s significance arithmetic, but
413 the semantics are isolated from the rest of the system.

414 11 OTHER PROJECTS THAT USE SYMPY

415 There are several projects that use SymPy as a library for implementing a part of their function-
416 ality. Some of them are listed below:

417 • Cadabra: Cadabra is a CAS designed specifically for the resolution of problems encountered
418 in field theory.

419 • Octave Symbolic: The Octave-Forge Symbolic package adds symbolic calculation features
420 to GNU Octave. These include common CAS tools such as algebraic operations, calculus,
421 equation solving, Fourier and Laplace transforms, variable precision arithmetic, and other
422 features.

423 • SymPy.jl: Provides a Julia interface to SymPy using PyCall.

424 • Mathics: Mathics is a free, general-purpose online CAS featuring Mathematica compatible
425 syntax and functions. It is backed by highly extensible Python code, relying on SymPy for
426 most mathematical tasks.

427 • Mathpix: An iOS App, that detects handwritten math as input, and uses SymPy Gamma
428 to evaluate the math input and generate the relevant steps to solve the problem.

429 • IKFast: IKFast is a robot kinematics compiler provided by OpenRAVE. It analytically

430 solves robot inverse kinematics equations and generates optimized C++ files. It uses
431 SymPy for its internal symbolic mathematics.

432 • Sage: A CAS, visioned to be a viable free open source alternative to Magma, Maple,
433 Mathematica and MATLAB. Sage includes many open source mathematical libraries,
434 including SymPy.

11/16
435 • SageMathCloud: SageMathCloud is a web-based cloud computing and course manage-
436 ment platform for computational mathematics.

437 • PyDy: Multibody Dynamics with Python.

438 • galgebra: Geometric algebra (previously sympy.galgebra).

439 • yt: Python package for analyzing and visualizing volumetric data (yt.units uses SymPy).

440 • SfePy: Simple finite elements in Python, see section 11.1.

441 • Quameon: Quantum Monte Carlo in Python.

442 • Lcapy: Experimental Python package for teaching linear circuit analysis.

443 • Quantum Programming in Python: Quantum 1D Simple Harmonic Oscillator and

444 Quantum Mapping Gate.

445 • LaTeX Expression project: Easy LATEX typesetting of algebraic expressions in symbolic
446 form with automatic substitution and result computation.

447 • Symbolic statistical modeling: Adding statistical operations to complex physical

448 models.

449 11.1 SfePy

450 SfePy (Simple finite elements in Python), cf. [3], is a Python package for solving partial
451 differential equations (PDEs) in 1D, 2D and 3D by the finite element (FE) method [17]. SymPy
452 is used within this package mostly for code generation and testing, namely:

453 • generation of the hierarchical FE basis module, involving generation and symbolic differenti-
454 ation of 1D Legendre and Lobatto polynomials, constructing the FE basis polynomials [15]
455 and generating the C code;

456 • generation of symbolic conversion formulas for various groups of elastic constants [6]:
457 provide any two of the Young’s modulus, Poisson’s ratio, bulk modulus, Lamé’s first
458 parameter, shear modulus (Lamé’s second parameter) or longitudinal wave modulus and
459 get the other ones;

460 • simple physical unit conversions, generation of consistent unit sets;

461 • testing FE solutions using method of manufactured (analytical) solutions: the differential
462 operator of a PDE is symbolically applied and a symbolic right-hand side is created,
463 evaluated in quadrature points, and subsequently used to obtain a numerical solution that
464 is then compared to the analytical one;

465 • testing accuracy of 1D, 2D and 3D numerical quadrature formulas (cf. [1]) by generating
466 polynomials of suitable orders, integrating them, and comparing the results with those
467 obtained by the numerical quadrature.

468 12 TENSORS
469 Ongoing work to provide the capabilities of tensor computer algebra has so far produced the
470 tensor module. It comprises three submodules whose purposes are quite different: sympy.
471 tensor.indexed and sympy.tensor.indexed_methods support indexed symbols, sympy.tensor.
472 array contains facilities to operate on symbolic N -dimensional arrays, and finally sympy.tensor.
473 tensor is used to define abstract tensors. The abstract tensors submodule is inspired by xAct [10]
474 and Cadabra [12]. Canonicalization based on the Butler-Portugal [9] algorithm is supported in
475 SymPy. Tensor support in SymPy is currently limited to polynomial tensor expressions.

12/16
476 13 NUMERICAL SIMPLIFICATION
477 The nsimplify function in SymPy (a wrapper of identify in mpmath) attempts to find a simple
478 symbolic expression that evaluates to the same numerical value as the given input. It works
479 by applying a few simple transformations (including square roots, reciprocals, logarithms and
480 exponentials) to the input and, for each transformed value, using the PSLQ algorithm [5] to
481 search for a matching algebraic number or optionally a linear combination of user-provided base
482 constants (such as π).

483 >>> t = 1 / (sin(pi/5)+sin(2*pi/5)+sin(3*pi/5)+sin(4*pi/5))**2

484 >>> nsimplify(t)
485 -2*sqrt(5)/5 + 1
486 >>> nsimplify(pi, tolerance=0.01)
487 22/7
488 >>> nsimplify(1.783919626661888, [pi], tolerance=1e-12)
489 pi/(-1/3 + 2*pi/3)

490 14 EXAMPLES
491 14.1 Simplification
492 • expand:

493 >>> expand((x + y)**3)

494 x**3 + 3*x**2*y + 3*x*y**2 + y**3

495 • factor:

496 >>> factor(x**3 + 3*x**2*y + 3*x*y**2 + y**3)

497 (x + y)**3

498 • collect:

499 >>> collect(y*x**2 + 3*x**2 - x*y + x - 1, x)

500 x**2*(y + 3) + x*(-y + 1) - 1

501 • cancel:

502 >>> cancel((x**2 + 2*x + 1)/(x**2 - 1))

503 (x + 1)/(x - 1)

504 • apart:

505 >>> apart((x**3 + 4*x - 1)/(x**2 - 1))

506 x + 3/(x + 1) + 2/(x - 1)

507 • trigsimp:

508 >>> trigsimp(cos(x)**2*tan(x) - sin(2*x))

509 -sin(2*x)/2

510 14.2 Polynomials

511 • Factorization:

512 >>> t = symbols('t')

513 >>> f = (2115*x**4*y + 45*x**3*z**3*t**2 - 45*x**3*t**2 -
514 ... 423*x*y**4 - 47*x*y**3 + 141*x*y*z**3 + 94*x*y*z*t -
515 ... 9*y**3*z**3*t**2 + 9*y**3*t**2 - y**2*z**3*t**2 +

13/16
516 ... y**2*t**2 + 3*z**6*t**2 + 2*z**4*t**3 - 3*z**3*t**2 -
517 ... 2*z*t**3)
518 >>> factor(f)
519 (t**2*z**3 - t**2 + 47*x*y)*(2*t*z + 45*x**3 - 9*y**3 - y**2 +
520 3*z**3)

521 • Gröbner bases:

522 >>> x0, x1, x2 = symbols('x:3')

523 >>> I = [x0 + 2*x1 + 2*x2 - 1,
524 ... x0**2 + 2*x1**2 + 2*x2**2 - x0,
525 ... 2*x0*x1 + 2*x1*x2 - x1]
526 >>> groebner(I, order='lex')
527 GroebnerBasis([7*x0 - 420*x2**3 + 158*x2**2 + 8*x2 - 7,
528 7*x1 + 210*x2**3 - 79*x2**2 + 3*x2,
529 84*x2**4 - 40*x2**3 + x2**2 + x2], x0, x1, x2, domain='ZZ',
530 order='lex')

531 • Root isolation:

532 >>> f = 7*z**4 - 19*z**3 + 20*z**2 + 17*z + 20

533 >>> intervals(f, all=True, eps=0.001)
534 ([],
535 [((-425/1024 - 625*I/1024, -1485/3584 - 2185*I/3584), 1),
536 ((-425/1024 + 2185*I/3584, -1485/3584 + 625*I/1024), 1),
537 ((3175/1792 - 2605*I/1792, 1815/1024 - 10415*I/7168), 1),
538 ((3175/1792 + 10415*I/7168, 1815/1024 + 2605*I/1792), 1)])

539 14.3 Solvers

540 • Single solution:

541 >>> solveset(x - 1, x)

542 {1}

543 • Finite solution set, quadratic equation:

544 >>> solveset(x**2 - pi**2, x)

545 {-pi, pi}

546 • No solution:

547 >>> solveset(1, x)

548 EmptySet()

549 • Interval solution:

550 >>> solveset(x**2 - 3 > 0, x, domain=S.Reals)

551 (-oo, -sqrt(3)) U (sqrt(3), oo)

552 • Infinitely many solutions:

553 >>> solveset(x - x, x, domain=S.Reals)

554 (-oo, oo)
555 >>> solveset(x - x, x, domain=S.Complexes)
556 S.Complexes

14/16
557 • Linear systems (linsolve)

558 >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])

559 >>> b = Matrix([3, 6, 9])
560 >>> linsolve((A, b), x, y, z)
561 {(-1, 2, 0)}
562 >>> linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z))
563 {(-y - 1, y, 2)}

564 Below are examples of solve applied to problems not yet handled by solveset.

565 • Nonlinear (multivariate) system of equations (the intersection of a circle and a parabola):

566 >>> solve([x**2 + y**2 - 16, 4*x - y**2 + 6], x, y)

567 [(-2 + sqrt(14), -sqrt(-2 + 4*sqrt(14))),
568 (-2 + sqrt(14), sqrt(-2 + 4*sqrt(14))),
569 (-sqrt(14) - 2, -I*sqrt(2 + 4*sqrt(14))),
570 (-sqrt(14) - 2, I*sqrt(2 + 4*sqrt(14)))]

571 • Transcendental equations:

572 >>> solve((x + log(x))**2 - 5*(x + log(x)) + 6, x)

573 [LambertW(exp(2)), LambertW(exp(3))]
574 >>> solve(x**3 + exp(x))
575 [-3*LambertW((-1)**(2/3)/3)]

576 14.4 Matrices

577 • Matrix expressions

578 >>> m, n, p = symbols('m n p', integer=True)

579 >>> R = MatrixSymbol('R', m, n)
580 >>> S = MatrixSymbol('S', n, p)
581 >>> T = MatrixSymbol('T', m, p)
582 >>> U = R*S + 2*T
583 >>> U.shape
584 (m, p)
585 >>> U[0, 1]
586 2*T[0, 1] + Sum(R[0, _k]*S[_k, 1], (_k, 0, n - 1))

587 • Block Matrices

588 >>> n, m, l = symbols('n m l')

589 >>> X = MatrixSymbol('X', n, n)
590 >>> Y = MatrixSymbol('Y', m ,m)
591 >>> Z = MatrixSymbol('Z', n, m)
592 >>> B = BlockMatrix([[X, Z], [ZeroMatrix(m, n), Y]])
593 >>> B
594 Matrix([
595 [X, Z],
596 [0, Y]])
597 >>> B[0, 0]
598 X[0, 0]
599 >>> B.shape
600 (m + n, m + n)

15/16
601 15 REFERENCES
602 REFERENCES
603
[1] Abramowitz, M. and Stegun, I. A. (1964). Handbook of Mathematical Functions with
604 Formulas, Graphs, and Mathematical Tables. Dover Publications, New York, NY, USA, ninth
605 printing edition.
606
[Brent and Zimmermann] Brent, R. P. and Zimmermann, P. Modern Computer Arithmetic.

607 Cambridge University Press, version 0.5.1 edition.

608
[3] Cimrman, R. (2014). SfePy - write your own FE application. In de Buyl, P. and Varoquaux,

609 N., editors, Proceedings of the 6th European Conference on Python in Science (EuroSciPy
610 2013), pages 65–70. http://arxiv.org/abs/1404.6391.
611
[4] Fateman, R. J. (1992). A review of Mathematica. Journal of Symbolic Computation,

612 13(5):545–579.
613
[5] Ferguson, H. R. P., Bailey, D. H., and Arno, S. (1999). Analysis of PSLQ, an integer relation

614 finding algorithm. Mathematics of Computation, 68(225):351–369.

615
[6] Fung, Y. C. (1993). A first course in continuum mechanics. Pearson, third edition edition.

616
[7] Gruntz, D. (1996). On Computing Limits in a Symbolic Manipulation System. PhD thesis,

617 Swiss Federal Institute of Technology, Zürich, Switzerland.

618
[8] Gruntz, D. and Koepf, W. (1993). Formal power series.

619
[9] Manssur, L. R. U., Portugal, R., and Svaiter, B. F. (2002). Group-theoretic approach for

620 symbolic tensor manipulation. International Journal of Modern Physics C, 13.

621
[10] Martín-García, J. (2002-2016). xAct, efficient tensor computer algebra.

622
[11] Moskewicz, M., Madigan, C., and Malik, S. (2008). Method and system for efficient

623 implementation of boolean satisfiability. US Patent 7,418,369.

624
[12] Peeters, K. (2007). Cadabra: a field-theory motivated symbolic computer algebra system.

625 Computer Physics Communications.

626
[13] Rocklin, M. and Terrel, A. R. (2012). Symbolic statistics with SymPy. Computing in Science

627 and Engineering, 14.

628
[14] Sofroniou, M. and Spaletta, G. (2005). Precise numerical computation. Journal of Logic

629 and Algebraic Programming, 64(1):113–134.

630
[15] Solin, P., Segeth, K., and Dolezel, I. (2003). Higher-Order Finite Element Methods. Chapman

631 & Hall / CRC Press.

632
[16] Wolfram, S. (2003). The Mathematica Book. Wolfram Media, Champaign, IL, USA, fifth

633 edition.
634
[17] Zienkiewicz, O., Taylor, R., and Zhu, J. (2013). The Finite Element Method: Its Basis and

635 Fundamentals. Butterworth-Heinemann, seventh edition edition.

16/16