Introduction to exact real and complex arithmetic with Calcium

Author: Fredrik Johansson

This notebook gives an introduction to Calcium (https://fredrikj.net/calcium/), using the builtin ctypes-based Python wrapper.

The source notebook file: https://github.com/fredrik-johansson/calcium/blob/master/doc/introduction.ipynb

To run this Notebook locally, you need to build Calcium and add calcium/pycalcium to your Python path.

Caution: the ctypes wrapper is incomplete, slow, and has some problems (e.g. memory leaks). It is mainly intended for testing the library. A more robust Cython interface should eventually be developed for "production use", along with a Julia interface. Calcium itself is also still experimental.

Imports

In [1]:
# For convenience
from IPython.display import display

# Import the Calcium ctypes wrapper
# Don't "import *" at home -- this is for demonstration use only
from pyca import *

Quick example: Euler's identity $e^{\pi i} + 1 = 0$

In [2]:
exp(pi * i) + 1
Out[2]:
$$0$$

Basic types

The core types implemented in Calcium are the following:

  • Symbolic expressions - pyca.fexpr (fexpr_t in C)
  • Algebraic numbers - pyca.qqbar (qqbar_t in C)
  • Calcium numbers / field elements - pyca.ca (ca_t in C)

The types have different roles, and it is important to understand their differences.

Symbolic expressions are the most flexible representation: they preserve the exact form of the input without performing any computation or simplification. For example, $\sqrt{2} / 2$, $1 / \sqrt{2}$ and $2^{1/2} / 2$ are all different as expressions (although they represent the same number). Symbolic expressions are not very useful on their own for computations, but they are convenient as input and output for other, more "computational" types.

Algebraic numbers represent elements of the field $\overline{\mathbb{Q}}$ in canonical form, consisting of the number's minimal polynomial together with an isolating complex interval for a unique root. Thus one number can only have one representation: $\sqrt{2} / 2$ and $1 / \sqrt{2}$ will evaluate to exactly the same qqbar algebraic number. This is useful because the results are predictable and in a sense maximally simplified (contrary to the other types, it is impossible to construct a complicated qqbar instance that represents 0 without trivially being 0). The downsides are that field operations are much more expensive than in the ca representation, and in passing from a symbolic expression to a qqbar, structural information (such as a simple closed form) may be lost.

Calcium numbers represent numbers as elements of fields $\mathbb{Q}(a_1,\ldots,a_n)$. Calcium constructs extension numbers $a_k$ and fields $\mathbb{Q}(a_1,\ldots,a_n)$ automatically and lazily, so from the point of view of the user Calcium numbers behave just like numbers. The extension numbers $a_k$ can be "absolute" algebraic numbers (qqbar instances), but they can also be transcendental numbers like $\pi$ or $\exp(\sqrt{2} i)$. This representation is highly efficient for arithmetic, but in general does not guarantee a canonical form. Relations between extension elements (e.g. $\log(4) / \log(2) = 2$) are simplified automatically using ideal reduction, but Calcium will not be able to prove all relations, so it is possible to have a ca instance that represents 0 without trivially being 0 (such instances are nevertheless handled in a mathematically rigorous way as discussed below). The ca type can also represent special values and metavalues (signed and unsigned infinities, undefined, unknown).

Symbolic expressions

In the pyca namespace, lowercase names (example: sqrt) denote ca functions while while uppercase names (example: Sqrt) denote fexpr symbolic expressions. It is easy to construct symbolic expressions:

In [3]:
# Create a symbolic sum of three terms
expr = Add(Sqrt(2) / 2, 1 / Sqrt(2), 2 ** Div(1, 2) / 2)
# Display the formula
expr
Out[3]:
$$\frac{\sqrt{2}}{2} + \frac{1}{\sqrt{2}} + \frac{{2}^{1 / 2}}{2}$$
In [4]:
type(expr)
Out[4]:
pyca.fexpr

Evaluating expressions and displaying objects

Constant expressions can be passed as input to other types, resulting in evaluation:

In [5]:
ca(expr)     # evaluate expr using ca arithmetic, producing a ca
Out[5]:
$$\frac{3 a_{1}}{2}\; \text{ where } a_{1} = \sqrt{2}$$
In [6]:
type(_)
Out[6]:
pyca.ca
In [7]:
qqbar(expr)     # evaluate expr using qqbar arithmetic, producing a qqbar
Out[7]:
$$\left(\text{Root }\, x \approx {2.12132} \;\text{ of } \;{2 x^{2}-9}\right)$$
In [8]:
type(_)
Out[8]:
pyca.qqbar

Symbolic expressions are also generated automatically behind the scenes in order to display objects in LaTeX in the notebook (this is done in the examples above). You can of course also print objects in text form:

In [9]:
print(expr)
print(ca(expr))
print(qqbar(expr))

Add(Div(Sqrt(2), 2), Div(1, Sqrt(2)), Div(Pow(2, Div(1, 2)), 2))
2.12132 {(3*a)/2 where a = 1.41421 [a^2-2=0]}
2.12132 (deg 2)

By default, qqbar objects display as polynomial roots. To produce a closed-form expression (if possible), the qqbar.fexpr() method can be called:

In [10]:
qqbar(expr).fexpr()
Out[10]:
$$\frac{3 \sqrt{2}}{2}$$

Manipulating symbolic expressions

The fexpr type provides methods for rudimentary manipulation (accessing subexpressions and so on).

In [11]:
expr = Add(Sqrt(2) / 2, 1 / Sqrt(2), 2 ** Div(1, 2) / 2)
print(expr.head())
print(expr.nargs())
print(expr.args())
print(expr.is_atom())
print(expr.args()[0].args()[1].is_atom())
print(expr.num_leaves())
print(expr.nwords())
print(expr.replace(2, Pi))

Add
3
(Div(Sqrt(2), 2), Div(1, Sqrt(2)), Div(Pow(2, Div(1, 2)), 2))
False
True
16
24
Add(Div(Sqrt(Pi), Pi), Div(1, Sqrt(Pi)), Div(Pow(Pi, Div(1, Pi)), Pi))

The .nstr() method computes a numerical approximation using Arb, returning a decimal string:

In [12]:
print((Sqrt(2) / 2).nstr(30))
print((Exp(1 + 2j).nstr(10)))

# No symbolic simplification is done - Arb can generally not detect
# exact zeros, and zeros will be output in the form 0e-n
print((Exp(Pi*NumberI)).nstr(30))
print((Sqrt(2)/2 - 1/Sqrt(2)).nstr())

0.707106781186547524400844362105
-1.131204384 + 2.471726672*I
-1.00000000000000000000000000000 + 0e-37*I
0e-731

The .expanded_normal_form() method puts the given formula in a canonical form as a formal rational expression.

In [13]:
x = fexpr("x"); y = fexpr("y")
A = (x+y)**5 * (x-y) * (x + 1)
B = (x**2 - y**2) * (x**2 - 1)
(A / B).expanded_normal_form()
Out[13]:
$$\frac{{x}^{4} + 4 {x}^{3} y + 6 {x}^{2} {y}^{2} + 4 x {y}^{3} + {y}^{4}}{x-1}$$

Please note that .expanded_normal_form() only simplifies rational arithmetic operations, treating anything non-arithmetical as an atomic node. For example, square roots are treated as atomic. It also does not simplify nodes recursively.

In [14]:
display((Sqrt(x) / Sqrt(x)).expanded_normal_form())
display((Sqrt(x) * Sqrt(x)).expanded_normal_form())
display((Sqrt(2*x) / Sqrt(x*2)).expanded_normal_form())
$$1$$
$${\left(\sqrt{x}\right)}^{2}$$
$$\frac{\sqrt{2 x}}{\sqrt{x \cdot 2}}$$

We will not do anything more sophisticated with symbolic expressions in this notebook; we will instead move on to describing exact numerical calculations using the ca and qqbar types.

Calcium numbers

Calcium numbers encompass rational numbers, of course:

In [15]:
ca(1) / 3
Out[15]:
$$\frac{1}{3}$$
In [16]:
(ca(1) / 3) * 6
Out[16]:
$$2$$

Irrational numbers result in extensions of $\mathbb{Q}$:

In [17]:
# Alternative syntax: ca(2).sqrt()
sqrt(2)
Out[17]:
$$a_{1}\; \text{ where } a_{1} = \sqrt{2}$$
In [18]:
sqrt(2) ** 2
Out[18]:
$$2$$
In [19]:
2 * pi
Out[19]:
$$2 a_{1}\; \text{ where } a_{1} = \pi$$

Field arithmetic produces numbers represented as formal fraction field elements:

In [20]:
pi * i * sqrt(2)
Out[20]:
$$a_{1} a_{2} a_{3}\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2},\;a_{3} = i$$
In [21]:
(pi * i * sqrt(2)) ** 2   # note the simplifications
Out[21]:
$$-2 a^{2}_{1}\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2},\;a_{3} = i$$
In [22]:
((pi + i + sqrt(2)) / (pi + sqrt(2)))**3
Out[22]:
$$\frac{a^{3}_{1} + 3 a^{2}_{1} a_{2} + 3 a^{2}_{1} a_{3} + 6 a_{1} a_{2} a_{3} + 3 a_{1}- a_{2} + 5 a_{3}}{a^{3}_{1} + 3 a^{2}_{1} a_{2} + 6 a_{1} + 2 a_{2}}\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2},\;a_{3} = i$$

