All predicates on real or complex numbers that check whether these numbers are equal to an integer do so within the given precision of the parent field. Thus IsOne(c) for an element of a complex domain of precision 20 returns true if and only if the real part equals one and the imaginary part equals 0 up to 20 decimals.
Note that IsZero will match either the positive or negative zero in MPFR.
The (in)equality test on real numbers of only test for equality up to the given precision. Equality testing on complex numbers is done by testing the real and imaginary parts.
Again the positive and negative zero in MPFR are equal.
The comparison functions gt, ge, lt, le are not defined for complex numbers.
Returns true if and only if the real or complex number c is a rational integer.
Returns true if the complex number c is real, false otherwise. This checks whether the digits of the imaginary part of c are 0 up to the precision of the parent complex field.
The binary operations +, -, *, / allow combinations of arguments from the integers, the rationals, and real and complex fields; automatic coercion is applied where necessary (see the Introduction).
Here we list various ways to convert between integers, reals of fixed precision, complexes and their various representations, other than by the creation functions and !. See also the rounding functions in a later section.
Given a real number r, this function returns an integer m (the mantissa of r) and an integer e (the exponent of r) such that r approx m x 2e, with the size of m being the same as that of the internal representation of r.
Given a complex number c, return the modulus m≥0 and the argument a (with -π≤a≤π) of c as real numbers to the same precision as c.
Given real numbers m and a, construct the complex number meia. The result will have the smaller of the precisions of m and a; each of m and a is allowed to be an integer or rational number; if both are integral or rational then the result will have the default precision, otherwise the result will be of the same precision as the real argument.
Given a complex number c, return the real number (to the same precision) that is the argument (in radians between -π and π) of c.
> 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!
Given a complex number c, return the real number (to the same precision as c) that is the modulus of c.
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).
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).
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.
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.
The ceiling of the real number r, i.e. the smallest integer greater than or equal to r.
The floor of the real number r, i.e. the greatest integer less than or equal to r.
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.
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.
Gives a sequence of real or complex numbers, return the precision p of their parent field.
Coerces the real (r) or complex (c) number into a field of precision n.
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.
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.
The value of Euler's constant computed to the precision of R.
The value of π computed to the precision of R.
The absolute value of the real or complex number r.
Return one of the integer values +1, 0, -1 depending upon whether the real number r is positive, zero or negative, respectively.
The complex conjugate x - y i of a complex number x + y i.
The real norm of a real or complex number c; note that for complex c=x + y i this returns x2 + y2, while for elements of real domains it just returns the absolute value. The result lies in the same field as the argument.
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.
Given a real or complex number c, return the square root of r as an element of the same field to which r belongs.
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. miny∈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.
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.
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 = a0 + a1 z + ... + an zn ∈C[z], define the norm of p, |p|, by
|p| = |a0| + |a1| + ... + |an|.
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 Lj = uj z - vj (j = 1, ..., n = deg(p)) such that
|p - L1 ... Ln| < ε |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| < ε1 |p| for some ε1 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 u1, ..., uk) which lie inside Γ, and the other n - k roots (uk + 1, ..., un) 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 - u1) ... (z - uk) and G = (z - uk + 1) ... (z - un). Through shifts and scalings, we may assume that Γ = {c∈C: | z| = 1}.
For m in {1, ..., k}, let sm denote the m-th power sum of the roots of p which lie inside the splitting circle. That is, sm = u1m + ... + ukm. The residue theorem can then be used to calculate sm (1 ≤m ≤k): sm = (1 /2π i) intΓ zm (p'(z) /p(z)) dz. where the integration can be computed to the required precision by the discrete sum sm approx (1 /N) ∑j=0N - 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 sm (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 F0 of F, and then use the following refining technique.
Define G0 (an approximation of G) by p = F0 G0 + r, where deg(r) < deg(F0). We want polynomials f and g such that F1 = F0 + f and G1 = G0 + g are better approximations of F and G. Now p - F1 G1 = p - F0 G0 - f G0 - g F0 - fg. Hence choosing f and g such that p - F0 G0 = f G0 + g F0 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 G0 + L F0, where deg(H)<deg(F0) and deg(L)<deg(G0).
The polynomial H can be calculated using the formula H(z) = (1/2π i) intΓ (1 /(F0 G0)(t)) (F0(z) - F0(t)/z - t) dt. Again, rather than computing the integral to the required precision directly, we find only an approximation H0 and then refine it using Newton iteration: Hm + 1 ≡ Hm (2 - Hm G0) mod (F0). Assuming that |H - H0| is small, the sequence (Hm) converges quadratically to H.
Once H is known to a large enough precision, f can be computed by f ≡ H (p - F0 G0) mod (F0).
This gives us the new approximation F1 of F, and G1 is computed by division of p by F1. We repeat this process until |p - FG| < ε1 |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 r1(p) ≤r2(p) ≤ ... ≤rn(p) denote the moduli of the roots of p in ascending order. For each k, the computation of rk(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 rn(p) allows us to scale p so that its maximum root modulus is now close to 1. For j = 0, 1, 2, 3, set qj(z)=p(z + 2 ij). Then amongst these four polynomials there exists q such that (rn(q) /r1(q)) = exp(Δ), with Δ > 0.3. A dichotomic process from the computation of some rj(q) can then be applied to find k (1 ≤k ≤n - 1) such that
(rk + 1(q) /rk(q)) > exp((Δ /n - 1)).
Then the circle {c: |c|=Sqrt(rk(q) rk + 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.
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> ]
Given a polynomial p of degree n defined over a real or complex field, returns a sequence [v1, ..., vn] of complex numbers such that |p - a(z - v1) ... (z - vn)| < 10 - d |p|, where a is the leading coefficient of p, and d is the precision of the field.A second sequence [e1, ..., en] of (free) real numbers may also be returned. Given any polynomial hat p such that |p - hat p| < 10 - d |p|, we can write hat p = a(z - u1) ... (z - un) with |vi - ui| < ei. In some cases, such error bounds cannot be derived, because the value d of Digits is too small for the given polynomial. In these cases, this second sequence is not returned.
This function acknowledges the fact that the polynomial p may not be the exact polynomial wanted, but only an approximation (to a certain number of decimal places), and so the roots of the true polynomial can only be found to a limited number of decimal places. Increasing the precision will decrease the errors on the `roots'.
> 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 ]
Let f be a real or complex polynomial and x an approximation to a single zero of f. This function will apply the Newton-iteration to improve the accuracy of the root to the precision indicated by k.
The following functions use the continued fraction expansion of real numbers to get Diophantine approximations. They were obtained from corresponding Pari implementations.
Bound: RngIntElt Default: -1
Given an element r from a real field, return a sequence of integers s that form the partial quotient for the (regular) continued fraction expansion for r, so r is approximately equal tos_1 + 1/( s_2 + 1/( s_3 + ... + 1/(s_n)))The length n of the sequence is determined in such a way that the last significant partial quotient is obtained (determined by the precision with which r is known), unless the optional integer argument Bound is used to limit the length.
Given an element r from a real field and a positive integer n, this function determines a rational approximation to r with denominator not exceeding n. The approximation is at least as close as the best continued fraction convergent with denominator not exceeding n. Pari is used here.
Given a sequence s of n non-negative integers (forming the partial fractions of a real number r, say), this function returns a 2 x 2 matrix with integer coefficients pmatrix(pn&pn - 1
qn&qn - 1
); the quotients pn - 1/qn - 1 and pn/qn form the last two convergents for r as provided by s.
Al: MonStgElt Default: "LLL"
(Deprecated). Given a sequence q or a vector v with entries from a complex field, return an integer sequence or vector forming the coefficients for a (small) linear dependency among the entries. The algorithm now only uses "LLL" (not "Hastad"), and this function is deprecated in any case. The new version of this function is IntegerRelation.
Given a sequence q with entries from a real or complex field, return the lattice of all (small) integer linear dependencies among the entries. The precision p is used both in the relation-detection, and to define "small" dependencies.
Delta: FldReElt Default: 0.75
Given a sequence q=(q1, ..., qd) with entries from a real or complex field, along possibly with an integer N>0, return a sequence (n1, ..., nd) of integers such that: any ni has magnitude less or equal to N and the quantity |∑ni qi| is small. The routine is based on the LLL lattice reduction algorithm [LLL82] (see Chapter Reduction of Matrices and Lattices). The second return value is a measure of goodness-of-fit, and the second input value will be determined from the input precision of real/complex numbers if it is omitted.
Al: MonStgElt Default: "LLL"
Precision: RngIntElt Default:
(Deprecated). Given an element r from a real or complex field, and an integer k>0, return a univariate integer polynomial of degree at most k having r as an approximate root. The parameters here have the same usage and meaning as for LinearRelation. Pari is used here. The new version of this function is MinimalPolynomial.
Delta: FldReElt Default: 0.75
ExactDegree: BoolElt Default: false
Squarefree: BoolElt Default: true
Given an element r from a real or complex field and two integers d≥1 and N≥1, tries to compute a univariate integer polynomial of degree at most d and coefficients of absolute values no greater than N which is the minimal polynomial of a number close to r. The routine relies on the LLL lattice reduction algorithm [LLL82] (see Chapter Reduction of Matrices and Lattices). When N is not given, it is determined from the precision of r. The ExactDegree vararg requires that the output degree be exactly as stated, and Squarefree removes any repeated factors. The second return value is a measure of goodness-of-fit.
Notice that the polynomial factorisation routine Factorization (see Section Factorization) relies upon Van Hoeij's algorithm [vH02], [vH01]. Though both approaches have polynomial time complexities, Van Hoeij's algorithm is significantly faster.
> 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