Computer Algebra I: Mathematica, Maple, Maxima, Sage, SymPy

a side-by-side reference sheet

grammar and invocation | variables and expressions | arithmetic and logic | strings | arrays | sets | arithmetic sequences | dictionaries | functions | execution control | exceptions | streams | files | directories | libraries and namespaces | reflection

symbolic expressions | calculus | equations and unknowns | optimization | vectors | matrices | combinatorics | number theory | polynomials | trigonometry | special functions | permutations | descriptive statistics | distributions | statistical tests

bar charts | scatter plots | line charts | surface charts | chart options

mathematicamaplemaximasagesympy
version used
 
10.020165.376.10Python 2.7; SymPy 0.7
show version
 
select About Mathematica in Mathematica menuselect About Maple in Maple menu$ maxima --version$ sage --version

also displayed on worksheet
sympy.__version__
implicit prologue# unknowns other than x must be declared:
y = var('y')
from sympy import *

# enable LaTeX rendering in Jupyter notebook:
init_printing()

# unknown variables must be declared:
x, y = symbols('x y')
grammar and invocation
mathematicamaplemaximasagesympy
interpreter
 
$ cat > hello.m
Print["Hello, World!"]

$ MathKernel -script hello.m
$ cat >> hello.max
print("Hello, world!");

$ maxima -b hello.maxima
$ cat > hello.sage
print("Hello, World!")

$ sage hello.sage
if foo.py imports sympy:
$ python foo.py
repl
 
$ MathKernel$ maple$ maxima$ sage$ python
>>> from sympy import *
block delimiters
 
( stmt; )block([x: 3, y: 4], x + y);

/* Multiple stmts are separated by commas; a list of assignments can be used to set variables local to the block. */
: and offside rule: and offside rule
statement separator; or sometimes newline

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

The colon : suppresses output.

In the worksheet, typing RETURN causes a cell to be evaluated. The expression in the cell does not have to end with a semicolon. To enter a newline in a cell, use SHIFT + RETURN.

In the command line REPL, a statement is terminated by a semicolon.
; or $

The dollar sign $ suppresses output.
newline or ;

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

newlines not separators inside (), [], {}, triple quote literals, or after backslash: \
end-of-line comment
 
none1 + 1; # additionnone1 + 1 # addition1 + 1 # addition
multiple line comment
 
1 + (* addition *) 11 + (* addition *) 1;1 + /* addition */ 1;nonenone
variables and expressions
mathematicamaplemaximasagesympy
assignmenta = 3
Set[a, 3]

(* rhs evaluated each time a is accessed: *)
a := x + 3
SetDelayed[a, x + 3]
a := 3;

# rhs evaluated each time a is accessed:
a := 'x + 3';
a: 3;a = 3a = 3
parallel assignment{a, b} = {3, 4}
Set[{a, b}, {3, 4}]
a, b := 3, 4;[a, b]: [3, 4]a, b = 3, 4a, b = 3, 4
compound assignment+= -= *= /=
corresponding functions:
AddTo SubtractFrom TimeBy DivideBy
nonenone+= -= *= /= //= %= **=+= -= *= /= //= %= ^= **=
increment and decrement++x --x
PreIncrement[x] PreDecrement[x]
x++ x--
Increment[x] Decrement[x]
nonenonenonenone
non-referential identifierany unassigned identifier is non-referentialany unassigned identifier is non-referentialany unassigned identifier is non-referentialy, z, w = var('y z w')

# x is non-referential unless assigned a value
x, y, z, w = symbols('x y z w')
identifier as valuex = 3
y = HoldForm[x]
x: 3;
y: 'x;
global variablevariables are global by defaultvariables are global by defaultg1, g2 = 7, 8

def swap_globals():
  global g1, g2
  g1, g2 = g2, g1
g1, g2 = 7, 8

def swap_globals():
  global g1, g2
  g1, g2 = g2, g1
local variableModule[{x = 3, y = 4}, Print[x + y]]

(* makes x and y read-only: *)
With[{x = 3, y = 4}, Print[x + y]]

(* Block[ ] declares dynamic scope *)
block([x: 3, y: 4], print(x + y));assignments inside functions are to local variables by defaultassignments inside functions are to local variables by default
null
 
Nullno null valueNoneNone
null test
 
x == Nullno null valuex is Nonex is None
undefined variable access
 
treated as an unknown numbertreated as an unknown numbertreated as an unknown numberraises NameErrorraises NameError
remove variable bindingClear[x]
Remove[x]
x := 'x'kill(x);del xdel x
conditional expression
 
If[x > 0, x, -x]if x < 0 then -x else x end if;if x < 0 then -x else x;x if x > 0 else -xx if x > 0 else -x
arithmetic and logic
mathematicamaplemaximasagesympy
true and false
 
True Falsetrue falsetrue falseTrue FalseTrue False
falsehoods
 
FalseFalse None 0 0.0 '' [] {}False 0 0.0
logical operators! True || (True && False)
Or[Not[True], And[True, False]]
and or notand or notOr(Not(True), And(True, False))

# when arguments are symbols:
~ x | (y & z)
relational expression1 < 2evalb(1 < 2);is(1 < 2);
relational operators== != > < >= <=
corresponding functions:
Equal Unequal Greater Less GreaterEqual LessEqual
= != > < >= <== # > < >= <=== != > < >= <=Eq Ne Gt Lt Ge 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
+ - * / iquo mod

iquo and mod are functions
+ - * / quotitent() mod()

quotient and mod are functions, not binary infix operators.
+ - * / // %+ - * / ?? %

if an expression contains a symbol, then the above operators are rewritten using the following classes:
Add Mul Pow Mod
integer division
 
Quotient[a, b]iquo(7, 3);quotient(7, 3);7 // 3
integer division by zerodividend is zero:
Indeterminate
otherwise:
ComplexInfinity
error, numeric exceptionerrorraises ZeroDivisionError
float divisionexact division:
a / b
# exact division:
a / b;

# float division:
a * 1.0 / b;
a / b
float division by zerodividend is zero:
Indeterminate
otherwise:
ComplexInfinity
float(infinity)
float(undefined)
-float(infinity)
error
power2 ^ 32
Power[2, 32]
2 ^ 32

The worksheet rewrites the expression using a superscript. Do not put spaces around the operator.
2 ^ 32;
2 ** 32;
2 ^ 32
2 ** 32
2 ** 32
Pow(2, 32)
sqrtreturns symbolic expression:
Sqrt[2]
sqrt(2)sqrt(2);sqrt(2)sqrt(2)
sqrt -1
 
II%iII
transcendental functionsExp Log
Sin Cos Tan
ArcSin ArcCos ArcTan
ArcTan
ArcTan accepts 1 or 2 arguments
exp log
sin cos tan
arcsin arccos
arctan(y, x)
exp log
sin cos tan
asin acos atan
atan2
exp log
sin cos tan
asin acos atan
atan2
exp log
sin cos tan
asin acos atan
atan2
transcendental constants
π and Euler's number
Pi E EulerGammaPi exp(1) gamma%pi %e %gammapi e euler_gammapi E
float truncation
round towards zero, round to nearest integer, round down, round up
IntegerPart Round Floor Ceilingtrunc round floor ceiltruncate
round
floor
ceiling
int
round
floor
ceil
floor
ceiling
absolute value
and signum
Abs Signabs signabs sign

sign returns pos, neg, or zero
abs signAbs sign
integer overflow
 
none, has arbitrary length integer typenone, has arbitrary length integer typenone, has arbitrary length integer typenone, has arbitrary length integer typenone, has arbitrary length integer type
float overflow
 
nonenone
rational construction2 / 72 / 72 / 72 / 7Mul(2, Pow(7, -1))
Rational(2, 7)
rational decomposition
 
Numerator[2 / 7]
Denominator[2 / 7]
numer(2 / 7)
denom(2 / 7)
num(2 / 7);
denom(2 / 7);
numerator(2 / 7)
denominator(2 / 7)
numer, denom = fraction(Rational(2, 7))
decimal approximationN[2 / 7]
2 / 7 + 0.
2 / 7 // N
N[2 / 7, 100]
evalf(2 / 7)
2 / 7 + 0.0
evalf[100](2 / 7)

# use hardware floats:
evalhf(2 / 7);
2 / 7, numer;n(2 / 7)
n(2 / 7, 100)

# synonyms for n:
N(2 / 7)
numerical_approx(2 / 7)
N(Rational(2, 7))
N(Rational(2, 7), 100)
complex construction
 
1 + 3I1 + 3I1 + 3 * %i;1 + 3 * I1 + 3 * I
complex decomposition
real and imaginary part, argument and modulus, conjugate
Re Im
Arg Abs
Conjugate
Re Im
argument abs
conjugate
realpart imagpart
cabs carg
conjugate
(3 + I).real()
(3 + I).imag()
abs(3 + I)
arg(3 + I)
(3 + I).conjugate()
re im
Abs arg
conjugate
random number
uniform integer, uniform float
RandomInteger[{0, 99}]
RandomReal[]
rgen := rand(0..99);
rgen();

rgen := rand(0.0 .. 1.0);
rgen();
random(100);
random(1.0);
random seed
set, get
SeedRandom[17]
??
set_random_state(make_random_state(17));
??
bit operatorsBitAnd[5, 1]
BitOr[5, 1]
BitXor[5, 1]
BitNot[5]
BitShiftLeft[5, 1]
BitShiftRight[5, 1]
with(Bits):

And(5, 1)
Or(5, 1)
Xor(5, 1)
Not(1, bits=32)
Join([0, op(Split(5))])
Join(subsop(1 = NULL, Split(5)))
none
binary, octal, and hex literals2^^101010
8^^52
16^^2a
ibase: 2;
101010;

ibase: 8;
52;

/* If first hex digit is a letter, prefix a zero: */
ibase: 16;
2a;
radixBaseForm[42, 7]
BaseForm[7^^60, 10]
# [0, 6]:
convert(42, base, 7)
obase: 7;
42;
to array of digits(* base 10: *)
IntegerDigits[1234]
(* base 2: *)
IntegerDigits[1234, 2]
convert(1234, base, 10)
convert(1234, base, 2)
strings
mathematicamaplemaximasagesympy
string literal
 
"don't say \"no\"""don't say \"no\"""don't say \"no\""use Python stringsuse Python strings
newline in literal
 
yesYes. String literals separated only by spaces or linebreaks are merged into a single literal.Newlines are inserted into strings by continuing the string on the next line. However, if the last character on a line inside a string is a backslash, the backslash and the following newline are omitted.
literal escapes\\ \" \b \f \n \r \t \ooo\" \\
concatenate
 
"one " <> "two " <> "three""one "||"two "||"three";
cat("one ", "two ", "three");
concat("one ", "two ", "three");
translate caseToUpperCase["foo"]
ToLowerCase["FOO"]
with(StringTools):

UpperCase("foo");
LowerCase("FOO");
supcase("foo");
sdowncase("FOO");
trim
 
StringTrim[" foo "]with(StringTools):

Trim(" foo ");
strim(" ", " foo ");
number to string
 
"value: " <> ToString[8]cat("value: ", 8);concat("value: ", 8);
string to number7 + ToExpression["12"]
73.9 + ToExpression[".037"]
7 + parse("12");
73.9 + parse(".037");
7 + parse_string("12");
73.9 + parse_string(".037");

/* parse_string raises error if the string does
not contain valid Maxima code. Use numberp
predicate to verify that the return value is
numeric. */
string joinStringJoin[Riffle[{"foo", "bar", "baz"}, ","]]with(StringTools):

Join(["foo", "bar", "baz"], ",");
simplode(["foo", "bar", "baz"], ",");
split
 
StringSplit["foo,bar,baz", ","]with(StringTools):

Split("foo,bar,baz", ",");
split("foo,bar,baz", ",");
substitute

first occurrence, all occurences
s = "do re mi mi"
re = RegularExpression["mi"]

StringReplace[s, re -> "ma", 1]
StringReplace[s, re -> "ma"]
with(StringTools):

Substitute("do re mi mi mi", "mi", "ma");
SubstituteAll("do re mi mi mi", "mi", "ma");

# supports regexes:
RegSubs("mi"="ma", "do re mi mi mi");
ssubst("mi", "ma", "do re mi mi mi");
ssubstfirst("mi", "ma", "do re mi mi mi");
length
 
StringLength["hello"]length("hello");slength("hello");
index of substringStringPosition["hello", "el"][[1]][[1]]

(* The index of the first character is 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. *)
searchtext("el", "hello");

# Index of first character is 1.
# Returns 0 if not found.
ssearch("el", "hello");

/* 1 is index of first character;
returns false if substring not found */
_
extract substring(* "el": *)
StringTake["hello", {2, 3}]
substring("hello", 2..3);substring("hello", 2, 4);
character literal
 
nonenonenone
character lookupCharacters["hello"][[1]]# string of length 1:
"hello"[1]
chr and ordFromCharacterCode[{65}]
ToCharacterCode["A"][[1]]
with(StringTools):

Char(65);
Ord("A")
ascii(65);
cint("A");
delete charactersrules = {"a" -> "", "e" -> "", "i" -> "",
  "o" -> "", "u" -> ""}
StringReplace["disemvowel me", rules]
with(StringTools):

Remove(IsVowel, "disemvowel me");
arrays
mathematicamaplemaximasagesympy
literal{1, 2, 3}

List[1, 2, 3]
a := [1, 2, 3];[1, 2, 3];use Python listsuse Python lists
size
 
Length[{1, 2, 3}]numelems([1, 2, 3]);length([1, 2, 3]);
lookup(* access time is O(1) *)
(* indices start at one: *)
{1, 2, 3}[[1]]

Part[{1, 2, 3}, 1]
a := [6, 7, 8];
a[1];
a: [6, 7, 8];
a[1];
update
 
a[[1]] = 7a[1] := 7;a[1]: 7;
out-of-bounds behaviorleft as unevaluated Part[] expressionError for both lookup and update.Error for both lookup and update.
element index(* Position returns list of all positions: *)
First /@ Position[{7, 8, 9, 9}, 9]
with(ListTools):

Search(9, [7, 8, 9]);

# SearchAll returns all positions
a: [7, 8, 9, 9];
first(sublist_indices(a, lambda([x], x = 9)));
slice
 
{1, 2, 3}[[1 ;; 2]]
array of integers as index(* evaluates to {7, 9, 9} *)
{7, 8, 9}[[{1, 3, 3}]]
manipulate backa = {6,7,8}
AppendTo[a, 9]
elem = a[[Length[a]]]
a = Delete[a, Length[a]]
elem
manipulate fronta = {6, 7, 8}
PrependTo[a, 5]
elem = a[[1]]
a = Delete[a, 1]
elem
a: [6, 7, 8];
push(5, a);
elem: pop(a);
head
 
First[{1, 2, 3}]first([1, 2, 3]);
tail
 
Rest[{1, 2, 3}]rest([1, 2, 3]);
cons(* first arg must be an array *)
Prepend[{2, 3}, 1]
cons(1, [2, 3]);
concatenate
 
Join[{1, 2, 3}, {4, 5, 6}][op([1, 2, 3]), op([4, 5, 6])]append([1, 2, 3], [4, 5, 6]);
replicate
 
ten_zeros = Table[0, {i, 0, 9}]ten_zeros: makelist(0, 10);
copy
 
a2 = aa2: copylist(a);
iterate
 
Function[x, Print[x]] /@ {1, 2, 3}for i in [1, 2, 3] do
  print(i);
end do;
for i in [1, 2, 3] do print(i);
reverse
 
Reverse[{1, 2, 3}]with(ListTools):

Reverse([1, 2, 3]);
reverse([1, 2, 3]);
sort(* original list not modified: *)
a = Sort[{3, 1, 4, 2}]
# original list not modified:
a := sort([3, 1, 4, 2]);
sort([3, 1, 4, 2]);
dedupe
 
DeleteDuplicates[{1, 2, 2, 3}]with(ListTools):

# original list not modified:
a := MakeUnique([1, 1, 2, 3]);
unique([1, 2, 2, 3]);
membership
 
MemberQ[{1, 2, 3}, 2]evalb(2 in [1, 2, 3]);
member(2, [1, 2, 3]);
member(7, {1, 2, 3});
evalb(7 in {1, 2, 3});
mapMap[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}]
map(x -> x * x, [1, 2, 3]);map(lambda([x], x * x), [1, 2, 3]);
filter
 
