On this page
Mathematics
Mathematical Operators
(x)

Unary minus operator.
+(x, y...)

Addition operator.
x+y+z+...
calls this function with all arguments, i.e.+(x, y, z, ...)
.
(x, y)

Subtraction operator.
*(x, y...)

Multiplication operator.
x*y*z*...
calls this function with all arguments, i.e.*(x, y, z, ...)
.
/(x, y)

Right division operator: multiplication of
x
by the inverse ofy
on the right. Gives floatingpoint results for integer arguments.
\(x, y)

Left division operator: multiplication of
y
by the inverse ofx
on the left. Gives floatingpoint results for integer arguments.
^(x, y)

Exponentiation operator.
.+(x, y)

Elementwise addition operator.
.(x, y)

Elementwise subtraction operator.
.*(x, y)

Elementwise multiplication operator.
./(x, y)

Elementwise right division operator.
.\(x, y)

Elementwise left division operator.
.^(x, y)

Elementwise exponentiation operator.
fma(x, y, z)

Computes
x*y+z
without rounding the intermediate resultx*y
. On some systems this is significantly more expensive thanx*y+z
.fma
is used to improve accuracy in certain algorithms. Seemuladd
.
muladd(x, y, z)

Combined multiplyadd, computes
x*y+z
in an efficient manner. This may on some systems be equivalent tox*y+z
, or tofma(x,y,z)
.muladd
is used to improve performance. Seefma
.
div(x, y)
÷(x, y)

The quotient from Euclidean division. Computes
x/y
, truncated to an integer.
fld(x, y)

Largest integer less than or equal to
x/y
.
cld(x, y)

Smallest integer larger than or equal to
x/y
.
mod(x, y)

Modulus after flooring division, returning in the range \([0,y)\), if
y
is positive, or \((y,0]\) ify
is negative.x == fld(x,y)*y + mod(x,y)
mod2pi(x)

Modulus after division by
2π
, returning in the range \([0,2π)\).This function computes a floating point representation of the modulus after division by numerically exact
2π
, and is therefore not exactly the same asmod(x,2π)
, which would compute the modulus ofx
relative to division by the floatingpoint number2π
.
rem(x, y)
%(x, y)

Remainder from Euclidean division, returning a value of the same sign as
x
, and smaller in magnitude thany
. This value is always exact.x == div(x,y)*y + rem(x,y)
divrem(x, y)

The quotient and remainder from Euclidean division. Equivalent to
(div(x,y), rem(x,y))
or(x÷y, x%y)
.
fldmod(x, y)

The floored quotient and modulus after division. Equivalent to
(fld(x,y), mod(x,y))
.
fld1(x, y)

Flooring division, returning a value consistent with
mod1(x,y)
x == fld(x,y)*y + mod(x,y) x == (fld1(x,y)1)*y + mod1(x,y)
mod1(x, y)

Modulus after flooring division, returning a value
r
such thatmod(r, y) == mod(x, y)
in the range \((0, y]\) for positivey
and in the range \([y,0)\) for negativey
.
fldmod1(x, y)

Return
(fld1(x,y), mod1(x,y))
.
//(num, den)

Divide two integers or rational numbers, giving a
Rational
result.
rationalize([T<:Integer=Int, ]x; tol::Real=eps(x))

Approximate floating point number
x
as aRational
number with components of the given integer type. The result will differ fromx
by no more thantol
. IfT
is not provided, it defaults toInt
.julia> rationalize(5.6) 28//5 julia> a = rationalize(BigInt, 10.3) 103//10 julia> typeof(num(a)) BigInt
num(x)

Numerator of the rational representation of
x
.
den(x)

Denominator of the rational representation of
x
.
<<(x, n)

Left bit shift operator,
x << n
. Forn >= 0
, the result isx
shifted left byn
bits, filling with0
s. This is equivalent tox * 2^n
. Forn < 0
, this is equivalent tox >> n
.julia> Int8(3) << 2 12 julia> bits(Int8(3)) "00000011" julia> bits(Int8(12)) "00001100"
>>(x, n)

Right bit shift operator,
x >> n
. Forn >= 0
, the result isx
shifted right byn
bits, wheren >= 0
, filling with0
s ifx >= 0
,1
s ifx < 0
, preserving the sign ofx
. This is equivalent tofld(x, 2^n)
. Forn < 0
, this is equivalent tox << n
.julia> Int8(13) >> 2 3 julia> bits(Int8(13)) "00001101" julia> bits(Int8(3)) "00000011" julia> Int8(14) >> 2 4 julia> bits(Int8(14)) "11110010" julia> bits(Int8(4)) "11111100"
>>>(x, n)

Unsigned right bit shift operator,
x >>> n
. Forn >= 0
, the result isx
shifted right byn
bits, wheren >= 0
, filling with0
s. Forn < 0
, this is equivalent tox << n
.For
Unsigned
integer types, this is equivalent to>>()
. ForSigned
integer types, this is equivalent tosigned(unsigned(x) >> n)
.julia> Int8(14) >>> 2 60 julia> bits(Int8(14)) "11110010" julia> bits(Int8(60)) "00111100"
BigInt
s are treated as if having infinite size, so no filling is required and this is equivalent to>>()
.
:(start, [step, ]stop)

Range operator.
a:b
constructs a range froma
tob
with a step size of 1, anda:s:b
is similar but uses a step size ofs
. These syntaxes call the functioncolon
. The colon is also used in indexing to select whole dimensions.
colon(start, [step, ]stop)

Called by
:
syntax for constructing ranges.
range(start, [step, ]length)

Construct a range by length, given a starting value and optional step (defaults to 1).
Base.OneTo(n)

Define an
AbstractUnitRange
that behaves like1:n
, with the added distinction that the lower limit is guaranteed (by the type system) to be 1.
==(x, y)

Generic equality operator, giving a single
Bool
result. Falls back to===
. Should be implemented for all types with a notion of equality, based on the abstract value that an instance represents. For example, all numeric types are compared by numeric value, ignoring type. Strings are compared as sequences of characters, ignoring encoding.Follows IEEE semantics for floatingpoint numbers.
Collections should generally implement
==
by calling==
recursively on all contents.New numeric types should implement this function for two arguments of the new type, and handle comparison to other types via promotion rules where possible.
!=(x, y)
≠(x, y)