Some more number field arithmetic

Let us construct the golden ratio:

In [23]:
phi = (sqrt(5) + 1) / 2
phi
Out[23]:
$$\frac{a_{1} + 1}{2}\; \text{ where } a_{1} = \sqrt{5}$$

We compute the 200th Fibonacci number using Binet's formula:

In [24]:
(phi**200 - (1-phi)**200) / sqrt(5)
Out[24]:
$$280571172992510140037611932413038677189525$$

Depending on the operations, ca arithmetic may result in different field representations of the same number:

In [25]:
display(sqrt(2)*sqrt(3))
display(sqrt(6))
$$a_{1} a_{2}\; \text{ where } a_{1} = \sqrt{3},\;a_{2} = \sqrt{2}$$
$$a_{1}\; \text{ where } a_{1} = \sqrt{6}$$

The difference simplifies to zero:

In [26]:
display(sqrt(2)*sqrt(3) - sqrt(6))
display(sqrt(2)*sqrt(3) == sqrt(6))
$$0$$ 
True

Calcium will attempt to eliminate more complex extension numbers when it encounters extension numbers that are algebraically dependent. Here, it eliminates $\sqrt{6}$ from the expression, writing the result in terms of $\sqrt{2}$ and $\sqrt{3}$

In [27]:
sqrt(2)*sqrt(3) + sqrt(6)
Out[27]:
$$2 a_{2} a_{3}\; \text{ where } a_{1} = \sqrt{6},\;a_{2} = \sqrt{3},\;a_{3} = \sqrt{2}$$

