Elementary Math, Polynomials and Interpolation

This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production.

The reference material consulted while preparing these pages is listed in the section bibliography.

Written against MATLAB R2026b. The sessions in the "Octave 11.3.0 output" blocks were run in GNU Octave 11.3.0. Octave rejects the digits argument of round, and it has no makima, griddedInterpolant or scatteredInterpolant, so those examples are checked against the R2026b reference pages and labelled "not run here". Everything on this page is base MATLAB; no toolbox is needed (see Toolboxes and Add-On Products for the boundary). Floating-point background is on Numeric Types and Floating Point, and the heatflow scenario is described on the section home page.

The documentation groups these functions under "Elementary Math" (arithmetic, exponents, discrete math, polynomials, special functions) and "Interpolation".

Rounding

Function Rounds toward

round(x)

The nearest integer; halves go away from zero.

floor(x)

Minus infinity.

ceil(x)

Plus infinity.

fix(x)

Zero.

x = [-7.5 -2.5 2.5 7.5];
disp(round(x)); disp(floor(x)); disp(ceil(x)); disp(fix(x))

Octave 11.3.0 output (MATLAB layout may differ):

  -8  -3   3   8
  -8  -3   2   7
  -7  -2   3   8
  -7  -2   2   7

Because halves go away from zero, round(2.5) is 3 and round(-2.5) is -3, unlike the "round half to even" rule of some other languages.

Rounding to digits (R2014b). round(x,n) rounds to n digits: to the right of the decimal point for a positive n, to the nearest integer for n = 0, and to the left of the decimal point (tens, hundreds and so on) for a negative n. round(x,n,"significant") rounds to n significant digits instead, and n must then be a positive integer; "decimals" is the default type. Since R2022a, TieBreaker chooses what happens to a half: "fromzero" (the default), "tozero", "even", "odd", "plusinf" or "minusinf".

round(pi,3)                        % 3.1420 (displayed with four decimals)
round(863178137,-2)                % 863178100
round(3132,2,"significant")        % 3100
round(2.5,TieBreaker="even")       % 2

Not run here (Octave 11.3.0 accepts only round(x)); the first three results are the examples of the reference page. A decimal value such as 2.015 is not exact in binary, so round(2.015,2) follows the stored value; round for display with fprintf or compose and keep the full value for computation.

mod and rem

mod(a,b) takes the sign of the divisor b; rem(a,b) takes the sign of the dividend a. For positive operands they agree.

disp(mod(-7,3)); disp(rem(-7,3)); disp(mod(7,-3)); disp(rem(7,-3))
disp(mod(5,0)); disp(rem(5,0))
disp(mod(10,2.5)); fprintf('%.17g\n', mod(5.1,0.1))

Octave 11.3.0 output:

2
-1
-2
1
5
NaN
0
0.099999999999999728

mod(5,0) is 5 and rem(5,0) is NaN, which is what the MATLAB reference pages state as the convention: mod(a,0) returns a, and rem(a,0) is NaN. The last line shows the floating-point trap: 5.1 and 0.1 are not exactly representable, so the remainder is almost the divisor instead of 0. Use a tolerance, or work in integers, when testing divisibility of non-integer values.

Number theory and counting

[g,c,d] = gcd(48,18);      % g = c*48 + d*18
disp([g c d])
disp(lcm(4,6))
disp(primes(30))
disp(isprime([1 2 9 11 97]))
disp(factor(84))
disp(factorial(5))
disp(nchoosek(5,2))
disp(nchoosek(1:4,2))      % all 2-element combinations of 1:4, one per row

Octave 11.3.0 output:

   6  -1   3
12
    2    3    5    7   11   13   17   19   23   29
  0  1  0  1  1
   2   2   3   7
120
10
   1   2
   1   3
   1   4
   2   3
   2   4
   3   4

gcd has a second form returning the Bezout coefficients, so 6 = (-1)*48 + 3*18. isprime accepts real, nonnegative integers of the integer and floating-point types. factorial(n) is a double and exact only while the result fits 53 bits: factorial(20) is 2.43290200817664e+18 and factorial(171) overflows to Inf (Octave 11.3.0 run; the double limit is the same in MATLAB). For the count form nchoosek(n,k) the reference page warns that a sufficiently large result is accurate only to 15 digits for double inputs (8 for single) and that MATLAB then displays a warning; the combinations form nchoosek(v,k) is practical only while length(v) is less than about 15.

