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 of y on the right. Gives floating-point results for integer arguments.

\(x, y)

Left division operator: multiplication of y by the inverse of x on the left. Gives floating-point results for integer arguments.

^(x, y)

Exponentiation operator.

.+(x, y)

Element-wise addition operator.

.-(x, y)

Element-wise subtraction operator.

.*(x, y)

Element-wise multiplication operator.

./(x, y)

Element-wise right division operator.

.\(x, y)

Element-wise left division operator.

.^(x, y)

Element-wise exponentiation operator.

fma(x, y, z)

Computes x*y+z without rounding the intermediate result x*y. On some systems this is significantly more expensive than x*y+z. fma is used to improve accuracy in certain algorithms. See muladd.

muladd(x, y, z)

Combined multiply-add, computes x*y+z in an efficient manner. This may on some systems be equivalent to x*y+z, or to fma(x,y,z). muladd is used to improve performance. See fma.

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 division, returning in the range [0,``y``), if y is positive, or (y,0] if y is negative.

mod2pi(x)

Modulus after division by 2pi, returning in the range [0,2pi).

This function computes a floating point representation of the modulus after division by numerically exact 2pi, and is therefore not exactly the same as mod(x,2pi), which would compute the modulus of x relative to division by the floating-point number 2pi.

rem(x, y)
%(x, y)

Remainder from Euclidean division, returning a value of the same sign as x, and smaller in magnitude than y. This value is always exact.

divrem(x, y)

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

fldmod(x, y)

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

mod1(x, m)

Modulus after division, returning in the range (0,m]

rem1(x, m)

Remainder after division, returning in the range (0,m]

//(num, den)

Divide two integers or rational numbers, giving a Rational result.

rationalize([Type=Int, ]x; tol=eps(x))

Approximate floating point number x as a Rational number with components of the given integer type. The result will differ from x by no more than tol.

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)

Right bit shift operator, preserving the sign of x.

>>>(x, n)

Unsigned right bit shift operator.

:(start, [step, ]stop)

Range operator. a:b constructs a range from a to b with a step size of 1, and a:s:b is similar but uses a step size of s. These syntaxes call the function colon. 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).

==(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 floating-point 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)

Not-equals 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)

Less-than comparison operator. New numeric types should implement this function for two arguments of the new type. Because of the behavior of floating-point NaN values, < implements a partial order. Types with a canonical partial order should implement <, and types with a canonical total order should implement isless.

<=(x, y)
≤(x, y)

Less-than-or-equals comparison operator.

>(x, y)

Greater-than 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)

Greater-than-or-equals comparison operator.

.==(x, y)

Element-wise equality comparison operator.

.!=(x, y)
.≠(x, y)

Element-wise not-equals comparison operator.

.<(x, y)

Element-wise less-than comparison operator.

.<=(x, y)
.≤(x, y)

Element-wise less-than-or-equals comparison operator.

.>(x, y)

Element-wise greater-than comparison operator.

.>=(x, y)
.≥(x, y)

Element-wise greater-than-or-equals comparison operator.

cmp(x, y)

Return -1, 0, or 1 depending on whether x is less than, equal to, or greater than y, respectively. Uses the total order implemented by isless. For floating-point 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

Short-circuiting boolean AND

x || y

Short-circuiting boolean OR

