julstrat / damath Goto Github PK

The DAMath archive contains units for double precision accurate mathematical
methods without using multi precision arithmetic or assembler. The main
purpose is to make the AMath functions available for 64-bit systems without
Extended Precision or 387-FPU, but they can be used with 32-bit systems
and on the Raspberry Pi 3.

The source can be compiled with the usual Pascal versions that allow const
parameters (tested with BP 7.0, VP 2.1, FPC 1.0/2.0/2.2/2.4/2.6/3x, and
Delphi versions 2-7/9-10/12/17-18/25/26).

The units assume IEEE-754 53-bit double precision (binary64) and rounding to
nearest; since Aug. 2017 there is the separate unit DFPU with 64-bit/ARM
compatible rounding / precision / exception control functions based on the
compiler's math unit.


---------------------------------------------------------------------------

Last changes (Nov. 2018)

  * Unit DAMath V1.17 with the new functions: population covariance pcov
    and generic aliases Infinity, NegInfinity, NaN;

  * Special Functions V1.38.02 with the new functions: derivatives of the
    zero order Kelvin functions (kelvin_berp / beip / kerp / keip);

  * Unit DAMCmplx V0.70 with the new functions: zeta, sinpi, Ei, E1, li;
    changed functions: Set agm(x,-x) = 0, nroot special cases n=1,2; fixed:
    expm1 and ln1p for small pure imaginary arguments;

  * Updated internal archive complexdemo.zip;

  * AMath and DAMath document Special Functions Reference Manual with
    Implementation Notes Version 2.27 (including Complex Functions)



---------------------------------------------------------------------------
DAMath V1.17: This unit implements the double precision accurate elementary
mathematical functions: exponential, logarithmic, trigonometric, hyperbolic,
inverse circular / hyperbolic functions; and there are the polynomial,
vector, statistic operations and floating point functions. The file
damath_info.txt contains the complete list of DAMath functions with short
descriptions.

DAMath uses the system functions abs, arctan, frac, int, ln, and trunc for
64-bit (frac, int, ln, and trunc are bug-fixed for some 16/32-bit compilers).
Because FPC-64 (version <= 2.6.4) looses up to 13 of the 53 bits for exp,
DAMath implements its own exp function. System sin(x) and cos(x) are used for
|x| <= Pi/4, Payne/Hanek range reduction is performed if |x| > 2^30.

On Win7/64 the 64-bit DAMath's one argument elementary transcendental
functions (except cosd, sind, ..) and power have peak relative errors
< 2*eps_d (about 4.4E-16), the RMS values are < 0.6*eps_d, a complete list
and the Delphi, FPC, ARM figures can be found in the log files *.cmp.

The dbl2 (double-double) routines operate on pairs of 53-bit double
floating point numbers, which represent the unevaluated sum of the high and
low parts: a = (a.h, a.l) = a.h + a.l, normally with |a.l| < eps_d*|a.h|.

---------------------------------------------------------------------------
The DAMTools V0.47 unit provides accurate and reliable tools for finding
zeros and local minima of functions, numerical integration of one-dimensional
functions,convergence acceleration of sums and sequences, and solving
quadratic /cubic / polynomial equations:

The functions localmin, mbrent, and fmin (differing in parameter count/ease of
use) use Brent's algorithm with guaranteed convergence for finding a local
minimum of a function f in an interval (a,b). The algorithm combines golden
section search and successive parabolic interpolation using only function (not
derivative) evaluations. The functions zbrent and zeroin use the Brent/Dekker
algorithm with guaranteed convergence for finding a zero of a continuous
function f in the interval [a,b], when f(a) and f(b) have different signs;
zbrenty handles the zeros of f(x)-y. The algorithm is based on a combination
of successive interpolations and bisection. zridders implements Ridders' root
finding method based on regula falsi and exponential factors.

The qag* procedures are Pascal translations of Quadpack algorithms by R.
Piessens, E. de Doncker-Kapenga, C.W. Ueberhuber, D. Kahaner. These routines
perform global adaptive quadrature of functions over finite or infinite
intervals based on Gauss-Kronrod rules for the subintervals and acceleration
by Wynn's epsilon algorithm, they can handle rather difficult integrals
including integrand functions with local difficulties such as a disconti-
nuities and integrable singularities. quagk is a simple general purpose
shell for the qag* routines. The Quadpack algorithm qawc computes the Cauchy
principal value of f(x)/(x-c) using a globally adaptive strategy and
modified Clenshaw-Curtis integration on the subintervals containing the
point x=c.