(Implementation detail: the output shows that $\sqrt{6}$ is still kept as part of the field structure, to aid future simplifications, although it is unused in this particular field element.)

In many cases, it would be desirable to write the above result as an element of the simple algebraic number field $\mathbb{Q}(\sqrt{6})$ or $\mathbb{Q}(\sqrt{2} + \sqrt{3})$ instead. Calcium will always stay within a univariate field when performing field operations starting from a single extension number, but it will not automatically reduce multivariate number fields to univariate fields. In the future, Calcium will offer more customizability so that the user can choose between different behaviors concerning simplification of extension numbers and fields.

Predicates; contexts objects

Predicates in Calcium have mathematically semantics: they return True or False only if Calcium can prove the result. When the truth value is unknown, Calcium says so explicitly; in pyca, this is done by raising an exception (NotImplementedError). Consider, as an example, testing whether $\exp(\varepsilon) = 1$ where $\varepsilon$ is a small number. With $\varepsilon = 10^{-1000}$, Calcium finds that the numbers are not equal:

In [28]:
eps = ca(10) ** (-1000)
exp(eps) == 1
Out[28]:
False

With $\varepsilon = 10^{-10000}$, Calcium fails:

In [29]:
eps = ca(10) ** (-10000)
try:
    exp(eps) == 1
except Exception as e:
    print(e)

unable to decide predicate: equal

The comparison fails because the internal precision limit for numerical evaluation has been exceeded. The precision limit and many other settings are stored in a context object, which also serves as a cache of computed data (such as extension numbers and extension fields). The context object has type ca_ctx. There is a default context called ctx_default, with the following settings:

In [30]:
ctx_default
Out[30]:
ca_ctx(verbose=0, print_flags=3, mpoly_ord=0, prec_limit=4096, qqbar_deg_limit=120, low_prec=64, smooth_limit=32, lll_prec=128, pow_limit=20, use_gb=1, gb_length_limit=100, gb_poly_length_limit=1000, gb_poly_bits_limit=10000, vieta_limit=6)

If we create a new context object with higher prec_limit than the default value of 4096 bits, the computation succeeds:

In [31]:
ctx = ca_ctx(prec_limit=65536)
eps = ca(10, context=ctx) ** (-10000)
exp(eps) == 1
Out[31]:
False

The intention is that the user will be able to create multiple "Calcium fields" for different purposes. Right now, context objects support only limited configurability.

Algebraic number identities

Calcium implements a complete decision procedure for testing equality (or inequality) of algebraic numbers. This functionality is accessible using either qqbar or ca operations.

Rob Corless proposed checking the following identity from Bill Gosper. We will first construct the LHS and RHS as symbolic expressions. We could evaluate them directly with qqbar or ca operations, which would be more efficient, but this way we can print the input.

In [32]:
I = NumberI

lhs = Sqrt(36 + 3*(-54+35*I*Sqrt(3))**Div(1,3)*3**Div(1,3) + \
            117/(-162+105*I*Sqrt(3))**Div(1,3))/3 + \
            Sqrt(5)*(1296*I+840*Sqrt(3)-35*3**Div(5,6)*(-54+35*I*Sqrt(3))**Div(1,3)-\
            54*I*(-162+105*I*Sqrt(3))**Div(1,3)+13*I*(-162+105*I*Sqrt(3))**Div(2,3))/(5*(162*I+105*Sqrt(3)))

rhs = Sqrt(5) + Sqrt(7)

display(Equal(lhs, rhs))
$$\frac{\sqrt{36 + 3 {\left(-54 + 35 i \sqrt{3}\right)}^{1 / 3} \cdot {3}^{1 / 3} + \frac{117}{{\left(-162 + 105 i \sqrt{3}\right)}^{1 / 3}}}}{3} + \frac{\sqrt{5} \left(1296 i + 840 \sqrt{3} - 35 \cdot {3}^{5 / 6} {\left(-54 + 35 i \sqrt{3}\right)}^{1 / 3} - 54 i {\left(-162 + 105 i \sqrt{3}\right)}^{1 / 3} + 13 i {\left(-162 + 105 i \sqrt{3}\right)}^{2 / 3}\right)}{5 \left(162 i + 105 \sqrt{3}\right)} = \sqrt{5} + \sqrt{7}$$

We can check numerically that the expressions agree:

In [33]:
print(lhs.nstr()); print(rhs.nstr())

4.881819288564380 - 0e-20*I
4.881819288564380

Evaluating the expressions as qqbars gives the same result, proving the identity:

In [34]:
display(qqbar(lhs)); display(qqbar(rhs))
$$\left(\text{Root }\, x \approx {4.88182} \;\text{ of } \;{ x^{4}-24 x^{2}+4}\right)$$
$$\left(\text{Root }\, x \approx {4.88182} \;\text{ of } \;{ x^{4}-24 x^{2}+4}\right)$$

Explicitly:

In [35]:
qqbar(lhs) == qqbar(rhs)
Out[35]:
True

We can also perform the verification using ca arithmetic. The equality is not immediately apparent since the two expressions evaluate to elements of different extension fields of $\mathbb{Q}$:

In [36]:
display(ca(lhs))
$$\frac{5 a_{1}- a^{2}_{4} a^{2}_{6} a_{7} + 24 a_{4} a_{6} a_{7}-39 a_{7}}{15 a_{4} a_{6}}\; \text{ where } a_{1} = \sqrt{9 a^{3}_{4} + 36 a^{2}_{4} a^{2}_{6} + 117 a_{4} a_{6}},\;a_{2} = {\left(105 a_{8} a_{9}-162\right)}^{2 / 3},\;a_{3} = {\left(105 a_{8} a_{9}-162\right)}^{1 / 3},\;a_{4} = {\left(35 a_{8} a_{9}-54\right)}^{1 / 3},\;a_{5} = {3}^{5 / 6},\;a_{6} = {3}^{1 / 3},\;a_{7} = \sqrt{5},\;a_{8} = \sqrt{3},\;a_{9} = i$$
In [37]:
display(ca(rhs))
$$a_{1} + a_{2}\; \text{ where } a_{1} = \sqrt{7},\;a_{2} = \sqrt{5}$$

Fortunately, the equality operator manages to prove that the values are really equal:

In [38]:
ca(lhs) == ca(rhs)
Out[38]:
True

Indeed, in this case ca arithmetic simplifies the difference of the expressions to 0 automatically:

In [39]:
ca(lhs) - ca(rhs)
Out[39]:
$$0$$

If this example did not look challenging enough, you may try proving equality of the two huge expressions (involving 7000 operations) given in https://ask.sagemath.org/question/52653/equality-of-algebraic-numbers-given-by-huge-symbolic-expressions/. You will need a bit of string preprocessing to convert the given Sage expressions to fexpr expressions or ordinary Python syntax for use with the ca or qqbar types.

This test problem is implemented in the huge_expr.c program in the Calcium examples directory. It should take a few seconds to run.

Transcendental number identities

Calcium is capable of manipulating numbers involving some transcendental functions such as exp and log. Contrary to the case of algebraic numbers, it does not have a complete decision procedure for transcendental numbers, but it has some decent heuristics.

Exp and log numbers

Basic logarithmic and exponential simplifications work as expected:

In [40]:
log(10**20) / log(100)
Out[40]:
$$10$$
In [41]:
exp(pi) * exp(-pi + log(2))
Out[41]:
$$2$$

Different formulas for the same number may result in different internal field representations, in which case Calcium may yet be able to recognize relations when the numbers are subtracted or divided or when evaluating a predicate:

In [42]:
display(log(sqrt(2)))
display(log(2)/2)
display(log(sqrt(2)) - log(2)/2)
display(log(sqrt(2)) / (log(2)/2))
display(log(sqrt(2)) * (log(2)/2))
log(sqrt(2)) == log(2)/2
$$a_{1}\; \text{ where } a_{1} = \log(a_{2}),\;a_{2} = \sqrt{2}$$
$$\frac{a_{1}}{2}\; \text{ where } a_{1} = \log(2)$$
$$0$$
$$1$$
$$\frac{a^{2}_{2}}{4}\; \text{ where } a_{1} = \log(a_{3}),\;a_{2} = \log(2),\;a_{3} = \sqrt{2}$$
Out[42]:
True

Calcium is aware of branch cuts:

In [43]:
# exp(log(x)) == x
display(exp(log(1 + 10*i)))
# log(exp(x)) != x in general; here we get the correct branch!
display(log(exp(1 + 10*i)))
$$10 a_{1} + 1\; \text{ where } a_{1} = i$$
$$-4 a_{1} a_{2} + 10 a_{2} + 1\; \text{ where } a_{1} = \pi,\;a_{2} = i$$

Let us check the formulas given on the Calcium documentation front page:

In [44]:
log(sqrt(2)+sqrt(3)) / log(5 + 2*sqrt(6))
Out[44]:
$$\frac{1}{2}$$
In [45]:
i**i - exp(pi / ((sqrt(-2)**sqrt(2)) ** sqrt(2)))
Out[45]:
$$0$$
In [46]:
ca(10)**-30 < (640320**3 + 744)/exp(pi*sqrt(163)) - 1 < ca(10)**-29
Out[46]:
True

Trigonometric functions

Calcium does not yet have trigonometric or inverse trigonometric funtions builtin at the C level, but we can "fake" these functions using complex exponentials and logarithms. The functions pyca.sin, pyca.cos, pyca.tan and pyca.atan do precisely this.

In [47]:
sin(3)
Out[47]:
$$\frac{- a^{2}_{1} a_{2} + a_{2}}{2 a_{1}}\; \text{ where } a_{1} = e^{3 a_{2}},\;a_{2} = i$$

Calcium is capable of proving some trigonometric identities:

In [48]:
sin(sqrt(2)/2)**2 + cos(1/sqrt(2))**2
Out[48]:
$$1$$
In [49]:
sin(3 + pi) + sin(3)
Out[49]:
$$0$$
In [50]:
tan(atan(5)) == 5
Out[50]:
True
In [51]:
try:
    atan(tan(1)) == 1
except NotImplementedError:
    print("This simplification does not work yet")

This simplification does not work yet

We test some simplifications involving the Gudermannian function $\operatorname{gd}(x) = 2 \operatorname{atan}(e^x) - \pi/2$:

In [52]:
def gd(x):
    return 2*atan(exp(x))-pi/2