A_ldiv_Bc(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(A\) \ \(Bᴴ\)

A_ldiv_Bt(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(A\) \ \(Bᵀ\)

A_mul_B!(Y, A, B) → Y

Calculates the matrix-matrix or matrix-vector product \(A⋅B\) and stores the result in Y, overwriting the existing value of Y. Note that Y must not be aliased with either A or B.

julia>A=[1.02.0;3.04.0];B=[1.01.0;1.01.0];Y=similar(B);A_mul_B!(Y,A,B);julia>Y2x2Array{Float64,2}:3.03.07.07.0
A_mul_Bc(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(A⋅Bᴴ\)

A_mul_Bt(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(A⋅Bᵀ\)

A_rdiv_Bc(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(A / Bᴴ\)

A_rdiv_Bt(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(A / Bᵀ\)

Ac_ldiv_B(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᴴ\) \ \(B\)

Ac_ldiv_Bc(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᴴ\) \ \(Bᴴ\)

Ac_mul_B(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᴴ⋅B\)

Ac_mul_Bc(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᴴ Bᴴ\)

Ac_rdiv_B(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᴴ / B\)

Ac_rdiv_Bc(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᴴ / Bᴴ\)

At_ldiv_B(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᵀ\) \ \(B\)

At_ldiv_Bt(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᵀ\) \ \(Bᵀ\)

At_mul_B(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᵀ⋅B\)

At_mul_Bt(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᵀ⋅Bᵀ\)

At_rdiv_B(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᵀ / B\)

At_rdiv_Bt(A, B)

For matrices or vectors \(A\) and \(B\), calculates \(Aᵀ / Bᵀ\)

Mathematical Functions

isapprox(x, y; rtol::Real=sqrt(eps), atol::Real=0)

Inexact equality comparison: true if norm(x-y)<=atol+rtol*max(norm(x),norm(y)). The default atol is zero and the default rtol depends on the types of x and y.

For real or complex floating-point values, rtol defaults to sqrt(eps(typeof(real(x-y)))). This corresponds to requiring equality of about half of the significand digits. For other types, rtol defaults to zero.

x and y may also be arrays of numbers, in which case norm defaults to vecnorm but may be changed by passing a norm::Function keyword argument. (For numbers, norm is the same thing as abs.)

The binary operator ≈ is equivalent to isapprox with the default arguments, and x≉y is equivalent to !isapprox(x,y).

sin(x)

Compute sine of x, where x is in radians

cos(x)

Compute cosine of x, where x is in radians

tan(x)

Compute tangent of x, where x is in radians

sind(x)

Compute sine of x, where x is in degrees.

cosd(x)

Compute cosine of x, where x is in degrees

tand(x)

Compute tangent of x, where x is in degrees

sinpi(x)

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

cospi(x)

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

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 both x and y 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, where x is in radians

csc(x)

Compute the cosecant of x, where x is in radians

cot(x)

Compute the cotangent of x, where x is in radians

secd(x)

Compute the secant of x, where x is in degrees

cscd(x)

Compute the cosecant of x, where x is in degrees

cotd(x)

Compute the cotangent of x, where x 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 \(\sqrt{x^2+y^2}\) avoiding overflow and underflow

log(x)

Compute the natural logarithm of x. Throws DomainError for negative Real 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 of x. Throws DomainError for negative Real arguments.

log2(x)

Compute the logarithm of x to base 2. Throws DomainError for negative Real arguments.

log10(x)

Compute the logarithm of x to base 10. Throws DomainError for negative Real arguments.

log1p(x)

Accurate natural logarithm of 1+x. Throws DomainError for Real 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 that x 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^x-1\).

round([T, ]x[, digits[, base]][, r::RoundingMode])

round(x) rounds x to an integer value according to the default rounding mode (see get_rounding()), returning a value of the same type as x. By default (RoundNearest), this will round to the nearest integer, with ties (fractional values of 0.5) being rounded to the even integer.

julia>round(1.7)2.0julia>round(1.5)2.0julia>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 type T, throwing an InexactError 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.14julia>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 by 1.15 is actually less than 1.15, yet will be rounded to 1.2.

julia>x=1.151.15julia>@sprintf"%.20f"x"1.14999999999999991118"julia>x<115//100truejulia>round(x,1)1.2
RoundingMode

A type which controls rounding behavior. Currently supported rounding modes are:

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).

RoundToZero

round() using this rounding mode is an alias for trunc().

RoundUp

round() using this rounding mode is an alias for ceil().

RoundDown

round() using this rounding mode is an alias for floor().

round(z, RoundingModeReal, RoundingModeImaginary)

Returns the nearest integral value of the same type as the complex-valued z to z, breaking ties using the specified RoundingModes. The first RoundingMode 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 as x that is greater than or equal to x.

ceil(T,x) converts the result to type T, throwing an InexactError if the value is not representable.

digits and base work as for round().

floor([T, ]x[, digits[, base]])

floor(x) returns the nearest integral value of the same type as x that is less than or equal to x.

floor(T,x) converts the result to type T, throwing an InexactError if the value is not representable.

digits and base work as for round().

trunc([T, ]x[, digits[, base]])

trunc(x) returns the nearest integral value of the same type as x whose absolute value is less than or equal to x.

trunc(T,x) converts the result to type T, throwing an InexactError if the value is not representable.

digits and base work as for round().

unsafe_trunc(T, x)

unsafe_trunc(T,x) returns the nearest integral value of type T whose absolute value is less than or equal to x. If the value is not representable by T, an arbitrary value will be returned.

signif(x, digits[, base])

Rounds (in the sense of round) x so that there are digits significant digits, under a base base representation, default 10. E.g., signif(123.456,2) is 120.0, and signif(357.913,4,2) is 352.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))

clamp(x, lo, hi)

Return x if lo<=x<=hi. If x<lo, return lo. If x>hi, return hi. Arguments are promoted to a common type. Operates elementwise over x if it is an array.

abs(x)

Absolute value of x

abs2(x)

Squared absolute value of x

copysign(x, y)

Return x such that it has the same sign as y

sign(x)

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

signbit(x)

Returns true if the value of the sign of x is negative, otherwise false.

flipsign(x, y)

Return x with its sign flipped if y is negative. For example abs(x)=flipsign(x,x).

sqrt(x)

Return \(\sqrt{x}\). Throws DomainError for negative Real arguments. Use complex negative arguments instead. The prefix operator √ is equivalent to sqrt.

isqrt(n)

Integer square root: the largest integer m such that m*m<=n.

cbrt(x)

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

erf(x)

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

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 of n items

factorial(n)

Factorial of n. If n is an Integer, the factorial is computed as an integer (promoted to at least 64 bits). Note that this may overflow if n is not small, but you can use factorial(big(n)) to compute the result exactly in arbitrary precision. If n is not an Integer, factorial(n) is equivalent to gamma(n+1).

factorial(n, k)

Compute factorial(n)/factorial(k)

factor(n) → Dict

Compute the prime factorization of an integer n. Returns a dictionary. The keys of the dictionary correspond to the factors, and hence are of the same type as n. The value associated with each key indicates the number of times the factor appears in the factorization.

julia>factor(100)# == 2*2*5*5Dict{Int64,Int64}with2entries:2=>25=>2
gcd(x, y)

Greatest common (positive) divisor (or zero if x and y are both zero).

lcm(x, y)

Least common (non-negative) multiple.

gcdx(x, y)

Computes the greatest common (positive) divisor of x and y and their Bézout coefficients, i.e. the integer coefficients u and v that satisfy \(ux+vy = d = gcd(x,y)\).

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 Euclid algorithm. (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.) These coefficients u and v are minimal in the sense that \(|u| < |\frac y d\) and \(|v| < |\frac x d\). Furthermore, the signs of u and v are chosen so that d is positive.

ispow2(n) → Bool

Test whether n is a power of two

nextpow2(n)

The smallest power of two not less than n. Returns 0 for n==0, and returns -nextpow2(-n) for negative arguments.

prevpow2(n)

The largest power of two not greater than n. Returns 0 for n==0, and returns -prevpow2(-n) for negative arguments.

nextpow(a, x)

The smallest a^n not less than x, where n is a non-negative integer. a must be greater than 1, and x must be greater than 0.

prevpow(a, x)

The largest a^n not greater than x, where n is a non-negative integer. a must be greater than 1, and x 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.

prevprod([k_1, k_2, ..., ]n)

Previous integer not greater 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 modulo m: y such that \(xy = 1 \pmod m\).

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() for Realx, while for Complexx it computes the logarithm of gamma(x).

lfact(x)

Compute the logarithmic factorial of x

digamma(x)

Compute the digamma function of x (the logarithmic derivative of gamma(x))

invdigamma(x)

Compute the inverse digamma function of x.

trigamma(x)

Compute the trigamma function of x (the logarithmic second derivative of gamma(x))

polygamma(m, x)

Compute the polygamma function of order m of argument x (the (m+1)th derivative of the logarithm of gamma(x))

airy(k, x)

The kth 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 kth derivative of the Airy function, return \(\operatorname{Ai}(x) e^{\frac{2}{3} x \sqrt{x}}\) for k==0||k==1, and \(\operatorname{Ai}(x) e^{- \left| \operatorname{Re} \left( \frac{2}{3} x \sqrt{x} \right) \right|}\) for k==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, x)

Bessel function of the third kind of order nu (Hankel function). k is either 1 or 2, selecting hankelh1 or hankelh2, respectively.

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}(-)^{n-1}/n^{s}\).

