OCaml Scientific Computing

1st Edition (in progress)
Table of Contents

Mathematical Functions


Basic Functions

Note that functions in this chapter works on scalar values. The N-dimensional array module introduced in later chapters contains these basic functions that work on n-dimensional arrays, including vectors and matrices.

Basic Unary Math Functions

Many basic math functions takes one float number as input and returns one float number. We call them unary functions. You can use these unary functions easily from the Maths module. For example:

Maths.sqrt 2. ;;
>- : float = 1.41421356237309515
Table 1: Basic unary math functions
Function Explanation
abs |x|
neg -x
reci 1/x
floor the largest integer that is smaller than x
ceil the smallest integer that is larger than x
round rounds x towards the bigger integer when on the fence
trunc integer part of x
sqr \(x^2\)
sqrt \(\sqrt{x}\)

Basic Binary Functions

Binary functions takes two floats as inputs and returns one float as return. The most common arithmetic functions belong to this category.

Table 2: Binary math functions
Function Explanation
add x + y
sub x - y
mul x * y
div x / y
fmod x % y
pow \(x^y\)
hypot \(\sqrt{x^2 + y^2}\)
atan2 returns \(\arctan(y/x)\), accounting for the sign of the arguments; this is the angle to the vector \((x, y)\) counting from the x-axis.

Exponential and Logarithmic Functions

The constant \(e = \sum_{n=0}^{\infty}\frac{1}{n!}\) is what we called the “natural constant”. It is called this way because the exponential function and it inverse function logarithm are so frequently used in nature and our daily life: logarithmic spiral, population growth, carbon date ancient artifacts, computing bank investments, etc.

We also have this beautiful Euler’s formula that connects the two most frequently used constants and the base of complex number and natural numbers:

\[e^{i\pi}+ 1=0.\]

As an example, in a scientific experiment about bacteria, we can assume the number of bacterial follows an exponential function \(n(t) = Ce^rt\) where \(C\) is the initial population and \(r\) is the daily increase rate. With this model, we can predict how the population of bacterial grows within certain time.

The full list of exponential and logarithmic functions, together with some variants, are presented in tbl. 3.

Table 3: Exponential and logarithmic math functions
Function Explanation
exp exponential \(e^x\)
exp2 \(2^x\)
exp10 \(10^x\)
expm1 returns \(\exp(x) - 1\) but more accurate for \(x \sim 0\)
log \(log_e~x\)
log2 \(log_2~x\)
log10 \(log_10~x\)
logn \(log_n~x\)
log1p Inverse of expm1
logabs \(\log(|x|)\)
xlogy \(x \log(y)\)
xlog1py \(x \log(y+1)\)
logit \(\log(p/(1-p))\)
expit \(1/(1+\exp(-x))\)
log1mexp \(log(1-exp(x))\)
log1pexp \(log(1+exp(x))\)

Trigonometric Functions

In mathematics, the trigonometric functions are real functions which relate an angle of a right-angled triangle to ratios of two side lengths. They are widely used in all sciences that are related to geometry, such as navigation, solid mechanics, celestial mechanics, geodesy, and many others. They are among the simplest periodic functions, and as such are also widely used for studying periodic phenomena, through Fourier analysis. The most widely used trigonometric functions are the sine, the cosine, and the tangent. Their reciprocals are respectively the cosecant, the secant, and the cotangent, which are less used in modern mathematics. (COPY)

The triangular functions are all unary functions, for example:

Maths.sin (Owl_const.pi /. 2.);;
>- : float = 1.

And they are all included in the math module in Owl, as shown in tbl. 4.

Table 4: Trigonometric math functions
Function Explanation Derivatives Taylor Expansion
sin \(\sin(x)\) \(\cos(x)\) \(\sum_{n=1}(-1)^{n+1}\frac{x^{2n+1}}{(2n+1)!}\)
cos \(\cos(x)\) \(-\sin(x)\) \(\sum_{n=1}(-1)^n\frac{x^{2n}}{(2n)!}\)
tan \(\tan(x)\) \(1 + \tan^2(x)\) \(\sum_{n=1}\frac{4^n(4^n-1)B_n~x^{2n-1}}{(2n)!}\)
cot \(1/\tan(x)\) \(-(1 + \textrm{cot}^2(x))\) \(\sum_{n=0}\frac{E_n~x^{2n}}{(2n)!}\)
sec \(1/\cos(x)\) \(\textrm{sec}(x)\tan(x)\) \(\sum_{n=0}\frac{2(2^{2n-1})B_n~x^{2n-1}}{(2n)!}\)
csc \(1/\sin(x)\) \(-\textrm{csc}(x)\textrm{cot}(x)\) \(\frac{1}{x}-\sum_{n=1}\frac{4^n~B_n~x^{2n-1}}{(2n)!}\)

