NUMLIBC
Welcome to the NUMLIBC library of numerical routines written in the programming language C. These routines are free, but for any damages caused by their use we do not have any responsibility.

The library has been organized into the following categories. Please see below for how to obtain a single file or the whole library.

For each routine there are one file for the code with the extension C and another for the documentation with the extension DOC. You can get one single file by clicking on it or download a packed version of the whole library by clicking on numlibc.zip. If you fetch single files be sure also to get the header files you need.

Here follows a few examples of how the routines can be used:

Good luck with NUMLIBC.

Header files for the NUMLIBC library.
numglob.h
Contains global definitions. Include this file only once if you have more than one file making up your program system.
numlibc.h
Contains definitions of constants and types and also prototypes for functions. Include this file where you use elemens of NUMLIBC.

A1 - Real arithmetic.
mach_err
Computes the machine error in double precision.
(mach_err.c, mach_err.doc)

B1 - Trig. and inverse trig. functions.
dtan
Computes tan(x) for a given argument.
(dtan.c, dtan.doc)

B3 - Exponential and logarithmic functions.
fac
Computes n! for an integer n > 0.
(fac.c, fac.doc)
power
Raises t to power n.
(power.c, power.doc)

C1 - Operations on polynomials and power series.
hor_ner
Computes the value of a polynomial for a given argument using Horner's algorithm.
(hor_ner.c, hor_ner.doc)
syn_div
Syntetic division on polynomial for a given argument.
(syn_div.c, syn_div.doc)

C2 - Zeros of polynomials.
pol_rot
Calculation of a single root in a polynomial equation with real coefficients.
(pol_rot.c, pol_rot.doc)

C3 - Calculation of special functions.
sign
Computes the sign of a double precision number.
(sign.c, sign.doc)
ber_num
Calculation of the Bernoulli numbers and the sums of reciprocal powers.
(ber_num.c, ber_num.doc)
ber_func
Calculation of the Bernoulli functions and the related scaled functions.
(ber_func.c, ber_func.doc)
d2q_cof
Calculation of some coefficients related to Euler-Maclaurin type expansions for oscillatory integrals.
(d2q_cof.c, d2q_cof.doc)
d2q_func
Calculation of some functions related to Euler-Maclaurin type expansions for oscillatory integrals.
(d2q_func.c, d2q_func.doc)
I0
Bessel function I_0(x), x >= -8.
(i0.c, i0.doc)
I1
Bessel function I_1(x), x >= -8.
(i1.c, i1.doc)
J0
Bessel function J_0(x), x >= -8
(j0.c, j0.doc)
J1
Bessel function J_1(x), x >= -8.
(j1.c, j1.doc)
K0
Bessel function K_0(x), x > 0.
(k0.c, k0.doc)
K1
Bessel function K_1(x), x > 0.
(k1.c, k1.doc)
Y0
Bessel function Y_0(x), x > 0.
(y0.c, y0.doc)
Y1
Bessel function Y_1(x), x > 0.
(y1.c, y1.doc)
factorial
Computes the factorial function x! for real values of x.
(factoria.c, factoria.doc)
ErrorFunction
Computes the error function erf(x) for values of x >= -4.
(errorfun.c, errorfun.doc)
ExponentialIntegral
Computes the exponential integral E1(x) for values of x >= -4.
(expint.c, expint.doc)

C5 - Zeros of one or more trancendental equations.
int_halv
Calculation of a single root in the equation f(x)=0, using the interval halving method.
(int_halv.c, int_halv.doc)
reg_fals
Calculation of a single root in the equation f(x)=0, using Regula-falsi type methods.
(reg_fals.c, reg_fals.doc)
inv_int
Calculation of a single root in the equation f(x)=0, using inverse interpolation.
(inv_int.c, inv_int.doc)
root_func
A set of test functions for computing a root in the equation f(x)=0.
(root_fun.c, root_fun.doc)
ite_rate
Calculation of a single root in the equation x=f(x) using either ordinary fixpoint iteration or Steffensen's method.
(ite_rate.c, ite_rate.doc)
se_cant
Calculation of a single root in the equation f(x)=0 using the secant method.
(se_cant.c, se_cant.doc)

