Parent and element

When using Sage, it is important to understand where your objects live. Or more mathematically in which structure they belong to. The number 1 as an integer will behave differently as the number 1 as a rational number...

In Sage this ambient space is called parent:

sage: a = 4
sage: parent(a)

sage: a = 2/3
sage: parent(a)

sage: a = 2.5
sage: parent(a)

sage: a = pi
sage: parent(pi)

The parent determines the behavior of many operations:

sage: x = SR.var('x')
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()

sage: x = polygen(ZZ)
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()

sage: x = polygen(QQ)
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()

sage: x = polygen(QQbar)
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()

Even for real numbers there is a huge variety of them

• exact numbers (integer, rational, algebraic numbers)
• floating point numbers (real or complex)
• intervals and balls (real or complex)
• symbolic (the "big garbage")

Floating point numbers are inexact by nature... it is particularly dangerous in an unstable situation. They should be used with a lot of care.

Exercise 1: stability matters

Look at the following four cells and try to guess what will be the output before executing them.:

sage: x0 = 1/3
sage: print 4*x0 - 1 == x0

sage: x = x0
sage: for i in range(100):
....:     x = 4*x - 1
sage: print x

sage: x0 = RR(1/3)
sage: print 4*x0 - 1 == x0

sage: x = x0
sage: for i in range(100):
....:      x = 4*x - 1
sage: print x

sage: #edit

Exercise 2: parent matters

Compare the following two cells. Which computations is the most accurate? the fastest?:

sage: a = sqrt(2) + sqrt(3) + sqrt(5) + sqrt(7) - 6
sage: for i in range(8):
....:    a = a^4 - 4*a^3 + 4*a^2
sage: print a.numerical_approx()

sage: a = sqrt(2.) + sqrt(3.) + sqrt(5.) + sqrt(7.) - 6
sage: for i in range(8):
....:    a = a^4 - 4*a^3 + 4*a^2
sage: print a.numerical_approx()

sage: #edit

In both cases, add some code to see what is the parent of a.

Coercion

Coercion is what happens when you mix numbers (or more generally objects) of different kinds.

For example the addition of an integer and a rational will be a rational:

sage: parent(5 + 2/3)

But what about:

• the addition of a rational with a polynomial with integer coefficients?
• division of two integers?
• division of two polynomials?

sage: # edit
sage: # edit
sage: # edit

Interval and ball arithmetic

Interval and balls are two kinds of floating point arithmetic that take care of error propagation. These types correspond to pairs of floating point numbers that represent respectively:

• the left and right endpoints of an interval
• the center and the radius of a ball

The parents are RIF (for real interval field) and RBF (for real ball field):

sage: RIF(sqrt(2))

sage: RBF(sqrt(2))

The arithmetic is done in such way that for any operation f(X) applied to an interval or a ball X contains all of {f(x): x in X}:

sage: a = RIF(1/2,1)   # the interval [1/2,1]
sage: print a.endpoints()
sage: b = a^2
sage: print b.endpoints()

sage: a = RBF(RIF(1/2,1))  # the same seen as a ball
sage: print a.endpoints()
sage: b = a^2
sage: print b.endpoints()

You can prove theorems using balls or intervals (which is much harder using floating point)!

You can apply most of standard functions on intervals and balls:

sage: a = RIF(1/3)
sage: print a.cos()
sage: print a.tan()

sage: # etc

Exercise:

Find an example of a singleton interval X = [x,x] and a function f where the f(X) is not a singleton:

sage: # edit
sage: # edit
sage: # edit

Precision

Most floating point parents admit a precision attribute:

sage: R = RealField(64)
sage: R(1/3)

sage: R = RealField(256)
sage: R(1/3)

sage: R = RealIntervalField(64)
sage: R(1/3)

sage: R = RealIntervalField(256)
sage: R(1/3)

Playing with algebraic numbers

Since sage-7.1 it is possible to work with embedded number fields in the set of real numbers:

sage: R.<x> = PolynomialRing(ZZ)
sage: K.<cbrt2> = NumberField(x^3 - 2, embedding=AA(2)**(1/3))

sage: 1 < cbrt2 < 2

sage: continued_fraction(cbrt2)

These numbers are exact by nature, but comparing them can lead to expensive computations. Moreover the complexity of addition/multiplication increase (linearly) with the degree of the number field.

Integrating Abelian differentials

In this section, you are invited to play with Abelian differential in the complex plane. The aim is to:

• compute holonomy of paths $$ int_gamma f(x) dx = int_0^1 f(gamma(t)) gamma'(t) dt.$$
• integrate $f(x) dx$ in other words find numerica solutions to the one parameter family of ODEs $$f left(gamma(t) right)gamma'(t) = e^{i theta}.$$

Two basic functions in Sage to perform this task are:

• numerical_integral
• ode_solver

sage: class DifferentialCaller:
....:       def __init__(self, f):
....:         vars = f.variables()
....:         assert len(vars) == 1
....:         var = vars[0]
....:
....:         x = SR.var('x')
....:         y = SR.var('y')
....:         theta = SR.var('theta')
....:         I = SR('I')
....:         g = exp(I * theta) / f(x + I*y)
....:         g1 = g.real()
....:         self.g1 = fast_callable(g1, vars=[x,y,theta], domain=float)
....:         g2 = g.imag()
....:         self.g2 = fast_callable(g2, vars=[x,y,theta], domain=float)
....:         dF = exp(I * theta) * f.derivative(var)(x + I*y) / (f(x+I*y) * f(x+I*y))
....:         dF = dF.full_simplify()
....:         print dF
....:         self.dg1x = fast_callable(dF.real(), vars=[x,y,theta], domain=float)
....:         self.dg1y = fast_callable((I * dF).real(), vars=[x,y,theta], domain=float)
....:         self.dg2x = fast_callable(dF.imag(), vars=[x,y,theta], domain=float)
....:         self.dg2y = fast_callable((I * dF).imag(), vars=[x,y,theta], domain=float)
....:     def function(self, t, z, params):
....:         return (self.g1(z[0], z[1], params[0]), self.g2(z[0], z[1], params[0]))
....:
....:     def jacobian(self, t, z, params):
....:         return ((self.dg1x(z[0], z[1], params[0]), self.dg2x(z[0],z[1],params[0])),
....:                 (self.dg1y(z[0], z[1], params[0]), self.dg2y(z[0],z[1],params[0])))

sage: f(z) = exp(1/z)
sage: D = DifferentialCaller(f)
-e^(I*theta - 1/(x + I*y))/(x^2 + 2*I*x*y - y^2)
sage: T = ode_solver()
sage: T.function = D.function
sage: #T.jacobian = D.jacobian
sage: G = Graphics()
sage: Z0 = [(-1,0.3)]
sage: theta_step = 0.5
sage: for z0 in Z0:
....:     theta = 0.0
....:     while theta < 2*RR.pi()+theta_step/2:
....:         T.ode_solve(y_0=z0, t_span=[0,1], params=[theta], num_points=100)
....:         G += line2d([z for t,z in T.solution], color='blue')
....:         theta += theta_step
sage: G.show()