Here the \(B_n\) is the \(n\)th Bernoulli number, and \(E_n\) is the \(n\)th Euler number. The fig. 1 shows the relationship between these trigonometric functions (figure src). These functions also have corresponding inverse functions: asin, acos, atan, acot, asec, acsc. For example, if \(\sin(a) = b\), then \(\textrm{asin}(b) = a\).

Figure 1: Relationship between different trigonometric functions

Another related idea is the Hyperbolic functions. Hyperbolic functions are analogous of the ordinary trigonometric functions defined for the hyperbola rather than on the circle: just as the points (cos t, sin t) form a circle with a unit radius, the points (cosh(x), sinh(x)) form the right half of the equilateral hyperbola. Hyperbolic functions occur in the solutions of many linear differential equations, calculations of angles and distances in hyperbolic geometry, and of Laplace’s equation in Cartesian coordinates. (COPY)

These functions in Owl is shown below:

  • sinh: \(\frac{e^x - e^{-x}}{2}\), derivative is \(\cosh(x)\), and taylor expansion is \(\sum_{n=0}\frac{x^{2n+1}}{(2n+1)!}\).
  • cosh: \(\frac{e^x + e^{-x}}{2}\), derivative is \(\sinh(x)\), and taylor expansion is \(\sum_{n=0}\frac{x^{2n+1}}{(2n+1)!}\).
  • tanh: \(\frac{\sinh{x}}{\cosh{x}}\), derivative is \(1-\tanh^2(x)\), and taylor expansion is \(\sum_{n=1}\frac{4^n(4^n-1)B_{2n}~x^{2n-1}}{(2n)!}\).
  • coth: \(\frac{\cosh{x}}{\sinh{x}}\), derivative is \(1-\coth^2(x)\), and taylor expansion is \(\frac{1}{x}-\sum_{n=1}\frac{4^n~B_{2n}~x^{2n-1}}{(2n)!}\).
  • sech: \(1/\cosh(x)\), derivative is \(-\tanh(x)/\cosh(x)\), and taylor expansion is \(\sum_{n=0}\frac{E_{2n}~x^{2n}}{(2n)!}\).
  • csch:\(1/\sinh(x)\), derivative is \(-\coth(x)/\sinh(x)\), and taylor expansion is \(\frac{1}{x}+\sum_{n=1}\frac{2(1-2^{2n-1})B_{2n}~x^{2n-1}}{(2n)!}\).

(TODO: Change these information to table; beware that this table would lead to page overflow.)

Similarly, each of these functions has corresponding inverse functions: asinh, acosh, atanh, acoth, asech, acsch. The relationship between these hyperbolic trigonometric functions are clearly depicted in fig. 2.

Figure 2: Relationship between different hyperbolic trigonometric functions

Besides these functions, there are also some related functions. sinc returns \(\sin(x)/x\) and 1 for \(x=0\). logsinh returns \(\log(\sinh(x))\) but handles large \(|x|\). logcosh returns \(\log(\cosh(x))\) but handles large \(|x|\). sindg/cosdg/tandg/cotdg are the sine/cosine/tangent/cotangent of the angle given in degrees.

Other Math Functions

There are some other function that may not be very common in traditional math. Functions such as sigmoid and relu are frequently used in Deep Learning as the activation functions in a neural network. The activation functions are crucial to the neural network regarding various aspects, including output result, accuracy, convergence speed, etc.

  • sigmoid x: \(1 / (1 + \exp(-x))\)
  • signum x: returns the sign of x: -1, 0, or 1
  • softsign x: smoothed sign function
  • relu x: \(\max(0, x)\)

Special Functions