The procedures intdeo and intdei use the Double Exponential (DE) trans-
formation (developed by M. Mori, T. Ooura, and others) for automatic
quadrature of f(x) over the infinite interval (a,+INF) for functions with
and without oscillatory factors resp. intde integrates over finite
intervals (a,b) and intde_p, intdei_p are the corresponding procudures for
functions f(x,p) with parameters.

The adaptive quanc8 algorithm by G.E. Forsythe, M.A. Malcolm, C.B. Moler
estimates the integral of a smooth function over a finite interval using a
Newton-Cotes rule.

The procedures levinu1 and wynneps1 are stand-alone versions of convergence
acceleration methods, they perform one step of the Levin u-transformation for
sums or one step of Wynn's epsilon algorithm for the sequences or sums
(original customized versions are used in LerchPhi and the Quadpack routines,
respectively). The calling driver routines have analyze the convergence of
the processes.

The squad functions accurately solve quadratic equations with double
coefficients; they implement ideas of G.E. Forsythe, W. Kahan, P.H. Sterbenz
(high precision calculation of discriminant, scaling by powers of two etc).
The cubsolve procedure computes the solutions of real cubic equations with
double coefficients; it is based on lecture notes by W. Kahan. The PolyRoots
procedures computes the n (complex) roots x[k] + i*y[k] of the polynomial
p(z) = p[0] + p[1]*z + ... p[n]*z^n   using a companion matrix method,
balancing, and the QR algorithm for the eigenvalues of an upper Hessenberg
matrix.


---------------------------------------------------------------------------
The unit SpecFunD V1.38.02 implements many special functions for double and
extended precision. The DAMath based double precision special functions are
derived from the AMath implementations (i.e. the descriptions in the AMath
introduction and the specialfunctions.pdf manual are essentially valid), the
interfaced functions have the same names, etc.

In this release the following function groups are available:

 * Bessel functions and related
 * Elliptic integrals, elliptic / theta functions
 * Error function and related
 * Exponential integrals and related
 * Gamma function and related
 * Orthogonal polynomials, Legendre functions, and related
 * Hypergeometric functions and related
 * Statistical distributions
 * Zeta functions, polylogarithms, and related
 * Other special functions

The relative errors of the DAMath special functions are usually larger
(especially on 64-bit systems) than those of the corresponding AMath double
functions (which are often correctly rounded due to the internal extended
precision calculations).


---------------------------------------------------------------------------
The unit DAMCmplx V0.70 provides DAMath based complex arithmetic, basic and
transcendental functions. The complex data type is a record with real and
imaginary parts of type double. Most routines are procedures with const input
record(s) and a var output record.

Please note that the exponential, trigonometric, or hyperbolic functions
may overflow or return INFs or NaNs for inputs with real or imaginary parts
of order ln_MaxDbl or greater, this will be handled more systematically in
future versions.

Internal archive complexdemo with the complex expression parser / interpreter
Damccalc.pas and the demo programs PlotPhaseD (a GUI application for the
visualization of complex functions) and T_CCalcD (a console program for the
evaluation of complex expressions);


---------------------------------------------------------------------------
The unit DAMQuat V0.20 units provides DAMath based quaternionic arithmetic and
elementary transcendental functions (here is the complete list with short
descriptions) The quaternion data type is a record with real (or scalar) and
imaginary (or vector) parts using the base type double.

  type
    Quaternion = record
                   r:     double; {real or scalar part     }
                   x,y,z: double; {imaginary or vector part}
                 end;

Most routines are procedures with const input record(s) and a var output
record. A quaternion is often written as a = r*1 + x*i + y*j + z*k, where 1,
i, j, k are the fundamental quaternion units (as 4-dimensional vectors they
are the standard unit vectors).

Addition and subtraction are defined component-wise, i.e. the standard
vector addition in R4. Multiplication of quaternions is associative and
distributive, but it is not commutative, it is determined by the relations
i^2 = j^2 = k^2 = ijk = -1. The norm of a quaternion is r^2 + x^2 + y^2 +z^,
the absolute value |a| is the square root of the norm, and conj(a) = r*1 -
x*i - y*j - z*k is the conjugate. The multiplicative inverse is (1/a) =
conj(a)/norm(a), division is defined as a/b = a*(1/b).

Most (inverse) trigonometric / hyperbolic quaternion functions b=f(a) are
basically computed with the corresponding complex function w=F(z) with z =
Re(a) + i*abs(Im(a)) and the mapping to quaternions Re(b) = Re(w), Im(b) =
Im(w)*Im(a)/abs(Im(a)). Note that DAMQuat has no quaternionic power function
a^b, if b=x is a real number, a^x is defined as exp(x*ln(a)).


---------------------------------------------------------------------------
http://wolfgang-ehrhardt.de/