Computer Algebra Software: Mathematica, SymPy, GAP, Pari/GP

a side-by-side reference sheet

mathematica sympy gap pari/gp
version used

10.0 Python 2.7; SymPy 0.7 4.7 2.7
show version

select About Mathematica in Mathematica menu sympy.__version__ $gap.sh -h$ gp --version
implicit prologue import sympy

sympy.init_printing()
grammar and invocation
mathematica sympy gap pari/gp
interpreter

if foo.py imports sympy:
$python foo.py$ cat hello.gp
print("Hello, World!")
quit

$gp -q hello.gp Hello, World! repl$ MathKernel $python >>> import sympy$ gap.sh \$ gp
block delimiters

( stmt; ) : and offside rule function( ) end
if then elif then else fi
while do od
for do od
{ }

braces cannot be nested
statement separator ; or sometimes newline

A semicolon suppresses echoing value of previous expression.
newline or ;

newlines not separators inside (), [], {}, triple quote literals, or after backslash: \
;

Two trailing semicolons ;; suppress echoing value of previous expression.
newline or ;

Newlines don't separate statements inside braces.

A semicolon suppresses echoing value of previous expression.
end-of-line comment

none 1 + 1 # addition 1 + 1; # addition 1 + 1 \\ addition
multiple line comment

1 + (* addition *) 1 none none 1 + /* addition */ 1
variables and expressions
mathematica sympy gap pari/gp
assignment a = 3
Set[a, 3]

(* delayed assignment: *)
a := x + 3
SetDelayed[a, x + 3]
a = 3 a := 3; x = 3.14
parallel assignment {a, b} = {3, 4}
Set[{a, b}, {3, 4}]
a, b = 3, 4 none [a, b] = [3, 4]
compound assignment += -= *= /=
corresponding functions:
AddTo SubtractFrom TimeBy DivideBy
+= -= *= /= //= %= **= none += -= *= /= %=
increment and decrement ++x --x
PreIncrement[x] PreDecrement[x]
x++ x--
Increment[x] Decrement[x]
none none postmodifiers:
x++ x--
non-referential identifier any unassigned identifier is non-referential x, y, z, w = sympy.symbols('x y z w') none any unassigned identifier is non-referential
identifier as value x = 3
x // HoldForm
x = 3
'x
global variable variables are global by default g1, g2 = 7, 8

def swap_globals():
global g1, g2
g1, g2 = g2, g1
variables are global by default
local variable Module[{x = 3, y = 4}, Print[x + y]]

(* makes x and y read-only: *)
With[{x = 3, y = 4}, Print[x + y]]
assignments inside functions are to local variables by default tmp = 19

add(x, y, z) = {
\\ don't overwrite global tmp:
my(tmp = x + y);
tmp + z
}
null

Null None none
null test

x == Null x is None none
undefined variable access

treated as an unknown number raises NameError error treated as an unknown number
remove variable binding Clear[x]
Remove[x]
del x kill(x)
conditional expression

If[x > 0, x, -x] x if x > 0 else -x if(x > 0, x, -x)
arithmetic and logic
mathematica sympy gap pari/gp
true and false

True False True False true false 1 0
falsehoods

False False 0 0.0 false 0
logical operators ! True || (True && False)
Or[Not[True], And[True, False]]
sympy.Or(sympy.Not(True), sympy.And(True, False))

# when arguments are symbols:
~ x | (y & z)
not true or (true and false) && || !
relational operators == != > < >= <=
corresponding functions:
Equal Unequal Greater Less GreaterEqual LessEqual
sympy.Eq sympy.Ne sympy.Gt sympy.Lt sympy.Ge sympy.Le

# when arguments are symbols:
== != > < >= <=
= <> < > <= >= == != > < >= <=
arithmetic operators + - * / Quotient Mod
adjacent terms are multiplied, so * is not necessary. Quotient and Mod are functions, not binary infix operators. These functions are also available:
Plus Subtract Times Divide
+ - * / ?? %

if an expression contains a symbol, then the above operators are rewritten using the following classes:
sympy.Add sympy.Mul sympy.Pow sympy.Mod
+ - * / mod

the operators + - * / are overloaded for integers, rationals, and floats; other arithmetic functions aren't and there are no implicit conversions; use constructors to convert:
Rat(3.1)
Float(3)
Float(31/10)
+ - * / %
integer division

Quotient[a, b] QuoInt(a, b); a \ b
divrem(a, b)[1]
integer division by zero dividend is zero:
Indeterminate
otherwise:
ComplexInfinity
error error
float division exact division:
a / b
depending upon the types of a and b, the value can be an exact rational, a machine float, or an arbitrary precision float:
a / b
7 /3
float division by zero dividend is zero:
Indeterminate
otherwise:
ComplexInfinity
error error
power 2 ^ 32
Power[2, 32]
2 ** 32
sympy.Pow(2, 32)
2 ^ 32 2 ^ 32
sqrt returns symbolic expression:
Sqrt[2]
sympy.sqrt(2) 2.0 ^ 0.5 sqrt(2)
sqrt -1

I sympy.I -1.0 ^ 0.5 evaluates to -1. 1.000 * I
transcendental functions Exp Log
Sin Cos Tan
ArcSin ArcCos ArcTan
ArcTan
ArcTan accepts 1 or 2 arguments
symp.exp sympy.log
sympy.sin sympy.cos sympy.tan
sympy.asin sympy.acos sympy.atan
sympy.atan2
arguments must be floats; no implicit conversion of integers to floats:
Exp Log
Sin Cos Tan
Asin Acos Atan
Atan2(y, x)
exp log none
sin cos tan
asin acos atan
none
transcendental constants
π and Euler's number
Pi E sympy.pi sympy.E FLOAT.PI FLOAT.E Pi exp(1)
float truncation
round towards zero, round to nearest integer, round down, round up
IntegerPart Round Floor Ceiling sympy.floor
sympy.ceiling
Trunc Round Floor Ceil truncate(x)
round(x)
floor(x)
ceil(x)
absolute value
and signum
Abs Sign sympy.Abs sympy.sign AbsInt
no absolute value for floats?
SignInt
SignFloat
abs(x)
sign(x)
integer overflow

none, has arbitrary length integer type none, has arbitrary length integer type none, has arbitrary length integer type none, has arbitrary length integer type
float overflow

none # prints as inf:
FLOAT.INFINTY
error
rational construction 2 / 7 sympy.Mul(2, sympy.Pow(7, -1)) 2 / 7 2 / 7
rational decomposition

Numerator[x/y]
Denominator[x/y]
numer, denom = sympy.fraction(x, y) x := 2 / 7;
NumeratorRat(x);
DenominatorRat(x);
x = 2 / 7
numerator(x)
denominator(x)
decimal approximation N[2 / 7]
2 / 7 + 0.
2 / 7 // N
N[2 / 7, 100]
sympy.N(sympy.Rational(2, 7))
sympy.N(sympy.Rational(2, 7), 100)
2 / 7 + 0.

\\ change precision to 100:
\p 100
2 / 7 + 0.
complex construction

1 + 3I 1 + 3 * sympy.I none 1 + 3 * I
complex decomposition
real and imaginary part, argument and modulus, conjugate
Re Im
Arg Abs
Conjugate
sympy.re sympy.im
sympy.Abs sympy.arg
sympy.conjugate
none real(z) imag(z)
arg(z) abs(z)
conj(z)
random number
uniform integer, uniform float
RandomInteger[{0, 99}]
RandomReal[]
rs := RandomSource(IsMersenneTwister);
Random(rs, 0, 99);
??
random(100)
??
random seed
set, get
SeedRandom[17]
??
rs := RandomSource(IsMersenneTwister, 17);
State(rs);
bit operators BitAnd[5, 1]
BitOr[5, 1]
BitXor[5, 1]
BitNot[5]
BitShiftLeft[5, 1]
BitShiftRight[5, 1]
none setrand(17)
getrand()
binary, octal, and hex literals 2^^101010
8^^52
16^^2a
none
BaseForm[7^^60, 10]
none \\ 42 as powers of 7 up to 9th power:
42 + O(7^10)
strings
mathematica sympy gap pari/gp
string literals "don't say \"no\"" use Python strings "don't say \"no\"" "don't say \"no\""
newline in literal yes no no; use \n escape
string literal escapes \\ \" \b \f \n \r \t \ooo \b \c \n \r \" \' \\ \ooo

when writing to a buffered output stream, encountering a \c causes a flush of output.
\n \t \" \\
character access Characters["hello"][[1]] s := "hello";
# the character 'h':
s[1];

# cannot use index notation on string literal
length StringLength["hello"] Length("hello"); length("hello")
concatenate "one " <> "two " <> "three" Concatenation("one", "two", "three"); Str("one", "two", "three")
index of substring StringPosition["hello", "el"][[1]][[1]]
StringPosition returns an array of pairs, one for each occurrence of the substring. Each pair contains the index of the first and last character of the occurrence.
extract substring StringTake["hello", {1, 4}] s := "hello";
s{[2..3]};
split StringSplit["foo,bar,baz", ","] SplitString("foo,bar,baz", ",");
join StringJoin[Riffle[{"foo", "bar", "baz"}, ","]] a := ["foo", "bar", "baz"];
JoinStringsWithSeparator(a, ",");
substitution s = "do re mi mi"
re = RegularExpression["mi"]