The definition of numerous special functions of mathematical physics. Special functions are particular mathematical functions that have more or less established names and notations due to their importance in mathematical analysis, functional analysis, physics, or other applications.(COPY) It interfaces to the Cephes Mathematical Functions Library.

Airy Functions

In the physical sciences, the Airy function (or Airy function of the first kind) Ai(x) is a special function named after the British astronomer George Biddell Airy. (COPY) It is the solution of differential equation:

\[y''(x) = xy(x).\]

This differential equation has two linearly independent solutions Ai and Bi. Owl provides the airy to do that:

val airy : float -> float * float * float * float

The four returned numbers are Ai, its derivative Ai', Bi, and its derivative Bi'. Let’s look at an example.

let x = Mat.linspace (-15.) 5. 200

let y0 = Mat.map (fun x -> 
    let ai, _, _, _ = Maths.airy x in ai
) x 

let y1 = Mat.map (fun x -> 
    let _, _, bi, _ = Maths.airy x in bi
) x 

let _ = 
  let h = Plot.create "special_airy.png" in
  Plot.(plot ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2. ] x y0);
  Plot.(plot ~h ~spec:[ RGB (219, 68,  55); LineStyle 2; LineWidth 2. ] x y1);
  Plot.(set_yrange h (-0.5) 1.);
  Plot.(legend_on h ~position:SouthEast [|"Ai"; "Bi"|]);
  Plot.output h
Figure 3: Examples of the two solutions of an Airy equation

APPLICATION description

Bessel Functions

Bessel functions, first defined by the mathematician Daniel Bernoulli and then generalized by Friedrich Bessel, are canonical solutions y(x) of Bessel’s differential equation:

\[x^2y''+xy'+(x^2 - \alpha^2)y = 0.\]

The complex number \(\alpha\) is called the order of the bessel function.

The Bessel functions can be divided into two “kinds”. Bessel functions of the first kind \(J\) are solutions of Bessel’s differential equation that are finite at \(x=0\) integer or positive order and diverge as \(x\) approaches zero for negative non-integer order. Bessel functions of the second kind are solutions of the Bessel differential equation that have a singularity at \(x=0\) and are multivalued. (COPY ALERT)

A special case is when \(x\) is purely imaginary. In this case, the solutions to the Bessel equation are called the modified Bessel functions. These modified Bessel functions can also be categorised as first kind and second kind.

Based on these category, Owl provides these functions.

Table 5: Bessel functions
Function Explanation
j0 x Bessel function of the first kind of order 0
j1 x Bessel function of the first kind of order 1
jv x y Bessel function of the first kind of real order
y0 x Bessel function of the second kind of order 0
y1 x Bessel function of the second kind of order 1
yv x y Bessel function of the second kind of real order
yn a x Bessel function of the second kind of integer order
i0 x Modified Bessel function of order 0
i1 x Modified Bessel function of order 1
iv x y Modified Bessel function of real order
i0e x Exponentially scaled modified Bessel function of order 0
i1e x Exponentially scaled modified Bessel function of order 1
k0 x Modified Bessel function of the second kind of order 0
k1 Modified Bessel function of the second kind of order 1
k0e Exponentially scaled modified Bessel function of the second kind of order 0
k1e Exponentially scaled modified Bessel function of the second kind of order 1

Let’s look at one example.

let x = Mat.linspace (0.) 20. 200

let y0 = Mat.map Maths.j0 x

let y1 = Mat.map Maths.j1 x

let y2 = Mat.map (Maths.jv 2.) x

let _ =
  let h = Plot.create "example_bessel.png" in
  Plot.(plot ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2. ] x y0);
  Plot.(plot ~h ~spec:[ RGB (219, 68,  55); LineStyle 2; LineWidth 2. ] x y1);
  Plot.(plot ~h ~spec:[ RGB (244, 180,  0); LineStyle 3; LineWidth 2. ] x y2);
  Plot.(legend_on h ~position:NorthEast [|"j0"; "j1"; "j2"|]);
  Plot.output h
Figure 4: Examples of Bessel function of the first kind, with different order

(More examples can be added if we want to expand)