Notequals comparison operator. Always gives the opposite answer as
==
. New types should generally not implement this, and rely on the fallback definition!=(x,y) = !(x==y)
instead.
===(x, y)
≡(x, y)

See the
is()
operator.
!==(x, y)
≢(x, y)

Equivalent to
!is(x, y)
.
<(x, y)

Lessthan comparison operator. New numeric types should implement this function for two arguments of the new type. Because of the behavior of floatingpoint NaN values,
<
implements a partial order. Types with a canonical partial order should implement<
, and types with a canonical total order should implementisless
.
<=(x, y)
≤(x, y)

Lessthanorequals comparison operator.
>(x, y)

Greaterthan comparison operator. Generally, new types should implement
<
instead of this function, and rely on the fallback definition>(x,y) = y<x
.
>=(x, y)
≥(x, y)

Greaterthanorequals comparison operator.
.==(x, y)

Elementwise equality comparison operator.
.!=(x, y)
.≠(x, y)

Elementwise notequals comparison operator.
.<(x, y)

Elementwise lessthan comparison operator.
.<=(x, y)
.≤(x, y)

Elementwise lessthanorequals comparison operator.
.>(x, y)

Elementwise greaterthan comparison operator.
.>=(x, y)
.≥(x, y)

Elementwise greaterthanorequals comparison operator.
cmp(x, y)

Return 1, 0, or 1 depending on whether
x
is less than, equal to, or greater thany
, respectively. Uses the total order implemented byisless
. For floatingpoint numbers, uses<
but throws an error for unordered arguments.
~(x)

Bitwise not.
&(x, y)

Bitwise and.
(x, y)

Bitwise or.
$(x, y)

Bitwise exclusive or.
!(x)

Boolean not.
x && y

Shortcircuiting boolean AND.
x  y

Shortcircuiting boolean OR.
Mathematical Functions
isapprox(x, y; rtol::Real=sqrt(eps), atol::Real=0)

Inexact equality comparison:
true
ifnorm(xy) <= atol + rtol*max(norm(x), norm(y))
. The defaultatol
is zero and the defaultrtol
depends on the types ofx
andy
.For real or complex floatingpoint values,
rtol
defaults tosqrt(eps(typeof(real(xy))))
. This corresponds to requiring equality of about half of the significand digits. For other types,rtol
defaults to zero.x
andy
may also be arrays of numbers, in which casenorm
defaults tovecnorm
but may be changed by passing anorm::Function
keyword argument. (For numbers,norm
is the same thing asabs
.) Whenx
andy
are arrays, ifnorm(xy)
is not finite (i.e.±Inf
orNaN
), the comparison falls back to checking whether all elements ofx
andy
are approximately equal componentwise.The binary operator
≈
is equivalent toisapprox
with the default arguments, andx ≉ y
is equivalent to!isapprox(x,y)
.
sin(x)

Compute sine of
x
, wherex
is in radians.
cos(x)

Compute cosine of
x
, wherex
is in radians.
tan(x)

Compute tangent of
x
, wherex
is in radians.
sind(x)

Compute sine of
x
, wherex
is in degrees.
cosd(x)

Compute cosine of
x
, wherex
is in degrees.
tand(x)

Compute tangent of
x
, wherex
is in degrees.
sinpi(x)

Compute \(\sin(\pi x)\) more accurately than
sin(pi*x)
, especially for largex
.
cospi(x)

Compute \(\cos(\pi x)\) more accurately than
cos(pi*x)
, especially for largex
.
sinh(x)

Compute hyperbolic sine of
x
.
cosh(x)

Compute hyperbolic cosine of
x
.
tanh(x)

Compute hyperbolic tangent of
x
.
asin(x)

Compute the inverse sine of
x
, where the output is in radians.
acos(x)

Compute the inverse cosine of
x
, where the output is in radians
atan(x)

Compute the inverse tangent of
x
, where the output is in radians.
atan2(y, x)

Compute the inverse tangent of
y/x
, using the signs of bothx
andy
to determine the quadrant of the return value.
asind(x)

Compute the inverse sine of
x
, where the output is in degrees.
acosd(x)

Compute the inverse cosine of
x
, where the output is in degrees.
atand(x)

Compute the inverse tangent of
x
, where the output is in degrees.
sec(x)

Compute the secant of
x
, wherex
is in radians.
csc(x)

Compute the cosecant of
x
, wherex
is in radians.
cot(x)

Compute the cotangent of
x
, wherex
is in radians.
secd(x)

Compute the secant of
x
, wherex
is in degrees.
cscd(x)

Compute the cosecant of
x
, wherex
is in degrees.
cotd(x)

Compute the cotangent of
x
, wherex
is in degrees.
asec(x)

Compute the inverse secant of
x
, where the output is in radians.
acsc(x)

Compute the inverse cosecant of
x
, where the output is in radians.
acot(x)

Compute the inverse cotangent of
x
, where the output is in radians.
asecd(x)

Compute the inverse secant of
x
, where the output is in degrees.
acscd(x)

Compute the inverse cosecant of
x
, where the output is in degrees.
acotd(x)

Compute the inverse cotangent of
x
, where the output is in degrees.
sech(x)

Compute the hyperbolic secant of
x
csch(x)

Compute the hyperbolic cosecant of
x
.
coth(x)

Compute the hyperbolic cotangent of
x
.
asinh(x)

Compute the inverse hyperbolic sine of
x
.
acosh(x)

Compute the inverse hyperbolic cosine of
x
.
atanh(x)

Compute the inverse hyperbolic tangent of
x
.
asech(x)

Compute the inverse hyperbolic secant of
x
.
acsch(x)

Compute the inverse hyperbolic cosecant of
x
.
acoth(x)

Compute the inverse hyperbolic cotangent of
x
.
sinc(x)

Compute \(\sin(\pi x) / (\pi x)\) if \(x \neq 0\), and \(1\) if \(x = 0\).
cosc(x)

Compute \(\cos(\pi x) / x  \sin(\pi x) / (\pi x^2)\) if \(x \neq 0\), and \(0\) if \(x = 0\). This is the derivative of
sinc(x)
.
deg2rad(x)

Convert
x
from degrees to radians.
rad2deg(x)

Convert
x
from radians to degrees.
hypot(x, y)