display(sin(gd(1)) - tanh(1))
display(tan(gd(1)) - sinh(1))
display(sin(gd(sqrt(2))) - tanh(sqrt(2)))
display(tan(gd(1)/2) - tanh(ca(1)/2))
$$0$$
$$0$$
$$0$$
$$0$$

Let us try to prove a famous identity; Machin's formula for $\pi$:

In [53]:
lhs = 4*Atan(Div(1,5)) - Atan(Div(1,239))
rhs = Pi / 4
Equal(lhs, rhs)
Out[53]:
$$4 \operatorname{atan}\!\left(\frac{1}{5}\right) - \operatorname{atan}\!\left(\frac{1}{239}\right) = \frac{\pi}{4}$$
In [54]:
print(lhs.nstr())
print(rhs.nstr())

0.7853981633974483
0.7853981633974483

Evaluating the left-hand side using ca arithmetic gives us nothing as simple as $\pi / 4$:

In [55]:
4*atan(ca(1)/5) - atan(ca(1)/239)
Out[55]:
$$\frac{a_{1} a_{3}-4 a_{2} a_{3}}{2}\; \text{ where } a_{1} = \log\!\left(\frac{239 a_{3} + 28560}{28561}\right),\;a_{2} = \log\!\left(\frac{5 a_{3} + 12}{13}\right),\;a_{3} = i$$

Nevertheless, Calcium finds the identity when $\pi$ is given as part of the input:

In [56]:
4*atan(ca(1)/5) - atan(ca(1)/239) - pi/4
Out[56]:
$$0$$
In [57]:
4*atan(ca(1)/5) - atan(ca(1)/239) == pi/4
Out[57]:
True

Here is a more complicated formula:

In [58]:
12*atan(ca(1)/49) + 32*atan(ca(1)/57) - 5*atan(ca(1)/239) + 12*atan(ca(1)/110443)
Out[58]:
$$\frac{-12 a_{1} a_{5} + 5 a_{2} a_{5}-32 a_{3} a_{5}-12 a_{4} a_{5}}{2}\; \text{ where } a_{1} = \log\!\left(\frac{110443 a_{5} + 6098828124}{6098828125}\right),\;a_{2} = \log\!\left(\frac{239 a_{5} + 28560}{28561}\right),\;a_{3} = \log\!\left(\frac{57 a_{5} + 1624}{1625}\right),\;a_{4} = \log\!\left(\frac{49 a_{5} + 1200}{1201}\right),\;a_{5} = i$$
In [59]:
12*atan(ca(1)/49) + 32*atan(ca(1)/57) - 5*atan(ca(1)/239) + 12*atan(ca(1)/110443) - pi/4
Out[59]:
$$0$$

Hyperbolic formulas also work:

In [60]:
atanh = lambda x: -i*atan(i*x)
32*atanh(ca(1)/31) + 24*atanh(ca(1)/49) + 14*atanh(ca(1)/161)
Out[60]:
$$-7 a_{1}-12 a_{2}-16 a_{3}\; \text{ where } a_{1} = \log\!\left(\frac{80}{81}\right),\;a_{2} = \log\!\left(\frac{24}{25}\right),\;a_{3} = \log\!\left(\frac{15}{16}\right),\;a_{4} = i$$
In [61]:
32*atanh(ca(1)/31) + 24*atanh(ca(1)/49) + 14*atanh(ca(1)/161) - log(5)
Out[61]:
$$0$$

Matrices

The ca_mat type provides matrices with ca entries. We look at some examples of basic manipulation:

In [62]:
A = ca_mat([[1, i*pi], [-i*pi, 2]])
A
Out[62]:
$$\displaystyle{\begin{pmatrix}1 & a_{1} a_{2} \\- a_{1} a_{2} & 2\end{pmatrix}}\; \text{ where } a_{1} = \pi,\;a_{2} = i$$
In [63]:
A * A
Out[63]:
$$\displaystyle{\begin{pmatrix}a^{2}_{1} + 1 & 3 a_{1} a_{2} \\-3 a_{1} a_{2} & a^{2}_{1} + 4\end{pmatrix}}\; \text{ where } a_{1} = \pi,\;a_{2} = i$$
In [64]:
display(A.det())
display(A.trace())
display(A.rank())
$$- a^{2}_{1} + 2\; \text{ where } a_{1} = \pi,\;a_{2} = i$$
$$3$$ 
2

Solving linear systems:

In [65]:
B = ca_mat([[1], [2]])
X = A.solve(B)
display(X)
display(A * X)
$$\displaystyle{\begin{pmatrix}\frac{2 a_{1} a_{2}-2}{a^{2}_{1}-2} \\\frac{- a_{1} a_{2}-2}{a^{2}_{1}-2}\end{pmatrix}}\; \text{ where } a_{1} = \pi,\;a_{2} = i$$
$$\displaystyle{\begin{pmatrix}1 \\2\end{pmatrix}}$$

Computing row echelon forms and characteristic polynomials:

In [66]:
ca_mat([[1,pi,2,pi],[1/pi,3,1/(pi+1),4],[1,1,1,1]]).rref()
Out[66]:
$$\displaystyle{\begin{pmatrix}1 & 0 & 0 & \frac{a^{3}_{1}- a^{2}_{1}-2 a_{1}}{3 a^{2}_{1} + 3 a_{1}-2} \\0 & 1 & 0 & \frac{4 a^{2}_{1} + 4 a_{1}-2}{3 a^{2}_{1} + 3 a_{1}-2} \\0 & 0 & 1 & \frac{- a^{3}_{1} + a_{1}}{3 a^{2}_{1} + 3 a_{1}-2}\end{pmatrix}}\; \text{ where } a_{1} = \pi$$
In [67]:
A = ca_mat([[5, pi], [1, -1]])**4
display(A)
display(A.charpoly())
display(A.charpoly()(A))
$$\displaystyle{\begin{pmatrix}a^{2}_{1} + 66 a_{1} + 625 & 8 a^{2}_{1} + 104 a_{1} \\8 a_{1} + 104 & a^{2}_{1} + 18 a_{1} + 1\end{pmatrix}}\; \text{ where } a_{1} = \pi$$
$$ x^{2}+ \left(-2 a^{2}_{1}-84 a_{1}-626\right) x+ \left(a^{4}_{1} + 20 a^{3}_{1} + 150 a^{2}_{1} + 500 a_{1} + 625\right)\; \text{ where } a_{1} = \pi$$
$$\displaystyle{\begin{pmatrix}0 & 0 \\0 & 0\end{pmatrix}}$$

Calcium correctly recognizes singular matrices, even matrices that are nontrivially singular.

In [68]:
A = ca_mat([[pi, pi**2], [pi**3, pi**4]])
A
Out[68]:
$$\displaystyle{\begin{pmatrix}a_{1} & a^{2}_{1} \\a^{3}_{1} & a^{4}_{1}\end{pmatrix}}\; \text{ where } a_{1} = \pi$$
In [69]:
try:
    A.solve(ca_mat([[1], [2]]))
except ZeroDivisionError:
    print("matrix is singular!")

matrix is singular!
In [70]:
A.rref()
Out[70]:
$$\displaystyle{\begin{pmatrix}1 & a_{1} \\0 & 0\end{pmatrix}}\; \text{ where } a_{1} = \pi$$
In [71]:
A.rank()
Out[71]:
1

Exact matrix operations easily result in large expressions:

In [72]:
A = ca_mat([[sqrt(i+j+1) for i in range(6)] for j in range(6)])
A
Out[72]:
$$\displaystyle{\begin{pmatrix}1 & a_{7} & a_{6} & 2 & a_{5} & a_{4} \\a_{7} & a_{6} & 2 & a_{5} & a_{4} & a_{3} \\a_{6} & 2 & a_{5} & a_{4} & a_{3} & 2 a_{7} \\2 & a_{5} & a_{4} & a_{3} & 2 a_{7} & 3 \\a_{5} & a_{4} & a_{3} & 2 a_{7} & 3 & a_{2} \\a_{4} & a_{3} & 2 a_{7} & 3 & a_{2} & a_{1}\end{pmatrix}}\; \text{ where } a_{1} = \sqrt{11},\;a_{2} = \sqrt{10},\;a_{3} = \sqrt{7},\;a_{4} = \sqrt{6},\;a_{5} = \sqrt{5},\;a_{6} = \sqrt{3},\;a_{7} = \sqrt{2}$$
In [73]:
A.det()
Out[73]:
$$-4 a_{1} a_{3} a_{5} a_{6}-20 a_{1} a_{3} a_{5} a_{7}-24 a_{1} a_{3} a_{5}-4 a_{1} a_{3} a_{6} a_{7} + 8 a_{1} a_{3} a_{6} + 136 a_{1} a_{3}-28 a_{1} a_{5} a_{6} a_{7}-116 a_{1} a_{5} a_{6}-88 a_{1} a_{5} a_{7} + 64 a_{1} a_{5} + 112 a_{1} a_{6} a_{7} + 164 a_{1} a_{6}-60 a_{1} a_{7} + 244 a_{1} + 204 a_{3} a_{5} a_{6} a_{7}-96 a_{3} a_{5} a_{6} + 8 a_{3} a_{5} a_{7} + 144 a_{3} a_{5}-152 a_{3} a_{6} a_{7}-240 a_{3} a_{6} + 500 a_{3} a_{7} + 548 a_{3}-216 a_{5} a_{6} a_{7} + 116 a_{5} a_{6} + 628 a_{5} a_{7} + 764 a_{5}-500 a_{6} a_{7} + 24 a_{6}-1440 a_{7}-3868\; \text{ where } a_{1} = \sqrt{11},\;a_{2} = \sqrt{10},\;a_{3} = \sqrt{7},\;a_{4} = \sqrt{6},\;a_{5} = \sqrt{5},\;a_{6} = \sqrt{3},\;a_{7} = \sqrt{2}$$

Eigenvalues and matrix functions

Calcium can calculate exact eigenvalues of matrices, with correct multiplicities. This is accomplished by factoring the characteristic polynomial. Currently, computing polynomial roots will only work for polynomials with very simple structure or with rational entries, so do not expect too much!

In [74]:
A = ca_mat([[1,pi],[-pi,1]])
display(A)

for lamda, mult in A.eigenvalues():
    print("Multiplicity", mult)
    display(lamda)
$$\displaystyle{\begin{pmatrix}1 & a_{1} \\- a_{1} & 1\end{pmatrix}}\; \text{ where } a_{1} = \pi$$ 
Multiplicity 1
$$a_{1} a_{2} + 1\; \text{ where } a_{1} = \pi,\;a_{2} = i$$ 
Multiplicity 1
$$- a_{1} a_{2} + 1\; \text{ where } a_{1} = \pi,\;a_{2} = i$$