Bessel’s equation arises when finding separable solutions to Laplace’s equation and the Helmholtz equation in cylindrical or spherical coordinates. Bessel functions are therefore especially important for many problems of wave propagation and static potentials. In solving problems in cylindrical coordinate systems, one obtains Bessel functions of integer order or half integer order. For example, electromagnetic waves in a cylindrical waveguide, pressure amplitudes of inviscid rotational flows, heat conduction in a cylindrical object, etc. (COPY ALERT)

Elliptic Functions

Table 6: Elliptic functions
Function Explanation
ellipj u m Jacobian elliptic functions of parameter m between 0 and 1, and real argument u
ellipk m Complete elliptic integral of the first kind
ellipkm1 p Complete elliptic integral of the first kind around m = 1
ellipkinc phi m Incomplete elliptic integral of the first kind
ellipe m Complete elliptic integral of the second kind
ellipeinc phi m Incomplete elliptic integral of the second kind

The Jacobian elliptic functions are found in the description of the motion of a pendulum, as well as in the design of the electronic elliptic filters. These functions are periodic, with quarter-period on the real axis equal to the complete elliptic integral. There are twelve Jacobi elliptic functions and ellipj returns three of them: sn, cn, dn. And the fourth result phi is called the amplitude of u. (COPY)

Elliptic integrals arose from the attempts to find the perimeter of an ellipse. elliptic integral. A Elliptic integral function can be expressed in the form of: \[f(x)=\int_c^xR(t, \sqrt(P(t)))dt,\] where \(R\) is a rational function of its two arguments, \(P\) is a polynomial of degree 3 or 4 with no repeated roots, and \(c\) is a constant. Incomplete elliptic integrals are functions of two arguments; complete elliptic integrals are functions of a single argument. (COPY)

In general, integrals in this form cannot be expressed in terms of elementary functions. Exceptions to this general rule are when P has repeated roots, or when \(R(x,y)\) contains no odd powers of y. However, with the appropriate reduction formula, every elliptic integral can be brought into a form that involves integrals over rational functions and the three Legendre canonical forms (i.e. the elliptic integrals of the first, second and third kind). (COPY)

We can use ellipe to compute the circumference of an ellipse. To compute that requires calculus, and the elliptic functions provides a solution. Suppose an ellipse has semi-major axis \(a=4\) and semi-minor axis \(b=3\). We an compute its circumference using \(4a\textrm{ellipe}(1 - \frac{b^2}{a^2})\).

let a = 4. ;;
>val a : float = 4.
let b = 3.;;
>val b : float = 3.
let c = 4. *. a *. Maths.(ellipe (1. -. pow (b /. a) 2.));;
>val c : float = 22.1034921607095072

Gamma Functions

For a positive integer n, the Gamma function is the factorial function.

\[\Gamma(n) = (n-1)!\]

For a complex numbers \(z\) with a positive real part,

\[\Gamma(z) = \int_0^{\infty}x^{z-1}e^{-x}dx.\]

The Gamma function is widely used in a range of areas such as fluid dynamics, geometry, astrophysics, etc. It is especially suitable for describing a common pattern of processes that decay exponentially in time or space. The Gamma function and related function provided in Owl are list in tbl. 7.

Table 7: Gamma functions
Function Explanation
gamma z Returns the value of the Gamma function
rgamma z Reciprocal of the Gamma function
loggamma z Principal branch of the logarithm of the Gamma function
gammainc a x Regularized lower incomplete gamma function
gammaincinv a y Inverse function of gammainc
gammaincc a x Complemented incomplete gamma integral
gammainccinv a y Inverse function of gammaincc
psi z The digamma function

The incomplete gamma functions are similarly to the gamma function but with different or “incomplete” integral limits. The gamma function is defined as an integral from zero to infinity. This contrasts with the lower incomplete gamma function, which is defined as an integral from zero to a variable upper limit. Similarly, the upper incomplete gamma function is defined as an integral from a variable lower limit to infinity. The digamma function is defined as the logarithmic derivative of the gamma function. (COPY)

Here is an example of using gamma.

let x = Mat.linspace (-3.5) 5. 2000

let y = Mat.map Maths.gamma x

let _ =
  let h = Plot.create "example_gamma.png" in
  Plot.(plot ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2. ] x y);
  Plot.(set_yrange h (-10.) 20.);
  Plot.output h
Figure 5: Examples of Gamma function along part of the real axis

(TODO: this figure should not have the vertical lines)