zeta(s)

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

zeta(s, z)

Hurwitz zeta function \(\zeta(s, z)\). (This is equivalent to the Riemann zeta function \(\zeta(s)\) for the case of z=1.)

ndigits(n, b)

Compute the number of digits in number n written in base b.

widemul(x, y)

Multiply x and y, giving the result as a larger type.

@evalpoly(z, c...)

Evaluate the polynomial \(\sum_k c[k] z^{k-1}\) for the coefficients c[1], c[2], ...; that is, the coefficients are given in ascending order by power of z. This macro expands to efficient inline code that uses either Horner’s method or, for complex z, a more efficient Goertzel-like algorithm.

Statistics

mean(v[, region])

Compute the mean of whole array v, or optionally along the dimensions in region. Note: Julia does not ignore NaN values in the computation. For applications requiring the handling of missing data, the DataArray package is recommended.

mean!(r, v)

Compute the mean of v over the singleton dimensions of r, and write results to r.

std(v[, region])

Compute the sample standard deviation of a vector or array v, optionally along dimensions in region. The algorithm returns an estimator of the generative distribution’s standard deviation under the assumption that each entry of v is an IID drawn from that generative distribution. This computation is equivalent to calculating sqrt(sum((v-mean(v)).^2)/(length(v)-1)). Note: Julia does not ignore NaN values in the computation. For applications requiring the handling of missing data, the DataArray package is recommended.