We demonstrate computing the Jordan decomposition of a matrix with nontrivial Jordan form:

In [75]:
A = ca_mat([[20,77,59,40], [0,-2,-3,-2], [-10,-35,-23,-15], [2,7,3,1]])
J, P = A.jordan_form(transform=True)
display(J)
display(P)
display(P * J * P.inv())
$$\displaystyle{\begin{pmatrix}a_{1} & 0 & 0 & 0 \\0 & - a_{1} & 0 & 0 \\0 & 0 & -2 & 1 \\0 & 0 & 0 & -2\end{pmatrix}}\; \text{ where } a_{1} = i$$
$$\displaystyle{\begin{pmatrix}\frac{3 a_{1}-11}{5} & \frac{-3 a_{1}-11}{5} & 7 & -\frac{3}{2} \\\frac{11 a_{1} + 68}{65} & \frac{-11 a_{1} + 68}{65} & -2 & 0 \\\frac{-6 a_{1}-17}{13} & \frac{6 a_{1}-17}{13} & 0 & 0 \\1 & 1 & 0 & 1\end{pmatrix}}\; \text{ where } a_{1} = i$$
$$\displaystyle{\begin{pmatrix}20 & 77 & 59 & 40 \\0 & -2 & -3 & -2 \\-10 & -35 & -23 & -15 \\2 & 7 & 3 & 1\end{pmatrix}}$$

We construct a simple matrix and compute its matrix logarithm, which is uses Jordan decomposition internally:

In [76]:
A = ca_mat([[-1, -2], [1, 1]])
display(A)
display(A.log())
$$\displaystyle{\begin{pmatrix}-1 & -2 \\1 & 1\end{pmatrix}}$$
$$\displaystyle{\begin{pmatrix}-\frac{ a_{1}}{2} & - a_{1} \\\frac{a_{1}}{2} & \frac{a_{1}}{2}\end{pmatrix}}\; \text{ where } a_{1} = \pi,\;a_{2} = i$$

We evaluate the exponential of the logarithm, recovering the original matrix. This will only work in very simple cases.

In [77]:
A.log().exp()
Out[77]:
$$\displaystyle{\begin{pmatrix}-1 & -2 \\1 & 1\end{pmatrix}}$$

Another nice example:

In [78]:
B = ca_mat([[0,0,1],[0,1,0],[1,0,0]])
display(B.log())
display(B.log().exp())
$$\displaystyle{\begin{pmatrix}\frac{a_{1} a_{2}}{2} & 0 & -\frac{ a_{1} a_{2}}{2} \\0 & 0 & 0 \\-\frac{ a_{1} a_{2}}{2} & 0 & \frac{a_{1} a_{2}}{2}\end{pmatrix}}\; \text{ where } a_{1} = \pi,\;a_{2} = i$$
$$\displaystyle{\begin{pmatrix}0 & 0 & 1 \\0 & 1 & 0 \\1 & 0 & 0\end{pmatrix}}$$

The following matrix is nilpotent; its exponential is a polynomial expression of the matrix:

In [79]:
A = ca_mat([[10,32,3,-13], [-8,-30,-4,11], [25,90,11,-34], [-6,-28,-5,9]])
display(A.exp())
display(A**0 + A**1 + A**2/2 + A**3/6)
$$\displaystyle{\begin{pmatrix}\frac{19}{2} & 29 & 3 & -\frac{23}{2} \\-\frac{21}{2} & -40 & -\frac{11}{2} & 15 \\\frac{57}{2} & 109 & 15 & -\frac{81}{2} \\-\frac{25}{2} & -53 & -8 & \frac{39}{2}\end{pmatrix}}$$
$$\displaystyle{\begin{pmatrix}\frac{19}{2} & 29 & 3 & -\frac{23}{2} \\-\frac{21}{2} & -40 & -\frac{11}{2} & 15 \\\frac{57}{2} & 109 & 15 & -\frac{81}{2} \\-\frac{25}{2} & -53 & -8 & \frac{39}{2}\end{pmatrix}}$$

We construct the 5x5 Hilbert matrix and check that its eigenvalues satisfy the expected determinant and trace relations:

In [80]:
H = ca_mat([[ca(1)/(i+j+1) for i in range(5)] for j in range(5)])
H
Out[80]:
$$\displaystyle{\begin{pmatrix}1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\\frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\\frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\\frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\\frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9}\end{pmatrix}}$$
In [81]:
eig = H.eigenvalues()
for c, mult in eig:
    display(c)
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {1.56705} \;\text{ of } \;{266716800000 x^{5}-476703360000 x^{4}+92708406000 x^{3}-1022881200 x^{2}+307505 x-1}\right)$$
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {0.208534} \;\text{ of } \;{266716800000 x^{5}-476703360000 x^{4}+92708406000 x^{3}-1022881200 x^{2}+307505 x-1}\right)$$
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {0.0114075} \;\text{ of } \;{266716800000 x^{5}-476703360000 x^{4}+92708406000 x^{3}-1022881200 x^{2}+307505 x-1}\right)$$
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {0.000305898} \;\text{ of } \;{266716800000 x^{5}-476703360000 x^{4}+92708406000 x^{3}-1022881200 x^{2}+307505 x-1}\right)$$
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {3.28793 \cdot 10^{-6}} \;\text{ of } \;{266716800000 x^{5}-476703360000 x^{4}+92708406000 x^{3}-1022881200 x^{2}+307505 x-1}\right)$$
In [82]:
display(sum(c * mult for (c, mult) in eig)); display(H.trace())
display(prod(c ** mult for (c, mult) in eig)); display(H.det())
$$\frac{563}{315}$$
$$\frac{563}{315}$$
$$\frac{1}{266716800000}$$
$$\frac{1}{266716800000}$$