(* replace first occurrence: *)
StringReplace[s, re -> "ma", 1]
(* replace all occurrences: *)
StringReplace[s, re -> "ma"]
# replace all occurrences:
ReplacedString("do re mi mi", "mi", "ma");
trim StringTrim[" foo "] none
number to string "value: " <> ToString[8] Concatenation("value: ", String(8));
string to number 7 + ToExpression["12"]
73.9 + ToExpression[".037"]
7 + Int("12");
73.9 + Float(".037");
case manipulation ToUpperCase["foo"]
ToLowerCase["FOO"]
UppercaseString("foo");
LowercaseString("FOO");
character literal 'h'
chr and ord FromCharacterCode[{65}]
ToCharacterCode["A"][[1]]
CharInt(65)
IntChar('A')
delete characters s := "disemvowel me";
# no retval; modifies s in place:
RemoveCharacters(s, "aeiou");
resizable arrays
mathematica sympy gap pari/gp
literal {1, 2, 3}

List[1, 2, 3]
use Python lists [1, 2, 3];

# creates array with gap at fourth index;
# reading a[4] causes an error:

a := [1, 2, 3, , 5];
\\ [1, 2, 3] is a vector literal:
List([1, 2, 3])
size

Length[{1, 2, 3}] Length([1, 2, 3]); length(List([1, 2, 3]))
#List([1, 2, 3])
lookup (* access time is O(1) *)
(* indices start at one: *)
{1, 2, 3}[[1]]

Part[{1, 2, 3}, 1]
# indices start at one:
a := [1, 2, 3];
a[1];
\\ access time is O(1).
\\ indices start at one:
List([1, 2, 3])[1]
update

a[[1]] = 7 a[1] := 7; listput(a, 7, 1)
out-of-bounds behavior left as unevaluated Part[] expression Lookups result in errors; arrays can have gaps which also cause lookup errors.

An update will expand the array, possibly creating gaps.
out of allowed range error
element index (* Position returns list of all positions: *)
First /@ Position[{7, 8, 9, 9}, 9]
# returns 3:
Position([7, 8, 9, 9], 9);

# returns [3, 4]:
Positions([7, 8, 9, 9], 9);
none
slice

{1, 2, 3}[[1 ;; 2]] none
array of integers as index (* evaluates to {7, 9, 9} *)
{7, 8, 9}[[{1, 3, 3}]]
none
manipulate back a = {6,7,8}
AppendTo[a, 9]
elem = a[[Length[a]]]
a = Delete[a, Length[a]]
elem
a = [6, 7, 8];
elem := Remove(a);
a = List([6, 7, 8])
listput(a, 9)
elem = listpop(a)
manipulate front a = {6,7,8}
PrependTo[a, 5]
elem = a[[1]]
a = Delete[a, 1]
elem
a = List([6, 7, 8]);
listinsert(a, 5, 1);
elem = a[1];
listpop(a, 1);

First[{1, 2, 3}] List([1, 2, 3])[1]
tail

Rest[{1, 2, 3}] none
cons (* first arg must be an array *)
Prepend[{2, 3}, 1]
a = List([1, 2, 3]);
listinsert(a, 1, 1);
concatenate

Join[{1, 2, 3}, {4, 5, 6}] Concatenation([1, 2, 3], [4, 5, 6]); concat(List([1, 2, 3]), List([4, 5, 6]))
replicate

ten_zeros = Table[0, {i, 0, 9}]
copy

a2 = a a2 = a
iterate

Function[x, Print[x]] /@ {1, 2, 3} Perform([1, 2, 3], function(x) Print(x); Print("\n"); end); a = List([1, 2, 3])

for(i=1, length(a), print(a[i]))
reverse

Reverse[{1, 2, 3}] Reversed([1, 2, 3]) a = List([1, 2, 3])
a2 = listcreate()
while(i > 0, listput(a2, a[i]); i—)
sort Sort[{3, 1, 4, 2}] A := [3, 1, 4, 2]
Sort(A);
a = List([3,1,4,2])
listsort(a)
a
dedupe

DeleteDuplicates[{1, 2, 2, 3}] Set([1, 2, 2, 3]);
Unique([1, 2, 2, 3]);
Set([1, 2, 2, 3])
membership

MemberQ[{1, 2, 3}, 2] 2 in [1, 2, 3] \\ returns 1-based index of first occurrence

setsearch([1, 2, 3], 2)
intersection Intersect[{1, 2}, {2, 3, 4}] Intersection(Set([1, 2]), Set([2, 3, 4])); setintersect([1, 2], [2, 3, 4])
union

Union[{1, 2}, {2, 3, 4}] Union(Set([1, 2]), Set([2, 3, 4])); setunion([1, 2], [2, 3, 4])
relative complement, symmetric difference Complement[{1, 2, 3}, {2}]
none
setminus([1, 2, 3], [2])
??
map Map[Function[x, x x], {1, 2, 3}]

Function[x, x x] /@ {1, 2, 3}

(* if function has Listable attribute, Map is unnecessary: *)
sqr[x_] := x * x
SetAttributes[sqr, Listable]
sqr[{1, 2, 3, 4}]
A := [1, 2, 3];

# modifies A:
Apply(A, x -> x * x);
filter

