Example FldRe_demonstrate-arg (H26E5)

> C := ComplexField(30);
> x := -One(C);
> Imaginary(x); // imaginary part is negative zero
-0.000000000000000000000000000000
> Arg(x); // approximately negative pi as Arg
-3.14159265358979323846264338328
> y := C!-1; // coerce, so imag part is positive zero
> Imaginary(y); // positive zero
0.000000000000000000000000000000
> Arg(y); // approximately positive pi as Arg
3.14159265358979323846264338328
> assert x eq y; // note that x and y are indeed equal!

Modulus(c) : FldComElt -> FldReElt

Given a complex number c, return the real number (to the same precision as c) that is the modulus of c.

Real(c) : FldComElt -> FldReElt

Re(c) : FldComElt -> FldReElt

Given a complex number c=x + y i, return the real part x of c (as a real number to the same precision as c).

Imaginary(c) : FldComElt -> FldReElt

Im(c) : FldComElt -> FldReElt

Given a complex number c=x + y i, return the imaginary part y of c (as a real number to the same precision as c).

Rounding

Round(r) : FldReElt -> FldReElt

Round(r) : FldComElt -> RngQuadElt

Given a real number r, return the integer i for which |r - i| is a minimum. i.e., the integer closest to r. If there are two such integers, the one of larger magnitude is chosen (rounding away from zero). Given a (non-real) complex number r, return the Gaussian integer i for which |r - i| is a minimum, i.e. the Gaussian integer closest to r.

Truncate(r) : FldReElt -> RngIntElt

Given a real number r, return ⌊r⌋ if r is positive, and return -⌊ - r⌋ + 1 if r is negative. Thus, the effect of this function is to round towards zero.

Ceiling(r) : Infty -> Infty

Ceiling(r) : FldReElt -> RngIntElt

The ceiling of the real number r, i.e. the smallest integer greater than or equal to r.

Floor(r) : Infty -> Infty

Floor(r) : FldReElt -> RngIntElt

The floor of the real number r, i.e. the greatest integer less than or equal to r.

Precision

Precision(r) : FldReElt -> RngIntElt

Precision(c) : FldComElt -> RngIntElt

Given a real or complex number c belonging to the real or complex field C, return the decimal precision p to which calculations are performed in C.

BitPrecision(r) : FldReElt -> RngIntElt

BitPrecision(c) : FldComElt -> RngIntElt

Given a real or complex number c belonging to the real or complex field C, return the (internally used) bit precision p to which calculations are performed in C.

Precision(L) : [FldReElt] -> RngIntElt

Precision(L) : [FldComElt] -> RngIntElt

Gives a sequence of real or complex numbers, return the precision p of their parent field.

ChangePrecision(r, n) : FldReElt, RngIntElt -> FldReElt

ChangePrecision(c, n) : FldComElt, RngIntElt -> FldComElt

Coerces the real (r) or complex (c) number into a field of precision n.

Constants

Let R denote a real or complex field. The functions described below will return an approximation of certain constants to the precision associated with a given real or complex field R. If R is real, a real number is returned; if R is complex, a complex number with imaginary part zero is returned.

Catalan(R) : FldRe -> FldReElt

Catalan(R) : FldCom -> FldComElt

The value of Catalan's constant computed to the accuracy associated with the real or complex field R. Catalan's constant is the sum from k equals 0 to infinity of ( - 1)^k by (2k + 1) ^{- 2}. MPFR calculates this constant using formula (31) of Victor Adamchik's document "33 representations for Catalan's constant"footnote{*} {http://www-2.cs.cmu.edu/~{adamchik/articles/catalan/catalan.htm}}, for more information see mpfr.org.

EulerGamma(R) : FldRe -> FldReElt

EulerGamma(R) : FldCom -> FldComElt

The value of Euler's constant computed to the precision of R.

Pi(R) : FldRe -> FldReElt

Pi(R) : FldCom -> FldComElt

The value of π computed to the precision of R.

Simple Element Functions

AbsoluteValue(r) : FldReElt-> FldReElt

AbsoluteValue(r) : FldComElt-> FldReElt

Abs(r) : FldReElt-> FldReElt

Abs(r) : FldComElt-> FldReElt

The absolute value of the real or complex number r.