Beta Functions

Beta function is defined as:

\[B(x,y) = \int_0^1t^{x-1}(1-t)^{y-1}dt = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\]

The incomplete beta function extends this definition to:

\[B(x, a, b) = \int_0^xt^{a-1}(1-t)^{b-1}dt.\]

They are both included in the special functions provided by Owl.

Table 8: Beta functions
Function Explanation
beta x y Beta function
betainc a b x Incomplete Beta integral
betaincinv a b y Inverse function of betainc

The Beta function has several properties:

let x = Maths.beta 3. 4. ;;
>val x : float = 0.0166666666666666664
let y = Maths.((gamma 3.) *. (gamma 4.) /. (gamma (7.))) ;;
>val y : float = 0.0166666666666666664

This validate the relationship between beta funtion and gamma function. Another property of beta function is it is symmetric, which means \(B(x,y) = B(y, x)\).

let x = Maths.beta 3. 4. ;;
>val x : float = 0.0166666666666666664
let y = Maths.beta 4. 3. ;;
>val y : float = 0.0166666666666666664

Beta function is the first known scattering amplitude in String theory in physics. It can also be used to model a preferential attachment process, which describes the distribution of resources among individuals based on the resource amount they already have. (COPY)

Struve Functions

The Struve function is defined as: \[H_v(x) = (z/2)^{v + 1} \sum_{n=0}^\infty \frac{(-1)^n (z/2)^{2n}}{\Gamma(n + \frac{3}{2}) \Gamma(n + v + \frac{3}{2})},\]

where \(\Gamma\) is the gamma funcction. \(x\) must be positive unless \(v\) is an integer. The function struve v x returns the value of Struve function. The paramter \(v\) is called the order of this Struve function. Here is an example.

let _ =
  let h = Plot.create "example_struve.png" in
  Plot.(plot_fun ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2.] (Maths.struve 0.) (-12.) 12.);
  Plot.(plot_fun ~h ~spec:[ RGB (219, 68,  55); LineStyle 2; LineWidth 2.] (Maths.struve 1.) (-12.) 12.);
  Plot.(plot_fun ~h ~spec:[ RGB (244, 180,  0); LineStyle 3; LineWidth 2.] (Maths.struve 2.) (-12.) 12.);
  Plot.(plot_fun ~h ~spec:[ RGB (77,  81,  57); LineStyle 1; LineWidth 2.] (Maths.struve 3.) (-12.) 12.);
  Plot.(plot_fun ~h ~spec:[ RGB (111, 51, 129); LineStyle 2; LineWidth 2.] (Maths.struve 4.) (-12.) 12.);
  Plot.(set_yrange h (-3.) 5.);
  Plot.(legend_on h ~position:SouthEast [|"H0"; "H1"; "H2"; "H3"; "H4"|]);
  Plot.output h
Figure 6: Examples of Struve function for different orders.

Struve functions have some specific uses across many different fields of physics in a wide variety of applications. For example, they can be found in water-wave and surface-wave problems (specifically flow of liquid near a turning ship) as well as calculations to do with the distribution of fluid pressure over a vibrating disk and other unsteady aerodynamics. They also crop up when considering aspects of optical diffraction, plasma stability (specifically resistive magnetohydrodynamics instability theory), quantum dynamical studies of spin decoherence and excitation in carbon nanotubes. (COPY)

Zeta Functions

The Hurwitz zeta function zeta x q returns the Hurwitz zeta function:

\[\zeta(x, q) = \sum_{k=0}^{\infty}\frac{1}{(k+q)^x}.\]

When \(q\) is set to 1, this function is reduced to Riemann zeta function. The function zetac x returns Riemann zeta function minus 1. We can evaluate the zeta function at certain points, for example:

Maths.zeta 4. 1. ;;
>- : float = 1.08232323371113837
(Maths.pow Owl_const.pi 4.) /. 90.;;
>- : float = 1.08232323371113792

The Riemann zeta function plays a pivotal role in analytic number theory and has applications in physics, probability theory, and applied statistics. Zeta function regularization is used as one possible means of regularization of divergent series and divergent integrals in quantum field theory. In one notable example, the Riemann zeta-function shows up explicitly in one method of calculating the Casimir effect. The zeta function is also useful for the analysis of dynamical systems. (COPY)

