Arbitrary precision real balls
Arbitrary precision real ball arithmetic is supplied by Arb which provides a ball representation which tracks error bounds rigorously. Real numbers are represented in mid-rad interval form $[m \pm r] = [m-r, m+r]$.
The types of real balls in Nemo are given in the following table, along with the libraries that provide them and the associated types of the parent objects.
Library | Field | Element type | Parent type |
---|---|---|---|
Arb | $\mathbb{R}$ (balls) | RealFieldElem | RealField |
The real field types belong to the Field
abstract type and the types of elements in this field, i.e. balls in this case, belong to the FieldElem
abstract type.
Real ball functionality
Real balls in Nemo provide all the field functionality described in AbstractAlgebra:
https://nemocas.github.io/AbstractAlgebra.jl/stable/field
Below, we document the additional functionality provided for real balls.
Precision management
Precision for ball arithmetic and creation of elements can be controlled using the functions:
precision
— Methodprecision(::Type{Balls})
Return the precision for ball arithmetic.
Examples
julia> set_precision!(Balls, 200); precision(Balls)
200
set_precision!
— Methodset_precision!(::Type{Balls}, n::Int)
Set the precision for all ball arithmetic to be n
.
Examples
julia> const_pi(RealField())
[3.141592653589793239 +/- 5.96e-19]
julia> set_precision!(Balls, 200); const_pi(RealField())
[3.14159265358979323846264338327950288419716939937510582097494 +/- 5.73e-60]
set_precision!
— Methodset_precision!(f, ::Type{Balls}, n::Int)
Change ball arithmetic precision to n
for the duration of f
..
Examples
julia> set_precision!(Balls, 4) do
const_pi(RealField())
end
[3e+0 +/- 0.376]
julia> set_precision!(Balls, 200) do
const_pi(RealField())
end
[3.1415926535897932385 +/- 3.74e-20]
This functions are not thread-safe.
Constructors
In order to construct real balls in Nemo, one must first construct the Arb real field itself. This is accomplished with the following constructor.
RealField()
Here is an example of creating the real field and using the resulting parent object to coerce values into the resulting field.
Examples
julia> RR = RealField()
Real field
julia> a = RR("0.25")
0.25000000000000000000
julia> b = RR("0.1 +/- 0.001")
[0.1 +/- 1.01e-3]
julia> c = RR(0.5)
0.50000000000000000000
julia> d = RR(12)
12.000000000000000000
Note that whilst one can coerce double precision floating point values into an Arb real field, unless those values can be represented exactly in double precision the resulting ball can't be any more precise than the double precision supplied.
If instead, values can be represented precisely using decimal arithmetic then one can supply them to Arb using a string. In this case, Arb will store them to the precision specified when creating the Arb field.
If the values can be stored precisely as a binary floating point number, Arb will store the values exactly. See the function is_exact
below for more information.
Real ball constructors
Using coercion into the real field, new elements can be created.
Examples
julia> c = RR(1)
1.0000000000000000000
julia> d = RR(1//2)
0.50000000000000000000
Note that for the construction, also the precision can be supplied:
julia> c = RR(1//3, precision=100)
[0.33333333333333333333 +/- 3.34e-21]
julia> d = RR(1//3, precision=4)
[0.3 +/- 0.0438]
Conversions
julia> convert(Float64, RR(1//3))
0.3333333333333333
Basic manipulation
is_nonzero
— Methodis_nonzero(x::RealFieldElem)
Return true
if $x$ is certainly not equal to zero, otherwise return false
.
isfinite
— Methodisfinite(x::RealFieldElem)
Return true
if $x$ is finite, i.e. having finite midpoint and radius, otherwise return false
.
is_exact
— Methodis_exact(x::RealFieldElem)
Return true
if $x$ is exact, i.e. has zero radius, otherwise return false
.
isinteger
— Methodisinteger(x::RealFieldElem)
Return true
if $x$ is an exact integer, otherwise return false
.
is_positive
— Methodis_positive(x::RealFieldElem)
Return true
if $x$ is certainly positive, otherwise return false
.
is_nonnegative
— Methodis_nonnegative(x::RealFieldElem)
Return true
if $x$ is certainly non-negative, otherwise return false
.
is_negative
— Methodis_negative(x::RealFieldElem)
Return true
if $x$ is certainly negative, otherwise return false
.
is_nonpositive
— Methodis_nonpositive(x::RealFieldElem)
Return true
if $x$ is certainly nonpositive, otherwise return false
.
midpoint
— Methodmidpoint(x::RealFieldElem)
Return the midpoint of the ball $x$ as an Arb ball.
radius
— Methodradius(x::RealFieldElem)
Return the radius of the ball $x$ as an Arb ball.
accuracy_bits
— Methodaccuracy_bits(x::RealFieldElem)
Return the relative accuracy of $x$ measured in bits, capped between typemax(Int)
and -typemax(Int)
.
Examples
julia> RR = RealField()
Real field
julia> a = RR("1.2 +/- 0.001")
[1.20 +/- 1.01e-3]
julia> b = RR(3)
3.0000000000000000000
julia> is_positive(a)
true
julia> isfinite(b)
true
julia> isinteger(b)
true
julia> is_negative(a)
false
julia> c = radius(a)
[0.0010000000038417056203 +/- 1.12e-23]
julia> d = midpoint(b)
3.0000000000000000000
julia> f = accuracy_bits(a)
9
Printing
Printing real balls can at first sight be confusing. Lets look at the following example:
julia> a = RR(1)
1.0000000000000000000
julia> b = RR(2)
2.0000000000000000000
julia> c = RR(12)
12.000000000000000000
julia> x = ball(a, b)
[+/- 3.01]
julia> y = ball(c, b)
[1e+1 +/- 4.01]
julia> mid = midpoint(x)
1.0000000000000000000
julia> rad = radius(x)
[2.0000000037252902985 +/- 3.81e-20]
julia> print(x, "\n", y, "\n", mid, "\n", rad)
[+/- 3.01]
[1e+1 +/- 4.01]
1.0000000000000000000
[2.0000000037252902985 +/- 3.81e-20]
The first reason that c
is not printed as [1 +/- 2]
is that the midpoint does not have a greater exponent than the radius in its scientific notation. For similar reasons y
is not printed as [12 +/- 2]
.
The second reason is that we get an additional error term after our addition. As we see, radius(c)
is not equal to $2$, which when printed rounds it up to a reasonable decimal place. This is because real balls keep track of rounding errors of basic arithmetic.
Containment
It is often necessary to determine whether a given exact value or ball is contained in a given real ball or whether two balls overlap. The following functions are provided for this purpose.
overlaps
— Methodoverlaps(x::RealFieldElem, y::RealFieldElem)
Returns true
if any part of the ball $x$ overlaps any part of the ball $y$, otherwise return false
.
contains
— Methodcontains(x::RealFieldElem, y::RealFieldElem)
Returns true
if the ball $x$ contains the ball $y$, otherwise return false
.
contains
— Methodcontains(x::RealFieldElem, y::Integer)
Returns true
if the ball $x$ contains the given integer value, otherwise return false
.
contains
— Methodcontains(x::RealFieldElem, y::ZZRingElem)
Returns true
if the ball $x$ contains the given integer value, otherwise return false
.
contains
— Methodcontains(x::RealFieldElem, y::QQFieldElem)
Returns true
if the ball $x$ contains the given rational value, otherwise return false
.
contains
— Methodcontains(x::RealFieldElem, y::Rational{T}) where {T <: Integer}
Returns true
if the ball $x$ contains the given rational value, otherwise return false
.
contains
— Methodcontains(x::RealFieldElem, y::BigFloat)
Returns true
if the ball $x$ contains the given floating point value, otherwise return false
.
The following functions are also provided for determining if a ball intersects a certain part of the real number line.
contains_zero
— Methodcontains_zero(x::RealFieldElem)
Returns true
if the ball $x$ contains zero, otherwise return false
.
contains_negative
— Methodcontains_negative(x::RealFieldElem)
Returns true
if the ball $x$ contains any negative value, otherwise return false
.
contains_positive
— Methodcontains_positive(x::RealFieldElem)
Returns true
if the ball $x$ contains any positive value, otherwise return false
.
contains_nonnegative
— Methodcontains_nonnegative(x::RealFieldElem)
Returns true
if the ball $x$ contains any non-negative value, otherwise return false
.
contains_nonpositive
— Methodcontains_nonpositive(x::RealFieldElem)
Returns true
if the ball $x$ contains any nonpositive value, otherwise return false
.
Examples
julia> x = RR("1 +/- 0.001")
[1.00 +/- 1.01e-3]
julia> y = RR("3")
3.0000000000000000000
julia> overlaps(x, y)
false
julia> contains(x, y)
false
julia> contains(y, 3)
true
julia> contains(x, ZZ(1)//2)
false
julia> contains_zero(x)
false
julia> contains_positive(y)
true
Comparison
Nemo provides a full range of comparison operations for Arb balls. Note that a ball is considered less than another ball if every value in the first ball is less than every value in the second ball, etc.
In addition to the standard comparison operators, we introduce an exact equality. This is distinct from arithmetic equality implemented by ==
, which merely compares up to the minimum of the precisions of its operands.
isequal
— Methodisequal(x::RealFieldElem, y::RealFieldElem)
Return true
if the balls $x$ and $y$ are precisely equal, i.e. have the same midpoints and radii.
We also provide a full range of ad hoc comparison operators. These are implemented directly in Julia, but we document them as though isless
and ==
were provided.
Function |
---|
==(x::RealFieldElem, y::Integer) |
==(x::Integer, y::RealFieldElem) |
==(x::RealFieldElem, y::ZZRingElem) |
==(x::ZZRingElem, y::RealFieldElem) |
==(x::RealFieldElem, y::Float64) |
==(x::Float64, y::RealFieldElem) |
isless(x::RealFieldElem, y::Integer) |
isless(x::Integer, y::RealFieldElem) |
isless(x::RealFieldElem, y::ZZRingElem) |
isless(x::ZZRingElem, y::RealFieldElem) |
isless(x::RealFieldElem, y::Float64) |
isless(x::Float64, y::RealFieldElem) |
isless(x::RealFieldElem, y::BigFloat) |
isless(x::BigFloat, y::RealFieldElem) |
isless(x::RealFieldElem, y::QQFieldElem) |
isless(x::QQFieldElem, y::RealFieldElem) |
Examples
julia> x = RR("1 +/- 0.001")
[1.00 +/- 1.01e-3]
julia> y = RR("3")
3.0000000000000000000
julia> z = RR("4")
4.0000000000000000000
julia> isequal(x, deepcopy(x))
true
julia> x == 3
false
julia> ZZ(3) < z
true
julia> x != 1.23
true
Absolute value
Examples
julia> x = RR("-1 +/- 0.001")
[-1.00 +/- 1.01e-3]
julia> a = abs(x)
[1.00 +/- 1.01e-3]
Shifting
Examples
julia> x = RR("-3 +/- 0.001")
[-3.00 +/- 1.01e-3]
julia> a = ldexp(x, 23)
[-2.52e+7 +/- 4.26e+4]
julia> b = ldexp(x, -ZZ(15))
[-9.16e-5 +/- 7.78e-8]
Miscellaneous operations
add_error!
— Methodadd_error!(x::RealFieldElem, y::RealFieldElem)
Adds the absolute values of the midpoint and radius of $y$ to the radius of $x$.
trim
— Methodtrim(x::RealFieldElem)
Return an RealFieldElem
interval containing $x$ but which may be more economical, by rounding off insignificant bits from the midpoint.
unique_integer
— Methodunique_integer(x::RealFieldElem)
Return a pair where the first value is a boolean and the second is an ZZRingElem
integer. The boolean indicates whether the interval $x$ contains a unique integer. If this is the case, the second return value is set to this unique integer.
setunion
— Methodsetunion(x::RealFieldElem, y::RealFieldElem)
Return an ArbFieldElem
containing the union of the intervals represented by $x$ and $y$.
Examples
julia> x = RR("-3 +/- 0.001")
[-3.00 +/- 1.01e-3]
julia> y = RR("2 +/- 0.5")
[2e+0 +/- 0.501]
julia> a = trim(x)
[-3.00 +/- 1.01e-3]
julia> b, c = unique_integer(x)
(true, -3)
julia> d = setunion(x, y)
[+/- 3.01]
Constants
const_pi
— Methodconst_pi(r::RealField)
Return $\pi = 3.14159\ldots$ as an element of $r$.
const_e
— Methodconst_e(r::RealField)
Return $e = 2.71828\ldots$ as an element of $r$.
const_log2
— Methodconst_log2(r::RealField)
Return $\log(2) = 0.69314\ldots$ as an element of $r$.
const_log10
— Methodconst_log10(r::RealField)
Return $\log(10) = 2.302585\ldots$ as an element of $r$.
const_euler
— Methodconst_euler(r::RealField)
Return Euler's constant $\gamma = 0.577215\ldots$ as an element of $r$.
const_catalan
— Methodconst_catalan(r::RealField)
Return Catalan's constant $C = 0.915965\ldots$ as an element of $r$.
const_khinchin
— Methodconst_khinchin(r::RealField)
Return Khinchin's constant $K = 2.685452\ldots$ as an element of $r$.
const_glaisher
— Methodconst_glaisher(r::RealField)
Return Glaisher's constant $A = 1.282427\ldots$ as an element of $r$.
Examples
julia> a = const_pi(RR)
[3.141592653589793239 +/- 5.96e-19]
julia> b = const_e(RR)
[2.718281828459045235 +/- 4.29e-19]
julia> c = const_euler(RR)
[0.5772156649015328606 +/- 4.35e-20]
julia> d = const_glaisher(RR)
[1.282427129100622637 +/- 3.01e-19]
Mathematical and special functions
rsqrt
— Methodrsqrt(x::RealFieldElem)
Return the reciprocal of the square root of $x$, i.e. $1/\sqrt{x}$.
sqrt1pm1
— Methodsqrt1pm1(x::RealFieldElem)
Return $\sqrt{1+x}-1$, evaluated accurately for small $x$.
sqrtpos
— Methodsqrtpos(x::RealFieldElem)
Return the sqrt root of $x$, assuming that $x$ represents a non-negative number. Thus any negative number in the input interval is discarded.
gamma
— Methodgamma(x::RealFieldElem)
Return the Gamma function evaluated at $x$.
lgamma
— Methodlgamma(x::RealFieldElem)
Return the logarithm of the Gamma function evaluated at $x$.
rgamma
— Methodrgamma(x::RealFieldElem)
Return the reciprocal of the Gamma function evaluated at $x$.
digamma
— Methoddigamma(x::RealFieldElem)
Return the logarithmic derivative of the gamma function evaluated at $x$, i.e. $\psi(x)$.
gamma
— Methodgamma(s::RealFieldElem, x::RealFieldElem)
Return the upper incomplete gamma function $\Gamma(s,x)$.
gamma_regularized
— Methodgamma_regularized(s::RealFieldElem, x::RealFieldElem)
Return the regularized upper incomplete gamma function $\Gamma(s,x) / \Gamma(s)$.
gamma_lower
— Methodgamma_lower(s::RealFieldElem, x::RealFieldElem)
Return the lower incomplete gamma function $\gamma(s,x) / \Gamma(s)$.
gamma_lower_regularized
— Methodgamma_lower_regularized(s::RealFieldElem, x::RealFieldElem)
Return the regularized lower incomplete gamma function $\gamma(s,x) / \Gamma(s)$.
zeta
— Methodzeta(x::RealFieldElem)
Return the Riemann zeta function evaluated at $x$.
atan2
— Methodatan2(y::RealFieldElem, x::RealFieldElem)
Return $\operatorname{atan2}(y,x) = \arg(x+yi)$. Same as atan(y, x)
.
agm
— Methodagm(x::RealFieldElem, y::RealFieldElem)
Return the arithmetic-geometric mean of $x$ and $y$
zeta
— Methodzeta(s::RealFieldElem, a::RealFieldElem)
Return the Hurwitz zeta function $\zeta(s,a)$.
root
— Methodroot(x::RealFieldElem, n::Int)
Return the $n$-th root of $x$. We require $x \geq 0$.
factorial
— Methodfactorial(x::RealFieldElem)
Return the factorial of $x$.
factorial
— Methodfactorial(n::Int, r::RealField)
Return the factorial of $n$ in the given Arb field.
binomial
— Methodbinomial(x::RealFieldElem, n::UInt)
Return the binomial coefficient ${x \choose n}$.
binomial
— Methodbinomial(n::UInt, k::UInt, r::RealField)
Return the binomial coefficient ${n \choose k}$ in the given Arb field.
fibonacci
— Methodfibonacci(n::ZZRingElem, r::RealField)
Return the $n$-th Fibonacci number in the given Arb field.
fibonacci
— Methodfibonacci(n::Int, r::RealField)
Return the $n$-th Fibonacci number in the given Arb field.
gamma
— Methodgamma(x::ZZRingElem, r::RealField)
Return the Gamma function evaluated at $x$ in the given Arb field.
gamma
— Methodgamma(x::QQFieldElem, r::RealField)
Return the Gamma function evaluated at $x$ in the given Arb field.
zeta
— Methodzeta(n::Int, r::RealField)
Return the Riemann zeta function $\zeta(n)$ as an element of the given Arb field.
bernoulli
— Methodbernoulli(n::Int, r::RealField)
Return the $n$-th Bernoulli number as an element of the given Arb field.
rising_factorial
— Methodrising_factorial(x::RealFieldElem, n::Int)
Return the rising factorial $x(x + 1)\ldots (x + n - 1)$ as an Arb.
rising_factorial(x::RingElement, n::Integer)
Return the rising factorial of $x$, i.e. $x(x + 1)(x + 2)\cdots (x + n - 1)$. If $n < 0$ we throw a DomainError()
.
Examples
julia> R, x = ZZ[:x];
julia> rising_factorial(x, 1)
x
julia> rising_factorial(x, 2)
x^2 + x
julia> rising_factorial(4, 2)
20
rising_factorial
— Methodrising_factorial(x::QQFieldElem, n::Int, r::RealField)
Return the rising factorial $x(x + 1)\ldots (x + n - 1)$ as an element of the given Arb field.
rising_factorial2
— Methodrising_factorial2(x::RealFieldElem, n::Int)
Return a tuple containing the rising factorial $x(x + 1)\ldots (x + n - 1)$ and its derivative.
rising_factorial2(x::RingElement, n::Integer)
Return a tuple containing the rising factorial $x(x + 1)\cdots (x + n - 1)$ and its derivative. If $n < 0$ we throw a DomainError()
.
Examples
julia> R, x = ZZ[:x];
julia> rising_factorial2(x, 1)
(x, 1)
julia> rising_factorial2(x, 2)
(x^2 + x, 2*x + 1)
julia> rising_factorial2(4, 2)
(20, 9)
polylog
— Methodpolylog(s::Union{RealFieldElem,Int}, a::RealFieldElem)
Return the polylogarithm Li$_s(a)$.
chebyshev_t
— Methodchebyshev_t(n::Int, x::RealFieldElem)
Return the value of the Chebyshev polynomial $T_n(x)$.
chebyshev_u
— Methodchebyshev_u(n::Int, x::RealFieldElem)
Return the value of the Chebyshev polynomial $U_n(x)$.
chebyshev_t2
— Methodchebyshev_t2(n::Int, x::RealFieldElem)
Return the tuple $(T_{n}(x), T_{n-1}(x))$.
chebyshev_u2
— Methodchebyshev_u2(n::Int, x::RealFieldElem)
Return the tuple $(U_{n}(x), U_{n-1}(x))$
bell
— Methodbell(n::ZZRingElem, r::RealField)
Return the Bell number $B_n$ as an element of $r$.
bell
— Methodbell(n::Int, r::RealField)
Return the Bell number $B_n$ as an element of $r$.
numpart
— Methodnumpart(n::ZZRingElem, r::RealField)
Return the number of partitions $p(n)$ as an element of $r$.
numpart
— Methodnumpart(n::Int, r::RealField)
Return the number of partitions $p(n)$ as an element of $r$.
airy_ai
— Methodairy_ai(x::RealFieldElem)
Return the Airy function $\operatorname{Ai}(x)$.
airy_ai_prime
— Methodairy_ai_prime(x::RealFieldElem)
Return the derivative of the Airy function $\operatorname{Ai}^\prime(x)$.
airy_bi
— Methodairy_bi(x::RealFieldElem)
Return the Airy function $\operatorname{Bi}(x)$.
airy_bi_prime
— Methodairy_bi_prime(x::RealFieldElem)
Return the derivative of the Airy function $\operatorname{Bi}^\prime(x)$.
Examples
julia> a = floor(exp(RR(1)))
2.0000000000000000000
julia> b = sinpi(QQ(5,6), RR)
0.50000000000000000000
julia> c = gamma(QQ(1,3), RR)
[2.678938534707747634 +/- 7.13e-19]
julia> d = bernoulli(1000, RR)
[-5.318704469415522036e+1769 +/- 6.61e+1750]
julia> f = polylog(3, RR(-10))
[-5.92106480375697 +/- 6.68e-15]
Linear dependence
lindep
— Methodlindep(A::Vector{RealFieldElem}, bits::Int)
Find a small linear combination of the entries of the array $A$ that is small (using LLL). The entries are first scaled by the given number of bits before truncating to integers for use in LLL. This function can be used to find linear dependence between a list of real numbers. The algorithm is heuristic only and returns an array of Nemo integers representing the linear combination.
Examples
julia> RR = RealField()
Real field
julia> a = RR(-0.33198902958450931620250069492231652319)
[-0.33198902958450932088 +/- 4.15e-22]
julia> V = [RR(1), a, a^2, a^3, a^4, a^5]
6-element Vector{RealFieldElem}:
1.0000000000000000000
[-0.33198902958450932088 +/- 4.15e-22]
[0.11021671576446420510 +/- 7.87e-21]
[-0.03659074051063616184 +/- 4.17e-21]
[0.012147724433904692427 +/- 4.99e-22]
[-0.004032911246472051677 +/- 6.25e-22]
julia> W = lindep(V, 20)
6-element Vector{ZZRingElem}:
1
3
0
0
0
1
simplest_rational_inside
— Method simplest_rational_inside(x::RealFieldElem)
Return the simplest fraction inside the ball $x$. A canonical fraction $a_1/b_1$ is defined to be simpler than $a_2/b_2$ iff $b_1 < b_2$ or $b_1 = b_2$ and $a_1 < a_2$.
Examples
julia> RR = RealField()
Real field
julia> simplest_rational_inside(const_pi(RR))
8717442233//2774848045
Random generation
rand
— Methodrand(r::RealField; randtype::Symbol=:urandom)
Return a random element in given Arb field.
The randtype
default is :urandom
which return an ArbFieldElem
contained in $[0,1]$.
The rest of the methods return non-uniformly distributed values in order to exercise corner cases. The option :randtest
will return a finite number, and :randtest_exact
the same but with a zero radius. The option :randtest_precise
return an ArbFieldElem
with a radius around $2^{-\mathrm{prec}}$ the magnitude of the midpoint, while :randtest_wide
return a radius that might be big relative to its midpoint. The :randtest_special
-option might return a midpoint and radius whose values are NaN
or inf
.
rand([rng=GLOBAL_RNG,] G::SymmetricGroup)
Return a random permutation from G
.
Examples
a = rand(RR)
b = rand(RR; randtype = :null_exact)
c = rand(RR; randtype = :exact)
d = rand(RR; randtype = :special)