stdm(v, m)

Compute the sample standard deviation of a vector v with known mean m. Note: Julia does not ignore NaN values in the computation.

var(v[, region])

Compute the sample variance of a vector or array v, optionally along dimensions in region. The algorithm will return an estimator of the generative distribution’s variance under the assumption that each entry of v is an IID drawn from that generative distribution. This computation is equivalent to calculating sumabs2(v-mean(v))/(length(v)-1). Note: Julia does not ignore NaN values in the computation. For applications requiring the handling of missing data, the DataArray package is recommended.

varm(v, m)

Compute the sample variance of a vector v with known mean m. Note: Julia does not ignore NaN values in the computation.

middle(x)

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

middle(x, y)

Compute the middle of two reals x and y, 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 in computing the mean of its extrema. Since a range is sorted, the mean is performed with the first and last element.

middle(array)

Compute the middle of an array, which consists in finding its extrema and then computing their mean.

median(v[, region])

Compute the median of whole array v, or optionally along the dimensions in region. For even number of elements no exact median element exists, so the result is equivalent to calculating mean of two median elements. NaN is returned if the data contains any NaN values. For applications requiring the handling of missing data, the DataArrays package is recommended.

median!(v)

Like median, but may overwrite the input vector.

hist(v[, n]) → e, counts

Compute the histogram of v, optionally using approximately n bins. The return values are a range e, which correspond to the edges of the bins, and counts containing the number of elements of v in each bin. Note: Julia does not ignore NaN values in the computation.

hist(v, e) → e, counts

Compute the histogram of v using a vector/range e as the edges for the bins. The result will be a vector of length length(e)-1, such that the element at location i satisfies sum(e[i].<v.<=e[i+1]). Note: Julia does not ignore NaN values in the computation.

hist!(counts, v, e) → e, counts

Compute the histogram of v, using a vector/range e as the edges for the bins. This function writes the resultant counts to a pre-allocated array counts.

hist2d(M, e1, e2) -> (edge1, edge2, counts)

Compute a “2d histogram” of a set of N points specified by N-by-2 matrix M. Arguments e1 and e2 are bins for each dimension, specified either as integer bin counts or vectors of bin edges. The result is a tuple of edge1 (the bin edges used in the first dimension), edge2 (the bin edges used in the second dimension), and counts, a histogram matrix of size (length(edge1)-1,length(edge2)-1). Note: Julia does not ignore NaN values in the computation.

