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.
The (in)equality test on real numbers of fixed precision only test for equality up to the given precision. In the free real field eqality is tested by comparing all available digits. Equality testing on complex numbers is done by testing the real and imaginary parts.
The comparison functions gt, ge, lt, le are not defined for complex numbers.
True if and only if the real or complex number c is a rational integer.
True if the complex number c is real, false otherwise. This checks whether the decimal digits of the imaginary part of c are 0 up to the precision of the complex parent 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).
If the first argument a for the exponentiation ^ is a complex number of fixed precision, then the exponent b must be an integer; if a is a fixed precision real however, then rational and real exponents are allowed as well. For free elements arbitrary real and complex exponents are allowed.
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.
MantissaExponent(r) : FldReElt -> FldReElt, RngIntElt
Given a real number r, this function returns a real number m (the mantissa of r) and an integer e (the exponent of r) such that 1 <= m< 10 and r=m x 10^e.
Given a complex number c, return the modulus m >= 0 and the argument a (with -pi <= a <= pi) of c as real numbers to the same precision as c.
Given free or fixed precision real numbers m and a, construct the complex number me^(ia). For fixed precision real arguments, the result will have the larger 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 -pi and pi) of c.
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 y 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. 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 Floor(r) if r is positive, and return - Floor( - 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 of fixed precision, return the precision p to which calculations are performed in C.
Given a free real or complex number s, this returns the relative precision to which s is accurately known; this is the number of `correct' decimal digits known for s. If s is known exactly (because it is an integer or rational number), -1 is returned to reflect infinite precision. For complex numbers the minimum of the relative precisions of real and imaginary parts is returned.
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 Catalan constant is not currently available in the free field.
Catalan(R) : FldRe -> FldReElt
The value of Catalan's constant computed to the accuracy associated with the real or complex field D of fixed precision. Catalan's constant is the sum from k equals 0 to infinity of ( - 1)^k by (2k + 1)^(-2).
The value of Euler's constantcomputed to the accuracy associated with R if R is a field of fixed precision, and to the default precision if R is the free real field. To allow the calculation of the constant to a precision different from the default in the latter case, a version of this function with the precision as a second argument is available.
The value of pi computed to the accuracy associated with R if R is a field of fixed precision, and to the default precision if R is the free real field. To allow the calculation of the constant to a precision different from the default in the latter case, a version of this function with the precision as a second argument is available.
The absolute value of the real number r; here r may be either free or of some fixed precision, and the result will be in the same field as r.
Return one of the integer values +1, 0, -1 depending upon whether the real number r is positive, zero or negative, respectively. Here r may be either free or of some fixed precision.
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 x^2 - y^2, while for elements of real domains it just returns the absolute value. The result lies in the same field as the argument, which may be free or of fixed precision.
Given a real number r of fixed precision 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. Here c may be free or of some fixed precision; if c is real of fixed precision it must be non-negative.
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_0 + a_1 z + ... + a_n z^n in C[z], define the norm of p, |p|, by
|p| = |a_0| + |a_1| + ... + |a_n|.
Schönhage's algorithm (given in his technical report of 1982 [A. Schönhage, The fundamental theorem of algebra in terms of computational complexity, Technical Report, Univ. Tübingen, 1982.]) takes as input a univariate polynomial p in C[z] and a positive real number varepsilon, and finds linear factors L_j = u_j z - v_j (j = 1, ..., n = deg(p)) such that
|p - L_1 ... L_n| < varepsilon |p|.
The parameter varepsilon may be chosen so as to find the roots of p to within a certain varepsilon', 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| < varepsilon_1 |p| for some varepsilon_1 depending on varepsilon.
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 Gamma such that, for some integer k with 0 < k < n, there are k roots of p (say u_1, ..., u_k) which lie inside Gamma, and the other n - k roots (u_(k + 1), ..., u_n) lie outside Gamma. Note that the circle Gamma is chosen so that the roots of p are not too close to it. Then we can write p = FG, where F = (z - u_1) ... (z - u_k) and G = (z - u_(k + 1)) ... (z - u_n). Through shifts and scalings, we may assume that Gamma = {c in 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_1^m + ... + u_k^m. The residue theorem can then be used to calculate s_m (1 <= m <= k): s_m = (1 /2pi i) int_(Gamma) 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) sum_(j=0)^(N - 1) (p'( omega ^j) /p( omega ^j)) omega ^((m + 1)j). for a large enough integer N, where omega = exp(2pi 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_0 of F, and then use the following refining technique.
Define G_0 (an approximation of G) by p = F_0 G_0 + r, where deg(r) < deg(F_0). We want polynomials f and g such that F_1 = F_0 + f and G_1 = G_0 + g are better approximations of F and G. Now p - F_1 G_1 = p - F_0 G_0 - f G_0 - g F_0 - fg. Hence choosing f and g such that p - F_0 G_0 = f G_0 + g F_0 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 auxilliary polynomial) and L such that 1 = H G_0 + L F_0, where deg(H)<deg(F_0) and deg(L)<deg(G_0).
The polynomial H can be calculated using the formula H(z) = (1/2pi i) int_(Gamma) (1 /(F_0 G_0)(t)) (F_0(z) - F_0(t)/z - t) dt. Again, rather than computing the integral to the required precision directly, we find only an approximation H_0 and then refine it using Newton iteration: H_(m + 1) = H_m (2 - H_m G_0) pmod(F_0). Assuming that |H - H_0| 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_0 G_0) pmod(F_0).
This gives us the new approximation F_1 of F, and G_1 is computed by division of p by F_1. We repeat this process until |p - FG| < varepsilon_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 r_1(p) <= r_2(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_1(q)) = exp(Delta), with Delta > 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((Delta /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.When working in the free real or complex field, 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.
> P<z> := PolynomialRing(ComplexField());
> p := (z-1.1)^6;
> p;
z^6 - 6.59999999999999999999999999987*z^5 +
18.149999999999999999999999998707*z^4 -
26.61999999999999999999999999767*z^3 +
21.96149999999999999999999999727*z^2 - 9.66305999999999999999999999861*z
+ 1.77156099999999999999999999973
> R := Roots(p);
> R;
[ <1.1000293845730291258969218315631755913348, 1>,
<1.1000146918029343139114323477003295381159 +
0.0000254477867168741508816635851578423578367*i, 1>,
<1.1000146918029343139114323477003295381159 -
0.0000254477867168741508816635851578423578367*i, 1>,
<1.0999706163941115999599722351289318378239, 1>,
<1.0999853077134953231601206188930322297719 +
0.0000254469491484368882655714786543303037228*i, 1>,
<1.0999853077134953231601206188930322297719 -
0.0000254469491484368882655714786543303037228*i, 1> ]
> q := (z-11/10)^6;
> Roots(q);
[ <1.0999999999999999999999999999999999999952, 6> ]
Digits: RngIntElt Default:
Given a polynomial p of degree n defined over the free real or complex field, returns a sequence [v_1, ..., v_n] of (free) complex numbers such that |p - a(z - v_1) ... (z - v_n)| < 10^(-d) |p|, where a is the leading coefficient of p, and d is the value of Digits. The default value for the optional argument Digits is the current precision of the free real field.A second sequence [e_1, ..., e_n] 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 - u_1) ... (z - u_n) with |v_i - u_i| < e_i. 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 parameter Digits will decrease the errors on the `roots'. However, to be meaningful, any non-exact coefficients of p must be given with at least Digits decimal places.
> P<z> := PolynomialRing(ComplexField()); > p := (z-1.1)^6; > R, E := RootsNonExact(p); > R; [ 1.1000293845730287249289598261713699697597 - 1.4992730221 E-34*i, 1.1000146918029341134742131109147227252400 + 0.0000254477867165269024409404759949678919020*i, 1.1000146918029341134742131109147225301256 - 0.0000254477867165269024409404759948179640380*i, 1.0999706163941120008344123617543512327655 + 1.3323163399 E-34*i, 1.0999853077134955236441007950618321355285 + 0.0000254469491480897208171711146641455283797*i, 1.0999853077134955236441007950618323714326 - 0.0000254469491480897208171711146642787575313*i ] > E; [ 0.00011963070943, 0.00011962911178, 0.00011962911178, 0.00011962431846, 0.00011962591606, 0.00011962591606 ]
The following functions use the continued fraction expansion of
real numbers to get Diophantine approximations. They were obtained
from corresponding Pari implementations.
ContinuedFraction(r) : FldPrElt -> [ RngIntElt ]
Bound: RngIntElt Default: -1
Given an element r from a free 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 to s_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 free real field and a positive integer n, this function determines the best rational approximation to r with denominator not exceeding n.
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(p_n&p_(n - 1)cr q_n&q_(n - 1)cr); the quotients p_(n - 1)/q_(n - 1) and p_(n)/q_(n) form the last two convergents for r as provided by s.
Al: MonStgElt Default: "Hastad"
Precision: RngIntElt Default:
Given a sequence q or a vector v with entries from the free real or complex field, return an integer sequence or vector forming the coefficients for a (small) linear dependency among the entries. The algorithm used may be specified by the optional parameter Al. The default is "Hastad", which uses a variation of the LLL algorithm due to Hastad, Lagarias and Schnorr; the alternative is "LLL", which uses a straight LLL algorithm. The second optional parameter "Precision" may be used to specify the precision used in the computation: the default is the current precision of the free real field. Note that to be meaningful this precision must be no larger than the precision of the input.
Al: MonStgElt Default: "Hastad"
Precision: RngIntElt Default:
Given an element r from the free 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.[Next] [Prev] [Right] [Left] [Up] [Index] [Root]