Error Functions

The error functions are not about error processing in programming. In mathematics, it is defined as: \[\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2})dt.\]

Table 9: Error functions
Function Explanation
erf x Error function
erfc x Complementary error function: \(1 - \textrm{erf}(x)\)
erfcx x Scaled complementary error function: \(\exp(x^2) \mathrm{erfc}(x)\)
erfinv x Inverse function of erf
erfcinv x Inverse function of erfc

The error function is a sigmoid function. We can observe its shape by the code below.

let _ =
  let h = Plot.create "example_erf.png" in
  Plot.(plot_fun ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2.] Maths.erf (-3.) 3.);
  Plot.output h
Figure 7: Plot of the Error function.

The error function occurs often in probability, statistics, and partial differential equations describing diffusion. In statistics, for nonnegative values of x, the error function has the following interpretation: for a random variable Y that is normally distributed with mean 0 and variance 0.5, then erf x is the probability that Y falls in the range [-x, x]. (COPY)

Integral Functions

Owl also provides several special integral functions. The Dawson function is defined as: \[D(x) = e^{-x^2}\int_0^x~e^{t^2}dt\]

And the Fresnel trigonometric integral returns a tuple that contains two parts: \[S(x) = \int_0^x~sin(t^2)dt, C(x) = \int_0^x~cos(t^2)dt.\]

We can observe the functions of these integrals with plots.

let _ =
  let h = Plot.create ~m:1 ~n:2 "example_integrals.png" in
  Plot.subplot h 0 0;
  Plot.(plot_fun ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2.] Maths.dawsn (-5.) 5.);
  Plot.set_ylabel h "dawsn(x)";
  Plot.subplot h 0 1;
  Plot.(plot_fun ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2.] (fun x -> let s, _ = Maths.fresnel x in s) 0. 5.);
  Plot.(plot_fun ~h ~spec:[ RGB (219, 68,  55); LineStyle 2; LineWidth 2.] (fun x -> let _, c = Maths.fresnel x in c) 0. 5.);
  Plot.(legend_on h ~position:SouthEast [|"S(x)"; "C(x)"|]);
  Plot.set_ylabel h "fresnel(x)";
  Plot.output h
Figure 8: Plot of the Dawson and Fresnel integral function.

Besides these two, other type of special integral functions are also provided, as shown in tbl. 10.

Table 10: Integral functions
Function Explanation
expn n x Generalized exponential integral \(E_n(x) = x^{n-1}\int_x^{\infty}\frac{e^{-t}}{t^n}dt\)
shi x Hyperbolic sine integral: \(\int_0^x~\frac{\sinh~t}{t}dt\)
chi x Hyperbolic cosine integral: \(\gamma + \log(x) + \int_0^x~\frac{\cosh~t -1}{t}dt\)
shichi x (shi x, chi x)
si x Sine integral: \(\int_0^x~\frac{\sin~t}{t}dt\)
ci x Cosine integral: \(\gamma + \log(x) + \int_0^x~\frac{\cos~t -1}{t}dt\)
sici x (si x, ci x)

Dawson integrals is motivated by research on the electromagnetic radiation propagation across the surface of earth. The Fresnel integrals were originally used in the calculation of the electromagnetic field intensity in an environment where light bends around opaque objects. More recently, they have been used in the design of highways and railways, specifically their curvature transition zones. Other applications are roller coasters or calculating the transitions on a velodrome track to allow rapid entry to the bends and gradual exit. (COPY)


The definition of factorials is simple:

\(F(n) = n! = n \times (n - 1) \times (n-2) \ldots \times 1\)

The factorial function, together with several variants, are contained in the math module.

Table 11: Factorial functions
Function Explanation
fact n Factorial function \(!n\)
log_fact n Logarithm of factorial function
doublefact n Double factorial function calculates \(n!! = n(n-2)(n-4)\dots 2\) (or 1)
log_doublefact n Logarithm of double factorial function

The factorial functions accepts integer as input, for example:

Maths.fact 5 ;;
>- : float = 120.