Sign(r) : FldReElt -> RngIntElt

Return one of the integer values +1, 0, -1 depending upon whether the real number r is positive, zero or negative, respectively.

ComplexConjugate(r) : FldReElt -> FldReElt

ComplexConjugate(c) : FldComElt -> FldComElt

Conjugate(c) : FldComElt -> FldComElt

The complex conjugate x - y i of a complex number x + y i.

Norm(c) : FldComElt -> FldReElt

Norm(r) : FldReElt -> FldReElt

The real norm of a real or complex number c; note that for complex c=x + y i this returns x² + y², while for elements of real domains it just returns the absolute value. The result lies in the same field as the argument.

Root(r, n) : FldReElt, RngIntElt -> FldReElt

Root(r, n) : FldComElt, RngIntElt -> FldComElt

Given a real number R and a positive integer n, calculate the n-th root of r (using Newton's method without divisions) with the same precision. If n is even then r must be non-negative.

SquareRoot(c) : FldComElt -> FldComElt

Sqrt(c) : FldComElt -> FldComElt

SquareRoot(r) : FldReElt -> FldReElt

Sqrt(r) : FldReElt -> FldReElt

Given a real or complex number c, return the square root of r as an element of the same field to which r belongs.

Distance(x, L) : FldReElt, [FldReElt] -> FldReElt, RngIntElt

Distance(x, L) : FldReElt, [FldComElt] -> FldReElt, RngIntElt

Distance(x, L) : FldComElt, [FldReElt] -> FldReElt, RngIntElt

Distance(x, L) : FldComElt, [FldComElt] -> FldReElt, RngIntElt

    Max: ExtReElt                       Default: ∞

Given a sequence L of real or complex numbers and an additional number x compute the distance between x and L, ie. min_y∈L |x - y|, that is the shortest distance between x and any element of L. Furthermore, the index in L of an element realising the distance is returned as a second argument.

Note: prior to V2.21, the Distance was only to be used for short distances, and thus a maximum of 1 was returned. Through the use of a new Max vararg (setting this equal to 1), this old behaviour can be obtained.

Diameter(L) : [FldReElt] -> FldReElt

Diameter(L) : [FldComElt] -> FldReElt

    Max: ExtReElt                       Default: ∞

Given a sequence L of real of complex numbers, compute the diameter of the set defined by L, ie. the smallest distance between distinct elements of L.

Note: prior to V2.21, the Diameter was only to be used for short distances, and thus a maximum of 1 was returned. Through the use of a new Max vararg (setting this equal to 1), this old behaviour can be obtained.

Roots

Magma contains a very powerful algorithm for finding highly accurate approximations to the complex roots of a polynomial; it is based on Xavier Gourdon's implementation of Schönhage's algorithm, which we will summarize below.

Given a polynomial p = a₀ + a₁ z + ... + a_n zⁿ ∈C[z], define the norm of p, |p|, by

|p| = |a₀| + |a₁| + ... + |a_n|.

Schönhage's algorithm (given in his technical report of 1982 [Sch82]) takes as input a univariate polynomial p in C[z] and a positive real number ε, and finds linear factors L_j = u_j z - v_j (j = 1, ..., n = deg(p)) such that

|p - L₁ ... L_n| < ε |p|.

The parameter ε may be chosen so as to find the roots of p to within a certain ε', and this is how the function Roots described below works (when run with Schönhage's algorithm).