Select[{1, 2, 3}, # > 2 &]select(x -> x > 2, [1, 2, 3]);sublist([1, 2, 3], lambda([x], x > 2));
reduce
 
Fold[Plus, 0, {1, 2, 3}]
universal and existential testsnone
min and max elementMin[{6, 7, 8}]
Max[{6, 7, 8}]
with(ListTools):

FindMinimalElement([6, 7, 8]);
FindMaximalElement([6, 7, 8]);
apply(min, [6, 7, 8]);
apply(max, [6, 7, 8]);
shuffle and samplex = {3, 7, 5, 12, 19, 8, 4}

RandomSample[x]
RandomSample[x, 3]
flatten
one level, completely
Flatten[{1, {2, {3, 4}}}, 1]
Flatten[{1, {2, {3, 4}}}]
with(ListTools):

FlattenOnce([1, [2, [3, 4]]]);
Flatten([1, [2, [3, 4]]]);
/* completely: */
flatten([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]
with(ListTools):

# list of six elements:
Interleave([1, 2, 3], ["a", "b", "c"]);

# list of lists with two elements:
zip(`[]`, [1, 2, 3], ["a", "b", "c"]);

# same as DotProduct([1, 2, 3], [2, 3, 4]):
add(zip(`*`, [1, 2, 3], [2, 3, 4]))
/* list of six elements: */
join([1, 2, 3], ["a", "b", "c"]);
cartesian productOuter[List, {1, 2, 3}, {"a", "b", "c"}]
sets
mathematicamaplemaximasagesympy
literal(* same as arrays: *)
{1, 2, 3}
{1, 2, 3}{1, 2, 3}{1, 2, 3}{1, 2, 3}
sizeLength[{1, 2, 3}]cardinality({1, 2, 3});len({1, 2, 3})len({1, 2, 3})
array to setDeleteDuplicates[{1, 2, 2, 3}]{op([1, 2, 3])};setify([1, 2, 3]);set([1, 2, 3])set([1, 2, 3])
set to arraynone; sets are arrays[op({1, 2, 3})];listify({1, 2, 3});list({1, 2, 3})list({1, 2, 3})
membership test
 
MemberQ[{1, 2, 3}, 7]evalb(7 in {1, 2, 3});
member(7, {1, 2, 3});
elementp(7, {1, 2, 3});7 in {1, 2, 3}7 in {1, 2, 3}
subset testSubsetQ[{1, 2, 3}, {1, 2}]subsetp({1, 2}, {1, 2, 3});{1, 2} <= {1, 2, 3}
{1, 2}.issubset({1, 2, 3})

{1, 2, 3} >= {1, 2}
{1, 2, 3}.issuperset({1, 2})
{1, 2} <= {1, 2, 3}
{1, 2}.issubset({1, 2, 3})

{1, 2, 3} >= {1, 2}
{1, 2, 3}.issuperset({1, 2})
universal and existential testsevery(lambda([x], x > 2), [1, 2, 3]);

some(lambda([x], x > 2), [1, 2, 3]);
unionUnion[{1, 2}, {2, 3, 4}]{1, 2, 3} union {2, 3, 4};union({1, 2, 3}, {2, 3, 4});{1, 2, 3} | {2, 3, 4}
{1, 2, 3}.union({2, 3, 4})
{1, 2, 3} | {2, 3, 4}
{1, 2, 3}.union({2, 3, 4})
intersection
 
Intersect[{1, 2}, {2, 3, 4}]{1, 2, 3} intersect {2, 3, 4};intersection({1, 2, 3}, {2, 3, 4});{1, 2, 3} & {2, 3, 4}
{1, 2, 3}.intersection({2, 3, 4})
{1, 2, 3} & {2, 3, 4}
{1, 2, 3}.intersection({2, 3, 4})
relative complementComplement[{1, 2, 3}, {2}]{1, 2, 3} minus {2, 3, 4};setdifference({1, 2, 3}, {2, 3, 4});{1, 2, 3} - {2, 3, 4}
{1, 2, 3}.difference({2, 3, 4})
{1, 2, 3} - {2, 3, 4}
{1, 2, 3}.difference({2, 3, 4})
powersetpowerset({1, 2, 3});set(Set({1, 2, 3}).subsets())
cartesian productOuter[List, {1, 2, 3}, {"a", "b", "c"}]cartesian_product({1, 2, 3}, {"a", "b", "c"});
arithmetic sequences
mathematicamaplemaximasagesympy
unit difference
 
Range[1, 100]i $ i = 1..100;makelist(i, i, 1, 100);range(1, 101)range(1, 101)
difference of 10
 
Range[1, 100, 10]10 * i + 1 $ i = 1..9;makelist(i, i, 1, 100, 10);range(1, 100, 10)range(1, 100, 10)
difference of 1/10
 
Range[1, 100, 1/10]1 + 1/10 * i $ i = 0..990makelist(i, i, 1, 100, 1/10);[1 + (1/10)*i for i in range(0, 991)][1 + Rational(1,10)*i for i in range(0, 991)]
dictionaries
mathematicamaplemaximasagesympy
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]
d := table(["t"=1, "f"=0]);d: [["t", 1], ["f", 0]];use Python dictionariesuse Python dictionaries
size
 
Length[Keys[d]]numelems(d);length(d);
lookup
 
d["t"]d["t"]assoc("t", d);
updated["f"] = -1d["f"] := -1;d2: cons(["f", -1],
  sublist(d, lambda([p], p[1] # "f")));
missing key behavior
 
Returns a symbolic expression with head "Missing". If the lookup key was "x", the expression is:

  Missing["KeyAbsent", "x"]
Returns unevaluated expression.assoc returns false
is key present
 
KeyExistsQ[d, "t"]
iterate
 
keys and values as arraysKeys[d]
Values[d]
with(ListTools):

FlattenOnce(indices(d));
FlattenOnce(entries(d));
map(lambda([p], p[1]), d);
map(lambda([p], p[2]), d);
sort by valuesSort[d]
functions
mathematicamaplemaximasagesympy
define functionAdd[a_, b_] := a + b

(* alternate syntax: *)
Add = Function[{a, b}, a + b]
add2 := (a, b) -> a + b

add2 := proc(a, b)
  a + b
end proc;
add(a, b) := a + b;

define(add(a, b), a + b);

/* block body: */
add(a, b) := block(print("adding", a, "and", b), a + b);

/* square bracket syntax: */
I[row, col] := if row = col then 1 else 0;
I[10, 10];
invoke functionAdd[3, 7]

Add @@ {3, 7}

(* syntax for unary functions: *)
2 // Log
add2(3, 7)add(3, 7);
boolean function attributes
list, set, clear
Attributes[add]
SetAttributes[add, {Orderless, Flat, Listable}]
ClearAtttibutes[add, Listable]
undefine functionClear[Add]add2 := 'add2';remfunction(add);
redefine functionAdd[a_, b_] := b + aadd2 := (a, b) -> b + a;add(a, b) := b + a;
missing function behaviorThe expression is left unevaluated. The head is the function name as a symbol, and the parts are the arguments.The expression is left unevaluated.The expression is left unevaluated.
missing argument behaviorThe expression is left unevaluated. The head is the function name as a symbol, and the parts are the arguments.Invalid input error.Too few arguments error.
extra argument behaviorThe expression is left unevaluated. The head is the function name as a symbol, and the parts are the arguments.Extra arguments are ignored.Too many arguments error.
default argumentOptions[myLog] = {base -> 10}
myLog[x_, OptionsPattern[]] :=
  N[Log[x]/Log[OptionValue[base]]]

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

(* override default: *)
myLog[100, base -> E]
return valuelast expression evaluated, or argument of Return[]last expression evaluated

Inside a block(), the last expression evaluated or the argument of return()
anonymous functionFunction[{a, b}, a + b]

(#1 + #2) &
f: lambda([x, y], x + y);

f(3, 7);
variable number of arguments(* one or more arguments: *)
add[a__] := Plus[a]

(* zero or more arguments: *)
add[a___] := Plus[a]
add([a]) := sum(a[i], i, 1, length(a));
pass array elements as separate argumentsApply[f, {a, b, c}]

f @@ {x, y, z}
add(a, b) := a + b;
apply(add, [3, 7]);
a = [x, y, z]
f(*a)
execution control
mathematicamaplemaximasagesympy
ifIf[x > 0,
  Print["positive"],
  If[x < 0,
    Print["negative"],
    Print["zero"]]]
if x > 0 then
  print("positive");
elif x < 0 then
  print("negative");
else
  print("zero");
end if
if x > 0
  then print("positive")
  else if x < 0
    then print("negative")
    else print("zero");
use Python execution controluse Python execution control
whilei = 0
While[i < 10, Print[i]; i++]
for i: 0 step 1 while i < 10 do print(i);
forFor[i = 0, i < 10, i++, Print[i]]for i from 1 by 1 while i < 10 do
  print(i);
end do;
for i: 1 step 1 thru 10 do print(i);
breakBreak[]
continueContinue[]
exceptions
mathematicamaplemaximasagesympy
raise exceptionThrow["failed"]error "failed";error("failed");use Python exceptionsuse Python exceptions
handle exceptionPrint[Catch[Throw["failed"]]]errcatch(error("failed"));
streams
mathematicamaplemaximasagesympy
standard file handlesStreams["stdout"]
Streams["stderr"]

(* all open file handles: *)
Streams[]
write line to stdoutPrint["hello"]
open file for readingf = OpenRead["/etc/hosts"]
open file for writingf = OpenWrite["/tmp/test"]
open file for appendingf = OpenAppend["/tmp/test"]
close fileClose[f]
read file into strings = ReadString[f]
write stringWriteString[f, "lorem ipsum"]
read file into array of stringss = Import["/etc/hosts"]
a = StringSplit[s, "\n"]
file handle position

get, set
f = StringToStream["foo bar baz"]

StreamPosition[f]

(* beginning of stream: *)
SetStreamPosition[f, 0]
(# end of stream: *)
SetStreamPosition[f, Infinity]
open temporary filef = OpenWrite[]
path = Part[f, 1]
files
mathematicamaplemaximasagesympy
file exists test, regular file testFileExistsQ["/etc/hosts"]
FileType["/etc/hosts"] == File
file sizeFileByteCount["/etc/hosts"]
is file readable, writable, executablewith(FileTools):

isReadable("/etc/hosts");
IsWritable("/etc/hosts");
IsExecutable("/etc/hosts");
last modification timeFileDate["/etc/hosts"]
copy file, remove file, rename fileCopyFile["/tmp/foo", "/tmp/bar"]
DeleteFile["/tmp/foo"]
RenameFile["/tmp/bar", "/tmp/foo"]
directories
mathematicamaplemaximasagesympy
working directorydir = Directory[]

SetDirectory["/tmp"]
dir := currentdir();

currentdir("/tmp");
build pathnameFileNameJoin[{"/etc", "hosts"}]with(FileTools):

JoinPath(["/etc", "hosts"]);
dirname and basenameDirectoryName["/etc/hosts"]
FileBaseName["/etc/hosts"]
with(FileTools):

ParentDirectory("/etc/hosts");
Filename("/etc/hosts");
absolute pathname(* file must exist;
   symbolic links are resolved: *)

AbsoluteFileName["foo"]
AbsoluteFileName["/foo"]
AbsoluteFileName["../foo"]
AbsoluteFileName["./foo"]
AbsoluteFileName["~/foo"]
glob pathsFunction[x, Print[x]] /@ FileNames["/tmp/*"]
make directoryCreateDirectory["/tmp/foo.d"]mkdir("/tmp/foo.d");
recursive copyCopyDirectory["/tmp/foo.d", "/tmp/baz.d"]
remove empty directoryDeleteDirectory["/tmp/foo.d"]rmdir("/tmp/foo.d");
remove directory and contentsDeleteDirectory["/tmp/foo.d",
  DeleteContents -> True]
directory testDirectoryQ["/etc"]isdir("/etc");
libraries and namespaces
mathematicamaplemaximasagesympy
load libraryGet["foo.m"]load(grobner);
reflection
mathematicamaplemaximasagesympy
get function documentation?Tan
Information[Tan]
Describe(solve);describe(solve);

? solve;
solve?print(solve.__doc__)

# in IJupyter:
solve?
help(solve)
function optionsOptions[Solve]
Options[Plot]
function sourceinterface(verboseproc = 2);
print(ScatterPlot);
import inspect

inspect.getsourcelines(integrate)
query data typeHead[x]whattype(x);symbolp(x);
numberp(7);
stringp("seven");
listp([1, 2, 3]);
type(x)
list variables in scopeNames[$Context <> "*"]/* user defined variables: */
values;

/* user defined functions: */
functions;
symbolic expressions
mathematicamaplemaximasagesympy
literalexpr = 1 + x + x^2expr := 1 + x + x^2;expr = 1 + x + x^2;expr = 1 + x + x^2x = symbols('x')

expr = 1 + x + x^2
prevent simplificationHoldForm[x + x]
x + x // HoldForm
variable updateexpr = 1 + x
x = 3
(* 4: *)
expr
expr := 1 + x;
x := 3;
# 4:
expr;
expr: 1 + x;
x: 3;
/* 1 + x: */
expr;
expr = 1 + x
x = 7
# 1 + x:
expr
x = symbols('x')
expr = 1 + x
x = 3
# 1 + x:
expr
substitute(* {3, 3}: *)
ReplaceAll[{x, x}, x -> 3]

(* {3, 3}: *)
{x, x} /. x -> 3

(* {3, 4}: *)
{x, y} /. {x -> 3, y -> 4}
# [3, 3]:
subs(x = 3, [x, x]);

# [3, 4]:
subs([x = 3, y = 4], [x, y])
/* [3, 3]: */
subst(3, x, [x, x]);
vector([x, x]).subs({x: 3})Matrix([x, x]).subs(x, 3)
piecewise-defined expressionPiecewise[{{x, x >= 0}, {-x, x < 0}}]

(* otherwise case: *)
Piecewise[{{-x, x < 0}}, x]
piecewise(x < 0, -x, x >= 0, x);

# otherwise case:
piecewise(x < 0, -x, x);
if x < 0 then -x else x;

/* integrating over piecewise-defined expression fails */
piecewise([
  ((-infinity,0), -x),
  ((0,infinity), x)])
Piecewise((-x, x < 0), (x, x >= 0))

# otherwise case:
Piecewise((-x, x < 0), (x, True))
simplifySimplify[Cos[x]^2 + Sin[x]^2]

(* perform more simplications: *)
FullSimplify[-(1/2) I E^(-I x) (-1 + E^(2 I x))]
simplify(cos(x)^2 + sin(x)^2);simplify(cos(x)**2 + sin(x)**2)
assumptionSimplify[Sqrt[x^2], Assumptions -> x >= 0]
Simplify[(-1)^(n * (n + 1)),
  Assumptions -> Element[n, Integers]]

(* perform fewer simplications: *)
Refine[Sqrt[x^2], Assumptions -> x >= 0]
Refine[(-1)^(n * (n + 1)), Element[n, Integers]]
assume(x > 0);
# Maple puts a tilde after an unknown with
# an assumption; e.g: x~:

sqrt(x^2);

# expression-local assumption:
sqrt(x^2) assuming(x > 0);

simplify((-1)^(n * (n + 1))) assuming(n, integer);
assume(x > 0);
sqrt(x^2);

/* There is no assumption predicate for
   integer variables. */
assume(x > 0)
sqrt(x^2)
x = symbols('x', positive=True)
sqrt(x ** 2)

n = symbols('n', integer=True)
(-1)**((n) * (n + 1))
assumption predicatesElement[x, Complexes]
Element[x, Reals]
Element[x, Algebraics]
Element[x, Rationals]
Element[x, Integers]
Element[x, Primes]
Element[x, Integers] && Mod[x, 5] == 0
Element[x, Booleans]

(* assumptions can use inequalities and logical operators: *)
x > 0 || x < 0
# a partial list:
complex
real
rational
integer
prime
odd
even
positive
nonnegative
negative
Assumptions can only be created using relational operators.assume(x, 'complex')
assume(x, 'real')
assume(x, 'rational')
assume(x, 'integer')
assume(x, 'odd')
assume(x, 'even')
# a partial list:
complex
real
algebraic
rational
integer
positive
nonpositive
negative
nonnegative
nonzero
prime
odd
even
list assumptionsNone. Assumptions are always local.getassumptions(x);facts(x);

# assumptions on all symbols:
facts();
assumptions()x.assumptions0
remove assumptionNone. Assumptions are always local.# removes all assumptions about x:
x := 'x';
forget(x > 0);forget(x > 0)

# rm all assumptions:
forget()
# removes all assumptions about x:
x = symbols('x')
calculus
mathematicamaplemaximasagesympy
limit
 
Limit[Sin[x]/x, x -> 0]limit(sin(x) / x, x = 0);limit(sin(x)/x, x, 0);limit(sin(x)/x, x=0)limit(sin(x)/x, x, 0)
limit at infinity
 
Limit[1/x, x -> Infinity]limit(1 / x, x = infinity);limit(1/x, x, inf);limit(1/x, x=infinity)limit(1/x, x, oo)
one-sided limit

from left, from right
Limit[1/x, x -> 0, Direction -> 1]
Limit[1/x, x -> 0, Direction -> -1]
limit( 1/ x, x = 0, left);
limit( 1/ x, x = 0, right);
limit(1/x, x, 0, minus);
limit(1/x, x, 0, plus);
limit(1/x, x=0, dir='-')
limit(1/x, x=0, dir='+')
limit(1/x, x, 0, '-')
limit(1/x, x, 0, '+')
derivativeD[x^3 + x + 3, x]

D[x^3 + x + 3, x] /. x -> 2
diff(x^3 + x + 3, x);diff(x^3 + x + 3, x);

at(diff(x^3 + x + 3, x), [x=2]);
diff(x^3 + x + 3, x)

diff(x^3 + x + 3, x).subs({x: 2})

# derivative is synonym of diff
diff(x**3 + x + 3, x)

diff(x**3 + x + 3, x).subs(x, 2)
derivative of a functionf[x_] = x^3 + x + 3

(* returns expression: *)
D[f[x, x]

(* return functions: *)
f'
Derivative[1][f]

(* evaluating derivative at a point: *)
f'[2]
Derivative[1][f][2]
f := (x) -> x^3 + x + 3;

# returns expression:
diff(f(x), x);

# returns function:
D(f);

# evaluating derivative at a point:
D(f)(2);
f(x) = x^3 + x + 3

diff(f)

diff(f)(2)
constants(* a depends on x; b does not: *)
D[a x + b, x, NonConstants -> {a}]

Dt[a x + b, x, Constants -> {b}]
/* symbols constant unless declared with depends: */
depends(a, x);
diff(a * x + b, x);

/* makes a constant again: */
remove(a, dependency);
higher order derivativeD[Log[x], {x, 3}]
Log'''[x]
Derivative[3][Log][x]
diff(log(x), [x$3]);diff(log(x), x, 3);diff(log(x), x, 3)diff(log(x), x, 3)
mixed partial derivativeD[x^9 * y^8, x, y, y]
D[x^9 * y^8, x, {y, 2}]
diff(x^9*y^8, x, y, y);
diff(x^9 * y^8, x, [y$2]);
diff(x^9 * y^8, x, 1, y, 2);diff(x^9 * y^8, x, 1).diff(y, 2)diff(x**9 * y**8, x, y, y)
div, grad, and curlDiv[{x^2, x * y, x * y * z}, {x, y, z}]

Grad[2 * x * y * z^2, {x, y, z}]

Curl[{x * y * z, y^2, 0}, {x, y, z}]
with(VectorCalculus):
SetCoordinates('cartesian'[x, y, z]);

Divergence(VectorField(<x^2, x*y, x*y*z>)));
Gradient(2*x*y*z^2, [x, y, z]);
Curl(VectorField(<x*y*z, y^2, 0>));
antiderivative
 
Integrate[x^3 + x + 3, x]integrate(x^3 + x + 3, x)integrate(x^3 + x + 3, x);integral(x^3 + x + 3, x)integrate(x**3 + x + 3, x)
definite integral
 
Integrate[x^3 + x + 3, {x, 0, 1}]integrate(x^3 + x + 3, [x = 0 .. 1]);integrate(x^3 + x + 3, x, 0, 1);integral(x^3 + x + 3, x, 0, 1)integrate(x**3 + x + 3, [x, 0, 1])
improper integral
 
Integrate[Exp[-x], {x, 0, Infinity}]integrate(exp(-x), [x = 0 .. infinity]);integrate(exp(-x), x, 0, inf);integral(exp(-x), x, 0, infinity)integrate(exp(-x), (x, 0, oo))
double integral(* integrates over y first: *)
Integrate[x^2 + y^2, {x, 0, 1}, {y, 0, x}]
# integrates over y first:
integrate(x^2 + y^2, [y = 0 .. x, x = 0 .. 1]);
integrate(
  integrate(x^2+y^2, y, 0, x), x, 0, 1);
integral(integral(x^2+y^2, y, 0, x), x, 0, 1)f = integrate(x**2 + y**2, (y, 0, x))
integrate(f, (x, 0, 1))
find poles
 
singular(1/(z-I));
residue
 
Residue[1/(z - I), {z, I}]residue(1/(z-I), z = I);residue(1 / (z - %i), z, %i);f(z) = 1/(z - I)
f.maxima_methods().residue(z, I)
residue(1/(z-I), z, I)
sum
 
Sum[2^i, {i, 1, 10}]sum(2^i, i = 1 .. 10);sum(2^i, i, 1, 10);sum(2^i for i in (1..10))Sum(2**i, (i, 1, 10)).doit()
series sum
 
Sum[2^-n, {n, 1, Infinity}]sum(2^(-n), n = 1 .. infinity);sum(2^-n, n, 1, inf), simpsum;sum(2^-n, n, 1, infinity)Sum(2**(-n), (n, 1, oo)).doit()
series expansion of functionSeries[Cos[x], {x, 0, 10}]series(cos(x), x = 0, 10);taylor(cos(x), [x, 0, 10]);taylor(cos(x), x, 0, 10)series(cos(x), x, n=11)
omitted order termexpr = 1 + x + x/2 + x^2/6 + O[x]^3

(* remove omitted order term: *)
Normal[expr]
product
 
Product[2*i + 1, {i, 0, 9}]product(2*i+1, i = 0 .. 9);product(2*i + 1, i, 0, 9);prod(2*i + 1 for i in (0..9))Product(2*i + 1, (i, 0, 9)).doit()
equations and unknowns
mathematicamaplemaximasagesympy
solve equation
 
Solve[x^3 + x + 3 == 0, x]solve(x^3 + x + 3 = 0, x);solve(x^3 + x + 3 = 0, x);solve(x^3 + x + 3 == 0, x)solve(x**3 + x + 3, x)
solve equationsSolve[x + y == 3 && x == 2 * y,
  {x, y}]

(* or: *)
Solve[{x + y == 3, x == 2 * y}, {x, y}]
solve({x = 2 * y, x + y = 3}, {x, y});solve([x + y = 3, x = 2*y], [x, y]);solve([x + y == 3, x == 2*y], x, y)solve([x + y - 3, 3*x - 2*y], [x, y])
differential equationDSolve[y'[x] == y[x], y[x], x]dsolve(diff(y(x), x) = y(x), y(x));desolve([diff(y(x), x) = y(x)], [y(x)]);y = function('y')(x)

desolve(diff(y, x) == y, y)
y = Function('y')

dsolve(Derivative(y(x), x) - y(x), y(x))
differential equation with boundary conditionDSolve[{y'[x] == y[x], y[0] == 1}, y[x], x]

DSolve[{y''[x] == y[x], y[0] == 1, y'[0] == 2},
  y[x], x]
dsolve({diff(y(x), x) = y(x), y(0) = 1}, y(x));

dsolve({diff(y(x), x, x) = y(x),
    y(0) = 1, (D(y))(0) = 2},
  y(x))
atvalue(y(x), x=0, 1);
desolve([diff(y(x), x) = y(x)], [y(x)]);
y = function('y')(x)

# y(0) = 1:
desolve(diff(y, x) == y, y, [0, 1])

# y(0) = 1 and y'(0) = 2:
desolve(diff(y, x, x) == y, y, [0, 1, 2])
support for boundary conditions is limited
differential equationseqn1 = x'[t] == x[t] - x[t] * y[t]
eqn2 = y'[t] == x[t] * y[t] - y[t]
DSolve[{eqn1, eqn2}, {x[t], y[t]}, t]
eqn1 := diff(x(t), t) = x(t)-x(t)*y(t);
eqn2 := diff(y(t), t) = x(t)*y(t)-y(t);
dsolve([eqn1, eqn2], [x(t), y(t)]);
eqn1: diff(x(t), t) = x(t) - x(t) * y(t);
eqn2: diff(y(t), t) = x(t) * y(t) - y(t);
desolve([eqn1, eqn2], [x(t), y(t)]);
recurrence equationeqns = {a[n + 2] == a[n + 1] + a[n],
  a[0] == 0,
  a[1] == 1}

RSolve[eqns, a[n], n]

(* remove Fibonacci[] from solution: *)
FunctionExpand[RSolve[eqns, a[n], n]]
eqns := {a(0) = 0,
  a(1) = 1,
  a(n+2) = a(n+1)+a(n)};

rsolve(eqns, a);
solve_rec(a[n]=a[n-1]+a[n-2], a[n], a[0] = 0, a[1] = 1);n = symbols('n')
a = Function('a')
eqn = a(n+2) - a(n+1) - a(n)

rsolve(eqn, a(n), {a(0): 0, a(1): 1})
optimization
mathematicamaplemaximasagesympy
minimize(* returns list of two items: min value and rule
   transforming x to argmin *)

Minimize[x^2 + 1, x]

(* 2 ways to get min value: *)
Minimize[x^2 + 1, x][[1]]
MinValue(x^2 + 1, x]

(* 2 ways to get argmin: *)
x /. Minimize[x^2 + 1, x][[2]]
ArgMin[x^2 + 1, x]
# 1:
minimize(x^2 + 1, x);

# 1, {[{x = 0}, 1]}:
minimize(x^2+1, x, location);
maximizeMaximize[-x^4 + 3 x^3, x]

Maxvalue[-x^4 + 3 x^3, x]
ArgMax[-x^4 + 3 x^3, x]
# 2187/256:
maximize(-x^4+3*x^3, x);

# 2187/256, {[{x = 9/4}, 2187/256]}:
maximize(-x^4+3*x^3, x, location);
objective with unknown parameter(* minval and argmin are expressions
   containing a: *)

Minimize[(x - a)^2 + x, x]
Doesn't work; returns expression unevaluated.
unbounded behavior(* MaxValue will be Infinity; MinValue will be
   -Infinity *)
# minimize will return -infinity;
# maximize will return infinity
multiple variables(* returns one solution: *)
Minimize[x^4 - 2 x^2 + 2 y^4 - 3 y^2, {x, y}]
expr := x^4+2*y^4-2*x^2-3*y^2;
# returns all four solutions:
minimize(expr, [x, y], location);
constraintsMinimize[{-x - 2 y^2, y^2 <= 17, 2 x + y <= 5},
  {x, y}]
none
infeasible behavior(* MaxValue will be -Infinity; MinValue will be
   Infinity; ArgMax or ArgMin will be
   Indeterminate *)
no infeasible expressions
integer variablesMaximize[{x^2 + 2*y,
    x >= 0, y >= 0,
    2 x + Pi * y <= 4},
  {x, y}, Integers]
vectors
mathematicamaplemaximasagesympy
vector literal(* row vector is same as array: *)
{1, 2, 3}
Vector(3, [1, 2, 3]);/* row vector is same as array: */
[1, 2, 3];
vector([1, 2, 3])# column vector:
Matrix([1, 2, 3])
constant vector

all zeros, all ones
Table[0, {i, 1, 100}]
Table[1, {i, 1, 100}]
Vector(100, fill=0);
Vector(100, fill=1);
makelist(0, 100);
makelist(1, 100);
vector([0] * 100)
vector([1] * 100)
Matrix([0] * 100)
Matrix([1] * 100)
vector coordinate(* indices start at one: *)
{6, 7, 8}[[1]]
Vector(3, [6, 7, 8])[1];[6, 7, 8][1];vector([6, 7, 8])[0]Matrix([6, 7, 8])[0]
vector dimension
 
Length[{1, 2, 3}]nops(Vector(3, [1, 2, 3]));length([1, 2, 3]);len(vector([1, 2, 3]))len(Matrix([6, 7, 8]))
Matrix([6, 7, 8]).shape[0]
element-wise arithmetic operators+ - * /
adjacent lists are multiplied element-wise
+ - *~ /~+ - * /+ -+ -

# element-wise multiplication:
A = Matrix([1, 2, 3])
B = Matrix([2, 3, 4])
A.multiply_elementwise(B)
vector length mismatch
 
errordimensions do not match errorerrorraises TypeErrorraises ShapeError
scalar multiplication3 {1, 2, 3}
{1, 2, 3} 3
* may also be used
3 * Vector(3, [1, 2, 3]);
Vector(3, [1, 2, 3]) * 3;
3 * [1, 2, 3];
[1, 2, 3] * 3;
3 * vector([1, 2, 3])
vector([1, 2, 3]) * 3
3 * Matrix([1, 2, 3])
Matrix([1, 2, 3]) * 3
dot product{1, 1, 1} . {2, 2, 2}
Dot[{1, 1, 1}, {2, 2, 2}]
Vector(3, [1, 1, 1]) . Vector(3, [2, 2, 2]);[1, 1, 1] . [2, 2, 2];vector([1, 1, 1]) * vector([2, 2, 2])
vector([1,1,1]).dot_product(vector([2,2,2]))
v1 = Matrix([1, 1, 1])
v2 = Matrix([2, 2, 2])
v1.dot(v2)
cross productCross[{1, 0, 0}, {0, 1, 0}]with(LinearAlgebra):

e1 := Vector(3, [1, 0, 0]);
e2 := Vector(3, [0, 1, 0]);
CrossProduct(e1, e2);
e1 &x e2;
e1 = vector([1, 0, 0])
e2 = vector([0, 1, 0])
e1.cross_product(e2)
e1 = Matrix([1, 0, 0])
e2 = Matrix([0, 1, 0])
e1.cross(e2)
normsNorm[{1, 2, 3}, 1]
Norm[{1, 2, 3}]
Norm[{1, 2, 3}, Infinity]
with(LinearAlgebra):

v := Vector(3, [1, 2, 3]);
VectorNorm(v, 1);
VectorNorm(v, 2);
VectorNorm(v, infinity);
vector([1, 2, 3]).norm(1)
vector([1, 2, 3]).norm()
vector([1, 2, 3]).norm(infinity)
vec = Matrix([1, 2, 3])

vec.norm(1)
vec.norm()
vec.norm(inf)
orthonormal basisOrthogonalize[{{1, 0, 1}, {1, 1, 1}}]load(eigen);

gramschmidt([[1, 0, 1], [1, 1, 1]]);
A = matrix([[1, 0, 1], [1, 1, 1]]

# Rows of B are orthogonal and span same
# space as rows of A. 2nd return value
# expresses rows of A as linear combos
# of rows of B.

B, _ = A.gram_schmidt()
matrices
mathematicamaplemaximasagesympy
literal or constructor(* used a nested array for each row: *)
{{1, 2}, {3, 4}}

(* display as grid with aligned columns: *)
MatrixForm[{{1, 2}, {3, 4}}]
A := <1, 2; 3, 4>;matrix([1, 2], [3, 4]);matrix([[1, 2], [3, 4]])Matrix([[1, 2], [3, 4]])
construct from sequenceArrayReshape[{1, 2, 3, 4, 5, 6}, {2, 3}]Matrix(2, 3, [1, 2, 3, 4, 5, 6])matrix([1, 2, 3, 4, 5, 6], nrows=2)Matrix(2, 3, [1, 2, 3, 4, 5, 6])
constant matrices

all zeros, all ones
Table[0, {i, 3}, {j, 3}]
Table[1, {i, 3}, {j, 3}]
Matrix(3, 3, fill=0);
Matrix(3, 3, fill=1);
zeromatrix(3, 3);

f[i, j] := 1;
genmatrix(f, 3, 3);
matrix([0] * 9, nrows=3)
matrix([1] * 9, nrows=3)
zeros(3, 3)
ones(3, 3)
diagonal matrices
and identity
DiagonalMatrix[{1, 2, 3}]
IdentityMatrix[3]
with(LinearAlgebra):

IdentityMatrix(3);
DiagonalMatrix([1, 2, 3]);
ident(3) * [1, 2, 3];
ident(3);
matrix.identity(3) * vector([1, 2, 3])
matrix.identity(3)
diag(*[1, 2, 3])
eye(3)
matrix by formulaTable[1/(i + j - 1), {i, 1, 3}, {j, 1, 3}]h := (i, j) -> 1/(i+j-1);
Matrix(3, 3, h);
h2[i, j] := 1/(i + j -1);
genmatrix(h2, 3, 3);
dimensions(* returns {3, 2}: *)
Dimensions[{{1, 2}, {3, 4}, {5, 6}}]
with(LinearAlgebra):

A := Matrix(3, 2, [1, 2, 3, 4, 5, 6]);
r, c := Dimension(A);
A: matrix([1, 2, 3], [4, 5, 6]);
matrix_size(A);
A = matrix([[1, 2], [3, 4], [5, 6]])
A.nrows()
A.ncols()
A = matrix([[1, 2], [3, 4], [5, 6]])

# returns (3, 2):
A.shape
element lookup(* top left corner: *)
{{1, 2}, {3, 4}}[[1, 1]]
<1, 2; 3, 4>[1][1];
<1, 2; 3, 4>[1, 1];
A: matrix([1, 2], [3, 4]);

A[1, 1];
A[1][1];
A = matrix([[1, 2], [3, 4]])
A[0, 0]
A[0][0]
A = Matrix([[1, 2], [3, 4]])

# top left corner:
A[0, 0]
extract row(* first row: *)
{{1, 2}, {3, 4}}[[1]]
<1, 2; 3, 4>[1];row(matrix([1, 2], [3, 4]), 1);
matrix([1, 2], [3, 4])[1];
# first row as vector:
A[0]
A.rows()[0]
# first row:
A[0, :]
extract column(* first column as array: *)
{{1, 2}, {3, 4}}[[All, 1]]
with(LinearAlgebra):

Column(<1, 2; 3, 4>, 1);
col(matrix([1, 2], [3, 4]), 1);# first column as vector:
A.columns()[0]
# first column as 1x2 matrix:
A[:, 0]
extract submatrixA = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}
A[[1;;2, 1;;2]]
A := <1, 2, 3; 4, 5, 6; 7, 8, 9>;
A[1..2, 1..2];
A = matrix(range(1, 10), nrows=3)

# takes two lists of indices:
A.matrix_from_rows_and_columns([0, 1], [0, 1])
rows = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
A = Matrix(rows)
A[0:2, 0:2]
scalar multiplication3 * {{1, 2}, {3, 4}}
{{1, 2}, {3, 4}} * 3
3 * <1, 2; 3, 4>;
<1, 2; 3, 4> * 3;
3 * matrix([1, 2], [3, 4]);
matrix([1, 2], [3, 4]) * 3;
3 * matrix([[1, 2], [3, 4]])
matrix([[1, 2], [3, 4]]) * 3
3 * Matrix([[1, 2], [3, 4]])
Matrix([[1, 2], [3, 4]]) * 3
element-wise operators+ - * /
adjacent matrices are multiplied element-wise
+ - *~+ - * /+ -+ -

# for Hadamard product:
A.multiply_elementwise(B)
productA = {{1, 2}, {3, 4}}
B = {{4, 3}, {2, 1}}
Dot[A, B]
(* or use period: *)
A . B
A := <1, 2; 3, 4>;
B := <4, 3; 2, 1>;
A . B;
A: matrix([1, 2], [3, 4]);
B: matrix([4, 3], [2, 1]);
A . B;
A = matrix([[1, 2], [3, 4]])
B = matrix([[4, 3], [2, 1]])
A * B
A = matrix([[1, 2], [3, 4]])
B = matrix([[4, 3], [2, 1]])
A * B
powerMatrixPower[{{1, 2}, {3, 4}}, 3]

(* element-wise operator: *)
A ^ 3
with(LinearAlgebra):

MatrixPower(<1, 2; 3, 4>, 3);
matrix([1, 2], [3, 4]) ^^ 3;A ^ 3
A ** 3
A ** 3
exponentialMatrixExp[{{1, 2}, {3, 4}}]with(LinearAlgebra):

MatrixExponential(<1, 2; 3, 4>);
exp(matrix([[1, 2], [3, 4]]))
logMatrixLog[{{1, 2}, {3, 4}}]
kronecker productA = {{1, 2}, {3, 4}}
B = {{4, 3}, {2, 1}}
KroneckerProduct[A, B]
with(LinearAlgebra):

A := <1, 2; 3, 4>;
B := <4, 3; 2, 1>;
KroneckerProduct(A, B);
A: matrix([1, 2], [3, 4]);
B: matrix([4, 3], [2, 1]);
kronecker_product(A, B);
A = matrix([[1, 2], [3, 4]])
B = matrix([[4, 3], [2, 1]])
A.tensor_product(B)
normsA = {{1, 2}, {3, 4}}

Norm[A, 1]
Norm[A, 2]
Norm[A, Infinity]
Norm[A, "Frobenius"]
with(LinearAlgebra):

A := <1, 2; 3, 4>;

Norm(A, 1);
Norm(A, 2);
Norm(A, infinity);
Norm(A, Frobenius);
A: matrix([1, 2], [3, 4]);

mat_norm(A, 1);
/* none */
mat_norm(A, inf);
mat_norm(A, frobenius);
A = matrix([[1, 2], [3, 4]])

# floating point values:
A.norm(1)
A.norm()
A.norm(infinity)
A.norm('frob')
transposeTranspose[{{1, 2}, {3, 4}}]

(* or ESC tr ESC for T exponent notation *)
with(LinearAlgebra):

Transpose(<1, 2; 3, 4>);
transpose(A);A.transpose()A.T
conjugate transposeA = {{1, I}, {2, -I}}
ConjugateTranspose[A]

(* or ESC ct ESC for dagger exponent notation *)
with(LinearAlgebra):

A := <1, I; 2, -I>;
HermitianTranspose(A);
ctranspose(matrix([1, %i], [2, -%i]));M = matrix([[1, I], [2, -I]])
M.conjugate_transpose()
M = Matrix([[1, I], [2, -I]])
M.adjoint()
inverseInverse[{{1, 2}, {3, 4}}]

(* expression left unevaluated: *)
Inverse[{{1, 0}, {0, 0}}]
with(LinearAlgebra):

MatrixInverse(<1, 2; 3, 4>);

# error:
MatrixInverse(<1, 0; 0, 0>);
invert(A);
A ^^ -1;

/* error: */
invert(matrix([1, 0], [0, 0]));
A.inverse()
A ^ -1
A ** -1
A.inv()

# raises ValueError:
Matrix([[1, 0], [0, 0]]).inv()
row echelon formRowReduce[{{1, 1}, {1, 1}}]with(LinearAlgebra):

ReducedRowEchelonForm(<1, 1; 1, 1>);
echelon(matrix([1, 1], [1, 1]));matrix([[1, 1], [1, 1]]).echelon_form()
pseudoinversePseudoInverse[{{1, 0}, {3, 0}}]
determinantDet[{{1, 2}, {3, 4}}]with(LinearAlgebra):

Determinant(<1, 2; 3, 4>);
determinant(A);A.determinant()A.det()
traceTr[{{1, 2}, {3, 4}}]with(LinearAlgebra):

Trace(<1, 2; 3, 4>);
load("nchrpl");

mattrace(matrix([1, 2], [3, 4]));
A.trace()
characteristic polynomialCharacteristicPolynomial[{{1, 2}, {3, 4}}, x]with(LinearAlgebra):

CharacteristicPolynomial(<1, 2; 3, 4>, x);
A: matrix([1, 2], [3, 4]);

charpoly(A, x);
matrix([[1, 2], [3, 4]]).charpoly('x')
minimal polynomialwith(LinearAlgebra):

MinimalPolynomial(IdentityMatrix(3), x);
load(diag);

minimalPoly(jordan(ident(3)));
matrix.identity(3).minpoly('x')
rankMatrixRank[{{1, 1}, {0, 0}}]with(LinearAlgebra):

Rank(<1, 1; 0, 0>);
rank(matrix([1, 1], [0, 0]));matrix([[1, 1], [0, 0]]).rank()
nullspace basisNullSpace[{{1, 1}, {0, 0}}]with(LinearAlgebra):

NullSpace(<1, 1; 0, 0>);
nullspace(matrix([1, 1], [0, 0]));
eigenvaluesEigenvalues[{{1, 2}, {3, 4}}]with(LinearAlgebra):

Eigenvalues(<1, 2; 3, 4>);
/* returns list of two lists:
   first is the eigenvalues,
   second is their multiplicities */

eigenvalues(A);
matrix([[1, 2], [3, 4]]).eigenvalues()A.eigenvals()
eigenvectorsEigenvectors[{{1, 2}, {3, 4}}]Eigenvectors(<1, 2; 3, 4>);/* returns list of two lists. The first item is the return value of eigenvalues(). The second item is a list containing a list of eigenvectors for each eigenvalue. */
eigenvectors(A);
A.eigenvects()
spectral decompositionA = {{1, 2}, {2, 1}}
z := Eigensystem[A]
d := DiagonalMatrix[z[[1]]]
P := Transpose[z[[2]]]

P . d . Inverse[P] == A
LUP decomposition{lu, p, c} = LUDecomposition[{{1, 2}, {3, 4}}]
L = LowerTriangularize[lu]
U = UpperTriangularize[lu]
P = Permute[IdentityMatrix[2], p]
L, U, P := LUDecomposition(<1, 2; 3, 4>);
QR decompositionA := {{1, 2}, {3, 4}}
{Q, R} = QRDecomposition[A]
A == Q . R
jordan decompositionA := 1, 2, 3}, {4, 5, 6}, {7, 8, 9
z := JordanDecomposition[A]
P := z[[1]]
J := z[[2]]
A . P == P . J
singular value decompositionA := {{1, 2}, {3, 4}}
z := SingularValueDecomposition[A]
U := z[[1]]
S := z[[2]]
V := z[[3]]

N[A] == N[U . S . ConjugateTranspose[V]]
polar decompositionA := {{1, 2}, {3, 4}}
{u, s, v} := SingularValueDecomposition[A]
vt := ConjugateTranspose[v]

U := u * vt
P = v * s * vt
combinatorics
mathematicamaplemaximasagesympy
factorial

and permutations
5!
Factorial[5]

Permutations[Range[1, 5]]
5!
factorial(5)
5!
factorial(5);
factorial(5)
5.factorial()
factorial(5)
binomial coefficient

and combinations
Binomial[10, 3]binomial(10, 3)binomial(10, 3);binomial(10, 3)binomial(10, 3)
multinomial coefficientMultinomial[3, 4, 5]with(combinat):

multinomial(12, 3, 4, 5);
multinomial(12, [3, 4, 5]);multinomial([3, 4, 5])
rising and falling factorialPochhammer[1/2, 3]

FactorialPower[1/2, 3]
pochhammer(1/2, 3);

pochhammer(1/2 - 2, 3);
pochhammer(1/2, 3);

pochhammer(1/2 - 2, 3);
rising_factorial(1/2, 3)

falling_factorial(1/2, 3)
subfactorial

and derangments
Needs["Combinatorica`"]

NumberOfDerangements[10]
subfactorial(10)subfactorial(10)
integer partitions(* number of partitions: *)
PartitionsP[10]

(* the partitions as an array: *)
IntegerPartitions[10]
with(combinat):

# number of partitions:
numbpart(10);

# the partitions as an array:
partition(10);
length(integer_partitions(10));

/* the partitions as an array: */
integer_partitions(10);

Partitions(10).cardinality()
Partitions(10).list()
from sympy.utilities.iterables \
  import partitions

len(list(partitions(10)))

[p.copy() for p in partitions(10)]
compositionsNeeds["Combinatorica`"]

(* weak compositions of size 3 is 66: *)
NumberOfCompositions[10, 3]

Compositions[10, 3]
with(combinat):

# compositions of size 3 is 36:
numbcomp(10, 3);

composition(10, 3);
# compositions of all lengths:
Compositions(10).cardinality()

Compositions(10).list()

# of length 3:
Compositions(10, min_length=3,
  max_length=3).list()
set partitionsStirlingS2[10, 3]

Needs["Combinatorica`"]

KSetPartitions[10, 3]
SetPartititions[10]
Stirling2(10, 3);stirling2(10, 3);stirling_number2(10, 3)
bell numberBellB[10]with(combinat):

bell(10);
belln(10);bell_number(10)bell(10)
permutations with k disjoint cyclesAbs[StirlingS1[10, 3]]abs(Stirling1(10, 3));abs(stirling1(10, 3));stirling_number1(10, 3)
fibonacci number

and lucas number
Fibonacci[10]
LucasL[10]
with(combinat):

fibonacci(10);
fib(10);
lucas(10);
fibonacci(10)
lucas_number2(10, 1, -1)
fibonacci(10)
lucas(10)
bernoulli numberBernoulliB[100]bernoulli(100);bern(100);bernoulli(100)bernoulli(100)
harmonic numberHarmonicNumber[100]sum(1/n, n = 1 .. 100);harmonic(100)
catalan numberCatalanNumber[10]catalan_number(10)catalan(10)
number theory
mathematicamaplemaximasagesympy
divisible test
 
Divisible[1001, 7]evalb(mod(1001, 7) = 0);is(mod(1001, 7) = 0);7.divides(1001)1001 % 7 == 0
divisors(* {1, 2, 4, 5, 10, 20, 25, 50, 100}: *)
Divisors[100]
with(NumberTheory):

Divisors(100);
divisors(100);divisors(100)ntheory.divisors(100)
pseudoprime testPrimeQ[7]isprime(7);primep(7);is_prime(7)
is_pseudoprime(7)
ntheory.primetest.isprime(7)
ntheory.primetest.mr(7, [2, 3])
prime factors(* returns {{2, 2}, {3, 1}, {7, 1}} *)
FactorInteger[84]
# [1, [[2, 2], [3, 1], [7, 1]]]:
ifactors(84);
/* 2^2 3 7: */
factor(84);

/* [[2,2],[3,1],[7,1]]: */
ifactors(84);
# 2^2 * 3 * 7:
factor(84)
# {2: 2, 3: 1, 7: 1}:
ntheory.factorint(84)
next prime

and preceding
NextPrime[1000]
NextPrime[1000, -1]
nextprime(1000);
prevprime(1000);
next_prime(1000);
prev_prime(1000);
next_prime(1000)
previous_prime(1000)
ntheory.generate.nextprime(1000)
ntheory.generate.prevprime(1000)
nth prime(* 541: *)
Prime[100]
ithprime(100);primes_first_n(100)[-1]ntheory.generate.prime(100)
prime counting function(* 25: *)
PrimePi[100]
with(NumberTheory):

pi(100);
prime_pi(100)ntheory.generate.primepi(100)
divmod
 
QuotientRemainder[7, 3]divide(7, 3);divmod(7, 3)divmod(7, 3)
greatest common divisor

and relatively prime test
GCD[14, 21]
GCD[14, 21, 777]

CoprimeQ[14, 21]
igcd(14, 21);
igcd(14, 21, 777);

with(NumberTheory):

AreCoprime(14, 21);
gcd(14, 21);
gcd(gcd(14, 21), 777);
gcd(14, 21)
gcd(gcd(14, 21), 777)
gcd(14, 21)
gcd(gcd(14, 21), 777)
extended euclidean algorithm(* {1, {2, -1}}: *)
ExtendedGCD[3, 5]
# Returns 1;
# sets a and b to 2 and -1:

igcdex(3, 5, 'a', 'b')
/* [2,-1,1]: */
gcdex(3, 5);
# (1, 2, -1):
xgcd(3, 5)
from sympy.ntheory.modular import igcdex

# (2, -1, 1):
igcdex(3, 5)
least common multipleLCM[14, 21]ilcm(14, 21);lcm(14, 21);lcm(14, 21)lcm(14, 21)
power modulusPowerMod[3, 212, 7]power_mod(3, 212, 7)
multiplicative inverse(* inverse of 2 mod 7: *)
PowerMod[2, -1, 7]

(* left unevaluated: *)
PowerMod[2, -1, 4]
r = Integers(7)
r(2)^-1

r2 = Integers(4)
# raises ZeroDivisionError:
r2(4)^-1
chinese remainder theorem(* returns 173, which is equal to 3 mod 17 and 8 mod 11: *)
ChineseRemainder[{3, 8}, {17, 11}]
# 173:
chrem([3, 8], [17, 11]);
/* 173: */
chinese([3, 8], [17, 11]);
crt(3, 8, 17, 11)
euler totient
 
EulerPhi[256]with(NumberTheory):

Totient(256);
totient(256);euler_phi(256)ntheory.totient(256)
carmichael functionCarmichaelLambda[561]with(NumberTheory):

CarmichaelLambda(561);
from sage.crypto.util import carmichael_lambda

carmichael_lambda(561)
multiplicative orderMultiplicativeOrder[7, 108]with(NumberTheory):

MultiplicativeOrder(7, 108);
Mod(7, 108).multiplicative_order()
primitive rootsPrimitiveRoot[11]

(* all primitive roots: *)
PrimitiveRootList[11]
with(NumberTheory):

PrimitiveRoot(11);
primitive_root(11)

# raises ValueError if none
discrete logarithm(* solves 10 = 2^x (mod 11): *)
MultiplicativeOrder[2, 11, 10]
with(NumberTheory):

ModularLog(10, 2, 11);
log(Mod(10, 11), Mod(2, 11))
quadratic residuesSelect[Range[0, 4], KroneckerSymbol[#, 5] == 1 &]quadratic_residues(5)
discrete square rootPowerMod[4, 1/2, 5]Mod(4, 5).sqrt()
kronecker symbol

and jacobi symbol
KroneckerSymbol[3, 5]
JacobiSymbol[3, 5]
with(NumberTheory):

KroneckerSymbol(3, 5);
JacobiSymbol(3, 5);
??
jacobi(3, 5);
kronecker_symbol(3, 5)
moebius functionMoebiusMu[11]with(NumberTheory):

Moebius(11);
moebius(11);moebius(11)
riemann zeta functionZeta[2]Zeta(2);zeta(2);zeta(2)mpmath.zeta(2)
continued fraction(* {0, 1, 1, 1, 5}: *)
ContinuedFraction[11/17]

(* arrray of first 100 digits for for pi: *)
ca= ContinuedFraction[Pi, 100]

(* rational approximation of pi: *)
FromContinuedFraction[a]
# [0, 1, 1, 1, 5]:
convert(11/17, confrac);

convert(Pi, confrac, 100);
/* [0,1,1,1,5]: */
cf(11/17);

float_pi: %pi, numer;
a = cf(float_pi);

/* as continued fraction: */
as_cf: cfdisrep(a);

/* as simple fraction: */
ratsimp(as_cf);
continued_fraction(11/17)

continued_fraction(pi, 100)
convergentsConvergents[11/17]

(* for continued fraction: *)
Convergents[{0, 1, 1, 1, 5}]

(* first 100 rational approximations: *)
Convergents[Pi, 100]
# no retval; convergents printed to stdout:
convergs(convert(11/17, confrac));

convergs(convert(Pi, confrac, 100));
# [0, 1, 1/2, 2/3, 11/17]:
continued_fraction(11/17).convergents()

# iterable infinite list:
continued_fraction(pi, 100).convergents()
polynomials
mathematicamaplemaximasagesympy
literalp = 2 -3 * x + 2* x^2
p2 = (1 + x)^10
p := x^2 - 3*x + 2;
p2 := (1 + x)^10;
p: x^2 - 3*x + 2;
p2: (x + 1)^10;
p = x^2 - 3*x + 2
p2 = (x + 1)^10
extract coefficientCoefficient[(1 + x)^10, x, 3]coeff((1 + x)^10, x, 3);coeff(expand((x + 1)^10), x^3);

coeff(expand((x + 1)^10), x, 3);
p = (1 + x)^10

# coefficients() returns (power, coeff) pairs:
[pair[0] for pair in p.coefficients()
 if pair[1] == 3][0]
extract coefficientsCoefficientList[(x + 1)^10, x]coeffs(collect((1 + x)^10, x), x);p: expand((x+1)^10);
makelist(coeff(p, x^i), i, 0, 10);
from array of coefficientsa = {2, -3, 1}
Sum[a[[i]] * x^i, {i, 1, 3}]
a := [2, -3, 1];
sum(a[i] * x^i, i = 1 .. 3);
a: [2, -3, 1];
sum(x^i * a[i + 1], i, 0, 2);
degreeExponent[(x + 1)^10, x]degree((1 + x)^10, x);hipow(expand((1 + x)^10), x);
expandExpand[(1 + x)^5]expand((1 + x)^5);expand((1 + x)^5);expand((1 + x)^5)expand((1 + x)**5)
factorFactor[3 + 10 x + 9 x^2 + 2 x^3]

Factor[x^10 - y^10]
factor(2*x^3 + 9*x^2 + 10*x + 3);

factor(x^10 - y^10);
factor(2*x^3 + 9*x^2 + 10*x + 3);factor(2*x^3 + 9*x^2 + 10*x + 3)factor(3 + 10*x + 9*x**2 + 2*x**3)
rootsSolve[x^3 + 3 x^2 + 2 x - 1 == 0, x]

(* just the 2nd root: *)
Root[x^3 + 3 x^2 + 2 x - 1, 2]
solve(x^3 + 3*x^2 + 2*x - 1 = 0, x);

# rational roots only:
roots(2*x^4 - 17*x^3 + 23*x^2 - 17*x + 21, x);
solve(x^3 + 3*x^2 + 2*x - 1 = 0, x);
quotient and remainderPolynomialReduce[x^10 - 1, x - 1, {x}]quo(x^10 - 1, x - 1, x);
rem(x^10 - 1, x - 1, x);
[q, r]: divide(x^10-1, x - 1);
greatest common divisorp1 = -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);
p1: -2 - x + 2 * x^2 + x^3;
p2: 6 - 7*x + x^3;
gcd(p1, p2);
extended euclidean algorithmp1 = -2 - x + 2 x^2 + x^3
p2 = 6 - 7 x + x^3

(* returns list; first element is GCD; 2nd element is list of two polynomials *)
PolynomialExtendedGCD[p1, p2, x]
p1 := x^3 + 2*x^2 - x - 2;
p2 := x^3 - 7*x + 6;

# returns gcd; sets a and b to polynomials:
gcdex(p1, p2, x, a, b);
resultantResultant[(x-1) * (x-2), (x-3) * (x-3), x]resultant((x-1)*(x-2), (x-3)*(x-3), x);resultant((x - 1)*(x - 2), (x - 3)*(x - 3), x);
discriminantDiscriminant[(x + 1) * (x - 2), x]discrim((x+1)*(x-2), x);
collect terms(* write as polynomial in x: *)
Collect[(1 + x + y)^3, x]
collect((1 + x + y)^3, x);load(facexp);

facsum(expand((1 + x + y)^5), x);
collect(expand((x+y+1)**3), x)
multivariate quotient and remainderPolynomialReduce[x^10 - y^10, x - y, {x, y}][q, r]: divide(x^10 - y^10, x - y);
groebner basisp1 = 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}]
with(Groebner):

p1 := x^2+y+z-1;
p2 := y^2+x+z-1;
p3 := z^2+x+y-1;

Basis([p1, p2, p3], plex(x, y, z));
specify orderingGroebnerBasis[{p1, p2, p3}, {x, y, z},
  MonomialOrder -> DegreeReverseLexicographic]

(* possible values for MonomialOrder:

   Lexicographic
   DegreeLexicographic
   EliminationOrder
   {1, 2, 3} *)
elementary symmetric polynomialSymmetricPolynomial[3, {x1, x2, x3, x4}]
symmetric reduction(* returns list of two elements; 2nd element is remainder if polynomial not symmetric: *)
SymmetricReduction[x^3 + y^3 + z^3, {x, y, z}]
# error if not symmetric:
convert(x^3 + y^3 + z^3, 'elsymfun');
cyclotomic polynomialCyclotomic[10, x]with(NumberTheory):

CyclotomicPolynomial(10, x);
hermite polynomialHermiteH[4, x]with(orthopoly):

H(4, x)
chebyshev polynomial

first and second kind
ChebyshevT[4, x]
ChebyshevU[4, x]
with(orthopoly):

T(4, x);
U(4, x);
interpolation polynomialpts = Inner[List, {1, 2, 3}, {2, 4, 7}, List]
InterpolatingPolynomial[pts, x]
spline
partial fraction decompositionApart[(3*x + 2)/ (x^2 + x)]

(* can handle multiple vars in denominator: *)
Apart[(b * c + a * d)/(b * d)]
convert((2 + 3*x)/(x^2 + x), 'parfrac');partfrac((3*x + 2) / (x^2 + x), x);apart((3*x+2) / (x*(x+1)))
add fractionsTogether[a/b + c/d]simplify(a/b + c/d);ratsimp(a/b + c/d);together(x/y + z/w)
pade approximantPadeApproximant[Log[x], {x, 1, {2, 3}}]with(numapprox):

pade(log(x), x = 1, [3, 2]);
p: taylor(log(x + 1), [x, 0, 5]);

pade(p, 3, 2);
trigonometry
mathematicamaplemaximasagesympy
eliminate powers and products of trigonometric functionsTrigReduce[Sin[x]^2 + Cos[x] Sin[x]]combine(sin(x)^2+cos(x)*sin(x));trigreduce(sin(x)^2 + cos(x) * sin(x));
eliminate sums and multiples inside trigonometric functionsTrigExpand[Sin[2 * x + 1]]expand(sin(2*x + 1));trigexpand(sin(2*x + 1));
trigonometric to exponentialTrigToExp[Cos[x]]convert(cos(x), exp);exponentialize(cos(x));cos(x).rewrite(cos, exp)
exponential to trigonometricExpToTrig[Exp[I x]]convert(exp(I*x), trig);demoivre(exp(%i * x));from sympy import exp, sin, I

exp(I * x).rewrite(exp, sin)
fourier expansion(* in sin and cos: *)
FourierTrigSeries[SquareWave[x / (2*Pi)],
  x, 10]

(* in complex exponentials: *)
FourierSeries[SquareWave[x / (2*Pi)], x, 10]
periodic functions on unit interval(*1: [0, 0.5); -1: [0.5, 1.0) *)
SquareWave[x]

(* 0 at 0 and 0.5; 1 at 0.25; -1 at 0.75 *)
TriangleWave[x]

(* x on [0, 1) *)
SawtoothWave[x]
fourier transformf[w_] = FourierTransform[ Sin[t], t, w]

InverseFourierTransform[f[w], w, t]
with(inttrans):

fourier(sin(t), t, w);
invfourier(f, w, t);
heaviside step functionHeaviside(1);   # 1
Heaviside(0);   # undefined
Heaviside(-1);  # 0

diff(Heaviside(x), x);  # Dirac(x)
dirac deltaDirac(1);   # 0
Dirac(0);   # Dirac(0)
Dirac(-1);  # 0

integrate(Dirac(x), x); # Heaviside(x)
integrate(Dirac(x), x=-1..1); # 1
special functions
mathematicamaplemaximasagesympy
gamma functionGamma[1/2]GAMMA(1/2);gamma(1/2);gamma(1/2)gamma(Rational(1, 2))
error functionErf[1/2] // N

Erfc Erfi
InverseErf InverseErfc
evalf(erf(1/2));

erfc erfi
erf(1/2), numer;

erfc erfi
n(erf(1/2))N(erf(Rational(1, 2)))

erfc erfi
hyperbolic functionsSinh Cosh Tanh
ArcSinh ArcCosh ArcTanh
sinh cosh tanh
arcsinh arccosh arctanh
sinh cosh tanh
asinh acosh atanh
sinh cosh tanh
asinh acosh atanh
sinh cosh tanh
asinh acosh atanh
elliptic integeralsEllipticK
EllipticF
EllipticE
EllipticPi
EllipticK
EllipticF
EllipticE
EllipticPi
elliptic_f
elliptic_e
elliptic_pi
bessel functionsBesselJ
BesselY
BesselI
BesselK
BesselJ
BesselY
BesselI
BesselK
bessel_j
bessel_y
bessel_i
bessel_k
permutations
mathematicamaplemaximasagesympy
from disjoint cyclesp = Cycles[{{1, 2}, {3, 4}}]Permutation([(1, 2), (3, 4)])import sympy.combinatorics as comb

p = combinatorics.Permutation(0, 1)(2, 3)
to disjoint cycles
from arrayp = PermutationCycles[{2, 1, 4, 3}]Permutation((2, 1, 4, 3))import sympy.combinatorics as comb

p = combinatorics.Permutation([1, 0, 3, 2])
from two arrays with same elementsFindPermutation[{a, b, c}, {b, c, a}]
size
support

and fixed points
PermutationSupport[Cycles[{{1, 3, 5}, {7, 8}}]]import sympy.combinatorics as comb

p = comb.Permutation(0, 2, 4)(6, 7)
p.support()
act on elementp = Cycles[{{1, 2}, {3, 4}}]

PermutationReplace[1, p]
Permutation((2, 1, 4, 3))(1)p(0)
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 = var('a b c d')

p = Permutation([(1, 2), (3, 4)])
p.action([a, b, c, d])
import sympy.combinatorics as comb

p = comb.Permutation(0, 1)(2, 3)
p([a, b, c, d])
composep1 = Cycles[{{1, 2}, {3, 4}}]
p2 = Cycles[{{1, 3}}]
PermutationProduct[p1, p2]
p1 = Permutation([(1, 2), (3, 4)])
p2 = Permutation((1, 3))
p1 * p2
import sympy.combinatorics as comb

p1 = comb.Permutation(0, 1)(2, 3)
p2 = comb.Permutation(0, 2)
p1 * p2
inverseInversePermutation[Cycles[{{1, 2, 3}}]]Permutation((1, 2, 3)).inverse()import sympy.combinatorics as comb

comb.Permutation(0, 1, 2) ** -1
powerPermutationPower[Cycles[{{1, 2, 3, 4, 5}}], 3]Permutation((1, 2, 3, 4, 5))^3import sympy.combinatorics as comb

comb.Permutation(0, 1, 2, 3, 4) ** 3
orderPermutationOrder[Cycles[{{1, 2, 3, 4, 5}}]]p = Permutation((1,2,3,4,5))
p.to_permutation_group_element().order()
import sympy.combinatorics as comb

comb.Permutation(0, 1, 2, 3, 4).order()
number of inversionsPermutation((1, 3, 2)).length()import sympy.combinatorics as comb

comb.Permutation(0, 2, 1).inversions()
parityPermutation((1, 3, 2)).is_even()import sympy.combinatorics as comb

comb.Permutation(0, 2, 1).parity()
to inversion vectorPermutation((1, 3, 2)).to_inversion_vector()import sympy.combinatorics as comb

comb.Permutation(0, 2, 1).inversion_vector()
from inversion vectorimport sympy.combinatorics as comb

comb.Permutation.from_inversion_vector([2, 0])
list permutationsGroupElements[SymmetricGroup[4]]

(* of a list: *)
Permutations[{a, b, c, d}]
list(SymmetricGroup(4))
random permutationRandomPermutation[10]Permutation.random(10)
descriptive statistics
mathematicamaplemaximasagesympy
first moment statisticsvals = {1, 2, 3, 8, 12, 19}
X = NormalDistribution[0, 1]

Mean[vals]
Total[vals]
Mean[X]
with(Statistics):

vals := [1, 2, 3, 8, 12, 19];
X := RandomVariable(Normal(0, 1));

Mean(vals);
Mean(X);
load(distrib);

/* Other distributions have similar functions: */
mean_normal(0, 1);
second moment statisticsVariance[X]
StandardDeviation[X]
with(Statistics):

Variance(X);
StandardDeviation(X);
load(distrib);

/* Other distributions have similar functions: */
var_normal(0, 1);
std_normal(0, 1);
second moment statistics for samplesVariance[vals]
StandardDeviation[vals]
with(Statistics):

Variance(vals);
StandardDeviation(vals);
skewnessSkewness[vals]
Skewness[X]
with(Statistics):

Skewness(vals);
Skewness(X);
load(distrib);

/* Other distributions have similar functions: */
skewness_normal(0, 1);
kurtosisKurtosis[vals]
Kurtosis[X]
with(Statistics):

Kurtosis(vals);
Kurtosis(X);
load(distrib);

/* Other distributions have similar functions: */
kurtosis_normal(0, 1);
nth moment and nth central momentMoment[vals, 5]
CentralMoment[vals, 5]
Moment[X, 5]
CentralMoment[X, 5]

MomentGeneratingFunction[X, t]
with(Statistics):

Moment(vals, 5);
CentralMoment(vals, 5);
Moment(X, 5);
CentralMoment(X, 5);

MGF(X, t);
cumulantCumulant[vals, 1]
Cumulant[X, 1]

CumulantGeneratingFunction[X, t]
with(Statistics):

Cumulant(vals, 1);
Cumulant(X, 1);

CGF(X, t);
entropy
 
Entropy[vals]
mode
 
Commonest[{1, 2, 2, 2, 3, 3, 8, 12}]with(Statistics):

Mode([1, 2, 2, 2, 3, 3, 8, 12]);
quantile statisticsMin[vals]
Median[vals]
Max[vals]
InterquartileRange[vals]
Quantile[vals, 9/10]
with(Statistics):

min(vals);
Median(vals);
max(vals);
InterquartileRange(vals);
Quantile(vals, 9/10);
load(distrib);

/* Other distributions have similar functions: */
quantile_normal(9/10, 0, 1);
bivariate statistiscs
correlation, covariance, Spearman's rank
Correlation[{1, 2, 3}, {2, 4, 7}]
Covariance[{1, 2, 3}, {2, 4, 7}]
SpearmanRho[{1, 2, 3}, {2, 4, 7}]
with(Statistics):

Correlation([1, 2, 3], [2, 4, 7]);
Covariance([1, 2, 3], [2, 4, 7]);
data set to frequency tabledata = {1, 2, 2, 2, 3, 3, 8, 12}
(* list of pairs: *)
tab = Tally[data]
(* dictionary: *)
dict = Counts[data]
with(Statistics):

data := [1, 2, 2, 2, 3, 3, 8, 12];
# array of eqns:
Tally(data);
# dictionary:
Tally(data, output=table);
frequency table to data setf = Function[a, Table[a[[1]], {i, 1, a[[2]]}]]
data = Flatten[Map[f, tab]]
bindata = {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
mathematicamaplemaximasagesympy
binomial

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

PDF[X][50]
CDF[X][50]
RandomVariate[X]
with(Statistics):

X := RandomVariable(Binomial(100, 1/2));
PDF(X, 50);
CDF(X, 50);
Sample(X, 1);
load(distrib);

pdf_binomial(x, 50, 1/2);
cdf_binomial(x, 50, 1/2);
random_binomial(50, 1/2);
from sympy.stats import *

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

density(Y).dict[Integer(50)]
P(X < 50)
sample(X)
poissonX = PoissonDistribution[1]with(Statistics):

X := RandomVariable(Poisson(1));
load(distrib);

pdf_poisson(x, 1);
cdf_poisson(x, 1);
random_poisson(1);
# P(X < 4) raises NotImplementedError:
X = Poisson('X', 1)
discrete uniformX = DiscreteUniformDistribution[{0, 99}]with(Statistics):

X := RandomVariable(DiscreteUniform(0, 99));
load(distrib);

/* {1, 2, …, 100}: */
pdf_discrete_uniform(x, 100);
cdf_discrete_uniform(x, 100);
random_discrete_uniform(100);
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, 10]
with(Statistics):

X := RandomVariable(Normal(0, 1));

PDF(X, 0);
CDF(X, 0);
# no inverse cdf
Sample(X, 10);
with(distrib);

pdf_normal(x, 0, 1);
cdf_normal(x, 0, 1);
/* no inverse cdf */
random_normal(0, 1);
X = RealDistribution('gaussian', 1)

X.distribution_function(0)
X.cum_distribution_function(0)
X.cum_distribution_function_inv(0.5)
X.get_random_element()
from sympy.stats import *

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

density(X)(0)
P(X < 0)
??
sample(X)
gammaX = GammaDistribution[1, 1]with(Statistics):

X := RandomVariable(Gamma(1, 1));
with(distrib);

pdf_gamma(x, 1, 1);
cdf_gamma(x, 1, 1);
random_gamma(1, 1);
X = Gamma('X', 1, 1)
exponentialX = ExponentialDistribution[1]with(Statistics):

X := RandomVariable(Exponential(1));
with(distrib);

pdf_exponential(x, );
cdf_exponential(x, 1);
random_exponential(1);
X = Exponential('X', 1)
chi-squaredX = ChiSquareDistribution[2]with(Statistics):

X := RandomVariable(ChiSquare(2));
with(distrib);

pdf_chi2(x, 2);
cdf_chi2(x, 2);
random_chi2(2);
X = RealDistribution('chisquared', 2)X = ChiSquared('X', 2)
betaX = BetaDistribution[10, 90]with(Statistics):

X := RandomVariable(Beta(10, 90));
with(distrib);

pdf_beta(x, 10, 90);
cdf_beta(x, 10, 90);
random_beta(10, 90);
X = RealDistribution('beta', [10, 90])X = Beta('X', 10, 90)
uniformX = UniformDistribution[{0, 1}]with(Statistics):

X := RandomVariable(Uniform(0, 1));
with(distrib);

pdf_continuous_uniform(x, 0, 1);
cdf_continuous_uniform(x, 0, 1);
random_continuous_uniform(0, 1);
X = RealDistribution('uniform', [0, 1])

X.distribution_function(0.5)
X.cum_distribution_function(0.5)
X.cum_distribution_function_inv(0.5)
X.get_random_element()
X = Uniform('X', 0, 1)
student's tX = StudentTDistribution[2]with(Statistics):

X := RandomVariable(StudentT(2));
with(distrib);

pdf_student_t(x, 2);
cdf_student_t(x, 2);
random_student_t(2);
X = RealDistribution('t', 2)X = StudentT('X', 2)
snedecor's FX = FRatioDistribution[1, 1]with(Statistics):

X := RandomVariable(FRatio(1, 1));
with(distrib);

pdf_f(x, 1, 1);
cdf_f(x, 1, 1);
random_f(1, 1);
X = RealDistribution('F', [1, 1])X = FDistribution('X', 1, 1)
empirical density functionX = NormalDistribution[0, 1]
data = Table[RandomVariate[X], {i, 1, 30}]
Y = EmpiricalDistribution[data]
PDF[Y]
with(Statistics):

X := RandomVariable(Normal(0, 1));
Y := EmpiricalDistribution(Sample(X, 30));
Mean(Y);
empirical cumulative distributionX = NormalDistribution[0, 1]
data = Table[RandomVariate[X], {i, 1, 30}]
Y = EmpiricalDistribution[data]
Plot[CDF[Y][x], {x, -4, 4}]
with(Statistics):

X := RandomVariable(Normal(0, 1));
Y := EmpiricalDistribution(Sample(X, 30));
plot(CDF(Y, x), x = -4 .. 4);
statistical tests
mathematicamaplemaximasagesympy
wilcoxon signed-rank test
variable is symmetric around zero
X = UniformDistribution[{-1/2, 1/2}]
data = RandomVariate[X, 100]

(* null hypothesis is true: *)
SignedRankTest[data]

(* alternative hypothesis is true: *)
SignedRankTest[data + 1/2]
load(distrib); load(stats);

data: makelist(
  random_continuous_uniform(-1/2, 1/2),
  i, 1, 100);

/* null hypothesis is true: */
test_signed_rank(data);

/* alternative hypothesis is true: */
test_signed_rank(data + 1/2);
kruskal-wallis rank sum test
variables have same location parameter
X = NormalDistribution[0, 1]
Y = UniformDistribution[{0, 1}]

(* null hypothesis is true: *)
LocationEquivalenceTest[
  {RandomVariate[X, 100],
   RandomVariate[X, 200]}]

(* alternative hypothesis is true: *)
LocationEquivalenceTest[
  {RandomVariate[X, 100],
   RandomVariate[Y, 200]}]
load(distrib); load(stats);

x1: makelist(
  random_normal(0, 1),
  i, 1, 100);
x2: makelist(
  random_normal(0, 1),
  i, 1, 100);
y: makelist(
  random_continuous_uniform(-1/2, 1/2),
  i, 1, 100);

/* null hypothesis is true: */
test_rank_sum(x1, x2);

/* alternative hypothesis is true: */
test_rank_sum(x1, y);
kolmogorov-smirnov test
variables have same distribution
X = NormalDistribution[0, 1]
Y = UniformDistribution[{-1/2, 1/2}]

(* null hypothesis is true: *)
KolmogorovSmirnovTest[RandomVariate[X, 200], X]

(* alternative hypothesis is true: *)
KolmogorovSmirnovTest[RandomVariate[X, 200], Y]
one-sample t-test
mean of normal variable with unknown variance is zero
X = NormalDistribution[0, 1]

(* null hypothesis is true: *)
TTest[RandomVariate[X, 200]]

(* alternative hypothesis is true: *)
TTest[RandomVariate[X, 200] + 1]
with(Statistics):

X := RandomVariable(Normal(0, 1));

# null hypothesis is true:
OneSampleTTest(Sample(X, 100), 0);

# alternative hypothesis is true:
OneSampleTTest(Sample(X, 100) +~ 1, 0);
independent two-sample t-test
two normal variables have same mean
X = NormalDistribution[0, 1]

(* null hypothesis is true: *)
TTest[
  {RandomVariate[X, 100],
   RandomVariate[X, 200]}]

(* alternative hypothesis is true: *)
TTest[
  {RandomVariate[X, 100],
   RandomVariate[X, 100] + 1}]
with(Statistics):

X := RandomVariable(Normal(0, 1));
x := Sample(X, 100);

# null hypothesis is true:
TwoSampleTTest(x, Sample(X, 200), 0);

# alternative hypothesis is false:
TwoSampleTTest(x, Sample(X, 200) +~ 1, 0);
paired sample t-test
population has same mean across measurements
one-sample binomial test
binomial variable parameter are as given
two-sample binomial test
parameters of two binomial variables are equal
chi-squared test
parameters of multinomial variable are all equal
poisson test
parameter of poisson variable is as given
F test
ratio of variance of normal variables are the same
X = NormalDistribution[0, 1]
Y = NormalDistribution[0, 2]

(* null hypothesis is true: *)
FisherRatioTest[
  {RandomVariate[X, 100],
   RandomVariate[X, 200]}]

(* alternative hypothesis is true: *)
FisherRatioTest[
  {RandomVariate[X, 100],
   RandomVariate[Y, 100]}]
with(Statistics);

X := RandomVariable(Normal(0, 1));
Y := RandomVariable(Normal(0, 2));
x := Sample(X, 100);

# null hypothesis is true:
TwoSampleFTest(x, Sample(X, 100), 1);

# alternative hypothesis is true:
TwoSampleFTest(x, Sample(Y, 100), 1);
pearson product moment test
normal variables are not correlated
X = NormalDistrubtion[0, 1]
x = RandomVariate[X, 100]
y = RandomVariate[X, 100]
x2 = x + RandomVariate[X, 100]
data1 = Inner[List, x, y, List]
data2 = Inner[List, x, x2, List]

(* null hypothesis is true: *)
CorrelationTest[data1, 0, "PearsonCorrelation"]

(* alternative hypothesis is true: *)
CorrelationTest[data2, 0, "PearsonCorrelation"]
spearman rank test
variables are not correlated
X = UniformDistribution[{0, 1}]
x = RandomVariate[X, 100]
y = RandomVariate[X, 100]
x2 = x + RandomVariate[X, 100]
data1 = Inner[List, x, y, List]
data2 = Inner[List, x, x2, List]

(* null hypothesis is true: *)
CorrelationTest[data1, 0, "SpearmanRank"]

(* alternative hypothesis is true: *)
CorrelationTest[data2, 0, "SpearmanRank"]
shapiro-wilk test
variable has normal distribution
X = NormalDistribution[0, 1]
Y = UniformDistribution[{0, 1}]

(* null hypothesis is true: *)
ShapiroWilkTest[RandomVariate[X, 100]]

(* alternative hypothesis is true: *)
ShapiroWilkTest[RandomVariate[Y, 100]]
with(Statistics);

X := RandomVariable(Normal(0, 1));
Y := RandomVariable(Uniform(0, 1));

# null hypothesis is true:
ShapiroWilkWTest(Sample(X, 100));

# alternative hypothesis is true:
ShapiroWilkWTest(Sample(Y, 100));
load(distrib); load(stats);

x: makelist(
  random_normal(0, 1),
  i, 1, 100);
y: makelist(
  random_continuous_uniform(-1/2, 1/2),
  i, 1, 100);

/* null hypothesis is true: */
test_normality(x);

/* alternative hypothesis is true: */
test_normality(y);
bartlett's test
two or more normal variables have same variance
levene's test
two or more variables have same variance
X = NormalDistribution[0, 1]
Y = NormalDistribution[0, 2]

(* null hypothesis is true: *)
LeveneTest[
  {RandomVariate[X, 100],
   RandomVariate[X, 200]}]

(* alternative hypothesis is true: *)
LeveneTest[
  {RandomVariate[X, 100],
   RandomVariate[Y, 100]}]
one-way anova
two or more normal variables have same mean
Needs["ANOVA‘"]

X = NormalDistribution[0, 1]
ones = Table[1, {i, 1, 100}]
x1 = Inner[
  List, ones, RandomVariate[X, 100], List]
x2 = Inner[
  List, 2 * ones, RandomVariate[X, 100], List]
x3 = Inner[
  List, 3 * ones, RandomVariate[X, 100], List]
y = Inner[
  List,
  3 * ones,
  RandomVariate[X, 100] + 0.5,
  List]

(* null hypothesis is true: *)
ANOVA[Join[x1, x2, x3]]

(* alternative hypothesis is true: *)
ANOVA[Join[x1, x2, y]]
two-way anova
bar charts
mathematicamaplemaximasagesympy
vertical-bar-chart.jpg
vertical bar chart
BarChart[{7, 3, 8, 5, 5},
  ChartLegends->
    {"a","b","c","d","e"}]
with(Statistics):

ColumnGraph(
  ["a" =7, "b"=3, "c"=8, "d"=5, "e"=5]);
x: [7, 3, 8, 5, 5];
labs: [a, b, c, d, e];
data: makelist(makelist(labs[i], j, x[i]), i, 5);
wxbarsplot(flatten(data));
horizontal-bar-chart.jpg
horizontal bar chart
BarChart[{7, 3, 8, 5, 5}, BarOrigin -> Left]with(Statistics):

BarChart(["a"=7, "b"=3, "c"=8, "d"=5, "e"=5]);
none
grouped-bar-chart.jpg
grouped bar chart
data = {{7, 1}, {3, 2}, {8, 1}, {5, 3}, {5, 1}}
BarChart[data]
with(Statistics):

ColumnGraph([[7, 3, 8, 5, 5], [1, 2, 1, 3, 1]]);
x: [7, 3, 8, 5, 5];
y: [1,2,1,3,1];
labs: [a, b, c, d, e];
d1: makelist(makelist(labs[i], j, x[i]), i, 5);
d2: makelist(makelist(labs[i], j, y[i]), i, 5);
wxbarsplot(flatten(d1), flatten(d2));
stacked-bar-chart.jpg
stacked bar chart
data = {{7, 1}, {3, 2}, {8, 1}, {5, 3}, {5, 1}}
BarChart[data, ChartLayout -> "Stacked"]
with(Statistics):

ColumnGraph([[7, 3, 8, 5, 5], [1, 2, 1, 3, 1]],
  format = stacked);
x: [7, 3, 8, 5, 5];
y: [1,2,1,3,1];
labs: [a, b, c, d, e];
d1: makelist(makelist(labs[i], j, x[i]), i, 5);
d2: makelist(makelist(labs[i], j, y[i]), i, 5);
wxbarsplot(flatten(d1), flatten(d2),
  grouping=stacked);
pie-chart.jpg
pie chart
PieChart[{7, 3, 8, 5, 5}]with(Statistics):

PieChart(["a"=7, "b"=3, "c"=8, "d"=5, "e"=5]);
x: [7, 3, 8, 5, 5];
labs: [a, b, c, d, e];
data: makelist(makelist(labs[i], j, x[i]), i, 5);
wxpiechart(flatten(data));
histogram.jpg
histogram
X = NormalDistribution[0, 1]
(* 2nd arg is approx number of bins: *)
Histogram[RandomReal[X, 100], 10]
with(Statistics):

X := RandomVariable(Normal(0, 1));
Histogram(Sample(X, 100));
load(distrib);

data: makelist(random_normal(0, 1), i, 1, 100);
wxhistogram(data);
box-plot.jpg
box plot
X = NormalDistribution[0, 1]
n100 = RandomVariate[X, 100]
BoxWhiskerChart[n100]

Y = ExponentialDistribution[1]
e100 = RandomVariate[Y, 100]
u100 = RandomReal[1, 100]
data = {n100, e100, u100}
BoxWhiskerChart[data]
with(Statistics):

X := RandomVariable(Normal(0, 1));
BoxPlot(Sample(X, 100));
load(distrib);

data: makelist(random_normal(0, 1), i, 1, 100);
wxboxplot(data);
scatter plots
mathematicamaplemaximasagesympy
strip-chart.jpg
strip chart
X = NormalDistribution[0, 1]
data = {RandomReal[X], 0} & /@ Range[1, 50]
ListPlot[data]
with(Statistics):

X := RandomVariable(Normal(0, 1));
ScatterPlot(Sample(X, 50), [0$i = 1 .. 50]);
strip-chart-jitter.jpg
strip chart with jitter
X = NormalDistribution[0, 1]
Y = UniformDistribution[{-0.05, 0.05}]
data = {RandomReal[X], RandomReal[Y]} & /@
  Range[1, 50]
ListPlot[data,
  PlotRange -> {Automatic, {-1, 1}}]
with(Statistics):

X := RandomVariable(Normal(0, 1));
ScatterPlot(Sample(X, 50),
  jitter = true,
  view = [-2 .. 2, -19 .. 20]);
scatter-plot.jpg
scatter plot
X = NormalDistribution[0, 1]
rand = Function[RandomReal[X]]
data = {rand[], rand[]} & /@ Range[1, 50]
ListPlot[data]
with(Statistics):

X := RandomVariable(Normal(0, 1));
ScatterPlot(Sample(X, 50), Sample(X, 50));
load(distrib);

x: makelist(random_normal(0, 1), i, 1, 50);
y: makelist(random_normal(0, 1), i, 1, 50);
wxplot2d([discrete, x, y], [style, points]);
additional-point-set.jpg
additional point set
X = NormalDistribution[0, 1]
rand = Function[RandomReal[X]]
data1 = {rand[], rand[]} & /@ Range[1, 50]
data2 = {rand[]+1, rand[]+1} & /@ Range[1, 50]
Show[ListPlot[data1],
  ListPlot[data2, PlotStyle -> Red]]
with(Statistics):
with(plots):

X := RandomVariable(Normal(0, 1));
p1 := plot(Sample(X, 50), Sample(X, 50),
  color = black, style = point);
p2 := plot(Sample(X, 50) +~ 1,
  Sample(X, 50) +~ 1, color = red,
  style = point);
display([p1, p2]);
load(distrib);

x1: makelist(random_normal(0, 1), i, 1, 50);
y1: makelist(random_normal(0, 1), i, 1, 50);
x2: makelist(random_normal(0, 1), i, 1, 50) + 1;
y2: makelist(random_normal(0, 1), i, 1, 50); + 1;
wxplot2d([[discrete, x1, y1],
    [discrete, x2, y2]],
  [style, points], [color, black, red]);
point typesListPlot[data, PlotMarkers -> {"*"}]

(* shows standard sequence of point types: *)
Graphics`PlotMarkers[]

(* The elements of the PlotMarkers array can be strings, symbols, expressions, or images. *)
ScatterPlot(Sample(X, 50), Sample(X, 50),
  symbol = asterisk);

(* possible symbol values:

  asterisk
  box
  circle
  cross
  diagonalcross
  diamond
  point
  solidbox
  solidcircle
  soliddiamond
*)
wxplot2d([discrete, x, y],
  [style, points],
  [point_type, asterisk]);

/* possible point_type values:

  asterisk
  box
  bullet
  circle
  diamond
  plus
  square
  times
  triangle

The bullet and box are filled versions of circle and square.
*/
point sizeX = NormalDistribution[0, 1]
rand = Function[RandomReal[X]]
data = {rand[], rand[]} & /@ Range[1, 50]
(* point size is fraction of plot width: *)
ListPlot[data, PlotStyle -> {PointSize[0.03]}]
scatter-plot-matrix.jpg
scatter plot matrix
Needs["StatisticalPlots`"]

X = NormalDistribution[0, 1]
x = RandomReal[X, 50]
y = RandomReal[X, 50]
z = x + 3 * y
w = y + RandomReal[X, 50]
PairwiseScatterPlot[Transpose[{x, y, z, w}]]
load(distrib);

x: makelist(random_normal(0, 1), i, 1, 50);
y: makelist(random_normal(0, 1), i, 1, 50);
z: x + 3 * y;
w: y + makelist(random_normal(0, 1), i, 1, 50);
wxscatterplot(transpose(matrix(x, y, z, w)));
scatter-plot-3d.jpg
3d scatter plot
X = NormalDistribution[0, 1]
data = RandomReal[X, {50, 3}]
ListPointPlot3D[data]
with(Statistics):

X := RandomVariable(Normal(0, 1));
ScatterPlot3D([[Sample(X, 50)], [Sample(X, 50)],
    [Sample(X, 50)]]);
none
bubble-chart.jpg
bubble chart
X = NormalDistribution[0, 1]
data = RandomReal[X, {50, 3}]
BubbleChart[data]
with(Statistics);

X := RandomVariable(Normal(0, 1));
Z := RandomVariable(DiscreteUniform(1, 10));
# appears to ignore the 3rd argument and
# make a scatter plot:

BubblePlot(Sample(X, 50), Sample(X, 50),
  Sample(Z, 50));
linear-regression-line.jpg
linear regression line
data = Table[{i, 2 * i + RandomReal[{-5, 5}]},
  {i, 0, 20}]
model = LinearModelFit[data, x, x]
Show[ListPlot[data],
  Plot[model["BestFit"],
    {x, 0, 20}]]
with(Statistics);

X := RandomVariable(Normal(0, 1));
x := [i$i = 1 .. 100];
y := x +~ 10 * Sample(X, 100);
ScatterPlot(x, y, fit = [a*t+b, t]);
load(distrib);
load(lsquares);

X: makelist(i, i, 50);
Y: makelist(X[i] + random_normal(0, 1), i, 50);
M: transpose(matrix(X, Y));
fit: lsquares_estimates(M, [x, y], y = A*x + B,
  [A, B]);
A: second(fit[1][1]), numer;
B: second(fit[1][2]), numer;
Xhat: makelist(A*X[i] + B, i, 50);
wxplot2d([[discrete, X, Y], [discrete, X, Xhat]],
  [style, points, lines], [color, black, red]);
q-q-plot.jpg
quantile-quantile plot
X = NormalDistribution[0, 1]
data1 = RandomReal[1, 50]
data2 = RandomReal[X, 50]
QuantilePlot[data1, data2]
with(Statistics):

X := RandomVariable(Normal(0, 1));
Y := RandomVariable(Uniform(0, 1));
QuantilePlot(Sample(X, 200), Sample(Y, 200));
load(distrib);

x: makelist(random_continuous_uniform(0, 1),
  i, 200);
y: makelist(random_normal(0, 1), i, 200);
wxplot2d([discrete, sort(x), sort(y)],
  [style, points]);
line charts
mathematicamaplemaximasagesympy
polygonal-line-plot.jpg
polygonal line plot
X = NormalDistribution[0, 1]
rand = Function[RandomReal[X]]
f = Function[i, {i, rand[]}]
data = f /@ Range[1, 20]
ListLinePlot[data]
with(Statistics):

X := RandomVariable(Normal(0, 1));
pts := Sample(X, 20);
plot([i$i = 1 .. 20], pts);
load(distrib);

x: makelist(random_normal(0, 1), i, 1, 20);
wxplot2d([discrete, makelist(i, i, 20), x]);
additional-line.jpg
additional line
X = NormalDistribution[0, 1]
data1 = RandomReal[X, 20]
data2 = RandomReal[X, 20]
ListLinePlot[{data1, data2}
  PlotStyle->{Black, Red}]
with(Statistics):
with(plots):

X := RandomVariable(Normal(0, 1));
p1 := plot([i$i = 1 .. 50], Sample(X, 50),
  color = black);
p2 := plot([i$i = 1 .. 50], Sample(X, 50),
  color = red);
display(p1, p2);
load(distrib);

x: makelist(random_normal(0, 1), i, 1, 20);
y: makelist(random_normal(0, 1), i, 1, 20);
wxplot2d([[discrete, makelist(i, i, 20), x],
    [discrete, makelist(i, i, 20), y]],
  [color, black, red]);
line typesListLinePlot[data, PlotStyle -> Dashed]

(* PlotStyle values:

  Dashed
  DotDashed
  Dotted
*)
with(Statistics):

X := RandomVariable(Normal(0, 1));
plot([i$i = 1 .. 50], Sample(X, 50),
  linestyle = dot);

(* linestyle values:

  solid
  dot
  dash
  dashdot
  longdash
  spacedash
  spacedot
*)
none
line thicknessX = NormalDistribution[0, 1]
data1 = RandomReal[X, 20]
data2 = RandomReal[X, 20]
(* thickness is fraction of plot width: *)
ListLinePlot[{data1, data2},
  PlotStyle -> {Thickness[0.01], Thickness[0.02]}]
with(Statistics);
with(plots);

X := RandomVariable(Normal(0, 1));
p1 := plot([i$i = 1 .. 50], Sample(X, 50),
  thickness = 1);
p2 := plot([i$i = 1 .. 50], Sample(X, 50),
  thickness = 3);
display(p1, p2)
function-plot.jpg
function plot
Plot[Sin[x], {x, -4, 4}]plot(sin(x), x = -4 .. 4);wxplot2d(sin(x), [x, -4, 4]);
parametric-plot.jpg
parametric plot
ParametricPlot[{Sin[u], Sin[2 * u]},
  {u, 0, 2 * Pi}]
plot([sin(t), sin(2*t), t = 0 .. 2*Pi]);wxplot2d([parametric, sin(t), sin(2*t),
  [t, 0, 2*%pi]]);
implicit-plot.jpg
implicit plot
ContourPlot[x^2 + y^2 == 1, {x, -1, 1},
  {y, -1, 1}]
with(plots, implicitplot);

implicitplot(x^2+y^2 = 1, x = -1 .. 1,
  y = -1 .. 1);
load(implicit_plot);

wximplicit_plot(x^2 + y^2 = 1, [x, -1, 1],
  [y, -1, 1]);
polar-plot.jpg
polar plot
PolarPlot[Sin[3 * t], {t, 0, Pi}]with(plots):

polarplot(sin(3*u), u = 0 .. Pi);
f(x) := sin(3 * x);
wxplot2d([parametric, cos(t)*f(t), sin(t)*f(t),
  [t, 0, %pi]]);
cubic-spline.jpg
cubic spline
X = NormalDistribution[0, 1]
data = Table[{i, RandomReal[X]},
  {i, 0, 20}]
f = Interpolation[data,
  InterpolationOrder -> 3]
Show[ListPlot[data],
  Plot[f[x], {x, 0, 20}]]
with(CurveFitting):
with(Statistics):

X := RandomVariable(Normal(0, 1));
x := [i$i = 1 .. 20];
y := convert(Sample(X, 20), list);
plot(Spline(zip(`[]`, x, y), z), z = 0 .. 21);
load(interpol);
load(distrib);
load(draw);

data: makelist([i, random_normal(0, 1)], i, 20);
cspline(data);
f(x):=''%;
wxdraw2d(explicit(f(x),x,0,20));
area-chart.jpg
area chart
data = {{7, 1, 3, 2, 8}, {1, 5, 3, 5, 1}}
stacked = {data[[1]], data[[1]] + data[[2]]}
ListLinePlot[stacked, Filling ->
  {1 -> {Axis, LightBlue},
   2 -> {{1}, LightRed}}]
with(Statistics):

AreaChart([[7, 1, 3, 2, 8], [1, 5, 3, 5, 1]],
  format = stacked);
surface charts
mathematicamaplemaximasagesympy
contour-plot.jpg
contour plot
(* of function: *)
ContourPlot[x * (y - 1), {x, 0, 10},
  {y, 0, 10}]

(* of data: *)
X = NormalDistribution[0, 1]
rand = Function[RandomReal[X]]
data = Table[x * (y - 1) + 5 * rand[],
  {x, 0, 10}, {y, 0, 10}]
ListContourPlot[data]
with(plots):

contourplot(x*(y-1), x = 0 .. 10, y = 0 .. 10);
wxcontour_plot(x * (y-1), [x, 0, 10],
  [y, 0, 10]);
heat-map.jpg
heat map
(* of function: *)
DensityPlot[Sin[x] * Sin[y],
  {x, -4, 4},
  {y, -4, 4}]

(* of data: *)
X = NormalDistribution[0, 1]
rand = Function[RandomReal[X]]
data = Table[x * y + 10 * rand[],
  {x, 1, 10},
  {y, 1, 10}]
ListDensityPlot[data]
with(plots):

densityplot(sin(x)*sin(y), x = -4..4, y = -4..4);
wxplot3d (sin(x) * sin(y), [x,-4,4], [y,-4,4],
  [elevation, 0], [azimuth, 0],
  [grid, 100, 100], [mesh_lines_color, false]);
shaded-surface-plot.jpg
shaded surface plot
Plot3D[Exp[-(x^2 + y^2)], {x, -2, 2},
  {y, -2, 2}, MeshStyle -> None]
plot3d(exp(-x^2-y^2), x = -2 .. 2,
  y = -2 .. 2, style = surface);
light sourcelot3D[Exp[-(x^2 + y^2)],
  {x, -2, 2}, {y, -2, 2},
  MeshStyle -> None,
  Lighting -> {{"Point", White, {5, -5, 5}}}]
mesh-surface-plot.jpg
mesh surface plot
Plot3D[Exp[-(x^2 + y^2)], {x, -2, 2},
  {y, -2, 2}, Lighting -> {White},
  PlotStyle -> White]
plot3d(exp(-x^2-y^2), x = -2 .. 2,
  y = -2 .. 2, style = wireframe)
wxplot3d(exp(-(x^2 + y^2)),
  [x, -2, 2], [y, -2, 2],
  [palette, false], [color, black]);
view point(* (x, y, z) coordinates;
   (0, 0, 3) is from above: *)

Plot3D[Exp[-(x^2 + y^2)],
  {x, -2, 2}, {y, -2, 2},
  MeshStyle -> None,
  ViewPoint -> {0, 0, 3}]
vector-field-plot.jpg
vector field plot
StreamPlot[{x^2 + y, 1 + x - y^2}, {x, -4, 4}, {y, -4, 4}]with(plots):

fieldplot([x^2 + y, y^2 + x + 1],
  x = -4 .. 4, y = -4 .. 4);
plotdf([x^2 + y, 1 + x - y^2], [x, -4, 4],
  [y, -4, 4]);
chart options
mathematicamaplemaximasagesympy
chart title(* title on top by default *)
Plot[Sin[x], {x, -4, 4},
  PlotLabel -> "title example"]
# title on top; caption on bottom

plot(sin(x), title = "title example",
  caption = "caption example")
wxplot2d(sin(x), [x, -4, 4],
  [title, "title example"]);

data: 1, 1, 2, 2, 3, 3, 3, 3, 4];
wxboxplot(data, title="title example");

/* title on top by default */
axis labeldata = Table[{i, i^2}, {i, 1, 20}]
ListLinePlot[data, AxesLabel -> {x, x^2}]
plot([i$i = 1 .. 20], [i*i$i = 1 .. 20],
  labels = [x, x^2]);
x: makelist(i, i, 20);
y: makelist(i^2, i, 20);
wxplot2d([discrete, x, y],
  [xlabel, "x"], [ylabel, "x^2"]);
legend.jpg
legend
X = NormalDistribution[0, 1]
data1 = RandomReal[X, 20]
data2 = RandomReal[X, 20]
ListLinePlot[{data1, data2},
    PlotLegends -> {"first", "second"}]
with(Statistics);
with(plots);

X := RandomVariable(Normal(0, 1));
p1 := plot([i$i = 1 .. 50], Sample(X, 50),
  color = black, legend = "first");
p2 := plot([i$i = 1 .. 50], Sample(X, 50),
  color = red, legend = "second");
display(p1, p2)
/* wxplot2d includes a legend by default.

Provide [legend, false] as argument to suppress it. */
data labeldata = {{313, 3.7}, {62, .094}, {138, 6.6},
  {113, 0.76}, {126, 0.15}}

(* The {0, -1} argument of Text[] centers the label above the data point. *)
Show[ListPlot[data,
    AxesLabel -> {"pop", "area"}],
  Graphics[Text["USA", data[[1]], {0, -1}]],
  Graphics[Text["UK", data[[2]], {0, -1}]],
  Graphics[Text["Russia", data[[3]], {0, -1}]],
  Graphics[Text["Mexico", data[[4]], {0, -1}]],
  Graphics[Text["Japan", data[[5]], {0, -1}]]]
named colorsWhite Gray Black Transparent

Blue Brown Cyan Green Magenta Orange Pink Purple Red Yellow

LightBlue LightBrown LightCyan LightGray LightGreen LightMagenta LightOrange LightPink LightPurple LightRed LightYellow
aquamarine black blue brown coral cyan gold gray green grey khaki magenta maroon navy orange pink plum red sienna tan turquoise violet wheat white yellow

White WhiteSmoke LightGray LightGrey DarkGray DarkGrey Gray Grey DimGray DimGrey Black

AliceBlue AntiqueWhite Aqua Aquamarine Azure Beige Bisque BlanchedAlmond Blue BlueViolet Brown Burlywood CadetBlue Chartreuse Chocolate Coral CornflowerBlue Cornsilk Crimson Cyan DeepPink DeepSkyBlue DodgerBlue Feldspar Firebrick FloralWhite ForestGreen Fuchsia Gainsboro GhostWhite Gold Goldenrod Green GreenYellow Honeydew HotPink IndianRed Indigo Ivory Khaki Lavender LavenderBlush LawnGreen LemonChiffon Lime LimeGreen Linen Magenta Maroon MediumAquamarine MediumBlue MediumOrchid MediumPurple MediumSeaGreen MediumSlateBlue MediumSpringGreen MediumTurquoise MediumVioletRed MidnightBlue MintCream MistyRose Moccasin NavajoWhite Navy NavyBlue OldLace Olive OliveDrab Orange OrangeRed Orchid PaleGoldenrod PaleGreen PaleTurquoise PaleVioletRed PapayaWhip PeachPuff Peru Pink Plum PowderBlue Purple Red RosyBrown RoyalBlue SaddleBrown Salmon SandyBrown SeaGreen Seashell Sienna Silver SkyBlue SlateBlue SlateGray SlateGrey Snow SpringGreen SteelBlue Tan Teal Thistle Tomato Turquoise Violet VioletRed Wheat Yellow YellowGreen

DarkBlue DarkCyan DarkGoldenrod DarkGray DarkGreen DarkGrey DarkKhaki DarkMagenta DarkOliveGreen DarkOrange DarkOrchid DarkRed DarkSalmon DarkSeaGreen DarkSlateBlue DarkSlateGray DarkSlateGrey DarkTurquoise DarkViolet

LightBlue LightCoral LightCyan LightGoldenrod LightGoldenrodYellow LightGray LightGreen LightGrey LightPink LightSalmon LightSeaGreen LightSkyBlue LightSlateBlue LightSlateGray LightSlateGrey LightSteelBlue LightYellow
white gray black

gray0 gray10 gray100 gray20 gray30 gray40 gray50 gray60 gray70 gray80 gray90 gray100 grey grey0 grey10 grey20 grey30 grey40 grey50 grey60 grey70 grey80 grey90 grey100

aquamarine beige blue brown coral cyan forest_green gold goldenrod green khaki magenta medium_blue midnight_blue navy orange orange_red pink plum purple red royalblue salmon sea_green skyblue spring_green turquoise violet yellow

dark_blue dark_cyan dark_goldenrod dark_gray dark_green dark_grey dark_khaki dark_magenta dark_orange dark_pink dark_red dark_salmon dark_turquoise dark_violet dark_yellow

light_blue light_coral light_cyan light_goldenrod light_gray light_green light_grey light_magenta light_pink light_red light_salmon light_turquoise light_yellow
rgb colorRGBColor[1, 0, 0]
(* with opacity: *)
RGBColor[1, 0, 0, 0.5]
ColorTools:-Color([1, 0, 0])[color, "#FF0000"]
background colorPlot[Sin[x], {x, 0, 2 Pi},
  Background -> Black,
  PlotStyle -> White,
  AxesStyle -> White,
  TicksStyle -> White,
  GridLines -> Automatic,
  GridLinesStyle -> White]
p1 := plot([1, -1], x = -4 .. 4,
  filled = true,
  color = black,
  axis = [color = white]);
p2 := plot(sin(x), x = -4 .. 4,
  color = white,
  axis = [color = white]);
display(p1, p2);
axis limitsPlot[x^2, {x, 0, 20},
  PlotRange -> {{0, 20}, {-200, 500}}]
plot([i$i = 1 .. 20], [i*i$i = 1 .. 20],
  x = 0 .. 20, y = -200 .. 500);
logarithmic y-axisLogPlot[{x^2, x^3, x^4, x^5},
  {x, 0, 20}]
with(plots):

logplot([x^2, x^3, x^4, x^5], x = 0 .. 20);
x: makelist(i, i, 20);
wxplot2d([
    [discrete, x, makelist(i^2, i, 20)],
    [discrete, x, makelist(i^3, i, 20)]],
  [logy, true]);
aspect ratio(* aspect ratio is height divided by width: *)
Plot[Sin[x], {x, 0, 2 Pi}, AspectRatio -> 0.25]

(* In the notebook, dragging the corner of an image increases or decreases the size, but aspect ratio is preserved. *)
# size is width and height in pixels:
plot(sin(x), size = [800, 200]);

# In the notebook, dragging the corner
# of an image can change width and height
# independently
wxplot2d(sin(x), [x, -4, 4],
  [yx_ratio, 0.25]);

/* Image size can’t be changed in notebook. */
ticksPlot[Sin[x], {x, 0, 2 Pi}, Ticks -> None]

Plot[Sin[x], {x, 0, 2 Pi},
  Ticks -> {{0, Pi, 2*Pi}, {-1, 0, 1}}]
# rm x and y ticks:
plot(sin(x), x = -4 .. 4, tickmarks = [0, 0]);

plot(sin(x), x = -4 .. 4,
  tickmarks = [[0, Pi, 2*Pi], [-1, 0, 1]])
wxplot2d(sin(x), [x, -4, 4],
  [xtics, -4, 2, 4],
  [ytics, -1, 0.5, 1]);
grid linesPlot[Sin[x], {x, 0, 2 Pi},
  GridLines -> Automatic]

Plot[Sin[x], {x, 0, 2 Pi},
  GridLines -> {{0, 1, 2, 3, 4, 5, 6},
    {-1, -0.5, 0, 0.5, 1}}]
plot(sin(x), x = -4 .. 4, gridlines);
subplot-grid.jpg
grid of subplots
GraphicsGrid[Table[Table[
      Histogram[RandomReal[X, 100], 10],
      {i, 1, 2}],
    {j, 1, 2}]]
with(plots): with(Statistics):

A := Array(1 .. 2, 1 .. 2);
X := RandomVariable(Normal(0, 1));
for i while i <= 2 do
  for j while j <= 2 do
    A[i, j] := Histogram(Sample(X, 100));
  end do;
end do;
Display(A);
load(distrib);

x: makelist(makelist(random_normal(0, 1), i, 50),
  j, 4);
p: makelist(histogram_description(x[i]), i, 4);
wxdraw(gr2d(p[1]), gr2d(p[2]), gr2d(p[3]),
  gr2d(p[4]), columns=2);
save plot as pngExport["hist.png",
  Histogram[RandomReal[X, 100], 10]]
with(plottools): with(Statistics):

X := RandomVariable(Normal(0, 1));
plt := BoxPlot(Sample(X, 100));
exportplot("boxplot.png", plt);
After creating a plot, run gnuplot on the .gnuplot file generated in the home directory:

$ gnuplot maxout.gnuplot
____________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

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.

Grammar and Invocation

interpreter

How to execute a script.

mathematica:

The full path to MathKernel on Mac OS X:

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

repl

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

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.

multiple line comment

The syntax for a delimited comment which can span lines.

Variables and Expressions

assignment

How to perform assignment.

Mathematica, Sympy, and Pari/GP support the chaining of assignments. For example, in Mathematica one can assign the value 3 to x and y with:

x = y = 3

In Mathematica and Pari/GP, assignments are expressions. In Mathematica, the following code is legal and evaluates to 7:

(x = 3) + 4

In 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.

increment and decrement

Increment and decrement operators which can be used in expressions.

non-referential identifier

An identifier which does not refer to a value.

A non-referential identifier will usually print as a string containing its name.

Expressions containing non-referential identifiers will not be evaluated, though they may be simplified.

Non-referential identifiers represent "unknowns" or "parameters" when performing algebraic derivations.

identifier as value

How to get a value referring to an identifier.

The identifier may be the name of a variable containing a value. But the value referring to the identifier is distinct from the value in the variable.

One may manipulate a value referring to an identifier even if it is not the name of a variable.

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

The null literal.

null test

How to test if a value is null.

undefined variable access

What happens when an undefined variable is used in an expression.

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.

pari/gp:

A vector or matrix evaluates to false if all components evaluate to false.

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.

pari/gp:

When the argument of random() is an integer n, it generates an integer in the range $$\{0, ..., n - 1\}$$ .

When the argument is a arbitrary precision float, it generates a value in the range [0.0, 1.0]. The precision of the argument determines the precision of the random number.

random seed

How to set or get the random seed.

mathematica:

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

bit operators

binary, octal, and hex literals

Binary, octal, and hex integer literals.

mathematica:

The notation works for any base from 2 to 36.

radix

Convert a number to a representation using a given radix.

to array of digits

Convert a number to an array of digits representing the number.

Strings

string literal

The syntax for a string literal.

newline in literal

Are newlines permitted in string literals.

literal escapes

Escape sequences for putting unusual characters in string literals.

concatenate

How to concatenate strings.

translate case

How to convert a string to all lower case letters or all upper case letters.

trim

How to remove whitespace from the beginning or the end of string.

number to string

How to convert a number to a string.

string to number

How to parse a number from a string.

string join

How to join an array of strings into a single string, possibly separated by a delimiter.

split

How to split a string in to an array of strings. How to specify the delimiter.

substitute

How to substitute one or all occurrences of substring with another.

length

How to get the length of a string in characters.

index of substring

How to get the index of the first occurrence of a substring.

extract substring

How to get a substring from a string using character indices.

character literal

The syntax for a character literal.

character lookup

How to get a character from a string by index.

chr and ord

Convert a character code point to a character or a single character string.

Get the character code point for a character or single character string.

delete characters

Delete all occurrences of a set of characters from a string.

Arrays

sectionmathematicamaplemaximasympy
arraysListlistlistlist
multidimensional arraysListArrayarraynone
vectorsListVectorlistMatrix
matricesListMatrixmatrixMatrix

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.

array of integers as index

manipulate back

manipulate front

head

tail

cons

concatenate

replicate

copy

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

iterate

reverse

sort

dedupe

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.

map

filter

reduce

universal and existential tests

min and max element

shuffle and sample

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

flatten

zip

How to interleave two arrays.

cartesian product

Sets

Arithmetic Sequences

Dictionaries

record literal

record member access

Functions

definition

invocation

function value

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.

Streams

Files

Directories

Libraries and Namespaces

Reflection

function documentation

How to get the documentation for a function.

Symbolic Expressions

A defining feature of computer algebra systems is symbolic expressions.

In most programming languages, evaluating an expression with an undefined variable results in an error. Some languages assign a default value to variables so that such expressions can be evaluated.

In a CAS, undefined variables are treated as unknowns; expressions which contains them are symbolic expressions. When evaluating them, if the unknowns cannot be eliminated, the expression cannot be reduced to a numeric value. The expression then evaluates to a possibly simplified or normalized version of itself. Symbolic expressions are first class values; they can be stored in variables or passed to functions. An application of symbolic expressions is a function which solves a system of equations. Without symbolic expressions, it would be awkward for the caller to specify the equations to be solved.

literal

How to create a symbolic expression.

In most CAS systems, any expression an undefined variables is automatically a a symbolic expression.

sympy:

In SymPy, unknowns must be declared. This is a consequence of SymPy being implemented as a library in a language which throws exceptions when undefined variables are encountered.

prevent simplification

variable update

Do symbolic expressions "see" changes to the unknown variables they contain.

substitute

piecewise-defined expression

simplify

assumption

assumption predicates

list assumptions

remove assumption

Calculus

limit

limit at infinity

one-sided limit

derivative

derivative of a function

constants

higher order derivative

mixed partial derivative

div, grad, and curl

antiderivative

definite integral

improper integral

double integral

find poles

residue

sum

series sum

series expansion of function

omitted order term

product

Equations and Unknowns

solve equation

solve equations

differential equation

differential equation with boundary condition

differential equations

recurrence equation

Optimization

An optimization problem consists of a real-valued function called the objective function.

The objective function takes one or more input variables. In the case of a maximization problem, the goal is to find the value for the input variables where the objective function achieves its maximum value. Similarly for a minimization function one looks for the values for which the objective function achieves its minimum value.

minimize

How to solve a minimization problem in one variable.

maximize

How to solve a maximization problem.

We can use a function which solves minimization problems to solve maximization problems by negating the objective function. The downside is we might forget the minimum value returned is the negation of the maximum value we seek.

objective with unknown parameter

How to solve an optimization when the objective function contains unknown parameters.

unbounded behavior

What happens when attempting to solve an unbounded optimization problem.

multiple variables

How to solve an optimization problem with more than one input variable.

constraints

How to solve an optimization with constraints on the input variable. The constrains are represented by inequalities.

infeasible behavior

What happens when attempting to solve an optimization problem when the solution set for the constraints is empty.

integer variables

How to solve an optimization problem when the input variables are constrained to linear values.

Vectors

vector literal

The notation for a vector literal.

constant vector

How to create a vector with components all the same.

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.

construct from sequence

constant matrices

diagonal matrices

matrix by formula

dimensions

How to get the number of rows and columns of a matrix.

element lookup

How to access an element of a matrix.

The anguages described here follow the mathematical convention of putting the row index before the column index.

extract row

How to access a row.

extract column

How to access a column.

extract submatrix

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.

product

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.

power

exponential

log

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.

norms

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

transpose

conjugate transpose

inverse

row echelon form

pseudoinverse

determinant

trace

characteristic polynomial

minimal polynomial

rank

nullspace basis

eigenvalues

eigenvectors

spectral decomposition

The spectral decomposition of a square matrix A is a factorization P ⋅D ⋅ P-1 where P is invertible and D is diagonal.

If a spectral decomposition exists, the matrix A is said to be diagonalizable.

If an invertiable matrix P exists such that A = P ⋅ B ⋅ P-1, then A and B are said to be similar.

According to the spectral theorem, and spectral decomposition exists when the matrix A is normal, which means it commutes with its conjugate transpose.

If a matrix A is symmetric, then a spectral decomposition P ⋅ D ⋅ P-1 exists, and moreover P and D are real matrices.

The spectral decomposition is also called the eigendecomposition. The values on the diagonal of D are eigenvalues of the original matrix A.

LUP decomposition

Factor a square matrix A into a lower triangular matrix L, an upper triangular matrix U, and a permutation matrix P.

An LUP factorization of a square matrix always exists.

QR decomposition

jordan decomposition

singular value decomposition

polar decomposition

Combinatorics

Enumerative combinatorics is the study of the size of finite sets. The sets are defined by some property, and we seek a formula for the size of the set so defined.

For some simple examples, let A and B be disjoint sets of size n and m respectively. The size of the union A ∪ B is n + m and the size of the Cartesian product A × B is nm. The size of the power set of A is 2n.

factorial

The factorial function n! is the product of the first n positive integers 1 × 2 × … × n.

It is also the number of permutations or bijective functions on a set of n elements. It is the number of orderings that can be given to n elements.

See the section on permutations below for how to iterate through all n! permutations on a set of n elements.

As the factorial function grows rapidly with n, it is useful to be aware of this approximation:

$$\begin{align} \ln n! \approx n \ln n - n + \frac{1}{2} \ln 2 \pi n \end{align}$$

binomial coefficient

The binomial coefficient can be defined using the factorial function:

$$\begin{align} {n \choose k} = \frac{n!}{(n-k)! k!} \end{align}$$

The binomial coefficient appears in the binomial theorem:

$$\begin{align} (x+y)^n = \sum_{k=0}^{n} {n \choose k} x^k y^{n-k} \end{align}$$

The binomial cofficient $${ n \choose k }$$ is the number of sets of size k which can be drawn from a set of size n without replacement.

multinomial coefficient

rising and falling factorial

subfactorial

A derangement is a permutation on a set of n elements where every element moves to a new location.

The number of derangements is thus less than the number of permutations, n!, and the function for the number of derangements is called the subfactorial function.

Using a exclamation point as a prefix to denote the subfactorial, the following equations hold:

$$\begin{align} !n = n \cdot [!(n-1)] + (-1)^n \end{align}$$
$$\begin{align} !n = n! \sum_{i=0}^n \frac{(-1)^i}{i!} \end{align}$$
$$\begin{align} lim_{n \rightarrow \infty} \frac{!n}{n!} = \frac{1}{e} \end{align}$$

integer partitions

The number of multisets of positive integers which sum to a integer.

There are 5 integer partitions of 4:

    4
    3 + 1
    2 + 2
    2 + 1 + 1
    1 + 1 + 1 + 1

compositions

The number of sequences of positive integers which sum to an integer.

There are 8 compositions of 4:

    4
    3 + 1
    1 + 3
    2 + 2
    2 + 1 + 1
    1 + 2 + 1
    1 + 1 + 2
    1 + 1 + 1 + 1

mathematica:

The NumberOfCompositions and Compositions functions use weak compositions, which include zero as a possible summation.

The number of weak compositions of an integer is infinite, since there is no limit on the number of times zero can appear as a summand. The number of weak compositions of a fixed size is finite, however.

set partitions

bell number

permutations with k disjoint cycles

fibonacci number

bernoulli number

harmonic number

catalan number

Number Theory

divisible test

A test whether an integer a is divisible by another integer b.

Equivalently, does there exists a third integer m such that a = mb.

divisors

The list of divisors for an integer.

pseudoprime test

A fast primality test.

An integer p is prime if for any factorization p = ab, where a and b are integers, either a or b are in the set {-1, 1}.

A number of primality tests exists which give occasional false positives. The simplest of these use Fermat's Little Theorem, in which for prime p and a in $$\{1, ..., p - 1\}$$ :

$$\begin{align} a^{p-1} \equiv 1 \;(\text{mod}\; p) \end{align}$$

The test for a candidate prime p is to randomly choose several values for a in $$\{1, ..., p - 1\}$$ and evaluate

$$\begin{align} a^{p-1} \;(\text{mod}\; p) \end{align}$$

If any of them are not equivalent to 1, then the test shows that p is not prime. Unfortunately, there are composite numbers n, the Carmichael numbers, for which

$$\begin{align} a^{n-1} \equiv 1 \;(\text{mod}\; n) \end{align}$$

holds for all a in $$\{1, ..., n - 1\}$$ .

A stronger test is the Miller-Rabin primality test. Given a candidate prime n, we factor n - 1 as 2rd where d is odd. If n is prime, then one of the following must be true:

$$\begin{align} a^d \equiv 1 \;(\text{mod}\;n) \end{align}$$
$$\begin{align} a^{2^r \cdot d} \equiv -1 \;(\text{mod}\;n) \end{align}$$

Thus, one checks the above two equations for a small number of primes. If we use all primes p ≤ 41, then it is known that there are no false positives for n ≤ 3 × 1024.

Since pseuodoprime tests are known which are correct for all integers up to a very large size, and since conclusively showing that a number is prime is a slow operation for larger integers, a true prime test is often not practical.

prime factors

The list of prime factors for an integer, with their multiplicities.

next prime

The smallest prime number greater than an integer. Also the greatest prime number smaller than an integer.

nth prime

The n-th prime number.

prime counting function

The number of primes less than or equal to a value.

According to the prime number theorem:

$$\begin{align} \lim_{n \rightarrow \infty} \frac{\pi(n)}{n/\log n} = 1 \end{align}$$

divmod

The quotient and remainder.

If the divisor is positive, then the remainder is non-negative.

greatest common divisor

The greatest common divisor of a pair of integers. The divisor is always positive.

Two integers are relatively prime if their greatest common divisor is one.

extended euclidean algorithm

How to express a greatest common divisor as a linear combination of the integers it is a GCD of.

The functions described return the GCD in addition to the coefficients.

least common multiple

The least common multiple of a pair of integers.

The LCM can be calculated from the GCD using this formula:

$$\begin{align} \text{lcm}(m, n) = \frac{|m\cdot n|}{\text{gcd}(m, n)} \end{align}$$

power modulus

Raise an integer to a integer power, modulo a third integer.

Euler's theorem can often be used to reduce the size of the exponent.

integer residues

The integer residues or integers modulo n are the equivalence classes formed by the relation

$$\begin{align} a\;(\text{mod}\;n) = b\; (\text{mod}\;n) \end{align}$$

An element in of these equivalence classes is called a representative. We can extend addition and multiplication to the residues by performing integer addition or multiplication on representatives. This is well-defined in the sense that it does not depend on the representatives chosen. Addition and multiplication defined this way turn the integer residues into commutative rings with identity.

multiplicative inverse

How to get the multiplicative inverse for a residue.

If the representative for a residue is relatively prime to the modulus, then the residue has a multiplicative inverse. For that matter, if the modulus n is a prime, then the ring of residues is a field.

Note that we cannot in general find the inverse using a representative, since the only units in the integers are -1 and 1.

By Euler's theorem, we can find a multiplicative inverse by raising it to the power $$\phi(n) - 1$$ :

$$\begin{align} a^{\phi(n) - 1} \cdot a = a^{\phi(n)} \equiv 1 \;(\text{mod}\;n) \end{align}$$

When a doesn't have a multiplicative inverse, then we cannot cancel it from both sides of a congruence. The following is true, however:

$$\begin{align} az \equiv az' \;(\text{mod}\; n) \iff z \equiv z' \;\left(\text{mod}\; \frac{n}{\text{gcd}(a, n)}\right) \end{align}$$

chinese remainder theorem

A function which finds a solution to a system of congruences.

The Chinese remainder theorem asserts that there is a solution x to the system of k congruences

$$\begin{align} x \equiv a_i \;(\text{mod}\;n_i) \end{align}$$

provided that the ni are pairwise relatively prime. In this case there are an infinite number of solutions, all which are equal modulo $$N = n_1 \cdots n_k$$ . For this reason the solution is returned as a residue modulo N.

lift integer residue

How to get a representative from the equivalence class of integers modulo n.

Typically an integer in $$\{0, ..., n - 1\}$$ is chosen. A centered lift chooses a representative x such that $$-n/2 < x \leq n/2$$ .

euler totient

The Euler totient function is defined for any positive integer n as:

$$\begin{align} \phi(n) = n \prod_{p | n} \frac{p - 1}{p} \end{align}$$

Note that the product is over the primes that divide n.

The Euler totient is the number of integers in $$\{1, ..., n - 1\}$$ which are relatively prime to n. It is thus the size of the multiplicative group of integers modulo n.

The Euler totient appears in Euler's theorem:

$$\begin{align} a^{\phi(n)} \equiv 1 \;(\text{mod}\;n) \end{align}$$

carmichael function

The smallest number k such that ak ≡ 1 (mod n) for all residues a.

By Euler's theorem, the Carmichael function λ(n) is less that or equal to the Euler totient function φ(n). The functions are equal when there are primitive roots modulo n.

multiplicative order

The multiplicative order of a residue a is the smallest exponent k such that

$$\begin{align} a^k \equiv 1\;(\text{mod}\;n) \end{align}$$

In older literature, it is sometimes said that a belongs to the exponent k modulo n.

primitive roots

A primitive root is a residue module n with multiplicative order φ(n).

The multiplicative group is not necessarily cyclic, though it is when n is prime. If it is not cyclic, then there are no primitive roots.

Any primitive root is a generator for the multiplicative group, so it can be used to find the other primitive roots.

discrete logarithm

For a residue x and a base residue b, find a positive integer such that:

$$\begin{align} b^k \equiv x\;(\text{mod}\; n) \end{align}$$

quadratic residues

A quadratic residue is a non-zero residue a which has a square root modulo p. That is, there is x such that

$$\begin{align} x^2 \equiv a \;(\text{mod}\;p) \end{align}$$

If a is non-zero and doesn't have a square root, then it is a quadratic non-residue.

discrete square root

How to find the square root of a quadratic residue.

kronecker symbol

The Legendre symbol is used to indicate whether a number is a quadratic residue and is defined as follows:

$$\begin{align} \left( \frac{a}{p} \right) = \begin{cases} \;\; 1 \;\;\; a \; \text{is a quadratic residue} \\ \;\; 0 \;\;\; p \mid a \\ -1 \;\;\; a \; \text{is a quadratic nonresidue} \end{cases} \end{align}$$

The Legendre symbol is only defined when p is an odd prime, but if n is an odd positive integer with prime factorization

$$\begin{align} p_1^{\alpha_1} \ldots p_n^{\alpha_n} \end{align}$$

then the Jacobi symbol is defined as

$$\begin{align} \left( \frac{a}{n} \right) = \left( \frac{a}{p_1} \right)^{\alpha_1} \ldots \left( \frac{a}{p_n} \right)^{\alpha_n} \end{align}$$

The Kronecker symbol is a generalization of the Jacobi symbol to all integers, but we omit the details.

moebius function

The Möbius function μ(n) is 1, -1, or 0 depending upon when n is a square-free integer with an even number of prime factors, a square-free integer with an odd number of prime factors, or an integer which is divisible by p2 for some prime p.

The Möbius function is multiplicative: when a and b are relatively prime, μ(a)μ(b) = μ(ab).

The Möbius function appears in the Möbius inversion formula. If g and f are possibly complex-valued functions defined on the natural numbers such that for all integers n ≥ 1:

$$\begin{align} g(n) = \sum_{d | n} f(d) \end{align}$$

then for all integers n ≥ 1:

$$\begin{align} f(n) = \sum_{d | n} \mu(d) g(n | d) \end{align}$$

riemann zeta function

The Riemann zeta function is a complex-valued function defined as the analytic continuation of this series:

$$\begin{align} \zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s} \end{align}$$

The function has zeros (called the trivial zeros) at -2, -4, …. All other zeros must lie in the strip 0 ≤ ℜ(z) ≤ 1. In 1859 Riemann conjectured that all non-trivial zeros are on the line ℜ(z) = 1/2.

continued fraction

Convert a real number to a continued fraction.

A continued fraction is a sequence of integers a0, a1, …, an representing the fraction:

$$\begin{align} a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{\ddots + {\frac{1}{a_n}}}}} \end{align}$$

The sequence can even be infinite, in which case the fraction is the limit of the rational numbers defined by taking the first n digits in the sequence.

A continued fraction for a real number can be computed using the Euclidean algorithm. In the case of a rational number, one starts with the numerator and the denominator. In the case of a rational number, one can start with the number itself and 1.

A continued fraction is finite if and only if the number is a rational.

A continued fraction repeats if and only if it is a quadratic irrational.

convergents

The first n digits of a continued fraction define a sequence of rational numbers called the convergents. The rational numbers converge to the number defined by the continued fraction.

Each convergent r/q is the closest rational number to the continued fraction with denominator of size q or smaller.

Polynomials

literal

extract coefficient

extract coefficients

from array of coefficients

degree

expand

factor

[[ #collect-terms-note]]

collect terms

roots

quotient and remainder

greatest common divisor

extended euclidean algorithm

resultant

discriminant

groebner basis

specify ordering

elementary symmetric polynomial

symmetric reduction

cyclotomic polynomial

hermite polynomial

chebyshev polynomial

interpolation polynomial

spline

add fractions

partial fraction decomposition

pade approximant

Trigonometry

eliminate powers and products of trigonometric functions

eliminate sums and multiples inside trigonometric functions

trigonometric to exponential

exponential to trigonometric

fourier expansion

periodic functions on unit interval

fourier transform

heaviside step function

dirac delta

Special Functions

gamma function

The gamma function is defined for all complex numbers except the non-positive integers.

For positive integers, the following equation holds:

$$\begin{align} \Gamma(n) = (n-1)! \end{align}$$

If the real part of t is positive, then

$$\begin{align} \Gamma(t) = \int_0^\infty x^{t-1} e^{-x} dx \end{align}$$

error function

The error function is function from ℝ to [-1, 1] defined by:

$$\begin{align} \mathrm{erf}(x) = \frac{2}{\sqrt(\pi)} \int_0^x e^{-t^2} dt \end{align}$$

The complementary error function is

$$\begin{align} \mathrm{erfc}(x) = 1 - erf(x) \end{align}$$

The cumulative distribution of the standard normal distribution is related to the error function by scaling:

$$\begin{align} \Phi(x) = \frac{1}{2} + \frac{1}{2} \mathrm{erf}(\frac{x}{\sqrt(2)}) = \frac{1}{2} \mathrm{erfc}(\frac{-x}{\sqrt(2)}) \end{align}$$

hyperbolic functions

Definitions of the hyperbolic functions:

$$\begin{align} \mathrm{sinh}\;x = \frac{e^x - e^{-x}}{2} \end{align}$$
$$\begin{align} \mathrm{cosh}\;x = \frac{e^x + e^{-x}}{2} \end{align}$$
$$\begin{align} \mathrm{tanh}\;x = \frac{\mathrm{sinh}\;x}{\mathrm{cosh}\;x} \end{align}$$

sinh and cosh are odd and even functions, respectively. Like ex and e-x, sinh and cosh span the linear space of solutions to y''(x) = y(x).

elliptic functions

bessel functions

Permutations

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

The notation that Mathematica uses 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

from disjoint cycles

to disjoint cycles

from array

from two arrays with same elements

size

support

act on element

act on list

compose

inverse

power

order

number of inversions

parity

Permutations are classified as even or odd based on the number of inversions.

The composition of two even permutations is even.

to inversion vector

from inversion vector

list permutations

random permutation

Descriptive Statistics

Distributions

Statistical Tests

A selection of statistical tests. For each test the null hypothesis of the test is stated in the left column.

In a null hypothesis test one considers the p-value, which is the chance of getting data which is as or more extreme than the observed data if the null hypothesis is true. The null hypothesis is usually a supposition that the data is drawn from a distribution with certain parameters.

The extremeness of the data is determined by comparing the expected value of a parameter according to the null hypothesis to the estimated value from the data. Usually the parameter is a mean or variance. In a one-tailed test the p-value is the chance the difference is greater than the observed amount; in a two-tailed test the p-value is the chance the absolute value of the difference is greater than the observed amount.

wilcoxon signed-rank test

A non-parametric est whether a variable is drawn from a distribution that is symmetric about zero.

Often this test is used to test that the mean of the distribution is zero.

kruskal-wallis rank sum test

A non-parametric test whether variables have the same mean.

For two variables, this test is the same as the Mann-Whitney test.

maxima:

The Maxima function only supports testing two variables.

kolmogorov-smirnov test

Test whether two samples are drawn from the same distribution.

one-sample t-test

Student's t-test determines whether a sample drawn from a normal distribution has mean zero.

The test can be used to test for a different mean value; just subtract the value from each value in the sample.

One may know in advance that the sample is drawn from a normal distribution. For example, if the values in the sample are each means of large samples, then the distribution is normal by the central limit theorem.

The Shapiro-Wilk test can be applied to determine if the values come from a normal distribution.

If the distribution is not known to be normal, the Wilcoxon signed-rank test can be used instead.

The Student's t-test used the sample to estimate the variance, and as a result the test statistic has a t-distribution.

By way of contrast, the z-test assumes that the variance is known in advance, and simply scales the data to get a z-score, which has standard normal distribution.

independent two-sample t-test

Test whether two normal variables have same mean.

paired sample t-test

A t-test used when the same individuals are measure twice.

one-sample binomial test

two-sample binomial test

chi-squared test

poisson test

F test

pearson product moment test

pearson spearman rank test

shapiro-wilk test

bartlett's test

A test whether variables are drawn from normal distributions with the same variance.

levene's test

A test whether variables are drawn from distributions with the same variance.

one-way anova

two-way anova

Bar Charts

vertical bar chart

A chart in which the height of bars is used to represent a list of numbers.

maxima:

Maxima plots the frequency of the values, and not the values themselves. Non-positive values cannot be represented.

horizontal bar chart

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

grouped bar chart

stacked bar chart

pie chart

maxima:

Note that Maxima plots the frequency of the values, and not the values themselves.

histogram

A histogram is a bar chart in which each bar represents the frequency of values in a data set within a range. The width of the bars can be used to indicate the ranges.

box plot

Scatter Plots

strip chart

A strip chart represents a list of values by points on a line. The values are converted to pairs by assigning the y-coordinate a constant value of zero. Pairs are then displayed with a scatter plot.

strip chart with jitter

A strip chart in which in which a random variable with small range is used to fill the y-coordinate. Jitter makes it easier to see how many values are in dense regions.

scatter plot

How to plot a list of pairs of numbers by representing the pairs as points in the (x, y) plane.

additional point set

How to add a second list of pairs of numbers to a scatter plot. Color can be used to distinguish the two data sets.

point types

How to select the symbols used to mark data points. Choice of symbols can be use to distinguish data sets.

point size

How to change the size of the symbols used to mark points.

scatter plot matrix

A scatter plot matrix is a way of displaying a multivariate data set by means of a grid of scatter plots. Off-diagonal plots are scatter plots of two of the variables. On-diagonal plots can be used to to display the name or a histogram of one of the variables.

3d scatter plot

How to represent a list of triples of numbers by points in (x, y, z) space.

bubble chart

How to represent a list of triples of numbers by position in the (x, y) plane and size of the point marker.

It is probably better to associate the 3rd component of each triple with the area, and not the diameter of the point marker, but in general bubble charts suffer from ambiguity.

linear regression line

How to add a linear regression line to a scatter plot.

quantile-quantile plot

Line Charts

polygonal line plot

additional line

line types

line thickness

function plot

parametric plot

implicit plot

polar plot

cubic spline

area chart

Surface Charts

contour plot

heat map

shaded surface plot

light source

mesh surface plot

view point

vector field plot

Chart Options

chart title

axis label

legend

data label

named colors

rgb color

background color

axis limits

logarithmic y-axis

aspect ratio

ticks

grid lines

grid of subplots

save plot as png

Mathematica

Mathematica Documentation Center
WolframAlpha

Maple

http://www.maplesoft.com/support/help/

Maxima

http://maxima.sourceforge.net/docs/manual/maxima.html

Sage

http://doc.sagemath.org/html/en/index.html

SymPy

Welcome to SymPy’s documentation!

issue tracker | content of this page licensed under creative commons attribution-sharealike 3.0