The factorials are applied in many areas of mathematics, most notably the combinatorics. The permutation and combination are both defined in factorials. The permutation returns the number \(n!/(n-k)!\) of ordered subsets of length \(k\), taken from a set of \(n\) elements. THe combination returns the number \({n\choose k} = n!/(k!(n-k)!)\) of subsets of \(k\) elements of a set of \(n\) elements. tbl. 12 provides the combinatorics functions you can use in the math module.

Table 12: Permutation and combination functions
Function Explanation
permutation n k Permutation number
permutation_float n k Similar to permutation but deals with larger range and returns float
combination n k Combination number
combination_float n k Similar to combination but deals with larger range and returns float
log_combination n k Returns the logarithm of \({n\choose k}\)

We can see a simple example.

let x = Maths.combination 10 2;;
>val x : int = 45
let y = Maths.combination_float 10 2;;
>val y : float = 45.

Interpolation and Extrapolation

Sometimes we don’t know the full description of a function \(f\), but only some points on it, and therefore we cannot calculate its value at an aribitrary point. The target is to esimate the \(f(x)\) for an arbitrary \(x\) by drawing a smooth curve through the given data. If \(x\) is within the range of the given data, this taks is called interpolation, otherwise it’s called extrapolation, which is much more difficult to do.

The Owl_maths_interpolate module provides an polint function for interpolation and extrapolation:

val polint : float array -> float array -> float -> float * float

polint xs ys x performs polynomial interpolation of the given arrays xs and ys. Given arrays \(xs[0 \ldots (n-1)]\) and \(ys[0\ldots~(n-1)]\), and a value x. The function returns a value y, and an error estimate dy. The paramter xs is an array of input x values of P(x), and ys is an array of corresponding y values of P(x). It returns (y', dy) wherein y' is the returned value y' = P(x), and dy is the estimated error.

As its name suggests, the polint approximate complicated curves with polynomial of lowest possible degree that passes the given points. We can show how this interplation method works for an example. In the previous chapter we have introduced that the Gamma function is actually a interpolation solution to the integer function \(y(x) = (n-1)!\). So we can specify five nodes on a plane that are generated from this factorial functions.

let x = [|2; 3; 4; 5; 6|];;
>val x : int array = [|2; 3; 4; 5; 6|]
let y = Array.map (fun x -> Maths.fact (x - 1)) x ;;
>val y : float array = [|1.; 2.; 6.; 24.; 120.|]
let x = Array.map float_of_int x;;
>val x : float array = [|2.; 3.; 4.; 5.; 6.|]

Now we can define the interpolation function f that accept on float number and returns another float number. Also we convert the given data \(x\) and \(y\) into matrix format for plotting purpose.

let f a = 
  let v, _ = Owl_maths_interpolate.polint x y a in

let xm = Mat.of_array x 1 5
let ym = Mat.of_array y 1 5

Now we can plot the interpolation function. We compare it to the Gamma function. As can be seen in fig. 9, both lines cross the given nodes. We can see that the interpolated line fits well with the “true interpolation”, i.e. the Gamma function. However, the extrapolation fitting where the x-value falls out of given data, is less than ideal.

let _ =
  let h = Plot.create "interp.png" in
  Plot.(plot_fun ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2.] f 2. 6.5);
  Plot.(plot_fun ~h ~spec:[ RGB (219, 68,  55); LineStyle 2; LineWidth 2.] Maths.gamma 2. 6.5);
  Plot.(scatter ~h ~spec:[ Marker "#[0x229a]"; MarkerSize 5. ] xm ym);
  Plot.(legend_on h ~position:NorthWest [|"Interpolation"; "Gamma function"; "Given values"|]);
  Plot.output h
Figure 9: Plot of interpolation and corresponding Gamma function.


Given a function \(f\) that accepts a real variable and an interval \([a, b]\) of the real line, the integral of this function


can be thought of as the sum of signed area of the region in the cartesian plane that is bounded by the curve of f, the x-axis within the x-axis range \([a, b]\). The area above the x-axis adds to the sum and that below the x-axis subtracts from the area sum.

Owl provides several neumerical routines to help you to do integrations in Owl_maths_quandrature module. For example, we can compute \(\int_1^4x^2\) with the code below:

Owl_maths_quadrature.trapz (fun x -> x ** 2.) 1. 4.;;
>- : float = 21.0001344681758439

We can verify this result using the fundamental theorem of calculus:

\[\int_1^4x^2 = (4^3 -1^3) / 3 = 21\].

So you might be thinking, what is this trapz? Why the result is not exactly 21?

Using numerical methods (or quadrature) to do integration dates back to the invention of calculus or even earlier. The basic idea is to use summation of small areas to approximate that of an integration, as shown in fig. 10 (src).

Figure 10: Basic method of numerical integration

There exists a lot of algorithms to do numerical integration, and using the trapezoial rule is one of them. This classical method divide a to b into \(N\) equally spaced abscissas: \(x_0, x_1, \ldots, x_N\). Each area between \(x_i\) and \(x_j\) is seen as an Trapezoid and the area formula is computed as:

\[\int_{x_0}^{x_1}f(x)dx = h(\frac{f(x_0)}{2} + \frac{f(x_1)}{2}) + O(h^3f'').\]

Here the error term \(O(h^3f'')\) indicated that the error of approximation is related with that of abscissas size \(h\) and second order derivative of the original function.

Function trapz implements this method. It’s interface is:

val trapz : ?n:int -> ?eps:float -> (float -> float) -> float -> float -> float

trapz ~n ~eps f a b computes the integral of f on the interval [a,b] using the trapezoidal rule. It works by iterating for several stages, each stage improving the accuracy by adding more interior points. The argument \(n\) specifies the maximum step which defaults to 20, and eps is the desired fractional accuracy threshold, which defaults to 1e-6.

The other methods are similar to trapz in interface, only different in implementation. For example, the simpson uses the Simpson formula:

\[\int_{x_0}^{x_2}f(x)dx = h(\frac{f(x_0)}{3} + \frac{4f(x_1)}{3} + \frac{f(x_2)}{3}) + O(h^5f(4)).\]

Then there is the Romberg integration (romberg) that can choose methods of different orders to give good accuracy, and the algorithms is normally much faster than the trapz and simpson methods. Moreover, if the abscissas can be varied, then there is the adaptive Gaussian quadrature of fixed tolerance gaussian and Gaussian quadrature of fixed order gaussian.

As an example, we can compute the special integral function \(Si(x)=\int_0^x\frac{sin(t)}{t}dt\) from previous section using the numerical integration method. Let’s set \(x=4\).

let f t = Maths.(div (sin t) t);;
>val f : float -> float = <fun>
Owl_maths_quadrature.gaussian f 0. 4.;;
>- : float = 1.75820313914469306
Owl_maths.si 4.;;
>- : float = 1.75820313894905289

We can see the numerical method gaussian works well to approximate this special integral function.

Utility Functions

Besides what we have mentioned, there are also some utitlity functions that worth mentioning.

A prime number is a natural number greater than 1 that cannot be formed by multiplying two smaller natural numbers. The is_prime checks if an integer is a prime number. This function is deterministic for all numbers representable by an int. It is implemented using the Miller-Rabin primality test method.

Maths.is_prime 997;;
>- : bool = true

Primes are used in several routines in information technology, such as public-key cryptography, which relies on the difficulty of factoring large numbers into their prime factors. In abstract algebra, objects that behave in a generalized way like prime numbers include prime elements and prime ideals.(COPY)

Another number theory related idea is the Fermat’s factorization, which represents an odd integer as the difference of two squares: \(N = a^2 - b^2\), and therefore N can be factorised as \((a+b)(a-b)\). The function fermat_fact performs Fermat factorisation over odd number N, i.e. into two roughly equal factors \(x\) and \(y\) so that \(N=x\times~y\).

Maths.fermat_fact 6557;;
>- : int * int = (83, 79)
83 * 79;;
>- : int = 6557

Next two functions concerns the precision of float numbers in computer.

TODO: Explain the mechansim of float number in a computer.

nextafter from to returns the next representable double precision value of from in the direction of to. If from equals to, this value is returned. The other is nextafterf. nextafter from to returns the next representable single precision value of from in the direction of to. If from equals to, this value is returned. For example:

Maths.nextafterf 1. 2.;;;;
>- : float = 1.00000011920928955
Maths.nextafter 1. 2.;;;;
>- : float = 1.00000000000000022
Maths.nextafter 1. 0.;;;;
>- : float = 0.999999999999999889

Next: Chapter 05Statistical Functions