D1 - Quadrature.
simp_int
Approximate calculation of an integral using trapezoidal-, midpoint or Simpson's formula.
(simp_int.c, simp_int.doc)
trap_int
Approximate calculation of an integral using the trapezoidal approximation and repeated halving of the integration steplength.
(trap_int.c, trap_int.doc)
rom_int
Approximate calculation of an integral using classical Romberg integration with repeated halving of the integration step length.
(rom_int.c, rom_int.doc)
rom_kah
Approximate calculation of an integral using generalized Romberg-Kahan integration.
(rom_kah.c, rom_kah.doc)
gauss_cof
Abscissae and weights for Gauss-Legendre integration.
(gauss_co.c, gauss_co.doc)
gauss_legendre
Aproximate calculation of an integral using Gauss-Legendre quadrature for n=2(1)16.
(gauss_le.c, gauss_le.doc)
gauss_quad
Approximate calculation of an integral using Gauss-Legendre quadrature where abscissae and weights have to be supplied using gauss_cof.
(gauss_qu.c, gauss_qu.doc)
herm_cof
Abscissae and weights for Gauss-Hermite integration.
(herm_cof.c, herm_cof.doc)
herm_quad
Approximate calculation of an integral using Gauss-Hermite quadrature where abscissae and weights have to be supplied using herm_cof.
(herm_qua.c, herm_qua.doc)
laguer_cof
Abscissae and weights for Gauss-Laguere integration.
(laguer_c.c, laguer_c.doc)
laguer_quad
Approximate calculation of an integral using Gauss-Laguer quadrature where abscissae and weights have to be supplied using laguer_cof.
(laguer_q.c, laguer_q.doc)
log_cof
Abscissae and weights for Gauss-quadrature with logarithmic weightfunction, ln(x).
(log_cof.c, log_cof.doc)
log_quad
Approximate calculation of an integral with logarithmic weight function, ln(x), using Gauss-quadrature where abscissae and weights have to be supplied using log_cof.
(log_quad.c, log_quad.doc)
test_func
A selected set of test functions for numerical quadrature.
(test_fun.c, test_fun.doc)
val_int
Integration boundaries and related "exact" values for the test functions given in test_func above.
(val_int.c, val_int.doc)
dat_int
Approximate calculation of an integral where the integrand is given in discrete points x0 < x1 ... < xn. An averaging process is used.
(dat_int.c, dat_int.doc)
pv_int0
Computes an approximation to the Cauchy principal value integral with the integrand f(x)/(x-c) over the integration range [a, b] with a < c < b (c may also be slightly outside the integration range. The method applied is to approximate f(x) with an expansion in terms of Chebyshev polynomials.
(pv_int0.c, pv_int0.doc)
pv_int1
Computes an approximation to the Cauchy principal value integral with the integrand f(x)/(x-c) over the integration range [a, b] with a < c < b (c may also be slightly outside the integration range. The method applied is a Romberg type extrapolation procedure applied to Euler's second expansion using the midpoint approximation combined with a repeated halving of the integration interval.
(pv_int1.c, pv_int1.doc)
int_func
Computes the integral function for a given on a given interval in equidistant points. Simpson's rule is used.
(int_func.c, int_func.doc)

D2 - Numerical solution of ordinary differential equations.
eu_ler
Numerical integration of a system of first order differential equations using Euler's first order method.
(eu_ler.c, eu_ler.doc)
imp_euler
Numerical integration of a system of first order differential equations using improved Euler's (Heun's) second order method.
(imp_eule.c, imp_eule.doc)
mod_euler
Numerical integration of a system of first order differential equations using modified Euler's (Runge's) second order method.
(mod_eule.c, mod_eule.doc)
run_kut4
Numerical integration of a system of first order differential equations using Runge-Kutta's classical 4th order method.
(run_kut4.c, run_kut4.doc)
pre_cor2
Numerical integration of a system of first order differential equations using a second order predictor-corrector method.
(pre_cor2.c, pre_cor2.doc)

D3 - Numerical solution of partial differential equations.
Crank_Nicolson
Numeriacal integration of the heat conduction equation one step using Crank Nicolsons method.
(crank_ni.c, crank_ni.doc)
wave_eq
Numerical integration of the wave equation (vibrating string) one step using finite differences.
(wave_eq.c, wave_eq.doc)
pois_eq
Numerical integrationm of Poisson's equation using finite differences and block iteration.
(pois_eq.c, pois_eq.doc)

E1 - Interpolation.
lag_cof
Calculation of coefficients in Lagrange interpolation formula.
(lag_cof.c, lag_cof.doc)
lag_val
Lagrange interpolation using the Barycentric representation.
(lag_val.c, lag_val.doc)
cheb_cof
Calculation of gridpoints and weights for Lagrange-Chebyshev interpolation formulae and Lagrange equidistant abscissae interpolation formulae (using Barycentric representation).
(cheb_cof.c, cheb_cof.doc)
cheb_val
Lagrange-Chebyshev endpoint interpolation using the Barycentric representation.
(cheb_val.c, cheb_val.doc)
new_val
Newton's interpolation formula with given wanted relative accuracy in the interpolation.
(new_val.c, new_val.doc)
index
Computes indices for array elements ordered after increasing distances from a given number. A service function for new_val and others.
(index.c, index.doc)
new_cof
Computes the coefficients in Newton's interpolation formula from a given dl of a function f(x).
(new_cof.c, new_cof.doc)
new_sum
For given coefficients in Newton's interpolation formula and a given argument the interpolated value is computed.
(new_sum.c, new_sum.doc)
ait_ken
Aitken's interpolation/extrapolation formula with given wanted relative accuracy in the interpolation/extrapolation.
(ait_ken.c, ait_ken.doc)
nev_ille
Neville's interpolation/extrapolation formula with given wanted relative accuracy in the interpolation/extrapolation.
(nev_ille.c, nev_ille.doc)