Polynomials

The ca_poly type represents univariate polynomials with ca coefficients. We can construct polynomials and do arithmetic:

In [83]:
x = ca_poly([0,1])
f = 1 + (2 + ca(2).sqrt() * ca.pi()) * x + 3*x**2
f
Out[83]:
$$3 x^{2}+ \left(a_{1} a_{2} + 2\right) x+1\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2}$$
In [84]:
f ** 5
Out[84]:
$$243 x^{10}+ \left(405 a_{1} a_{2} + 810\right) x^{9}+ \left(540 a^{2}_{1} + 1080 a_{1} a_{2} + 1485\right) x^{8}+ \left(180 a^{3}_{1} a_{2} + 1080 a^{2}_{1} + 1620 a_{1} a_{2} + 1800\right) x^{7}+ \left(60 a^{4}_{1} + 240 a^{3}_{1} a_{2} + 1260 a^{2}_{1} + 1560 a_{1} a_{2} + 1590\right) x^{6}+ \left(4 a^{5}_{1} a_{2} + 40 a^{4}_{1} + 200 a^{3}_{1} a_{2} + 880 a^{2}_{1} + 1070 a_{1} a_{2} + 1052\right) x^{5}+ \left(20 a^{4}_{1} + 80 a^{3}_{1} a_{2} + 420 a^{2}_{1} + 520 a_{1} a_{2} + 530\right) x^{4}+ \left(20 a^{3}_{1} a_{2} + 120 a^{2}_{1} + 180 a_{1} a_{2} + 200\right) x^{3}+ \left(20 a^{2}_{1} + 40 a_{1} a_{2} + 55\right) x^{2}+ \left(5 a_{1} a_{2} + 10\right) x+1\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2}$$
In [85]:
ca_poly([1,1,1,1,1,1,1,1,1]).integral()
Out[85]:
$$\frac{1}{9} x^{9} + \frac{1}{8} x^{8} + \frac{1}{7} x^{7} + \frac{1}{6} x^{6} + \frac{1}{5} x^{5} + \frac{1}{4} x^{4} + \frac{1}{3} x^{3} + \frac{1}{2} x^{2}+ x$$

Polynomial division, GCD and other operations work as expected:

In [86]:
display(f ** 5 // f**4)
$$3 x^{2}+ \left(a_{1} a_{2} + 2\right) x+1\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2}$$
In [87]:
display(f**5 % f**4)
$$0$$
In [88]:
(f * (x**2-1)).gcd(f**2 * (x-1))
Out[88]:
$$ x^{3} + \frac{2 a^{2}_{1}-3 a_{1} a_{2} + 2}{3 a_{1} a_{2}-6} x^{2} + \frac{-2 a^{2}_{1} + a_{1} a_{2} + 2}{3 a_{1} a_{2}-6} x-\frac{1}{3}\; \text{ where } a_{1} = \pi,\;a_{2} = \sqrt{2}$$
In [89]:
(x**2 + ca.pi()**2).gcd(x + ca.i() * ca.pi())
Out[89]:
$$ x + a_{1} a_{2}\; \text{ where } a_{1} = \pi,\;a_{2} = i$$
In [90]:
ca_poly([9,6,7,-28,12]).squarefree_part()
Out[90]:
$$ x^{3}-\frac{5}{6} x^{2}-\frac{2}{3} x-\frac{1}{2}$$
In [91]:
# Squarefree factorization
for (fac, mult) in ca_poly([9,6,7,-28,12]).factor_squarefree()[1]:
    print(f"Multiplicity {mult}:")
    display(fac)

Multiplicity 1:
$$ x^{2} + \frac{2}{3} x + \frac{1}{3}$$ 
Multiplicity 2:
$$ x-\frac{3}{2}$$

Finding the roots of a high-degree polynomial:

In [92]:
f = 4*x**7 + 4*x**6 - 11*x**5 - 16*x**4 + x**3 + 15*x**2 + 10*x + 2

for root, mult in f.roots():
    print(f"Multiplicity {mult}:")
    display(root)

Multiplicity 1:
$$a_{1}\; \text{ where } a_{1} = \sqrt{2}$$ 
Multiplicity 1:
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {1.32472} \;\text{ of } \;{ x^{3}- x-1}\right)$$ 
Multiplicity 1:
$$- a_{1}\; \text{ where } a_{1} = \sqrt{2}$$ 
Multiplicity 1:
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {-0.662359 + 0.562280 i} \;\text{ of } \;{ x^{3}- x-1}\right)$$ 
Multiplicity 1:
$$a_{1}\; \text{ where } a_{1} = \left(\text{Root }\, x \approx {-0.662359-0.562280 i} \;\text{ of } \;{ x^{3}- x-1}\right)$$ 
Multiplicity 2:
$$-\frac{1}{2}$$
In [93]:
f(sqrt(2))
Out[93]:
$$0$$