Compute the hypotenuse \(\sqrt{x^2+y^2}\) avoiding overflow and underflow.
hypot(x...)

Compute the hypotenuse \(\sqrt{\sum x_i^2}\) avoiding overflow and underflow.
log(x)

Compute the natural logarithm of
x
. ThrowsDomainError
for negativeReal
arguments. Use complex negative arguments to obtain complex results.There is an experimental variant in the
Base.Math.JuliaLibm
module, which is typically faster and more accurate.
log(b, x)

Compute the base
b
logarithm ofx
. ThrowsDomainError
for negativeReal
arguments.julia> log(4,8) 1.5 julia> log(4,2) 0.5
Note
If
b
is a power of 2 or 10,log2
orlog10
should be used, as these will typically be faster and more accurate. For example,julia> log(100,1000000) 2.9999999999999996 julia> log10(1000000)/2 3.0
log2(x)

Compute the logarithm of
x
to base 2. ThrowsDomainError
for negativeReal
arguments.
log10(x)

Compute the logarithm of
x
to base 10. ThrowsDomainError
for negativeReal
arguments.
log1p(x)

Accurate natural logarithm of
1+x
. ThrowsDomainError
forReal
arguments less than 1.There is an experimental variant in the
Base.Math.JuliaLibm
module, which is typically faster and more accurate.
frexp(val)

Return
(x,exp)
such thatx
has a magnitude in the interval \([1/2, 1)\) or 0, and val = \(x \times 2^{exp}\).
exp(x)

Compute \(e^x\).
exp2(x)

Compute \(2^x\).
exp10(x)

Compute \(10^x\).
ldexp(x, n)

Compute \(x \times 2^n\).
modf(x)

Return a tuple (fpart,ipart) of the fractional and integral parts of a number. Both parts have the same sign as the argument.
expm1(x)

Accurately compute \(e^x1\).
round([T, ]x[, digits[, base]][, r::RoundingMode])

Rounds
x
to an integer value according to the providedRoundingMode
, returning a value of the same type asx
. When not specifying a rounding mode the global mode will be used (seerounding()
), which by default is round to the nearest integer (RoundNearest
mode), with ties (fractional values of 0.5) being rounded to the nearest even integer.julia> round(1.7) 2.0 julia> round(1.5) 2.0 julia> round(2.5) 2.0
The optional
RoundingMode
argument will change how the number gets rounded.round(T, x, [r::RoundingMode])
converts the result to typeT
, throwing anInexactError
if the value is not representable.round(x, digits)
rounds to the specified number of digits after the decimal place (or before if negative).round(x, digits, base)
rounds using a base other than 10.julia> round(pi, 2) 3.14 julia> round(pi, 3, 2) 3.125
Note
Rounding to specified digits in bases other than 2 can be inexact when operating on binary floating point numbers. For example, the
Float64
value represented by1.15
is actually less than 1.15, yet will be rounded to 1.2.julia> x = 1.15 1.15 julia> @sprintf "%.20f" x "1.14999999999999991118" julia> x < 115//100 true julia> round(x, 1) 1.2
RoundingMode

A type used for controlling the rounding mode of floating point operations (via
rounding()
/setrounding()
functions), or as optional arguments for rounding to the nearest integer (via theround()
function).Currently supported rounding modes are:
RoundNearest
(default)RoundNearestTiesAway
RoundNearestTiesUp
RoundToZero
RoundFromZero
(BigFloat
only)RoundUp
RoundDown
RoundNearest

The default rounding mode. Rounds to the nearest integer, with ties (fractional values of 0.5) being rounded to the nearest even integer.
RoundNearestTiesAway

Rounds to nearest integer, with ties rounded away from zero (C/C++
round()
behaviour).
RoundNearestTiesUp

Rounds to nearest integer, with ties rounded toward positive infinity (Java/JavaScript
round()
behaviour).
round(z, RoundingModeReal, RoundingModeImaginary)

Returns the nearest integral value of the same type as the complexvalued
z
toz
, breaking ties using the specifiedRoundingMode
s. The firstRoundingMode
is used for rounding the real components while the second is used for rounding the imaginary components.
ceil([T, ]x[, digits[, base]])

ceil(x)
returns the nearest integral value of the same type asx
that is greater than or equal tox
.ceil(T, x)
converts the result to typeT
, throwing anInexactError
if the value is not representable.digits
andbase
work as forround()
.
floor([T, ]x[, digits[, base]])

floor(x)
returns the nearest integral value of the same type asx
that is less than or equal tox
.floor(T, x)
converts the result to typeT
, throwing anInexactError
if the value is not representable.digits
andbase
work as forround()
.
trunc([T, ]x[, digits[, base]])

trunc(x)
returns the nearest integral value of the same type asx
whose absolute value is less than or equal tox
.trunc(T, x)
converts the result to typeT
, throwing anInexactError
if the value is not representable.digits
andbase
work as forround()
.
unsafe_trunc(T, x)

unsafe_trunc(T, x)
returns the nearest integral value of typeT
whose absolute value is less than or equal tox
. If the value is not representable byT
, an arbitrary value will be returned.
signif(x, digits[, base])

Rounds (in the sense of
round
)x
so that there aredigits
significant digits, under a basebase
representation, default 10. E.g.,signif(123.456, 2)
is120.0
, andsignif(357.913, 4, 2)
is352.0
.
min(x, y, ...)

Return the minimum of the arguments. Operates elementwise over arrays.
max(x, y, ...)

Return the maximum of the arguments. Operates elementwise over arrays.
minmax(x, y)

Return
(min(x,y), max(x,y))
. See also:extrema()
that returns(minimum(x), maximum(x))
.julia> minmax('c','b') ('b','c')
clamp(x, lo, hi)

Return
x
iflo <= x <= hi
. Ifx < lo
, returnlo
. Ifx > hi
, returnhi
. Arguments are promoted to a common type. Operates elementwise overx
ifx
is an array.julia> clamp([pi, 1.0, big(10.)], 2., 9.) 3element Array{BigFloat,1}: 3.141592653589793238462643383279502884197169399375105820974944592307816406286198 2.000000000000000000000000000000000000000000000000000000000000000000000000000000 9.000000000000000000000000000000000000000000000000000000000000000000000000000000
clamp!(array::AbstractArray, lo, hi)

