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 EulerMaclaurin type expansions for oscillatory integrals.
(d2q_cof.c, d2q_cof.doc)
 d2q_func
 Calculation of some functions related to EulerMaclaurin 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 Regulafalsi 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 RombergKahan integration.
(rom_kah.c, rom_kah.doc)
 gauss_cof
 Abscissae and weights for GaussLegendre integration.
(gauss_co.c, gauss_co.doc)
 gauss_legendre
 Aproximate calculation of an integral using GaussLegendre quadrature for n=2(1)16.
(gauss_le.c, gauss_le.doc)
 gauss_quad
 Approximate calculation of an integral using GaussLegendre quadrature where abscissae and weights have to be supplied using gauss_cof.
(gauss_qu.c, gauss_qu.doc)
 herm_cof
 Abscissae and weights for GaussHermite integration.
(herm_cof.c, herm_cof.doc)
 herm_quad
 Approximate calculation of an integral using GaussHermite quadrature where abscissae and weights have to be supplied using herm_cof.
(herm_qua.c, herm_qua.doc)
 laguer_cof
 Abscissae and weights for GaussLaguere integration.
(laguer_c.c, laguer_c.doc)
 laguer_quad
 Approximate calculation of an integral using GaussLaguer quadrature where abscissae and weights have to be supplied using laguer_cof.
(laguer_q.c, laguer_q.doc)
 log_cof
 Abscissae and weights for Gaussquadrature 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 Gaussquadrature 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)/(xc) 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)/(xc) 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 RungeKutta'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 predictorcorrector 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 LagrangeChebyshev interpolation
formulae and Lagrange equidistant abscissae interpolation formulae
(using Barycentric representation).
(cheb_cof.c, cheb_cof.doc)
 cheb_val
 LagrangeChebyshev 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 Enorms 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 Evaluatefunctionsdescribed above to do the actual function calculations.
(compute.c, compute.doc)
Lars Erik Lund
<larserik[+]math.ntnu.no>
Last modified: Aug 30 1996