E2 - Curve and surface fitting.
cub_spl
Calculation of the coefficients in a cubic spline approximation for a function f(x) given at the points x0 < x1 < ... < xn. For applications see spl_tab and spl_int.
(cub_spl.c, cub_spl.doc)
spl_tab
Tabulation of a function using the cubic spline approximation computed by cub_spl.
(spl_tab.c, spl_tab.doc)
spl_int
Calculation of the integral of a function approximating the integrand with the cubic spline approximation computed by cub_spl.
(spl_int.c, spl_int.doc)
spl_cub
Calculating of the coefficients in a cubic spline approximation for the integral function F(x) for a function f(x) given at the points x0 < x1 < ... < xn. For applications see spl_use.
(spl_cub.c, spl_cub.doc)
spl_use
Using the spline approximation computed by spl_cub this routine may be used for tabulating the derivative function f'(x), the integrand function f(x) and the integral function F(x).
(spl_use.c, spl_use.doc)

E4 - Approximation of functions.
trig_cof
Calculation of the coefficients in an interpolating trigonometric polynimial.
(trig_cof.c, trig_cof.doc)
trig_sum
Summation of trigonometric series.
(trig_sum.c, trig_sum.doc)
cheby_chev
Calculation of the coefficients in an approxmation in terms of Chebyshev polynomials for a given function f(x).
(cheby_ch.c, cheby_ch.doc)
cheb_sum
Summation of a Chebyshev expansion.
(cheb_sum.c, cheb_sum.doc)
pow_to_cheb
Conversion of a power expansion to an expansion in terms of ordinary Chebyshev polynomials of first kind.
(pow_to_c.c, pow_to_c.doc)
pow_to_shift_cheb
Conversion of a power expansion to an expansion in terms of shifted Chebyshev polynomials of first kind.
(pow_to_s.c, pow_to_s.doc)

F1 - Matrix and vector operations.
vector_norm
Calculation of the 1-, 2- and maximum norms for a given vector.
(vector_n.c, vector_n.doc)
matrix_norm
Calculation of the 1-, maximum and E-norms for a given matrix.
(matrix_n.c, matrix_n.doc)

F2 - Eigenvalues and eigenvectors.
pow_met
Computes the eigenvalue of largest abs. value and the related eigenvector for a general n*n matrix (with real elements). The calculation may be combined with the use of a shift parameter.
(pow_met.c, pow_met.doc)
Hy_man
Computes a single eigenvalue and the related eigenvector for an upper Hessenberg matrix using Hyman's method.
(hy_man.c, hy_man.doc)
Gi_vens
Computes a single eigenvalue and the related eigenvector for a symmetric tridiagonal matrix using Givens' method.
(gi_vens.c, gi_vens.doc)

F3 - Linear equations.
gauss_elim
Gaussian elimination on a system of linear equations.
(gauss_el.c, gauss_el.doc)
gauss_solv
Computes the solution of a system of equation reduced to upper triangular form by gauss_elim.
(gauss_so.c, gauss_so.doc)
tri_diag
Computes the solution of a linear system of tridiagonal equations with a single righthand side. The coefficient matrix and the righthand side is not changed.
(tri_diag.c, tri_diag.doc)
tri_lu
Computes the solution of a linear system of tridiagonal equations with several righthand sides. The coefficient matrix is changed.
(tri_lu.c, tri_lu.doc)
sym_eqs
Computes the solution of a linear system of equations with a symmetric and positive definit coefficient matrix. Cholesky's method is used.
(sym_eqs.c, sym_eqs.doc)
suc_rel
Computes the solution of a linear system of equations using successive under/over relaxation iteration.
(suc_rel.c, suc_rel.doc)
gauss_reduk
Gauss elimination on a linear system of equations. Single precision calculation is used.
(gauss_re.c, gauss_re.doc)
gauss_back
Computes the solution of a system of equations reduced to upper triangular form by gauss_reduk. Single precision calculation is used.
(gauss_ba.c, gauss_ba.doc)
g_it_imp
Computes the solution of a system of linear equations in single precision combined with an iterative improvement. The functions gauss_reduk and gauss_back are used for computing the start vector for the iteration.
(g_it_imp.c, g_it_imp.doc)

Z2 - Input routines/parser routines.
Parse
Translates a functional expression into a stack of operations.
(parse.c, parse.doc)
Evaluate
Computes function values with one to three independent variables, x, y and z from the stack made in the function Parse(). Evaluate() uses the function Compute() for the actual calculation.
(evaluate.c, evaluate.doc)
Compute
Used by the Evaluate-functionsdescribed above to do the actual function calculations.
(compute.c, compute.doc)

Lars Erik Lund <larserik[+]math.ntnu.no>
Last modified: Aug 30 1996