Number bases

dec2bin(D,minDigits) returns a character array, one row per element, padded to equal length. Fractions are truncated down, negatives are two’s complement (R2020a), and minDigits never truncates. dec2hex and bin2dec/hex2dec work the same way on characters.

dec2bin(10)        % '1010'
dec2bin(5,8)       % '00000101'
bin2dec('1010')    % 10
dec2hex(255)       % 'FF'
hex2dec('FF')      % 255

The values in the comments were produced by Octave 11.3.0 (the negative-number behaviour is not run here).

Cumulative functions

cumsum, cumprod, cummax and cummin return running totals, products, maxima and minima along the first non-singleton dimension, or along dim.

cumsum([1 2 3 4])           % 1 3 6 10
cumprod([1 2 3 4])          % 1 2 6 24
cummax([1 3 2 5 4])         % 1 3 3 5 5
cummin([5 3 4 1 2])         % 5 3 3 1 1
cumsum([1 2; 3 4],2)        % [1 3; 3 7]

Octave 11.3.0 produced exactly these values. cumsum(A,"reverse") runs from the end of the dimension to the start, and a nanflag decides what a NaN does: "includemissing" (the default) propagates it to every later element, while "omitmissing" skips it (R2023a names; "includenan" and "omitnan" behave the same). cumprod, cummax and cummin take the same options. Missing values in general are on Descriptive Statistics and Preprocessing.

Polynomials

A polynomial is a row vector of coefficients, highest power first, with 0 for a missing power: is [1 -6 11 -6].

p = [1 -6 11 -6];
polyval(p,[0 1 2 3 4])       % -6 0 0 0 6
r = roots(p)                 % 1, 2, 3 (a column vector; order may differ)
poly([1 2 3])                % back from roots: 1 -6 11 -6
c = conv([1 2],[1 -1 3])     % product: 1 1 1 6
[q,rm] = deconv(c,[1 2])     % q = 1 -1 3, rm = 0
polyder(p)                   % 3 -12 11
polyint(p)                   % 0.25 -2 5.5 -6 0  (constant 0)
polyint(p,5)                 % same, constant of integration 5

Octave 11.3.0 produced these values. roots builds the companion matrix and returns its eigenvalues, so the roots of high-degree or ill-conditioned polynomials can be inaccurate; leading zeros are dropped and trailing zeros give zero roots. Octave returns roots([1 0 1]) as 0+1i and 0-1i. conv computes the convolution, which equals polynomial multiplication; see Fourier Analysis.

Fitting with polyfit

coeffs = polyfit(x,y,n)
[coeffs,fitStats] = polyfit(x,y,n)
[coeffs,fitStats,centerScale] = polyfit(x,y,n,Normalize=tf)

fitStats (the books and old code call it S) is a structure with R (triangular factor of the QR decomposition of the Vandermonde matrix), df (degrees of freedom), normr (norm of residuals) and, since R2024a, rsquared. centerScale (mu in older code) holds mean(x) and std(x) when normalisation is on, and [0 1] otherwise. Pass it to polyval(coeffs,x,[],centerScale) to evaluate a centred and scaled fit. Normalising matters for high degrees or large x, where the Vandermonde matrix is badly conditioned. Normalize=tf is the R2026b documented syntax; the Octave run below uses the three-output call.

t = 0:5; y = [1.1 2.9 5.2 6.8 9.1 11.0];     % temperature rise of a sensor
pf = polyfit(t,y,1)
[pf,S,mu] = polyfit(t,y,1);
yfit = polyval(pf,t,[],mu);
S.normr, S.df

Octave 11.3.0 output (slope 1.9914 and intercept 1.0381 for the plain fit; mu is 2.5000 and 1.8708; normr is 0.3272 and df is 4). Fitted values are 1.0381 3.0295 5.0210 7.0124 9.0038 10.9952. With three outputs Octave returned the normalised coefficients 3.7256 and 6.0167, which is why they are evaluated with mu as shown.

Interpolation

Interpolation estimates values between known samples. The documentation separates gridded data (interp1, interp2, interp3, interpn, griddedInterpolant) from scattered data (griddata, scatteredInterpolant), and offers piecewise-polynomial constructors (spline, pchip, makima, ppval, mkpp).

interp1