Select[{1, 2, 3}, # > 2 &]
reduce Fold[Plus, 0, {1, 2, 3}]
universal and existential tests none
min and max element Min[{1, 2, 3}]
Max[{1, 2, 3}]
Minimum([1, 2, 3])
Maximum([1, 2, 3])
shuffle and sample x = {3, 7, 5, 12, 19, 8, 4}

RandomSample[x]
RandomSample[x, 3]
Shuffle([1, 2, 3, 4])
flatten
one level, completely
Flatten[{1, {2, {3, 4}}}, 1]
Flatten[{1, {2, {3, 4}}}]
# completely:
Flat([1, [2, [3, 4]]])
zip (* list of six elements: *)
Riffle[{1, 2, 3}, {"a", "b", "c"}]

(* list of lists with two elements: *)
Inner[List, {1, 2, 3}, {"a", "b", "c"}, List]

(* same as Dot[{1, 2, 3}, {2, 3, 4}]: *)
Inner[Times, {1, 2, 3}, {2, 3, 4}, Plus]
cartesian product Outer[List, {1, 2, 3}, {"a", "b", "c"}] Cartesian([1, 2, 3], ["a", "b", "c"])
arithmetic sequences
mathematica sympy gap pari/gp
unit difference Range[1, 100] range(1, 101) [1 .. 100] vector(100, i, i)
difference of 10 Range[1, 100, 10] range(1, 100, 10) [1,11 .. 91] vector(10, i, 10 * i - 9)
difference of 0.1 Range[1, 100, .1] none vector(1000 - 9, i, i / 10 + 9 / 10)
dictionaries
mathematica sympy gap pari/gp
literal

d = <|"t" -> 1, "f" -> 0|>

(* or convert list of rules: *)
d = Association[{"t" -> 1, "f" -> 0}]
(* and back to list of rules: *)
Normal[d]
use Python dictionaries
size

Length[Keys[d]]
lookup

d["t"]
update d["f"] = -1
missing key behavior

Returns a symbolic expression with head "Missing". If the lookup key was "x", the expression is:

Missing["KeyAbsent", "x"]
is key present

KeyExistsQ[d, "t"]
delete
from array of pairs, from even length array
merge
invert
iterate

keys and values as arrays Keys[d]
Values[d]
sort by values Sort[d]
default value, computed value
functions
mathematica sympy gap pari/gp
define function Add[a_, b_] := a + b

(* alternate syntax: *)
Add = Function[{a, b}, a + b]
add := function(x, y)
return x + y;
end;
add(x, y) = x + y

\\ function body w/ sequence of statements:
say(s1, s2, s3) = print(s1); print(s2); print(s3)

\\ function body w/ newlines:
dire(s1, s2, s3) = {
print(s1);
print(s2);
print(s3);
}
invoke function Add[3, 7]

Add @@ {3, 7}

(* syntax for unary functions: *)
2 // Log
boolean function attributes
list, set, clear
SetAttributes[add, {Orderless, Flat, Listable}]
redefine function Add[a_, b_] := b + a add(x, y, z) = x + y + z
missing function behavior The expression is left unevaluated. The head is the function name as a symbol, and the parts are the arguments. "not a function" error
missing argument behavior The expression is left unevaluated. The head is the function name as a symbol, and the parts are the arguments. set to zero
extra argument behavior The expression is left unevaluated. The head is the function name as a symbol, and the parts are the arguments. "too many parameters" error
default argument Options[myLog] = {base -> 10}
myLog[x_, OptionsPattern[]] :=
N[Log[x]/Log[OptionValue[base]]]

(* call using default: *)
myLog[100]

(* override default: *)
myLog[100, base -> E]
mylog(x = 1, base = 10) = log(x) / log(base)

\\ log10(3):
mylog(3)
\\ ln(3):
mylog(3, exp(1))
\\ ln(1):
mylog(, exp(1))

\\ If neither caller nor function definition
\\ provide a value, zero is used.
return value last expression evaluated, or argument of Return[]
anonymous function Function[{a, b}, a + b]

(#1 + #2) &
# unary functions only?
f := x -> x * x;

f2 := function(x, y) return 2 * x + 3 * y; end;
f = (x, y) -> x + y

f(1, 2)
variable number of arguments (* one or more arguments: *)

(* zero or more arguments: *)
pass array elements as separate arguments Apply[f, {a, b, c}]

f @@ {x, y, z}
a = [x, y, z]
f(*a)
execution control
mathematica sympy gap pari/gp
if If[x > 0,
Print["positive"],
If[x < 0,
Print["negative"],
Print["zero"]]]
use Python execution control if x > 0 then
Print("positive\n");
elif x < 0 then
Print("negative\n");
else
Print("zero\n");
fi;
if(x > 0, \
print("positive"), \
if(x < 0, \
print("negative"), \
print("zero")))
while i = 0
While[i < 10, Print[i]; i++]
i := 0;
while i < 10 do
Print(i, "\n");
i := i + 1;
od;
i = 0
while(i < 10, print(i); i++)
for For[i = 0, i < 10, i++, Print[i]] for i in [0..9] do
Print(i, "\n");
od;
for(i = 0, 9, print(i))
break Break[] break break
continue Continue[] continue next
exceptions
mathematica sympy gap pari/gp
raise exception Throw["failed"] use Python exceptions Error("failed"); error("failed")
handle exception Print[Catch[Throw["failed"]]] iferr(error("failed"), E, \
print(errname(E), ": ", component(E, 1)))
uncaught exception behavior Error() invokes the GAP debugger. Type

Quit;

file handles
mathematica sympy gap pari/gp
write to stdout Print["hello"] Print("hello"); print("hello")
read entire file into string or array s = Import["/etc/hosts"]
a = StringSplit[s, "\n"]
libraries and namespaces
mathematica sympy gap pari/gp
define library
library path
reflection
mathematica sympy gap pari/gp
list function documentation ?
get function documentation ?Tan
Information[Tan]
print(sympy.solve.__doc__)

# in IPython:
sympy.solve?
help(sympy.solve)
? tan
list function options Options[Solve]
Options[Plot]
query data type Head[x] type(x)
list types \t
list variables in scope variable()
list built-in functions ?*
list metacommands ?\
search documentation ??DirectProduct
vectors
mathematica sympy gap pari/gp
vector literal (* row vector is same as array: *)
{1, 2, 3}
# column vector:
sympy.Matrix([1, 2, 3])
# row vector is same as array:
[1, 2, 3]
[1, 2, 3]
constant vector

all zeros, all ones
Table[0, {i, 1, 3}]
Table[1, {i, 1, 3}]
vector(3, i, 0)
vector(3, i, 1)
vector coordinate (* indices start at one: *)
{1,v2, 3}[[1]]
vec := [1, 2, 3];
# indices start at one:
v[1];
\\ indices start at one:
[1, 2, 3][1]
vector dimension

Length[{1, 2, 3}] Length([1, 2, 3]) length([1, 2, 3])
#[1, 2, 3]
element-wise arithmetic operators + - * /
adjacent lists are multiplied element-wise
+ - * / + -
vector length mismatch

error shorter vector is zero-padded error
scalar multiplication 3 {1, 2, 3}
{1, 2, 3} 3
* may also be used
3 * [1, 2, 3];
[1, 2, 3] * 3;
3 * [1, 2, 3]
[1, 2, 3] * 3
dot product {1, 1, 1} . {2, 2, 2}
Dot[{1, 1, 1}, {2, 2, 2}]
v1 = sympy.Matrix([1, 1, 1])
v2 = sympy.Matrix([2, 2, 2])
v1.dot(v2)
[1, 1, 1] * [2, 2, 2]
cross product Cross[{1, 0, 0}, {0, 1, 0}] e1 = sympy.Matrix([1, 0, 0])
e2 = sympy.Matrix([0, 1, 0])
e1.cross(e2)
norms Norm[{1, 2, 3}, 1]
Norm[{1, 2, 3}]
Norm[{1, 2, 3}, Infinity]
vec = sympy.Matrix([1, 2, 3])

vec.norm(1)
vec.norm()
vec.norm(inf)
vec = [1, 2, 3]

normlp(vec, 1)
normlp(vec, 2)
normlp(vec)
orthonormal basis Orthogonalize[{{1, 0, 1}, {1, 1, 1}}]
matrices
mathematica sympy gap pari/gp
literal (* used a nested array for each row: *)
{{1, 2}, {3, 4}}

(* display as grid with aligned columns: *)
MatrixForm[{{1, 2}, {3, 4}}]
sympy.Matrix([[1, 2], [3, 4]]) [[1, 2], [3, 4]] [1, 2; 3, 4]

\\ from rows:
row1 = [1, 2]
row2 = [3, 4]
matconcat([row1; row2])
construct from sequence ArrayReshape[{1, 2, 3, 4, 5, 6}, {2, 3}] sympy.Matrix(2, 3, [1, 2, 3, 4, 5, 6])
construct from columns col1 = [1, 3]~
col2 = [2, 4]~
matconcat([col1, col2])
construct from submatrices A = [1, 2; 3, 4]
B = [4, 3; 2, 1]
\\ 4x4 matrix:
C = matconcat([A, B; B, A])
constant matrices Table[0, {i, 3}, {j, 3}]
Table[1, {i, 3}, {j, 3}]
sympy.zeros(3, 3)
sympy.ones(3, 3)
matrix(3, 3, i, j, 0)
matrix(3, 3, i, j, 1)

\\ 3x3 Hilbert matrix:
matrix(3, 3, i, j, 1 / (i + j - 1))
diagonal matrices
and identity
DiagonalMatrix[{1, 2, 3}]
IdentityMatrix[3]
sympy.diag(*[1, 2, 3])
sympy.eye(3)
DiagonalMat([1, 2, 3])
IdentityMat(3)
matdiagonal([1, 2, 3])
matid(3)
dimensions (* returns {3, 2}: *)
Dimensions[{{1, 2}, {3, 4}, {5, 6}}]
A = sympy.matrix([[1, 2], [3, 4], [5, 6]])

# returns (3, 2):
A.shape
# returns [3, 2]:
DimensionsMat([[1, 2], [3, 4], [5, 6]])
\\ [3, 2]:
matsize([1, 2; 3, 4; 5, 6])
element lookup (* top left corner: *)
{{1, 2}, {3, 4}}[[1, 1]]
A = sympy.Matrix([[1, 2], [3, 4]])

# top left corner:
A[0, 0]
A := [[1, 2], [3, 4]];

# top left corner:
A[1][1]
\\ top left corner:
A[1, 1]
extract row (* first row: *)
{{1, 2}, {3, 4}}[[1]]
# first row:
A[0, :]
\\ first row:
[1, 2; 3, 4][1, ]
extract column (* first column as array: *)
{{1, 2}, {3, 4}}[[All, 1]]
# first column as 1x2 matrix:
A[:, 0]
\\ first column:
[1, 2; 3, 4][, 1]
extract submatrix A = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}
A[[1;;2, 1;;2]]
rows = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
A = sympy.Matrix(rows)
A[0:2, 0:2]
A = [1, 2, 3; 4, 5, 6; 7, 8, 9]

vecextract(A, "1..2", "1..2")
element-wise operators + - * /
adjacent matrices are multiplied element-wise
+ -

# for Hadamard product:
A.multiply_elementwise(B)
+ - + -
product A = {{1, 2}, {3, 4}}
B = {{4, 3}, {2, 1}}
Dot[A, B]
(* or use period: *)
A . B
A = sympy.matrix([[1, 2], [3, 4]])
B = sympy.matrix([[4, 3], [2, 1]])
A * B
A := [[1, 2], [3, 4]];
B := [[4, 3], [2, 1]];
A * B;
A = [1, 2; 3, 4]
B = [4, 3; 2, 1]
A * B
power MatrixPower[{{1, 2}, {3, 4}}, 3]

(* element-wise operator: *)
A ^ 3
A ** 3 [[1, 2], [3, 4]] ^ 3 [1, 2; 3, 4] ^ 3
exponent MatrixExp[{{1, 2}, {3, 4}}]
log MatrixLog[{{1, 2}, {3, 4}}]
kronecker product A = {{1, 2}, {3, 4}}
B = {{4, 3}, {2, 1}}
KroneckerProduct[A, B]
A := [[1, 2], [3, 4]];
B := [[4, 3], [2, 1]];
KroneckerProduct(A, B);
norms A = {{1, 2}, {3, 4}}

Norm[A, 1]
Norm[A, 2]
Norm[A, Infinity]
Norm[A, "Frobenius"]
transpose Transpose[{{1, 2}, {3, 4}}]

(* or ESC tr ESC for T exponent notation *)
A.T A~
mattranspose(A)
conjugate transpose A = {{1, I}, {2, -I}}
ConjugateTranspose[A]

(* or ESC ct ESC for dagger exponent notation *)
M = sympy.Matrix([[1, sympy.I], [2, -sympy.I]])
inverse Inverse[{{1, 2}, {3, 4}}] A.inv() Inverse([[1, 2], [3, 4]])
row echelon form RowReduce[{{1, 1}, {1, 1}}]
pseudoinverse PseudoInverse[{{1, 0}, {3, 0}}]
determinant Det[{{1, 2}, {3, 4}}] A.det() Determinant([[1, 2], [3, 4]]) matdet([1, 2; 3, 4])
trace Tr[{{1, 2}, {3, 4}}] Trace([[1, 2], [3, 4]]) trace([1, 2; 3, 4])
rank MatrixRank[{{1, 1}, {0, 0}}] RankMat([[1, 1], [0, 0]]) matrank([1, 1; 0, 0])
nullspace basis NullSpace[{{1, 1}, {0, 0}}] matker([1, 1; 0, 0])
range basis matimage([1, 1; 0, 0])
eigenvalues Eigenvalues[{{1, 2}, {3, 4}}] A.eigenvals()
eigenvectors Eigenvectors[{{1, 2}, {3, 4}}] A.eigenvects() mateigen([1, 2; 3, 4])
singular value decomposition SingularValueDecomposition[{{1, 1}, {1, 0}}]
qr decomposition QRDecomposition[{{1, 2}, {3, 4}}] matqr([1, 2; 3, 4])
solve system of equations A = [1, 2; 3, 4]
matsolve(A, [2, 3]~)
symbolic expressions
mathematica sympy gap pari/gp
replace head Apply[Times, Plus[x, 3]]
Times @@ Plus[x, 3]
Times @@ (x + 3)
prevent simplification HoldForm[1 + 2]
1 + 2 // HoldForm
list evaluation steps Trace[Apply[And, Map[EvenQ, {2, 3, 4}]]]
expand polynomial Expand[(1 + x)^5] sympy.expand((1+x)**5)
factor polynomial Factor[3 + 10 x + 9 x^2 + 2 x^3] sympy.factor(3 + 10*x + 9*x**2 + 2*x**3)
collect terms (* write as polynomial in x: *)
Collect[(1 + x + y)^3, x]
sympy.collect(sympy.expand((x+y+1)3), x)
add fractions Together[a/b + c/d] sympy.together(x/y + z/w)
partial fraction decomposition Apart[(b c + a d)/(b d)] # only one symbol allowed in denominator:
sympy.apart((3*x+2) / (x*(x+1)))
eliminate sums and multiples inside trig functions TrigExpand[Sin[2 x + y]]
eliminate powers of trig functions TrigReduce[Sin[x]^2]
trig to complex exponential
complex exponential to trig
quick simplify
slow simplify
simplify with assumption Assuming and Refine[]
calculus
mathematica sympy gap pari/gp
limit Limit[Sin[x]/x, x -> 0]
residue
differentiation D[x^3 + x + 3, x] sympy.diff(x**3 + x + 3, x) P = x^3 + x + 3
P'
sin(x)'
deriv(y^2 + 2, y)
higher order differentiation D[Log[x], {x, 3}] sympy.diff(sympy.log(x), x, 3)
unevaluated derivative Hold[D[x^2, x]] sympy.Derivative(x2, x)
mixed partial derivative
antiderivative Integrate[x^3 + x + 3, x] sympy.integrate(x**3 + x + 3, x)
definite integral Integrate[x^3 + x + 3, {x, 0, 1}] sympy.integrate(x**3 + x + 3, [x, 0, 1])
improper integral sympy.integrate(sympy.exp(-x), (x, 0, sympy.oo))
equations and unknowns
mathematica sympy gap pari/gp
solution to an equation Solve[x^3 + x + 3 == 0, x] solve(x**3 + x + 3, x)
solution to two equations Solve[x + y == 3 && x == 2y,
{x, y}]
solve([x + y - 3, 3*x - 2*y], [x, y])
solve diophantine equation Solve[a^2 + b^2 == c^2 &&
a > 0 && a < 10 &&
b > 0 && b < 10 &&
c > 0 && c < 10,
{a, b, c}, Integers]
optimization
mathematica sympy gap pari/gp
minimize Minimize[x^2 + 1, x]
minimize subject to constraint
combinatorics
mathematica sympy gap pari/gp
factorial 10!
Factorial[10]
factorial(10) 10!
falling factorial FactorialPower[10, 3] 10! / (10 - 3)!
binomial coefficient Binomial[10, 3] binomial(10, 3) binomial(10, 3)
multinomial coefficient Multinomial[3, 4, 5]
integer partitions IntegerPartitions[10] partitions(10)
set partitions

and Bell number
StirlingS2[10, 3]
BellB[10]
stirling(10, 3, 2)
sum(i=1, 10, stirling(10, i, 2))
permutations with k disjoint cycles Abs[StirlingS1[n, k]] abs(stirling(n, k, 1))
fibonacci number

and lucas number
Fibonacci[10]
LucasL[10]
fibonacci(10)
??
bernoulli number BernoulliB[100] bernfrac(100)
number theory
mathematica sympy gap pari/gp
prime test PrimeQ[7] sympy.ntheory.primetest.isprime(7) IsPrimeInt(7); isprime(7)
pseudoprime test ispseudoprime(7)
divisors Divisors[100] divisors(100)
prime factors returns {{2, 2}, {3, 1}, {7, 1}}
FactorInteger[84]
sympy.ntheory.factorint(84) \\ [2,2; 3,1; 7,1]:
factor(84)
next prime

and preceding
NextPrime[1000] nextprime(1000)
precprime(1000)
nth prime Prime[100] \\ first 100 primes: primes(100)
primes(100)[100]
prime counting function PrimePi[100]
divmod QuotientRemainder[7, 3] divrem(7, 3)
coprime test CoprimeQ[14, 45]
greatest common divisor GCD[14, 21] sympy.igcd(14, 21) gcd(14, 21)
least common multiple LCM[14, 21] lcm(14, 21)
integer residues \\ (2 + 3) % 5:
Mod(2, 5) + Mod(3, 5)
chinese remainder theorem chinese(Mod(3, 17), Mod(8, 11))
lift integer residue to integer \\ 7:
lift(-17, 12)
\\ -5:
centerlift(-17, 12)
p-adic number \\ p is 2 and precision in powers of 2 is 100:
1/2 + O(2^100)
lift p-adic to rational lift(1/2 + O(2^100)
gaussian integer norm norm(1 + I)
Euler totient EulerPhi[256] sympy.ntheory.totient(256) Phi(256); eulerphi(256)
jacobi symbol

and kronecker symbol
JacobiSymbol[3, 5]
KroneckerSymbol[3, 5]
??
kronecker(3, 5)
dirichlet character Table[DirichletCharacter[2, 1, i], {i, 100}]
moebius function MoebiusMu[11] moebius(11)
mangoldt lambda MangoldtLambda[11]
digits (* base 10: *)
IntegerDigits[1234]
(* base 2: *)
IntegerDigits[1234, 2]
\\ base 10:
digits(1234)
\\ base 2:
digits(1234, 2)
\\ number of digits in base 10:
sizedigits(1234)
to continued fraction cf = ContinuedFraction[Pi, 100] \p 100
contfrac(Pi)
from continued fraction (* rational approx. of π: *)
FromContinudFraction[cf]
elliptic curves
mathematica sympy gap pari/gp
elliptic curve from coefficients \\ ellinit([a, b, c, d, e]) where
\\
\\   y^2 + axy + by = x^3 + cx^2 + dx + e
\\

e0 = ellinit([0,0,1,-7,6])

\\ ellinit([a, b]) where
\\
\\   y^2 = x^3 + ax + b
\\

e1 = ellinit([-1, 0])
discriminant e0.disc
conductor ellglobalred(e0)[1]
singularity test e0.disc == 0
convert to minimal model e0 = ellinit([6, -3, 9, -16, -14])
e = ellminimalmodel(e0)
coordinate transformation on point e0 = ellinit([6, -3, 9, -16, -14])
e = ellminimalmodel(e0, &v)
\\ minimal to original:
ellchangepointinv([0, 0], v)
\\ original to minimal:
ellchangepoint([-2, 2], v)
coordinate transformation on curve: ellchangecurve e0 = ellinit([6, -3, 9, -16, -14])
e = ellminimalmodel(e0, &v)
\\ same as e0:
ellchangecurve(e, v)
point on curve test ellisoncurve(e, [0, 2])
abscissa to ordinates \\ vector of size 0, 1, or 2:
ellordinate(e, 0)
group identity [0]
group operation elladd(e, [0, 2], [1, -1])
group inverse ellneg(e, [0, 2])
group multiplication ellmul(e, [0, 2], 3)
canonical height of point ellheight(e, [0, -3])
order of point \\ returns 0 for infinite order:
ellorder(e, [0, 2])

ellorder(e1, [0, 0])
torsion subgroup e1 = ellinit([-1, 0])

\\ returns [t, v1, v2]:
\\
\\   t: order of torsion group
\\   v1: orders of component cyclic groups
\\   v2: generators of same cyclic groups
\\

elltors(e1)
analytic rank \\ first value is rank:
[a, b] = ellanalyticrank(e)

\\ recompute second value to higher precision:
\p 100
b = ellL1(e, a)
L-function value elllseries(e, 1 + I)
L-function coefficients \\ tenth coefficient:
ellak(e, 10)
\\ first ten coefficients:
ellan(e, 10)
algebraic numbers
mathematica sympy gap pari/gp
quadratic extension \\ make w equal to sqrt(D)/4:
D = -4
quadratic number (1 + w)^2
polynomials
mathematica sympy gap pari/gp
from expression with indeterminates (x - 1) (x - 2)

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

(1+x)^2 * (2+y)^3
from coefficient array coeff = {1, -3, 2}
Plus @@ Table[coeff[[i]] * x^(3 - i), {i, 1, 3}]
Pol([1, -3, 2])

@@\\@ zero-degree coefficient first:
Polrev([2, -3, 1])
to coefficient array CoefficientList[(x + 1)^10, x]
lookup coefficient Coefficient[(1 + x)^10, x, 3] polcoeff((x+1)^10, 3)
substitute indeterminate \\ replace x with 3:
subst((x-1)*(x-2), x, 3)
\\ replace x with (x-1):
subst((x-1)*(x-2), x, (x-1))
degree Exponent[(x + 1)^10, x] poldegree((x-1)^10)
operations + - * / + - * /
division and remainder PolynomialReduce[x^10 - 1, x - 1, {x}]
PolynomialReduce[x^10 - y^10, x - y, {x, y}]
factor Factor[x^10 - y^10] factor(x^2-1)
roots Solve[x^3 + 3 x^2 + 2 x - 1 == 0, x] polroots(x^3+3*x^2+2*x-1)
greatest common divisor p1 = -2 - x + 2 x^2 + x^3
p2 = 6 - 7 x + x^3
PolynomialGCD[p1, p2]
p1 = x^3 + 2*x^2 -x -2
p2 = x^3 -7*x + 6
gcd(p1, p2)
resultant Resultant[(x-1)(x-2), (x-3)(x-3), x] polresultant((x-1)*(x-2), (x-3)^2)
discriminant Discriminant[(x + 1) (x - 2), x] poldisc((x+1)*(x-2))
groebner basis p1 = x^2 + y + z - 1
p2 = x + y^2 + z - 1
p3 = x + y + z^2 - 1

(* uses lexographic order by default: *)
GroebnerBasis[{p1, p2, p3}, {x, y, z}]
none
specify ordering none
symmetric polynomial SymmetricPolynomial[3, {x1, x2, x3, x4}] none
cyclotomic polynomial Cyclotomic[10, x] polcyclo(10)
hermite polynomial HermiteH[4, x] polhermite(4)
chebyshev polynomial

first and second kind
ChebyshevT[4, x]
ChebyshevU[4, x]
polchebyshev(4, 1)
polychebyshev(4, 2)
interpolation polynomial pts = Inner[List, {1, 2, 3}, {2, 4, 7}, List]
InterpolatingPolynomial[pts, x]
polinterpolate([1, 2, 3], [2, 4, 7])
characteristic polynomial CharacteristicPolynomial[{{1, 2}, {3, 4}}, x] charpoly([1, 2; 3, 4])
minimal polynomial
piecewise polynomial
rational function (x - 1) / (x - 2)^2 (x - 1) / (x - 2)^2
power series
mathematica sympy gap pari/gp
power series of differentiable function Series[Cos[x], {x, 0, 10}] sympy.series(sympy.cos(x), x, n=11) \ps 10
Ser(cos(x))
power series by formula Plus @@ Table[x^i / i!, {i, 0, 10}] + O[x]^11
special functions
mathematica sympy gap pari/gp
gamma Gamma[1/2] gamma(1/2)
hyperpolic Sinh Cosh Tanh
ArcSinh ArcCosh ArcTanh
elliptic integerals EllipticK EllipticF
EllipticE
EllipticPi
bessel functions BesselJ
BesselY
BesselI
BesselK
Riemann zeta Zeta[2] zeta(2)
permutations
mathematica sympy gap pari/gp
permutation from disjoint cycles p = Cycles[{{1, 2}, {3, 4}}] import sympy.combinatorics as combinatorics

p = combinatorics.Permutation(0, 1)(2, 3)
p := (1, 2)(3, 4);
permutation from list p = PermutationCycles[{2, 1, 4, 3}] import sympy.combinatorics as combinatorics

p = combinatorics.Permutation([1, 0, 3, 2])
p2 := PermList([2, 1, 4, 3]);
permutation from two lists FindPermutation[{a, b, c}, {b, c, a}] # must be positive integers:
p := MappingPermListList([6, 8, 4, 2], [2, 4, 6, 8])
act on element p = Cycles[{{1, 2}, {3, 4}}]

PermutationReplace[1, p]
p(0) 1 ^ p;

# preimage of 1 under p:
1 / p;
act on list (* if list is too big, extra elements retain
their positions; if list is too small,
expression is left unevaluated. *)

Permute[{a, b, c, d}, p12n34]
a, b, c, d = sympy.symbols('a b c d')

p = combinatorics.Permutation(0, 1)(2, 3)
p([a, b, c, d])
compose p1 = Cycles[{{1, 2}, {3, 4}}]
p2 = Cycles[{{1, 3}}]
PermutationProduct[p1, p2]
p1 = combinatorics.Permutation(0, 1)(2, 3)
p2 = combinatorics.Permutation(0, 2)

p1 * p2
(1, 2)(3, 4) * (1, 3);
invert InversePermutation[Cycles[{{1, 2, 3}}]] p = combinatorics.Permutation(0, 1, 2)

p ** -1
(1, 2, 3) ^ -1;
power PermutationPower[Cycles[{{1, 2, 3, 4, 5}}], 3] p = combinatorics.Permutation(0, 1, 2, 3, 4)

p ** 3
(1, 2, 3, 4, 5) ^ 3;
order PermutationOrder[Cycles[{{1, 2, 3, 4, 5}}]] combinatorics.Permutation(0, 1, 2, 3, 4).order()
support PermutationSupport[Cycles[{{1, 3, 5}, {7, 8}}]] p = combinatorics.Permutation(0, 2, 4)(6, 7)

p.support()
MovedPoints((1, 3, 5)(7, 8));
number of inversions Permutation(0, 2, 1).inversions()
parity Permutation(0, 2, 1).parity()
to inversion vector Permutation(0, 2, 1).inversion_vector()
from inversion vector Permutation.from_inversion_vector([2, 0])
all permutations GroupElements[SymmetricGroup[4]]

(* of a list: *)
Permutations[{a, b, c, d}]
random permutation RandomPermutation[10] Permutation.random(10)
groups
mathematica sympy gap pari/gp
group from permutation generators e1 = Cycles[{{1, 3, 5, 2}}]
e2 = Cycles[{{1, 2}}]
g := PermutationGroup[{e1, e2}]
from sympy.combinatorics import *

p1 = Permutation(0, 2, 4, 1)
p2 = Permutation(0, 1)
g = PermutationGroup(p1, p2)
g := Group((1, 3, 5, 2), (1, 2));

# or
g := GroupWithGenerators([(1, 3, 5, 2), (1, 2)]);
named groups

symmetric, alternating, cyclic, dihedral
s4 = SymmetricGroup[4]
a4 = AlternatingGroup[4]
z5 = CyclicGroup[5]
d10 = DihedralGroup[10]
from sympy.combinatorics import *

s4 = SymmetricGroup(4)
a4 = AlternatingGroup(4)
z5 = CyclicGroup(5)
d10 = DihedralGroup(10)
s4 := SymmetricGroup(4);
a4 := AlternatingGroup(4);
z5 := CyclicGroup(5);
d10 := DihedralGroup(2 * 10);
groups by size AllSmallGroups(8);
conjugate group ConjugateGroup(SymmetricGroup(4), (4, 5));
direct product from sympy.combinatorics import *

z3 = CyclicGroup(3)
a4 = AlternatingGroup(4)
g = DirectProduct(z3, a4)
z3 := CyclicGroup(3);
a4 := AlternatingGroup(4);
g := DirectProduct(z3, a4);
free product f := FreeProduct(CyclicGroup(3), CyclicGroup(2));
free group # integers under addition:
z := FreeGroup(1);

# free group with 2 generators:
f := FreeGroup("a", "b");
all elements GroupElements[DihedralGroup[10]]
identity element Cycles[{}] Identity(g)
random element RandomSample[GroupElements[g], 1][[1]] g.random() Random(g)
group operation e1 := RandomSample[GroupElements[g], 1][[1]]
e2 := RandomSample[GroupElements[g], 1][[1]]
PermutationProduct[e1, e2]
e1 = g.random()
e2 = g.random()
e1 * e2
e1 := Random(g);
e2 := Random(g);
e1 * e2;
inverse element e1**-1 Inverse(e1);
# or:
e1^-1;
commutator # e2 ** -1 * e1 ** -1 * e2 * e1:
e1.commutator(e2)
# e1^-1 * e2^-1 * e1 * e2:
Comm(e1, e2);
generators g.generators s10 := SymmetricGroup(10);
# return generators in an array:
GeneratorsOfGroup(s10);

# notation for individual generators:
s10.1;
s10.2;
express element using generators s10 := SymmetricGroup(10);
Factorization(s10, (1,3,8,10,5,9,2,7));
number of elements by generator word length s6 := SymmetricGroup(6);
GrowthFunctionOfGroup(s6);
group from finite presentation f := FreeGroup( "a", "b" );
g := f / [ f.1^2, f.2^3, (f.1 * f.2)^5 ];
order of group element d10 := DihedralGroup(10);

Order(d10.1);
Order(d10.2);
order GroupOrder[g] g.order() Size(g)
cyclic test IsCyclic(AlternatingGroup(10));
abelian test g.is_abelian IsAbelian(CyclicGroup(10));
identify StructureDescription(g);
cosets RightCoset()
CanonicalRightCosetElement()
CosetDecomposition()
RightTraversal(G, U)
subgroups
mathematica sympy gap pari/gp
all subgroups AllSubgroups(SymmetricGroup(4));
subgroup lattice s4 := SymmetricGroup(4);
lat := LatticeSubgroups(s4);
DotFileLatticeSubgroups(lat, "lattice.dot");
# dot -Tpng < lattice.dot > lattice.png
maximal subgroups MaximalSubgroups(s4);
frattini subgroup FrattiniSubgroup(DihedralGroup(8));
subgroup from generators g := Group((1, 3, 5, 7), (2, 4));
h := Subgroup(g, [(2, 4)]);
normal subgroups NormalSubgroups(s4);
center g.center() g := DirectProduct(CyclicGroup(4), DihedralGroup(6));
Center(g);
centralizer g := SymmetricGroup(5);
h := Centralizer(g, (1, 3)(4, 5));
normalizer s4 := SymmetricGroup(4);
g := Group([(1,2)(3,4)]);
Normalizer(s4, g);
commutator subgroup g1 := Group((1,2,3),(1,2));
g2 := Group((2,3,4),(3,4));
CommutatorSubgroup(g1, g2);
subgroup test
subgroup index Index(g, h);
normal test IsNormal(g, h);
subnormal test IsSubnormal(g, h);
nonabelian simple groups # argument is list of orders:
AllSmallNonabelianSimpleGroups([1..10000]);
simple test IsSimple(SymmetricGroup(4));
solvable test g.is_solvable IsSolvable(SymmetricGroup(4));
derived series g.derived_series() DerivedSeriesOfGroup(SymmetricGroup(4));
characteristic test s4 := SymmetricGroup(4);
h := Subgroup(s4, [(1,4)(2,3), (1,3)(2,4), (2,4,3)]);
IsCharacteristicSubgroup(s4, h);
semidirect product
group homomorphisms
mathematica sympy gap pari/gp
all homomorphisms s4 := SymmetricGroup(4);
s3 := SymmetricGroup(3);
AllHomomorphisms(s3, s4);
all homomorphims classes AllHomomorphismClasses(s3, s4);
endomorphisms and automorphisms AllEndomorphisms(s4);
AllAutomorphisms(s4);
homomorphism from generator images hom := GroupHomomorphismByImages(s3, s4,
[(1,2,3), (1,2)],
[(2,3,4), (2,3)]);

# uses generators of s3:
hom := GroupHomomorphismByImages(s3, s4,
[(2,3,4), (2,3)]);
surjective test IsSurjective(hom);
injective test IsInjective(hom);
bijective test IsBijective(hom);
kernel Kernel(AllHomomorphisms(s3, s4)[1]);
image Image(AllHomomorphisms(s3, s4)[1]);
actions
mathematica sympy gap pari/gp
conjugate element # (1,2,3)^-1 * (1,2) * (1,2,3):
(1,2)^(1,2,3)
conjugate set s3 := SymmetricGroup(3);
s3^(3,4);
(3,4)^s3;
conjugacy class s4: SymmetricGroup(4);
AsList(ConjugacyClass(s4, (1,2,3)));
conjugacy classes ConjugacyClasses(SymmetricGroup(4));
stabilizer
orbit
transitive test
descriptive statistics
mathematica sympy gap pari/gp
first moment statistics vals = {1, 2, 3, 8, 12, 19}
X = NormalDistribution[0, 1]

Mean[vals]
Total[vals]
Mean[X]
second moment statistics Variance[X]
StandardDeviation[X]
second moment statistics for samples Variance[vals]
StandardDeviation[vals]
skewness Skewness[vals]
Skewness[X]
kurtosis Kurtosis[vals]
Kurtosis[X]
nth moment and nth central moment Moment[vals, 5]
CentralMoment[vals, 5]
Moment[X, 5]
CentralMoment[X, 5]

MomentGeneratingFunction[X, t]
cumulant Cumulant[vals, 1]
Cumulant[X, 1]

CumulantGeneratingFunction[X, t]
entropy Entropy[vals]
mode Commonest[{1, 2, 2, 2, 3, 3, 8, 12}]
quantile statistics Min[vals]
Median[vals]
Max[vals]
InterquartileRange[vals]
Quantile[vals, 9/10]
bivariate statistiscs
correlation, covariance
Correlation[{1, 2, 3}, {2, 4, 7}]
Covariance[{1, 2, 3}, {2, 4, 7}]
data set to frequency table data = {1, 2, 2, 2, 3, 3, 8, 12}
(* list of pairs: *)
tab = Tally[data]
(* dictionary: *)
dict = Counts[data]
frequency table to data set f = Function[a, Table[a[[1]], {i, 1, a[[2]]}]]
data = Flatten[Map[f, tab]]
bin data = {1.1, 3.7, 8.9, 1.2, 1.9, 4.1}
(* bins are [0, 3), [3, 6), and [6, 9): *)
bins = BinCounts[data, 0, 3, 6, 9]
distributions
mathematica sympy gap pari/gp
binomial

density, cumulative distribution, sample
X = BinomialDistribution[100, 1/2]

PDF[X][50]
CDF[X][50]
RandomVariate[X]
from sympy.stats import *

X = Binomial('X', 100, sympy.Rational(1, 2))

density(Y).dict[sympy.Integer(50)]
P(X < 50)
sample(X)
poisson X = PoissonDistribution[1] # P(X < 4) raises NotImplementedError:
X = Poisson('X', 1)
discrete uniform X = DiscreteUniformDistribution[{0, 99}] X = DiscreteUniform('X', list(range(0, 100)))
normal

density, cumulative distribution, quantile, sample
X = NormalDistribution[0, 1]

PDF[X][0]
CDF[X][0]
InverseFunction[CDF[X]][1/2]
RandomVariate[X]
from sympy.stats import *

X = Normal('X', 0, 1)

density(X)(0)
P(X < 0)
??
sample(X)
gamma X = GammaDistribution[1, 1] X = Gamma('X', 1, 1)
exponential X = ExponentialDistribution[1] X = Exponential('X', 1)
chi-squared X = ChiSquareDistribution[2] X = ChiSquared('X', 2)
beta X = BetaDistribution[10, 90] X = Beta('X', 10, 90)
uniform X = UniformDistribution[{0, 1}] X = Uniform('X', 0, 1)
student's t X = StudentTDistribution[2] X = StudentT('X', 2)
snedecor's F X = FRatioDistribution[1, 1] X = FDistribution('X', 1, 1)
empirical density function X = NormalDistribution[0, 1]
data = Table[RandomVariate[X], {i, 1, 30}]
Y = EmpiricalDistribution[data]
PDF[Y]
empirical cumulative distribution X = NormalDistribution[0, 1]
data = Table[RandomVariate[X], {i, 1, 30}]
Y = EmpiricalDistribution[data]
Plot[CDF[Y][x], {x, -4, 4}]
univariate charts
mathematica sympy gap pari/gp
vertical bar chart BarChart[{7, 3, 8, 5, 5},
ChartLegends->
{"a","b","c","d","e"}]

horizontal bar chart
BarChart[{7, 3, 8, 5, 5}, BarOrigin -> Left]
pie chart PieChart[{7, 3, 8, 5, 5}]

stem-and-leaf plot
Needs["StatisticalPlots"]
nd = NormalDistribution[0, 1]
n100 = RandomVariate[nd, 100]
StemLeafPlot[20 * n100]
histogram nd = NormalDistribution[0, 1]
Histogram[RandomReal[nd, 100], 10]
box-and-whisker plot nd = NormalDistribution[0, 1]
n100 = RandomVariate[nd, 100]
BoxWhiskerChart[d]

ed = ExponentialDistribution[1]
e100 = RandomVariate[ed, 100]
u100 = RandomReal[1, 100]
d = {n100, e100, u100}
BoxWhiskerChart[d]
set chart title BoxWhiskerChart[data,
PlotLabel -> "chart example"]
chart options PlotLabel -> "an example"

AxisLabel -> {"time", "distance"}
bivariate charts
mathematica sympy gap pari/gp

stacked bar chart
d = {{7, 1}, {3, 2}, {8, 1},
{5, 3}, {5, 1}}
BarChart[d, ChartLayout ->
"Stacked"]
scatter plot nd = NormalDistribution[0, 1]
rn = Function[RandomReal[nd]]
d = {rn[],rn[]} & /@ Range[1,50]
ListPlot[d]
linear regression line d = Table[{i,
2 i + RandomReal[{-5, 5}]},
{i, 0, 20}]
model = LinearModelFit[d, x, x]
Show[ListPlot[d],
Plot[model["BestFit"],
{x, 0, 20}]]
polygonal line plot f = Function[i, {i, rn[]}]
d = f /@ Range[1, 20]
ListLinePlot[d]
area chart d = {{7, 1, 3, 2, 8},
{1, 5, 3, 5, 1}}
sd = {d[[1]], d[[1]] + d[[2]]}
ListLinePlot[sd, Filling ->
{1 -> {Axis, LightBlue},
2 -> {{1}, LightRed}}]
cubic spline d = Table[{i, RandomReal[nd]},
{i, 0, 20}]
f = Interpolation[d,
InterpolationOrder -> 3]
Show[ListPlot[d],
Plot[f[x], {x, 0, 20}]]
function plot Plot[Sin[x], {x, -4, 4}]
quantile-quantile plot nd = NormalDistribution[0, 1]
d1 = RandomReal[1, 50]
d2 = RandomReal[nd, 50]
QuantilePlot[d1, d2]
axis label d = Table[{i, i^2}, {i, 1, 20}]
ListLinePlot[d,
AxesLabel -> {x, x^2}]
axis limits Plot[x^2, {x, 0, 20},
PlotRange -> {{0, 20}, {-200, 500}}]
logarithmic y-axis LogPlot[{x^2, x^3, x^4, x^5},
{x, 0, 20}]
trivariate charts
mathematica sympy gap pari/gp
3d scatter plot nd = NormalDistribution[0,1]
d = RandomReal[nd, {50, 3}]
ListPointPlot3D[d]
additional data set nd = NormalDistribution[0, 1]
x1 = RandomReal[nd, 20]
x2 = RandomReal[nd, 20]
ListLinePlot[{x1, x2}]
bubble chart nd = NormalDistribution[0,1]
d = RandomReal[nd, {50, 3}]
BubbleChart[d]
surface plot Plot3D[Sinc[Sqrt[x^2 + y^2]],
{x, -25, 25},
{y, -25, 25}]
_______________________________________________________ _______________________________________________________ _______________________________________________________ _______________________________________________________

## version used

The version of software used to check the examples in the reference sheet.

## show version

How to determine the version of an installation.

## implicit prologue

Code assumed to have been executed by the examples in the sheet.

sympy:

# Grammar and Invocation

## interpreter

How to execute a script.

## repl

How to launch a command line read-eval-print loop for the language.

mathematica:

The full path to MathKernel on Mac OS X:

/Applications/Mathematica.app/Contents/MacOS/MathKernel


## block delimiters

How blocks are delimited.

## statement separator

How statements are separated.

## end-of-line comment

Character used to start a comment that goes to the end of the line.

# Variables and Expressions

## assignment

How to perform assignment.

In all three languages an assignment is an expression that evaluates to the right side of the expression. Assignments can be chained to assign the same value to multiple variables.

mathematica:

The Set function behaves identically to assignment and can be nested:

Set[a, Set[b, 3]]


## delayed assignment

How to assign an expression to a variable name. The expression is re-evaluated each time the variable is used.

mathematica:

GNU make also supports assignment and delayed assignment, but = is used for delayed assignment and := is used for immediate assignment. This is the opposite of how Mathematica uses the symbols.

The POSIX standard for make only has = for delayed assignment.

## parallel assignment

How to assign values in parallel.

Parallel assignment can be used to swap the values held in two variables.

## compound assignment

The compound assignment operators.

## global variable

How to declare a global variable.

## local variable

How to declare a local variable.

pari/gp:

There is my for declaring a local variable with lexical scope and local for declaring a variable with dynamic scope.

local can be used to change the value of a global as seen by any functions which are called while the local scope is in effect.

## null test

How to test if a value is null.

## remove variable binding

How to remove a variable. Subsequent references to the variable will be treated as if the variable were undefined.

## conditional expression

A conditional expression.

# Arithmetic and Logic

## true and false

The boolean literals.

## falsehoods

Values which evaluate to false in a conditional test.

sympy:

Note that the logical operators Not, And and Or do not treat empty collections or None as false. This is different from the Python logical operators not, and, and or.

## logical operators

The Boolean operators.

sympy:

In Python, &, |, and & are bit operators. SymPy has defined __and__, __or__, and __invert__ methods to make them Boolean operators for symbols, however.

## relational operators

The relational operators.

sympy:

The full SymPy names for the relational operators are:

sympy.Equality             # ==
sympy.Unequality           # !=
sympy.GreaterThan          # >=
sympy.LessThan             # <=
sympy.StrictGreaterThan    # >
sympy.StrictLessThan       # <


The SymPy functions are attatched to the relational operators ==, !=, for symbols … using the methods __eq__, __ne__, __ge__, __le__, __gt__, __lt__. The behavior they provide is similar to the default Python behavior, but when one of the arguments is a SymPy expression, a simplification will be attempted before the comparison is made.

## arithmetic operators

The arithmetic operators.

## integer division

How to compute the quotient of two integers.

## integer division by zero

The result of dividing an integer by zero.

## float division

How to perform float division, even if the arguments are integers.

## float division by zero

The result of dividing a float by zero.

## power

How to compute exponentiation.

Note that zero to a negative power is equivalent to division by zero, and negative numbers to a fractional power may have multiple complex solutions.

## sqrt

The square root function.

For positive arguments the positive square root is returned.

## sqrt -1

How the square root function handles negative arguments.

mathematica:

An uppercase I is used to enter the imaginary unit, but Mathematica displays it as a lowercase i.

## transcendental functions

The standard transcendental functions such as one might find on a scientific calculator.

The functions are the exponential (not to be confused with exponentiation), natural logarithm, sine, cosine, tangent, arcsine, arccosine, arctangent, and the two argument arctangent.

## transcendental constants

The transcendental constants pi and e.

The transcendental functions can used to computed to compute the transcendental constants:

pi = acos(-1)
pi = 4 * atan(1)
e = exp(1)


## float truncation

Ways to convert a float to a nearby integer.

## absolute value

How to get the absolute value and signum of a number.

## integer overflow

What happens when the value of an integer expression cannot be stored in an integer.

The languages in this sheet all support arbitrary length integers so the situation does not happen.

## float overflow

What happens when the value of a floating point expression cannot be stored in a float.

## rational construction

How to construct a rational number.

## rational decomposition

How to extract the numerator and denominator from a rational number.

## decimal approximation

How to get a decimal approximation of an irrational number or repeating decimal rational.

## complex construction

How to construct a complex number.

## complex decomposition

How to extract the real and imaginary part from a complex number; how to extract the argument and modulus; how to get the complex conjugate.

## random number

How to generate a random integer or a random float.

## random seed

How to set or get the random seed.

mathematica:

The seed is not set to the same value at start up.

## binary, octal, and hex literals

Binary, octal, and hex integer literals.

mathematica:

The notation works for any base from 2 to 36.

# Strings

## convert from string, to string

How to convert strings to numbers and vice versa.

# Resizable Arrays

## literal

The notation for an array literal.

## size

The number of elements in the array.

## lookup

How to access an array element by its index.

## update

How to change the value stored at an array index.

## out-of-bounds behavior

What happens when an attempt is made to access an element at an out-of-bounds index.

## element index

How to get the index of an element in an array.

## slice

How to extract a subset of the elements. The indices for the elements must be contiguous.

## copy

How to copy an array. Updating the copy will not alter the original.

## membership

How to test whether a value is an element of a list.

## intersection

How to to find the intersection of two lists.

## union

How to find the union of two lists.

## relative complement, symmetric difference

How to find all elements in one list which are not in another; how to find all elements which are in one of two lists but not both.

## shuffle and sample

How to shuffle an array. How to extract a random sample from an array without replacement.

## zip

How to interleave two arrays.

# Execution Control

## if

How to write a branch statement.

mathematica:

The 3rd argument (the else clause) of an If expression is optional.

## while

How to write a conditional loop.

mathematica:

Do can be used for a finite unconditional loop:

Do[Print[foo], {10}]


## for

How to write a C-style for statement.

## break/continue

How to break out of a loop. How to jump to the next iteration of a loop.

# Exceptions

## raise exception

How to raise an exception.

## handle exception

How to handle an exception.

## uncaught exception behavior

gap:

Calling Error() invokes the GAP debugger, which is similar to a Lisp debugger. In particular, all the commands available in the GAP REPL are still available. Variables can be inspected and modified while in the debugger but any changes will be lost when the debugger is quitted.

One uses quit; or ^D to exit the debugger. These commands also cause the top-level GAP REPL exit if used while not in a debugger.

If Error() is invoked while in the GAP debugger, the debugger will be invoked recursively. One must use quit; for each level of debugger recursion to return to the top -level GAP REPL.

Use

brk> Where(4);


to print the top four functions on the stack when the error occurred. Use DownEnv() and UpEnv() to move down the stack—i.e. from callee to caller—and UpEnv() to move up the stack. The commands take the number of levels to move down or up:

brk> DownEnv(2);
brk> UpEnv(2);


When the debugger is invoked, it will print a message. It may give the user the option of providing a value with the return statement so that a computation can be continued:

brk> return 17;
`

## finally block

How to write code that executes even if an exception is raised.

# Reflection

## function documentation

How to get the documentation for a function.

# Vectors

## vector literal

The notation for a vector literal.

## vector coordinate

How to get one of the coordinates of a vector.

## vector dimension

How to get the number of coordinates of a vector.

## element-wise arithmetic operators

How to perform an element-wise arithmetic operation on vectors.

## vector length mismatch

What happens when an element-wise arithmetic operation is performed on vectors of different dimension.

## scalar multiplication

How to multiply a scalar with a vector.

## dot product

How to compute the dot product of two vectors.

## cross product

How to compute the cross product of two three-dimensional vectors.

## norms

How to compute the norm of a vector.

# Matrices

## literal or constructor

Literal syntax or constructor for creating a matrix.

mathematica:

Matrices are represented as lists of lists. No error is generated if one of the rows contains too many or two few elements. The MatrixQ predicate can be used to test whether a list of lists is matrix: i.e. all of the sublists contain numbers and are of the same length.

Matrices are displayed by Mathematica using list notation. To see a matrix as it would be displayed in mathematical notation, use the MatrixForm function.

## dimensions

How to get the dimensions of a matrix.

## element access

How to access an element of a matrix. All languages described here follow the convention from mathematics of specifying the row index before the column index.

## row access

How to access a row.

## column access

How to access a column.

## submatrix access

How to access a submatrix.

## scalar multiplication

How to multiply a matrix by a scalar.

## element-wise operators

Operators which act on two identically sized matrices element by element. Note that element-wise multiplication of two matrices is used less frequently in mathematics than matrix multiplication.

## multiplication

How to multiply matrices. Matrix multiplication should not be confused with element-wise multiplication of matrices. Matrix multiplication in non-commutative and only requires that the number of columns of the matrix on the left match the number of rows of the matrix. Element-wise multiplication, by contrast, is commutative and requires that the dimensions of the two matrices be equal.

## kronecker product

The Kronecker product is a non-commutative operation defined on any two matrices. If A is m x n and B is p x q, then the Kronecker product is a matrix with dimensions mp x nq.

## comparison

How to test two matrices for equality.

## norms

How to compute the 1-norm, the 2-norm, the infinity norm, and the frobenius norm.

# Elliptic Curves

The multidegree of a multivariate monomial is the sum of the exponents of the indeterminates. E.g. the multidegree of x²y⁵z³ is 10.

The multidegree of a multivariate polynomial is the maximum of the multidegrees of its terms.

A cubic multivariate polynomial is a multivariate polynomial of multidegree 3.

An elliptic curve can be represented by a cubic bivariate polynomial.

Geometrically, the zeros of a cubic bivariate polynomial form a curve in the coordinate plane. It is customary, however, to study elliptic curves in the projective plane. To get the projective plane, we take the set of triples (a, b, c), not all zero, and we define an equivalence relation on them where (a, b, c) = (ta, tb, tc) for any non-zero t. If c is non-zero, the (a/c, b/c, 1) is an element of this class, and by associating x with a/c and y with b/c this gives us a way of embedding the coordinate plane in the projective plane. Points in the projective plane with c = 0 can be thought of as points at infinity.

Two curves are birationally equivalent if there exist coordinate transformations going both directions which are rational functions and which map one curve to the other.

Every cubic bivariate polynomial is birationally equivalent to a polynomial in Weierstrass normal form:

(1)
$$y^2 = x^3 + a x + b$$

A consequence of Bezout's theorem is that any line will intersect a cubic bivariate polynomial in three places. The field of coefficients must be ℂ, and the coordinates must be projective to account for intersection at infinity. Also, one must allow for multiple intersections at a point.

The discriminant of an elliptic curve in Weierstrass normal form is:

(2)
$$-16(4a^3 + 27b^2)$$

An elliptic curve is said to be singular if the discriminant is non-zero. Geometrically, non-singular curves do not intersect themselves and don't have any cusps or isolated points.

A line will intersect a non-singular elliptic curve at three points in projective space. One can thus use the elliptic curve to define an abelian group in the following manner. Given points P and Q, define P * Q as the third point of intersection on the curve of the line defined by P and Q. Define P + Q as the third point of intersection of the line defined by P * Q and the point O at infinity. Note that + and not * is the group operator. To add a point P to itself, we use the tangent line of the elliptic curve at P. We define P * P as the third point of intersection of the tangent line with the curve, and P + P as the third point of intersection of the line defined by P * P and the point at infinity O.

# Permutations

A permutation is a bijection on a set of n elements.

The notation that Mathematica and GAP use assumes the set the permutation operates on is indexed by {1, .., n}. The notation that SymPy uses assumes the set is indexed by {0, …, n-1}.

Cayley two line notation

one line notation

cycle notation

inversions

# Groups

A group is a set G and a binary operator—written here as*—which takes elements of the set as operands and obeys the following axioms:

• Closure: For every g and h in G, g * h is also in G.
• Identity: There exists e in G such that for all g in G, e * g = g * e = g.
• Associativity: For all f, g and h in G, (f * g) * h = f * (g * h).
• Inverse: For every g in G, there exists g' in G such that g * g' = g' * g = e.

Abelian groups obey an additional axiom:

• Commutativity: For all g and h in G, g * h = h * g.

The order of a group is the number elements in the set. The smallest group is the trivial group, which contains a single element which is the identity. It has order 1.

The integers, rationals, real numbers, and complex numbers are Abelian groups under addition; zero is the identity element. The integers and rationals are countably infinite. The real numbers and complex numbers are uncountably infinite.

The integers modulo n are an Abelian group under addition; zero is the identity number. The group is finite and has order n.

The non-zero rationals, non-zero real numbers, and non-zero complex numbers are Abelian groups under multiplication; one is the identity element.

The set of permutations on a set of n elements are non-Abelian groups under composition; the identity permutation which maps each element of the set to itself is the identiy. The order of the group is n!.

## classical lie groups

The classical Lie groups provide examples of infinite, non-Abelian groups. In all cases the group operation is matrix multiplication:

 GL(n, ℝ) general linear group of degree n invertible n×n matrices SL(n, ℝ) special linear group of degree n n×n matrices with determinant one O(n, ℝ) orthogonal group of degree n n×n orthogonal matrices; i.e. MMT = I SO(n, ℝ) special orthogonal group of degree n n×n orthogonal matrices with determinant one U(n, ℂ) unitary group of degree n n×n unitary matrices; i.e. MM* = I SU(n, ℂ) special unitary group of degree n n×n unitary matrices with determinant one

## group from generators

gap:

When a group is created using GroupByGenerators, the generators returned by GeneratorsOfGroup will not necessarily be the same as the generators provided to the constructor.

If the group is created using GroupWithGenerators, then the generators returned by GeneratorsOfGroup will be the same.

# Subroups

A subgroup is a subset of a group which is itself a group.

testing whether a set is a subgroup

A nontrivial group always has at least two subgroups: the group itself and the trivial subgroup.

## subgroup from generators

gap:

ClosureGroup finds the smallest group containing a group and a list of elements, some of which might not be in the group which is the first argument.

# Group Homomorphisms

A homomorphism is a function φ from (G, *) to (H, *') such that

• φ(x * y) = φ(x) *' φ(y) for all x, y ∈ G.

An isomorphism is a bijective homomorphism.

If an isomorphism exists between two groups, they are said to be isomorphic. Isomorphic groups are in a sense the same, except that the elements and the operation are written differently. Isomorphism is an equivalence relation. Group theory is the study of properties which if they hold for one member of an isomorphism equivalence class, they hold for all members of the equivalence class.

# Actions

A group G is said to act on a set A if there is an operation ⋅: G × A → A such that

• g₁⋅(g₂⋅a) = (g₁*g₂)⋅a for all g₁, g₂ ∈ G and a ∈ A
• e⋅a = a for all a ∈ A where e is the identity in A

# Univariate Charts

A univariate chart can be used to display a list or array of numerical values. Univariate data can be displayed in a table with a single column or two columns if each numerical value is given a name. A multivariate chart, by contrast, is used to display a list or array of tuples of numerical values.

In order for a list of numerical values to be meaningfully displayed in a univariate chart, it must be meaningful to perform comparisons (<, >, =) on the values. Hence the values should have the same unit of measurement.

## vertical bar chart

A chart which represents values with rectangular bars which line up on the bottom. It is a deceptive practice for the bottom not to represent zero, even if a y-axis with labelled tick marks or grid lines is provided. A cut in the vertical axis and one of the bars may be desirable if the cut value is a large outlier. Putting such a cut all of the bars near the bottom is a deceptive practice similar not taking to the base of the bars to be zero, however.

Another bad practice is the 3D bar chart. In such a chart heights are represented by the height of what appear to be three dimensional blocks. Such charts impress an undiscriminating audience but make it more difficult to make a visual comparison of the charted quantities.

mathematica

## horizontal bar chart

A bar chart in which zero is the y-axis and the bars extend to the right.

## pie chart

A bar chart displays values using the areas of circular sectors or equivalently, the lengths of the arcs of those sectors. A pie chart implies that the values are percentages of a whole. The viewer is likely to make an assumption about what the whole circle represents. Thus, using a pie chart to show the revenue of some companies in a line of business could be regarded as deceptive if there are other companies in the same line of business which are left out. The viewer may mistakenly assume the whole circle represents the total market.

If two values are close in value, people cannot determine visually which of the corresponding sectors in a pie chart is larger without the aid of a protractor. For this reason many consider bar charts superior to pie charts.

Many software packages make 3D versions of pie charts which communicate no additional information and in fact make it harder to interpret the data.

# SymPy

Welcome to SymPy’s documentation!

# GAP

GAP - Reference Manual