hist2d!(counts, M, e1, e2) -> (e1, e2, counts)

Compute a “2d histogram” with respect to the bins delimited by the edges given in e1 and e2. This function writes the results to a pre-allocated array counts.

histrange(v, n)

Compute nice bin ranges for the edges of a histogram of v, using approximately n bins. The resulting step sizes will be 1, 2 or 5 multiplied by a power of 10. Note: Julia does not ignore NaN values in the computation.

midpoints(e)

Compute the midpoints of the bins with edges e. The result is a vector/range of length length(e)-1. Note: Julia does not ignore NaN values in the computation.

quantile(v, p; sorted=false)

Compute the quantile(s) of a vector v at a specified probability or vector p. The keyword argument sorted indicates whether v can be assumed to be sorted.

The p should be on the interval [0,1], and v should not have any NaN values.

Quantiles are computed via linear interpolation between the points ((k-1)/(n-1),v[k]), for k=1:n where n=length(v). This corresponds to Definition 7 of Hyndman and Fan (1996), and is the same as the R default.

  • Hyndman, R.J and Fan, Y. (1996) “Sample Quantiles in Statistical Packages”, The American Statistician, Vol. 50, No. 4, pp. 361-365
quantile!([q, ]v, p; sorted=false)

Compute the quantile(s) of a vector v at the probabilities p, with optional output into array q (if not provided, a new output array is created). The keyword argument sorted indicates whether v can be assumed to be sorted; if false (the default), then the elements of v may be partially sorted.

The elements of p should be on the interval [0,1], and v should not have any NaN values.

Quantiles are computed via linear interpolation between the points ((k-1)/(n-1),v[k]), for k=1:n where n=length(v). This corresponds to Definition 7 of Hyndman and Fan (1996), and is the same as the R default.

  • Hyndman, R.J and Fan, Y. (1996) “Sample Quantiles in Statistical Packages”, The American Statistician, Vol. 50, No. 4, pp. 361-365
cov(v1[, v2][, vardim=1, corrected=true, mean=nothing])

Compute the Pearson covariance between the vector(s) in v1 and v2. Here, v1 and v2 can be either vectors or matrices.

This function accepts three keyword arguments:

  • vardim: the dimension of variables. When vardim=1, variables are considered in columns while observations in rows; when vardim=2, variables are in rows while observations in columns. By default, it is set to 1.
  • corrected: whether to apply Bessel’s correction (divide by n-1 instead of n). By default, it is set to true.
  • mean: allow users to supply mean values that are known. By default, it is set to nothing, which indicates that the mean(s) are unknown, and the function will compute the mean. Users can use mean=0 to indicate that the input data are centered, and hence there’s no need to subtract the mean.

The size of the result depends on the size of v1 and v2. When both v1 and v2 are vectors, it returns the covariance between them as a scalar. When either one is a matrix, it returns a covariance matrix of size (n1,n2), where n1 and n2 are the numbers of slices in v1 and v2, which depend on the setting of vardim.

Note: v2 can be omitted, which indicates v2=v1.

cor(v1[, v2][, vardim=1, mean=nothing])

Compute the Pearson correlation between the vector(s) in v1 and v2.

Users can use the keyword argument vardim to specify the variable dimension, and mean to supply pre-computed mean values.

Signal Processing

Fast Fourier transform (FFT) functions in Julia are implemented by calling functions from FFTW. By default, Julia does not use multi-threaded FFTW. Higher performance may be obtained by experimenting with multi-threading. Use FFTW.set_num_threads(np) to use np threads.

fft(A[, dims])

Performs a multidimensional FFT of the array A. The optional dims argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most efficient if the size of A along the transformed dimensions is a product of small primes; see nextprod(). See also plan_fft() for even greater efficiency.

A one-dimensional FFT computes the one-dimensional discrete Fourier transform (DFT) as defined by

\[\operatorname{DFT}(A)[k] = \sum_{n=1}^{\operatorname{length}(A)} \exp\left(-i\frac{2\pi (n-1)(k-1)}{\operatorname{length}(A)} \right) A[n].\]

A multidimensional FFT simply performs this operation along each transformed dimension of A.

Higher performance is usually possible with multi-threading. Use FFTW.set_num_threads(np) to use np threads, if you have np processors.