vq = interp1(x,v,xq,method,extrapolation) interpolates a vector. The methods are "linear" (the default), "nearest", "next", "previous", "pchip", "cubic" ("v5cubic" is the same), "makima" and "spline". Query points outside the range of x give NaN by default, except with "pchip", "spline" and "makima", which extrapolate with the same method.

x = 0:4; v = [0 1 4 9 16];
xq = [0.5 1.5 2.5 3.5];
for m = {'linear','previous','next','pchip','spline'}
  fprintf('%-9s', m{1}); fprintf('%8.4f ', interp1(x,v,xq,m{1})); fprintf('\n');
end

Octave 11.3.0 output:

linear     0.5000   2.5000   6.5000  12.5000
previous   0.0000   1.0000   4.0000   9.0000
next       1.0000   4.0000   9.0000  16.0000
pchip      0.3125   2.2188   6.2396  12.2292
spline     0.2500   2.2500   6.2500  12.2500

The samples lie on , so spline (not-a-knot) reproduces the parabola exactly. Outside the sample range, interp1(x,v,5) returns NaN in MATLAB (Octave prints NA); 'linear','extrap' extrapolates (23 here) and a number as last argument fills (-1). A missing value in v yields NaN for the intervals around it.

Resampling a heatflow sensor

t = [0 1.5 2 4.5 6]; T = [20.0 21.2 21.9 24.0 24.6];   % irregular times
tu = 0:1:6;
Tu = interp1(t,T,tu,'linear')
Tp = interp1(t,T,tu,'pchip')
Ts = interp1(t,T,7,'linear','extrap')

Octave 11.3.0 output:

20.000 20.800 21.900 22.740 23.580 24.200 24.600
20.000 20.691 21.900 22.913 23.702 24.246 24.600
25.000

pchip is shape preserving (no overshoot), which is usually what you want for a physical signal. Do not extrapolate far beyond the sampled time.

spline, pchip and makima

spline(x,y,xq) uses not-a-knot end conditions; if y has two more entries than x, the first and last are taken as end slopes (clamped spline). pchip(x,y,xq) is a shape-preserving piecewise cubic Hermite. makima(x,y,xq) is the modified Akima: fewer undulations than spline, and, unlike pchip, still able to follow oscillatory data. Each also has a pp = f(x,y) form returning a piecewise polynomial for ppval. makima is not run here (not in Octave 11.3.0).

pp = spline(0:4,[0 1 4 9 16]);   class(pp)    % 'struct'
ppval(pp,2.5)                                 % 6.25
makima(0:4,[0 1 4 9 16],[0.5 1.5])            % not run here

interp2

Vq = interp2(X,Y,V,Xq,Yq,method,extrapval) interpolates on a meshgrid-format grid. Methods: linear (default), nearest, cubic (uniform spacing, at least three points), makima and spline. Queries outside the grid give NaN except for spline and makima, which extrapolate. The grid points must be unique and strictly monotonic, and integer V needs nearest.

[X,Y] = meshgrid(0:2,0:2); V = X.^2 + Y;
interp2(X,Y,V,0.5,1.5)     % 2

Octave 11.3.0 output: 2.

griddedInterpolant and scatteredInterpolant

F = griddedInterpolant(x,v,Method,ExtrapolationMethod) builds a reusable interpolant of 1-D to N-D gridded data, evaluated as F(xq). Methods are the interp1 ones; ExtrapolationMethod also takes "none". Evaluating one object at many query sets is faster than repeated interp1/interpn calls. scatteredInterpolant(x,y,v,Method,ExtrapolationMethod) interpolates unstructured 2-D or 3-D points through a Delaunay triangulation, with methods linear (default), nearest and natural, extrapolation linear, nearest, boundary (R2024a) or none. Differing coordinate scales hurt it (use normalize), and duplicate points are averaged with a warning.

F = griddedInterpolant(0:4,[0 1 4 9 16],'pchip');
F(2.5)
S = scatteredInterpolant(px,py,pT,'natural');   % px, py, pT: sensor positions and temperatures
S(0.5,0.5)

Neither function exists in Octave 11.3.0: not run here.

Where the books are outdated

Older books show interp1 with the 'cubic' method as a separate algorithm, and interp2/griddata as the only scattered tools. makima, griddedInterpolant and scatteredInterpolant supersede ad-hoc loops; polyfit now has rsquared (R2024a) and Normalize=tf. Rounding to digits (round(x,n)) replaces the round(x*10^n)/10^n idiom. Integration (quad) is covered on Integration, ODEs and Optimization.