Restrict values in
array
to the specified range, inplace. See alsoclamp()
.
abs(x)

The absolute value of
x
.When
abs
is applied to signed integers, overflow may occur, resulting in the return of a negative value. This overflow occurs only whenabs
is applied to the minimum representable value of a signed integer. That is, whenx == typemin(typeof(x))
,abs(x) == x < 0
, notx
as might be expected.
Base.checked_abs(x)

Calculates
abs(x)
, checking for overflow errors where applicable. For example, standard two’s complement signed integers (e.g.Int
) cannot representabs(typemin(Int))
, thus leading to an overflow.The overflow protection may impose a perceptible performance penalty.
Base.checked_neg(x)

Calculates
x
, checking for overflow errors where applicable. For example, standard two’s complement signed integers (e.g.Int
) cannot representtypemin(Int)
, thus leading to an overflow.The overflow protection may impose a perceptible performance penalty.
Base.checked_add(x, y)

Calculates
x+y
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_sub(x, y)

Calculates
xy
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_mul(x, y)

Calculates
x*y
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_div(x, y)

Calculates
div(x,y)
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_rem(x, y)

Calculates
x%y
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_fld(x, y)

Calculates
fld(x,y)
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_mod(x, y)

Calculates
mod(x,y)
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
Base.checked_cld(x, y)

Calculates
cld(x,y)
, checking for overflow errors where applicable.The overflow protection may impose a perceptible performance penalty.
abs2(x)

Squared absolute value of
x
.
copysign(x, y)

Return
x
such that it has the same sign asy
sign(x)

Return zero if
x==0
and \(x/x\) otherwise (i.e., ±1 for realx
).
signbit(x)

Returns
true
if the value of the sign ofx
is negative, otherwisefalse
.
flipsign(x, y)

Return
x
with its sign flipped ify
is negative. For exampleabs(x) = flipsign(x,x)
.
sqrt(x)

Return \(\sqrt{x}\). Throws
DomainError
for negativeReal
arguments. Use complex negative arguments instead. The prefix operator√
is equivalent tosqrt
.
isqrt(n)

Integer square root: the largest integer
m
such thatm*m <= n
.
cbrt(x)

Return \(x^{1/3}\). The prefix operator
∛
is equivalent tocbrt
.
erf(x)

Compute the error function of
x
, defined by \(\frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} dt\) for arbitrary complexx
.
erfc(x)

Compute the complementary error function of
x
, defined by \(1  \operatorname{erf}(x)\).
erfcx(x)

Compute the scaled complementary error function of
x
, defined by \(e^{x^2} \operatorname{erfc}(x)\). Note also that \(\operatorname{erfcx}(ix)\) computes the Faddeeva function \(w(x)\).
erfi(x)

Compute the imaginary error function of
x
, defined by \(i \operatorname{erf}(ix)\).
dawson(x)

Compute the Dawson function (scaled imaginary error function) of
x
, defined by \(\frac{\sqrt{\pi}}{2} e^{x^2} \operatorname{erfi}(x)\).
erfinv(x)

Compute the inverse error function of a real
x
, defined by \(\operatorname{erf}(\operatorname{erfinv}(x)) = x\).
erfcinv(x)

Compute the inverse error complementary function of a real
x
, defined by \(\operatorname{erfc}(\operatorname{erfcinv}(x)) = x\).
real(z)

Return the real part of the complex number
z
.
imag(z)

Return the imaginary part of the complex number
z
.
reim(z)

Return both the real and imaginary parts of the complex number
z
.
conj(z)

Compute the complex conjugate of a complex number
z
.
angle(z)

Compute the phase angle in radians of a complex number
z
.
cis(z)

Return \(\exp(iz)\).
binomial(n, k)

Number of ways to choose
k
out ofn
items.
factorial(n)

Factorial of
n
. Ifn
is anInteger
, the factorial is computed as an integer (promoted to at least 64 bits). Note that this may overflow ifn
is not small, but you can usefactorial(big(n))
to compute the result exactly in arbitrary precision. Ifn
is not anInteger
,factorial(n)
is equivalent togamma(n+1)
.
gcd(x, y)

Greatest common (positive) divisor (or zero if
x
andy
are both zero).
lcm(x, y)

Least common (nonnegative) multiple.
gcdx(x, y)

Computes the greatest common (positive) divisor of
x
andy
and their Bézout coefficients, i.e. the integer coefficientsu
andv
that satisfy \(ux+vy = d = gcd(x,y)\). \(gcdx(x,y)\) returns \((d,u,v)\).julia> gcdx(12, 42) (6,3,1)
julia> gcdx(240, 46) (2,9,47)
Note
Bézout coefficients are not uniquely defined.
gcdx
returns the minimal Bézout coefficients that are computed by the extended Euclidean algorithm. (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.) For signed integers, these coefficientsu
andv
are minimal in the sense that \(u < y/d\) and \(v < x/d\). Furthermore, the signs ofu
andv
are chosen so thatd
is positive. For unsigned integers, the coefficientsu
andv
might be near theirtypemax
, and the identity then holds only via the unsigned integers’ modulo arithmetic.
ispow2(n) → Bool

Test whether
n
is a power of two.
nextpow2(n)

The smallest power of two not less than
n
. Returns 0 forn==0
, and returnsnextpow2(n)
for negative arguments.
prevpow2(n)

The largest power of two not greater than
n
. Returns 0 forn==0
, and returnsprevpow2(n)
for negative arguments.
nextpow(a, x)

The smallest
a^n
not less thanx
, wheren
is a nonnegative integer.a
must be greater than 1, andx
must be greater than 0.
prevpow(a, x)

The largest
a^n
not greater thanx
, wheren
is a nonnegative integer.a
must be greater than 1, andx
must not be less than 1.
nextprod([k_1, k_2, ..., ]n)

Next integer not less than
n
that can be written as \(\prod k_i^{p_i}\) for integers \(p_1\), \(p_2\), etc.
invmod(x, m)

Take the inverse of
x
modulom
:y
such that \(x y = 1 \pmod m\), with \(div(x,y) = 0\). This is undefined for \(m = 0\), or if \(gcd(x,m) \neq 1\).
powermod(x, p, m)

Compute \(x^p \pmod m\).
gamma(x)

Compute the gamma function of
x
.
lgamma(x)