The algorithm uses the concept of a `splitting circle' to find polynomials F and G such that |p - FG| < ε₁ |p| for some ε₁ depending on ε.

This splitting circle method can then be applied recursively to F and G until we have only linear factors, as required.

The splitting circle method works as follows. For the purposes of this discussion assume that p is monic. Suppose we know a circle Γ such that, for some integer k with 0 < k < n, there are k roots of p (say u₁, ..., u_k) which lie inside Γ, and the other n - k roots (u_{k + 1}, ..., u_n) lie outside Γ. Note that the circle Γ is chosen so that the roots of p are not too close to it. Then we can write p = FG, where F = (z - u₁) ... (z - u_k) and G = (z - u_{k + 1}) ... (z - u_n). Through shifts and scalings, we may assume that Γ = {c∈C: | z| = 1}.

For m in {1, ..., k}, let s_m denote the m-th power sum of the roots of p which lie inside the splitting circle. That is, s_m = u₁^m + ... + u_k^m. The residue theorem can then be used to calculate s_m (1 ≤m ≤k): s_m = (1 /2π i) int_Γ z^m (p'(z) /p(z)) dz. where the integration can be computed to the required precision by the discrete sum s_m approx (1 /N) ∑_j=0^{N - 1} (p'(ω^j) /p(ω^j)) ω^{(m + 1)j}. for a large enough integer N, where ω = exp(2π i/N).

The coefficients of the polynomial F can then be computed from the Newton sums s_m (1 ≤m ≤k) using the classical Newton formulae. Then set G = p/F.

The integer N above needed to get F and G to the required precision can be quite large. It is more efficient to use a smaller value of N to give an approximation F₀ of F, and then use the following refining technique.

Define G₀ (an approximation of G) by p = F₀ G₀ + r, where deg(r) < deg(F₀). We want polynomials f and g such that F₁ = F₀ + f and G₁ = G₀ + g are better approximations of F and G. Now p - F₁ G₁ = p - F₀ G₀ - f G₀ - g F₀ - fg. Hence choosing f and g such that p - F₀ G₀ = f G₀ + g F₀ will lead to a second order error.

The Euclidean algorithm could be used to find f and g, but this is numerically unstable. It suffices to find polynomials H (called the auxiliary polynomial) and L such that 1 = H G₀ + L F₀, where deg(H)<deg(F₀) and deg(L)<deg(G₀).

The polynomial H can be calculated using the formula H(z) = (1/2π i) int_Γ (1 /(F₀ G₀)(t)) (F₀(z) - F₀(t)/z - t) dt. Again, rather than computing the integral to the required precision directly, we find only an approximation H₀ and then refine it using Newton iteration: H_{m + 1} ≡ H_m (2 - H_m G₀) mod (F₀). Assuming that |H - H₀| is small, the sequence (H_m) converges quadratically to H.

Once H is known to a large enough precision, f can be computed by f ≡ H (p - F₀ G₀) mod (F₀).

This gives us the new approximation F₁ of F, and G₁ is computed by division of p by F₁. We repeat this process until |p - FG| < ε₁ |p| and we are done.

The problem remains to find the splitting circle.

This relies mainly on the computation of the moduli of the roots of p. Let r₁(p) ≤r₂(p) ≤ ... ≤r_n(p) denote the moduli of the roots of p in ascending order. For each k, the computation of r_k(p) with a small number of digits can be achieved in a reliable way using the Graeffe process. The Graeffe process is a root squaring step transforming any given polynomial p into a polynomial q of the same degree whose roots are the square of the roots of p.

By the use of a suitable shift, we may assume that the sum of the roots of p is zero. If p(0)=0, then we have found a factorization p approx FG with F=z and G=p/z. If not, then the computation of the maximum root modulus r_n(p) allows us to scale p so that its maximum root modulus is now close to 1. For j = 0, 1, 2, 3, set q_j(z)=p(z + 2 i^j). Then amongst these four polynomials there exists q such that (r_n(q) /r₁(q)) = exp(Δ), with Δ > 0.3. A dichotomic process from the computation of some r_j(q) can then be applied to find k (1 ≤k ≤n - 1) such that

(r_{k + 1}(q) /r_k(q)) > exp((Δ /n - 1)).

Then the circle {c: |c|=Sqrt(r_k(q) r_{k + 1)(q)}} is a suitable splitting circle, with the roots not too close to it.

The function RootsNonExact (below) is more suitable for non-exact polynomials.

Roots(p) : RngUPolElt -> [ <FldComElt, RngIntElt> ]

    Al: MonStgElt                       Default: "Schonhage"

    Digits: RngIntElt                   Default: 

Given a univariate polynomial p over a real or complex field, this returns a sequence of complex approximations to the roots of p. The elements of this sequence are of the form <r, m>, where r is a root and m its multiplicity.

The algorithm used to find the roots of p may be specified by using the optional argument Al. This must be one of "Schonhage" (which is the default), "Laguerre", "NewtonRaphson" or "Combination" (a combination of Laguerre and Newton-Raphson). When using the (default) Schönhage algorithm, the roots given are correct to within an absolute error of 10 ^{- d}, where d is the value of Digits. This algorithm gives correct results in all cases. When using the other algorithms (for which correct answers are not guaranteed in all cases), the results are found with Digits significant figures. The default value for Digits is the current precision of the free real field.

Pari is used here for complex polynomials.

Warning: Beware of the problems of floating point numbers. Because real numbers are stored in the computer with finite precision, you may not be finding the roots of the polynomial you want. If you know the polynomial exactly, you should enter it with exact (that is, integer or rational) coefficients. This is illustrated in the following example.

> P<z> := PolynomialRing(ComplexField());
> p := (z-1.1)^6;
> p;
z^6 - 6.60000000000000000000000000001*z^5 +
    18.1500000000000000000000000000*z^4 -
    26.6200000000000000000000000000*z^3 +
    21.9615000000000000000000000000*z^2 -
    9.66306000000000000000000000003*z +
    1.77156100000000000000000000001
> R := Roots(p);
> R;
[ <1.10001330596590605421651999857, 1>,
<1.10000665289430860969298668917 +
1.15233044958179825651486651257E-5*i, 1>,
<1.10000665289430860969298668917 -
1.15233044958179825651486651257E-5*i, 1>,
<1.09998669421138030521834731234, 1>,
<1.09999334701704821058957965537 +
1.15231509613269858004669167711E-5*i, 1>,
<1.09999334701704821058957965537 -
1.15231509613269858004669167711E-5*i, 1> ]
> P<x> := PolynomialRing(Rationals());
> q := (x-11/10)^6;
> Roots(q);
[ <11/10, 6> ]

> P<z> := PolynomialRing(ComplexField());
> p := (z-1.1)^6;
> R, E := RootsNonExact(p);
> R;   
[ 1.10001483296913451410370191006 -
8.56404454142796527111307103383E-6*i,
1.10001483296913451410370191006 +
8.56404454142796527111304767692E-6*i,
1.09998516742199321904393975959 +
8.56336708832137724325720239457E-6*i,
1.09999999960887226685235833026 +
1.71274116266516824685125851069E-5*i,
1.09998516742199321904393771623 -
8.56336708832137723945115457519E-6*i,
1.09999999960887226685236037365 -
1.71274116266516824723187272572E-5*i ]
> E;
[ 0.00482314415421569719910621643066,
0.00482314415421569719910621643066,
0.00482301408919738605618476867676,
0.00482307911261159460991621017456,
0.00482301408919738605618476867676,
0.00482307911261159460991621017456 ]

    Bound: RngIntElt                    Default: -1

    Al: MonStgElt                       Default: "LLL"

    Delta: FldReElt                     Default: 0.75

    Al: MonStgElt                       Default: "LLL"

    Precision: RngIntElt                Default: 

    Delta: FldReElt                     Default: 0.75

    ExactDegree: BoolElt                Default: false

    Squarefree: BoolElt                 Default: true

> C<i> := ComplexField(1000);
> P<X> := PolynomialRing(Integers());
> P_C<Z> := PolynomialRing(C);
> p0 := 5*X^12 + 64*X^11 + 51*X^10 + 38*X^9 - 60*X^8 - 56*X^7 
>     - 51*X^6 - 14*X^5 - X^4 + 12*X^3 + 8*X^2 + 4*X;
> R := Roots (P_C!p0);
> root := R[1][1];
> res0 := MinimalPolynomial (root, 12, 100);
> p1 := P!(p0/res0);
> //
> R := Roots (P_C!p1);
> root := R[1][1];
> res1 := MinimalPolynomial (root, 11, 100);
> p2 := P!(p1/res1);
> //
> R := Roots (P_C!p2);
> root := R[1][1];
> res2 := MinimalPolynomial (root, 5, 100);
> p3 := P!(p2/res2);
> //
> R := Roots (P_C!p3);
> root := R[1][1];
> res3 := MinimalPolynomial (root, 4, 100);
> //
> assert p0 eq res0*res1*res2*res3;
> res0, res1, res2, res3;
X
X^6 + 13*X^5 + 13*X^4 + 13*X^3 - 4*X^2 - 4*X - 4
X - 1
5*X^4 + 4*X^3 + 3*X^2 + 2*X + 1