fft!(A[, dims])

Same as fft(), but operates in-place on A, which must be an array of complex floating-point numbers.

ifft(A[, dims])

Multidimensional inverse FFT.

A one-dimensional inverse FFT computes

\[\operatorname{IDFT}(A)[k] = \frac{1}{\operatorname{length}(A)} \sum_{n=1}^{\operatorname{length}(A)} \exp\left(+i\frac{2\pi (n-1)(k-1)} {\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 in-place on A.

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 than ifft() 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 in-place on A.

plan_fft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Pre-plan an optimized FFT along given dimensions (dims) of arrays matching the shape and type of A. (The first two arguments have the same meaning as for fft().) Returns an object P which represents the linear operator computed by the FFT, and which contains all of the information needed to compute fft(A,dims) quickly.

To apply P to an array A, use P*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 the A for which the plan was created.) You can also apply a plan with a preallocated output array  by calling A_mul_B!(Â,plan,A). You can compute the inverse-transform plan by inv(P) and apply the inverse plan with P\ (the inverse plan is cached and reused for subsequent calls to inv or \), and apply the inverse plan to a pre-allocated output array A with A_ldiv_B!(A,P,Â).

The flags argument is a bitwise-or of FFTW planner flags, defaulting to FFTW.ESTIMATE. e.g. passing FFTW.MEASURE or FFTW.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 optional timelimit argument specifies a rough upper bound on the allowed planning time, in seconds. Passing FFTW.MEASURE or FFTW.PATIENT may cause the input array A to be overwritten with zeros during plan creation.

plan_fft!() is the same as plan_fft() but creates a plan that operates in-place on its argument (which must be an array of complex floating-point numbers). plan_ifft() and so on are similar but produce plans that perform the equivalent of the inverse transforms ifft() and so on.

plan_ifft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as plan_fft(), but produces a plan that performs inverse transforms ifft().

plan_bfft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as plan_fft(), but produces a plan that performs an unnormalized backwards transform bfft().

plan_fft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as plan_fft(), but operates in-place on A.

plan_ifft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as plan_ifft(), but operates in-place on A.

plan_bfft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Same as plan_bfft(), but operates in-place on A.

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 with fft(). If A 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 of A to transform, similar to fft(). Instead of (roughly) halving the first dimension of A in the result, the dims[1] dimension is (roughly) halved in the same way.

irfft(A, d[, dims])

Inverse of rfft(): for a complex array A, gives the corresponding real array whose FFT yields A in the first half. As for rfft(), dims is an optional subset of dimensions to transform, defaulting to 1:ndims(A).

d is the length of the transformed real array along the dims[1] dimension, which must satisfy div(d,2)+1==size(A,dims[1]). (This parameter cannot be inferred from size(A) since both 2*size(A,dims[1])-2 as well as 2*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 to bfft()), 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)

Pre-plan an optimized real-input FFT, similar to plan_fft() except for rfft() instead of fft(). The first two arguments, and the size of the transformed result, are the same as for rfft().

plan_brfft(A, d [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Pre-plan an optimized real-input unnormalized transform, similar to plan_rfft() except for brfft() instead of rfft(). The first two arguments and the size of the transformed result, are the same as for brfft().

plan_irfft(A, d [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf)

Pre-plan an optimized inverse real-input FFT, similar to plan_rfft() except for irfft() and brfft(), respectively. The first three arguments have the same meaning as for irfft().

dct(A[, dims])

Performs a multidimensional type-II discrete cosine transform (DCT) of the array A, using the unitary normalization of the DCT. The optional dims argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most efficient if the size of A along the transformed dimensions is a product of small primes; see nextprod(). See also plan_dct() for even greater efficiency.

dct!(A[, dims])

Same as dct!(), except that it operates in-place on A, which must be an array of real or complex floating-point values.

idct(A[, dims])

Computes the multidimensional inverse discrete cosine transform (DCT) of the array A (technically, a type-III DCT with the unitary normalization). The optional dims argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most efficient if the size of A along the transformed dimensions is a product of small primes; see nextprod(). See also plan_idct() for even greater efficiency.

idct!(A[, dims])

Same as idct!(), but operates in-place on A.

plan_dct(A[, dims[, flags[, timelimit]]])

Pre-plan an optimized discrete cosine transform (DCT), similar to plan_fft() except producing a function that computes dct(). The first two arguments have the same meaning as for dct().

plan_dct!(A[, dims[, flags[, timelimit]]])

Same as plan_dct(), but operates in-place on A.

plan_idct(A[, dims[, flags[, timelimit]]])

Pre-plan an optimized inverse discrete cosine transform (DCT), similar to plan_fft() except producing a function that computes idct(). The first two arguments have the same meaning as for idct().

plan_idct!(A[, dims[, flags[, timelimit]]])

Same as plan_idct(), but operates in-place on A.

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 and b to vector x, with an optional initial filter state vector si (defaults to zeros).

filt!(out, b, a, x[, si])

Same as filt() but writes the result into the out argument, which may alias the input x to modify it in-place.

deconv(b, a)

Construct vector c such that b=conv(a,c)+r. Equivalent to polynomial division.

conv(u, v)

Convolution of two vectors. Uses FFT algorithm.

conv2(u, v, A)

2-D convolution of the matrix A with the 2-D separable kernel generated by the vectors u and v. Uses 2-D FFT algorithm

conv2(B, A)

2-D convolution of the matrix B with the matrix A. Uses 2-D FFT algorithm

xcorr(u, v)

Compute the cross-correlation of two vectors.

The following functions are defined within the Base.FFTW module.

r2r(A, kind[, dims])

Performs a multidimensional real-input/real-output (r2r) transform of type kind of the array A, as defined in the FFTW manual. kind specifies either a discrete cosine transform of various types (FFTW.REDFT00, FFTW.REDFT01, FFTW.REDFT10, or FFTW.REDFT11), a discrete sine transform of various types (FFTW.RODFT00, FFTW.RODFT01, FFTW.RODFT10, or FFTW.RODFT11), a real-input DFT with halfcomplex-format output (FFTW.R2HC and its inverse FFTW.HC2R), or a discrete Hartley transform (FFTW.DHT). The kind argument may be an array or tuple in order to specify different transform types along the different dimensions of A; 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 for dims[i], with kind[end] being used for i>length(kind).

See also plan_r2r() to pre-plan optimized r2r transforms.

r2r!(A, kind[, dims])

Same as r2r(), but operates in-place on A, which must be an array of real or complex floating-point numbers.

plan_r2r(A, kind[, dims[, flags[, timelimit]]])

Pre-plan an optimized r2r transform, similar to Base.plan_fft() except that the transforms (and the first three arguments) correspond to r2r() and r2r!(), respectively.

plan_r2r!(A, kind[, dims[, flags[, timelimit]]])

Similar to Base.plan_fft(), but corresponds to r2r!().

Numerical Integration

Although several external packages are available for numeric integration and solution of ordinary differential equations, we also provide some built-in 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) from a to b, and optionally over additional intervals b to c and so on. Keyword options include a relative error tolerance reltol (defaults to sqrt(eps) in the precision of the endpoints), an absolute error tolerance abstol (defaults to 0), a maximum number of function evaluations maxevals (defaults to 10^7), and the order of the integration rule (defaults to 7).

Returns a pair (I,E) of the estimated integral I and an estimated upper bound on the absolute error E. If maxevals is not exceeded then E<=max(abstol,reltol*norm(I)) will hold. (Note that it is useful to specify a positive abstol in cases where norm(I) may be zero.)

The endpoints a etcetera can also be complex (in which case the integral is performed over straight-line segments in the complex plane). If the endpoints are BigFloat, then the integration will be performed in BigFloat 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 floating-point 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 a norm (i.e., any normed vector space). Alternatively, a different norm can be specified by passing a norm-like function as the norm keyword argument (which defaults to vecnorm).

[Only one-dimensional integrals are provided by this function. For multi-dimensional integration (cubature), there are many different algorithms (often much better than simple nested 1d integrals) and the optimal choice tends to be very problem-dependent. See the Julia external-package 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 Gauss-Kronrod 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 at x=0.7 and you want to integrate from 0 to 1, you should use quadgk(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, a log(x) or 1/sqrt(x) singularity).

For real-valued 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.)