Compute the logarithm of the absolute value of
gamma()
forReal
x
, while forComplex
x
it computes the logarithm ofgamma(x)
.
lfact(x)

Compute the logarithmic factorial of
x
digamma(x)

Compute the digamma function of
x
(the logarithmic derivative ofgamma(x)
)
invdigamma(x)

Compute the inverse digamma function of
x
.
trigamma(x)

Compute the trigamma function of
x
(the logarithmic second derivative ofgamma(x)
).
polygamma(m, x)

Compute the polygamma function of order
m
of argumentx
(the(m+1)th
derivative of the logarithm ofgamma(x)
)
airy(k, x)

The
k
th derivative of the Airy function \(\operatorname{Ai}(x)\).
airyai(x)

Airy function \(\operatorname{Ai}(x)\).
airyprime(x)

Airy function derivative \(\operatorname{Ai}'(x)\).
airyaiprime(x)

Airy function derivative \(\operatorname{Ai}'(x)\).
airybi(x)

Airy function \(\operatorname{Bi}(x)\).
airybiprime(x)

Airy function derivative \(\operatorname{Bi}'(x)\).
airyx(k, x)

scaled
k
th derivative of the Airy function, return \(\operatorname{Ai}(x) e^{\frac{2}{3} x \sqrt{x}}\) fork == 0  k == 1
, and \(\operatorname{Ai}(x) e^{ \left \operatorname{Re} \left( \frac{2}{3} x \sqrt{x} \right) \right}\) fork == 2  k == 3
.
besselj0(x)

Bessel function of the first kind of order 0, \(J_0(x)\).
besselj1(x)

Bessel function of the first kind of order 1, \(J_1(x)\).
besselj(nu, x)

Bessel function of the first kind of order
nu
, \(J_\nu(x)\).
besseljx(nu, x)

Scaled Bessel function of the first kind of order
nu
, \(J_\nu(x) e^{  \operatorname{Im}(x) }\).
bessely0(x)

Bessel function of the second kind of order 0, \(Y_0(x)\).
bessely1(x)

Bessel function of the second kind of order 1, \(Y_1(x)\).
bessely(nu, x)

Bessel function of the second kind of order
nu
, \(Y_\nu(x)\).
besselyx(nu, x)

Scaled Bessel function of the second kind of order
nu
, \(Y_\nu(x) e^{  \operatorname{Im}(x) }\).
hankelh1(nu, x)

Bessel function of the third kind of order
nu
, \(H^{(1)}_\nu(x)\).
hankelh1x(nu, x)

Scaled Bessel function of the third kind of order
nu
, \(H^{(1)}_\nu(x) e^{x i}\).
hankelh2(nu, x)

Bessel function of the third kind of order
nu
, \(H^{(2)}_\nu(x)\).
hankelh2x(nu, x)

Scaled Bessel function of the third kind of order
nu
, \(H^{(2)}_\nu(x) e^{x i}\).
besselh(nu, [k=1, ]x)

Bessel function of the third kind of order
nu
(the Hankel function).k
is either 1 or 2, selectinghankelh1()
orhankelh2()
, respectively.k
defaults to 1 if it is omitted. (See alsobesselhx()
for an exponentially scaled variant.)
besselhx(nu, [k=1, ]z)

Compute the scaled Hankel function \(\exp(∓iz) H_ν^{(k)}(z)\), where \(k\) is 1 or 2, \(H_ν^{(k)}(z)\) is
besselh(nu, k, z)
, and \(∓\) is \(\) for \(k=1\) and \(+\) for \(k=2\).k
defaults to 1 if it is omitted.The reason for this function is that \(H_ν^{(k)}(z)\) is asymptotically proportional to \(\exp(∓iz)/\sqrt{z}\) for large \(z\), and so the
besselh()
function is susceptible to overflow or underflow whenz
has a large imaginary part. Thebesselhx
function cancels this exponential factor (analytically), so it avoids these problems.
besseli(nu, x)

Modified Bessel function of the first kind of order
nu
, \(I_\nu(x)\).
besselix(nu, x)

Scaled modified Bessel function of the first kind of order
nu
, \(I_\nu(x) e^{  \operatorname{Re}(x) }\).
besselk(nu, x)

Modified Bessel function of the second kind of order
nu
, \(K_\nu(x)\).
besselkx(nu, x)

Scaled modified Bessel function of the second kind of order
nu
, \(K_\nu(x) e^x\).
beta(x, y)

Euler integral of the first kind \(\operatorname{B}(x,y) = \Gamma(x)\Gamma(y)/\Gamma(x+y)\).
lbeta(x, y)

Natural logarithm of the absolute value of the beta function \(\log(\operatorname{B}(x,y))\).
eta(x)

Dirichlet eta function \(\eta(s) = \sum^\infty_{n=1}(1)^{n1}/n^{s}\).
zeta(s)

Riemann zeta function \(\zeta(s)\).
zeta(s, z)

Generalized zeta function \(\zeta(s, z)\), defined by the sum \(\sum_{k=0}^\infty ((k+z)^2)^{s/2}\), where any term with \(k+z=0\) is excluded. For \(\Re z > 0\), this definition is equivalent to the Hurwitz zeta function \(\sum_{k=0}^\infty (k+z)^{s}\). For \(z=1\), it yields the Riemann zeta function \(\zeta(s)\).
ndigits(n, b = 10)

Compute the number of digits in number
n
written in baseb
.
widemul(x, y)

Multiply
x
andy
, giving the result as a larger type.
@evalpoly(z, c...)

Evaluate the polynomial \(\sum_k c[k] z^{k1}\) for the coefficients
c[1]
,c[2]
, ...; that is, the coefficients are given in ascending order by power ofz
. This macro expands to efficient inline code that uses either Horner’s method or, for complexz
, a more efficient Goertzellike algorithm.
Statistics
mean(v[, region])

Compute the mean of whole array
v
, or optionally along the dimensions inregion
.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.
mean(f::Function, v)

Apply the function
f
to each element ofv
and take the mean.
mean!(r, v)

Compute the mean of
v
over the singleton dimensions ofr
, and write results tor
.
std(v[, region]; corrected::Bool=true, mean=nothing)

Compute the sample standard deviation of a vector or array
v
, optionally along dimensions inregion
. The algorithm returns an estimator of the generative distribution’s standard deviation under the assumption that each entry ofv
is an IID drawn from that generative distribution. This computation is equivalent to calculatingsqrt(sum((v  mean(v)).^2) / (length(v)  1))
. A precomputedmean
may be provided. Ifcorrected
istrue
, then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = length(x)
.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.
stdm(v, m::Number; corrected::Bool=true)

Compute the sample standard deviation of a vector
v
with known meanm
. Ifcorrected
istrue
, then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = length(x)
.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.
var(v[, region])

Compute the sample variance of a vector or array
v
, optionally along dimensions inregion
. The algorithm will return an estimator of the generative distribution’s variance under the assumption that each entry ofv
is an IID drawn from that generative distribution. This computation is equivalent to calculatingsumabs2(v  mean(v)) / (length(v)  1)
. Note: Julia does not ignoreNaN
values in the computation. For applications requiring the handling of missing data, theDataArray
package is recommended.
varm(v, m[, region]; corrected::Bool=true)

Compute the sample variance of a collection
v
with known mean(s)m
, optionally overregion
.m
may contain means for each dimension ofv
. Ifcorrected
istrue
, then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = length(x)
.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.
middle(x)

Compute the middle of a scalar value, which is equivalent to
x
itself, but of the type ofmiddle(x, x)
for consistency.
middle(x, y)

Compute the middle of two reals
x
andy
, which is equivalent in both value and type to computing their mean ((x + y) / 2
).
middle(range)

Compute the middle of a range, which consists of computing the mean of its extrema. Since a range is sorted, the mean is performed with the first and last element.
julia> middle(1:10) 5.5
middle(a)

Compute the middle of an array
a
, which consists of finding its extrema and then computing their mean.julia> a = [1,2,3.6,10.9] 4element Array{Float64,1}: 1.0 2.0 3.6 10.9 julia> middle(a) 5.95
median(v[, region])

Compute the median of an entire array
v
, or, optionally, along the dimensions inregion
. For an even number of elements no exact median element exists, so the result is equivalent to calculating mean of two median elements.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.
median!(v)

Like
median
, but may overwrite the input vector.
midpoints(e)

Compute the midpoints of the bins with edges
e
. The result is a vector/range of lengthlength(e)  1
. Note: Julia does not ignoreNaN
values in the computation.
quantile(v, p; sorted=false)

Compute the quantile(s) of a vector
v
at a specified probability or vectorp
. The keyword argumentsorted
indicates whetherv
can be assumed to be sorted.The
p
should be on the interval [0,1], andv
should not have anyNaN
values.Quantiles are computed via linear interpolation between the points
((k1)/(n1), v[k])
, fork = 1:n
wheren = length(v)
. This corresponds to Definition 7 of Hyndman and Fan (1996), and is the same as the R default.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.quantile
will throw anArgumentError
in the presence ofNaN
values in the data array. Hyndman, R.J and Fan, Y. (1996) “Sample Quantiles in Statistical Packages”, The American Statistician, Vol. 50, No. 4, pp. 361365
quantile!([q, ]v, p; sorted=false)

Compute the quantile(s) of a vector
v
at the probabilitiesp
, with optional output into arrayq
(if not provided, a new output array is created). The keyword argumentsorted
indicates whetherv
can be assumed to be sorted; iffalse
(the default), then the elements ofv
may be partially sorted.The elements of
p
should be on the interval [0,1], andv
should not have anyNaN
values.Quantiles are computed via linear interpolation between the points
((k1)/(n1), v[k])
, fork = 1:n
wheren = length(v)
. This corresponds to Definition 7 of Hyndman and Fan (1996), and is the same as the R default.Note
Julia does not ignore
NaN
values in the computation. For applications requiring the handling of missing data, theDataArrays.jl
package is recommended.quantile!
will throw anArgumentError
in the presence ofNaN
values in the data array. Hyndman, R.J and Fan, Y. (1996) “Sample Quantiles in Statistical Packages”, The American Statistician, Vol. 50, No. 4, pp. 361365
cov(x[, corrected=true])

Compute the variance of the vector
x
. Ifcorrected
istrue
(the default) then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = length(x)
.
cov(X[, vardim=1, corrected=true])

Compute the covariance matrix of the matrix
X
along the dimensionvardim
. Ifcorrected
istrue
(the default) then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = size(X, vardim)
.
cov(x, y[, corrected=true])

Compute the covariance between the vectors
x
andy
. Ifcorrected
istrue
(the default) then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = length(x) = length(y)
.
cov(X, Y[, vardim=1, corrected=true])

Compute the covariance between the vectors or matrices
X
andY
along the dimensionvardim
. Ifcorrected
istrue
(the default) then the sum is scaled withn1
, whereas the sum is scaled withn
ifcorrected
isfalse
wheren = size(X, vardim) = size(Y, vardim)
.
cor(x)

Return the number one.
cor(X[, vardim=1])

Compute the Pearson correlation matrix of the matrix
X
along the dimensionvardim
.
cor(x, y)

Compute the Pearson correlation between the vectors
x
andy
.
cor(X, Y[, vardim=1])

Compute the Pearson correlation between the vectors or matrices
X
andY
along the dimensionvardim
.
Signal Processing
Fast Fourier transform (FFT) functions in Julia are implemented by calling functions from FFTW.
fft(A[, dims])

Performs a multidimensional FFT of the array
A
. The optionaldims
argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most efficient if the size ofA
along the transformed dimensions is a product of small primes; seenextprod()
. See alsoplan_fft()
for even greater efficiency.A onedimensional FFT computes the onedimensional discrete Fourier transform (DFT) as defined by
\[\operatorname{DFT}(A)[k] = \sum_{n=1}^{\operatorname{length}(A)} \exp\left(i\frac{2\pi (n1)(k1)}{\operatorname{length}(A)} \right) A[n].\]A multidimensional FFT simply performs this operation along each transformed dimension of
A
.Note
 Julia starts FFTW up with 1 thread by default. Higher performance is usually possible by increasing number of threads. Use
FFTW.set_num_threads(Sys.CPU_CORES)
to use as many threads as cores on your system.  This performs a multidimensional FFT by default. FFT libraries in other languages such as Python and Octave perform a onedimensional FFT along the first nonsingleton dimension of the array. This is worth noting while performing comparisons. For more details, refer to the Noteworthy Differences from other Languages section of the manual.
 Julia starts FFTW up with 1 thread by default. Higher performance is usually possible by increasing number of threads. Use
fft!(A[, dims])

Same as
fft()
, but operates inplace onA
, which must be an array of complex floatingpoint numbers.
ifft(A[, dims])

Multidimensional inverse FFT.
A onedimensional inverse FFT computes
\[\operatorname{IDFT}(A)[k] = \frac{1}{\operatorname{length}(A)} \sum_{n=1}^{\operatorname{length}(A)} \exp\left(+i\frac{2\pi (n1)(k1)} {\operatorname{length}(A)} \right) A[n].\]A multidimensional inverse FFT simply performs this operation along each transformed dimension of
A
.
ifft!(A[, dims])

Same as
ifft()
, but operates inplace onA
.
bfft(A[, dims])

Similar to
ifft()
, but computes an unnormalized inverse (backward) transform, which must be divided by the product of the sizes of the transformed dimensions in order to obtain the inverse. (This is slightly more efficient thanifft()
because it omits a scaling step, which in some applications can be combined with other computational steps elsewhere.)\[\operatorname{BDFT}(A)[k] = \operatorname{length}(A) \operatorname{IDFT}(A)[k]\]
bfft!(A[, dims])

Same as
bfft()
, but operates inplace onA
.
plan_fft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Preplan an optimized FFT along given dimensions (
dims
) of arrays matching the shape and type ofA
. (The first two arguments have the same meaning as forfft()
.) Returns an objectP
which represents the linear operator computed by the FFT, and which contains all of the information needed to computefft(A, dims)
quickly.To apply
P
to an arrayA
, useP * A
; in general, the syntax for applying plans is much like that of matrices. (A plan can only be applied to arrays of the same size as theA
for which the plan was created.) You can also apply a plan with a preallocated output arrayÂ
by callingA_mul_B!(Â, plan, A)
. (ForA_mul_B!
, however, the input arrayA
must be a complex floatingpoint array like the outputÂ
.) You can compute the inversetransform plan byinv(P)
and apply the inverse plan withP \ Â
(the inverse plan is cached and reused for subsequent calls toinv
or\
), and apply the inverse plan to a preallocated output arrayA
withA_ldiv_B!(A, P, Â)
.The
flags
argument is a bitwiseor of FFTW planner flags, defaulting toFFTW.ESTIMATE
. e.g. passingFFTW.MEASURE
orFFTW.PATIENT
will instead spend several seconds (or more) benchmarking different possible FFT algorithms and picking the fastest one; see the FFTW manual for more information on planner flags. The optionaltimelimit
argument specifies a rough upper bound on the allowed planning time, in seconds. PassingFFTW.MEASURE
orFFTW.PATIENT
may cause the input arrayA
to be overwritten with zeros during plan creation.plan_fft!()
is the same asplan_fft()
but creates a plan that operates inplace on its argument (which must be an array of complex floatingpoint numbers).plan_ifft()
and so on are similar but produce plans that perform the equivalent of the inverse transformsifft()
and so on.
plan_ifft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as
plan_fft()
, but produces a plan that performs inverse transformsifft()
.
plan_bfft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as
plan_fft()
, but produces a plan that performs an unnormalized backwards transformbfft()
.
plan_fft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as
plan_fft()
, but operates inplace onA
.
plan_ifft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as
plan_ifft()
, but operates inplace onA
.
plan_bfft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as
plan_bfft()
, but operates inplace onA
.
rfft(A[, dims])

Multidimensional FFT of a real array
A
, exploiting the fact that the transform has conjugate symmetry in order to save roughly half the computational time and storage costs compared withfft()
. IfA
has size(n_1, ..., n_d)
, the result has size(div(n_1,2)+1, ..., n_d)
.The optional
dims
argument specifies an iterable subset of one or more dimensions ofA
to transform, similar tofft()
. Instead of (roughly) halving the first dimension ofA
in the result, thedims[1]
dimension is (roughly) halved in the same way.
irfft(A, d[, dims])

Inverse of
rfft()
: for a complex arrayA
, gives the corresponding real array whose FFT yieldsA
in the first half. As forrfft()
,dims
is an optional subset of dimensions to transform, defaulting to1:ndims(A)
.d
is the length of the transformed real array along thedims[1]
dimension, which must satisfydiv(d,2)+1 == size(A,dims[1])
. (This parameter cannot be inferred fromsize(A)
since both2*size(A,dims[1])2
as well as2*size(A,dims[1])1
are valid sizes for the transformed real array.)
brfft(A, d[, dims])

Similar to
irfft()
but computes an unnormalized inverse transform (similar tobfft()
), which must be divided by the product of the sizes of the transformed dimensions (of the real output array) in order to obtain the inverse transform.
plan_rfft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Preplan an optimized realinput FFT, similar to
plan_fft()
except forrfft()
instead offft()
. The first two arguments, and the size of the transformed result, are the same as forrfft()
.
plan_brfft(A, d [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Preplan an optimized realinput unnormalized transform, similar to
plan_rfft()
except forbrfft()
instead ofrfft()
. The first two arguments and the size of the transformed result, are the same as forbrfft()
.
plan_irfft(A, d [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Preplan an optimized inverse realinput FFT, similar to
plan_rfft()
except forirfft()
andbrfft()
, respectively. The first three arguments have the same meaning as forirfft()
.
dct(A[, dims])

Performs a multidimensional typeII discrete cosine transform (DCT) of the array
A
, using the unitary normalization of the DCT. The optionaldims
argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most efficient if the size ofA
along the transformed dimensions is a product of small primes; seenextprod()
. See alsoplan_dct()
for even greater efficiency.
dct!(A[, dims])

Same as
dct!()
, except that it operates inplace onA
, which must be an array of real or complex floatingpoint values.
idct(A[, dims])

Computes the multidimensional inverse discrete cosine transform (DCT) of the array
A
(technically, a typeIII DCT with the unitary normalization). The optionaldims
argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most efficient if the size ofA
along the transformed dimensions is a product of small primes; seenextprod()
. See alsoplan_idct()
for even greater efficiency.
idct!(A[, dims])

Same as
idct!()
, but operates inplace onA
.
plan_dct(A[, dims[, flags[, timelimit]]])

Preplan an optimized discrete cosine transform (DCT), similar to
plan_fft()
except producing a function that computesdct()
. The first two arguments have the same meaning as fordct()
.
plan_dct!(A[, dims[, flags[, timelimit]]])

Same as
plan_dct()
, but operates inplace onA
.
plan_idct(A[, dims[, flags[, timelimit]]])

Preplan an optimized inverse discrete cosine transform (DCT), similar to
plan_fft()
except producing a function that computesidct()
. The first two arguments have the same meaning as foridct()
.
plan_idct!(A[, dims[, flags[, timelimit]]])

Same as
plan_idct()
, but operates inplace onA
.
fftshift(x)

Swap the first and second halves of each dimension of
x
.
fftshift(x, dim)

Swap the first and second halves of the given dimension of array
x
.
ifftshift(x[, dim])

Undoes the effect of
fftshift
.
filt(b, a, x[, si])

Apply filter described by vectors
a
andb
to vectorx
, with an optional initial filter state vectorsi
(defaults to zeros).
filt!(out, b, a, x[, si])

Same as
filt()
but writes the result into theout
argument, which may alias the inputx
to modify it inplace.
deconv(b, a)

Construct vector
c
such thatb = conv(a,c) + r
. Equivalent to polynomial division.
conv(u, v)

Convolution of two vectors. Uses FFT algorithm.
conv2(u, v, A)

2D convolution of the matrix
A
with the 2D separable kernel generated by the vectorsu
andv
. Uses 2D FFT algorithm.
conv2(B, A)

2D convolution of the matrix
B
with the matrixA
. Uses 2D FFT algorithm.
xcorr(u, v)

Compute the crosscorrelation of two vectors.
The following functions are defined within the Base.FFTW
module.
r2r(A, kind[, dims])

Performs a multidimensional realinput/realoutput (r2r) transform of type
kind
of the arrayA
, as defined in the FFTW manual.kind
specifies either a discrete cosine transform of various types (FFTW.REDFT00
,FFTW.REDFT01
,FFTW.REDFT10
, orFFTW.REDFT11
), a discrete sine transform of various types (FFTW.RODFT00
,FFTW.RODFT01
,FFTW.RODFT10
, orFFTW.RODFT11
), a realinput DFT with halfcomplexformat output (FFTW.R2HC
and its inverseFFTW.HC2R
), or a discrete Hartley transform (FFTW.DHT
). Thekind
argument may be an array or tuple in order to specify different transform types along the different dimensions ofA
;kind[end]
is used for any unspecified dimensions. See the FFTW manual for precise definitions of these transform types, at http://www.fftw.org/doc.The optional
dims
argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along.kind[i]
is then the transform type fordims[i]
, withkind[end]
being used fori > length(kind)
.See also
plan_r2r()
to preplan optimized r2r transforms.
r2r!(A, kind[, dims])

Same as
r2r()
, but operates inplace onA
, which must be an array of real or complex floatingpoint numbers.
plan_r2r(A, kind[, dims[, flags[, timelimit]]])

Preplan an optimized r2r transform, similar to
Base.plan_fft()
except that the transforms (and the first three arguments) correspond tor2r()
andr2r!()
, respectively.
plan_r2r!(A, kind[, dims[, flags[, timelimit]]])

Similar to
Base.plan_fft()
, but corresponds tor2r!()
.
Numerical Integration
Although several external packages are available for numeric integration and solution of ordinary differential equations, we also provide some builtin integration support in Julia.
quadgk(f, a, b, c...; reltol=sqrt(eps), abstol=0, maxevals=10^7, order=7, norm=vecnorm)

Numerically integrate the function
f(x)
froma
tob
, and optionally over additional intervalsb
toc
and so on. Keyword options include a relative error tolerancereltol
(defaults tosqrt(eps)
in the precision of the endpoints), an absolute error toleranceabstol
(defaults to 0), a maximum number of function evaluationsmaxevals
(defaults to10^7
), and theorder
of the integration rule (defaults to 7).Returns a pair
(I,E)
of the estimated integralI
and an estimated upper bound on the absolute errorE
. Ifmaxevals
is not exceeded thenE <= max(abstol, reltol*norm(I))
will hold. (Note that it is useful to specify a positiveabstol
in cases wherenorm(I)
may be zero.)The endpoints
a
et cetera can also be complex (in which case the integral is performed over straightline segments in the complex plane). If the endpoints areBigFloat
, then the integration will be performed inBigFloat
precision as well.Note
It is advisable to increase the integration
order
in rough proportion to the precision, for smooth integrands.More generally, the precision is set by the precision of the integration endpoints (promoted to floatingpoint types).
The integrand
f(x)
can return any numeric scalar, vector, or matrix type, or in fact any type supporting+
,
, multiplication by real values, and anorm
(i.e., any normed vector space). Alternatively, a different norm can be specified by passing anorm
like function as thenorm
keyword argument (which defaults tovecnorm
).Note
Only onedimensional integrals are provided by this function. For multidimensional integration (cubature), there are many different algorithms (often much better than simple nested 1d integrals) and the optimal choice tends to be very problemdependent. See the Julia externalpackage listing for available algorithms for multidimensional integration or other specialized tasks (such as integrals of highly oscillatory or singular functions).
The algorithm is an adaptive GaussKronrod integration technique: the integral in each interval is estimated using a Kronrod rule (
2*order+1
points) and the error is estimated using an embedded Gauss rule (order
points). The interval with the largest error is then subdivided into two intervals and the process is repeated until the desired error tolerance is achieved.These quadrature rules work best for smooth functions within each interval, so if your function has a known discontinuity or other singularity, it is best to subdivide your interval to put the singularity at an endpoint. For example, if
f
has a discontinuity atx=0.7
and you want to integrate from 0 to 1, you should usequadgk(f, 0,0.7,1)
to subdivide the interval at the point of discontinuity. The integrand is never evaluated exactly at the endpoints of the intervals, so it is possible to integrate functions that diverge at the endpoints as long as the singularity is integrable (for example, alog(x)
or1/sqrt(x)
singularity).For realvalued endpoints, the starting and/or ending points may be infinite. (A coordinate transformation is performed internally to map the infinite interval to a finite one.)
© 2009–2016 Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors
Licensed under the MIT License.
https://docs.julialang.org/en/release0.5/stdlib/math/