delorie.com/archives/browse.cgi   search  
Mail Archives: djgpp/1994/09/28/20:59:09

Date: Wed, 28 Sep 94 18:33:06 +0100
From: buers AT dg1 DOT chemie DOT uni-konstanz DOT de (Dieter Buerssner)
To: djgpp AT sun DOT soe DOT clarkson DOT edu
Cc: djgpp-announce AT sun DOT soe DOT clarkson DOT edu
Subject: alternative math library for djgpp

Hello,

I have uploaded an alternative math library to

    omnigate.clarkson.edu:pub/msdos/djgpp/pub/dbmlib01.zip

It is packed like the usual djgpp contributions, so you
should unpack it in the djgpp directory with unzip -d dbmlib01.zip

Highlights: missing functions included, long double precision
versions of most functions, extensive test packages, some
functions do yield superior accuracy compared to libm.a.
I append the file readme.1st, which is the only documentation
available.

Dieter

--------- readme.1st

    libmn, an alternative math library for djgpp
    --------------------------------------------

This library can replace (or supplement) the original math library of
djgpp.  It includes all the functions of libm, plus long double
precison versions of almost all functions (and more, see below). Some
of the functions yield superior accuracy compared to the original
libm; they may need more time, though. Most of the functions, that are
not in the Borland C math library, will also work with Borland C.

I have only written a small part of the code myself, a lot of the
functions are collected from different sources.  Usually I have
'ANISfied' these sources and modified them to run faster (some inline
assembly, some modified algorithms) with djgpp. Also I may have
removed very few bugs in the originals.  Of course, I hope that I have
not introduced any bugs.  Some of the comments in the source might not
be appropriate anymore.

As far as I know, all sources are freely distributable in source and
in binary (no GNU public licence).  In most of the code written by me,
there is no copyright message, sometimes there will be a comment like
/* -buers */.  Everything I have written can be distributed freely.  I
have not always noted in the code, when I have modified anything.  You
can assume, that I modified all files slightly, excluding the
assembler files distributed with djgpp. So please don't complain to
the original authors, when anything fails.

The original authors are
    DJ Delorie
	Most of the functions written in assembly language
	(and some of the C code as well) are written
	by DJ Delorie. I have modified them slightly for long double
	precision. For copyright info see copying.dj.
    UCB
	Most of the hyperbolic functions, cbrt and probably more
	are derived from the BSD net2 distribution of the
	University of Berkley.
	Their copyright info is in the sources.
    Stephen L. Moshier
	powl, cbrtl and the extended precision argument reduction code
	of the long double precision trigonometric functions are from the
	brilliant cephes math library by Stephen L. Moshier. (See copying.sm)
	Also expl.S is almost a direct translation of expl in the cephes
	library into assembly.
    Richard Geers
	The ...87 functions were inspired by similar functions of
	the math library of Richard Geers distributed with some
	versions of BSD Unix.
    W.J. Cody
	The erf, gamma and some of the bessel functions were written
	by W.J. Cody, see the comments in the source code. These files
	were distributed with Richard Geers math library.

I wish to thank all the authors.
If I have forgotten anybody, I appologize. Please let me know.

Installation

To install the library, you can just move lib/libmn.a to a directory
where gcc will look for libraries and the include files to the usual
include directory. (There is a slightly modified version of math.h in
the include directory. The original of djgpp112 has a wrong prototype
for finite).  Then you have to link with -lmn instead of -lm. To
recompile the library use mk.bat

There are two header files in the include directory, math.h and mathext.h.
math.h contains the functions of the original math library of djgpp,
and mathext.h some extra functions, that, as far as I know, are
not in the ANSI C standard (this should avoid the pollution of the
name space). There are also some useful structures and constants
defined in mathext.h.

Differences to libm.a of djgpp

All of the functions in libm.a are included in this library.  Most of
them are unmodified. The hyperbolic functions of libm.a loose
precision for extreme arguments. These have been replaced by adapted
code from UCB. cosh is still the original, because it doesn't suffer
from this loss of precision.  Also frexp has been replaced (the
original will give wrong results for 0).

Almost all functions have a long double precision counterpart in this
library. Usually, the original name is just followed by 'l' (double
sin(double) -> long double sinl(long double)).  To use these
functions, you must #include <mathext.h> The functions sinl, cosl,
tanl and powl also have versions with an underscore (_sinl, ...). The
later are written in assembly an have a little bit less accuracy for
extreme arguments, but run considerably faster.

The following functions are not in libm:

1. functions for control of the 80x87 NPX

    void _clex87(void);
	clear 80x87 exceptions.
    unsigned short _set_cw87(unsigned short cw);
	sets the control word of the 80x87 to cw. Returns the
	previous set control word.
    unsigned short _get_sw87(void);
	Returns the status word of the 80x87.

2. IEEE recommended functions

    double copysign(double x, double y)
	returns x with the sign of y.
    double drem(double x, double y)
	returns IEEE 754 remainder r of x/y,
	r = x - n*y with integer n nearest to the exact value of x/y
	if |n-x/y| = 0.5, n will be even.
    double logb(double x);
	returns the unbiased exponent of its argument.
    double scalb(double x, int n);
	x*2^n.
    double rint(double x);
	returns x rounded to an int. With the normal setting
	of the 80x87 status word this means that some integer + 0.5
	will round to the next even integer. All other numbers (not
	ending exactly with .5) will round to the next integer.
    int isinf(double x)
	returns 1 if x is +/- Inf; 0 otherwise.
    int isnan(double x);
	returns 1 if x is a NaN; 0 otherwise
    int finite(double x);
	returns 1 if x is not a NaN and not +/- Inf; 0 otherwise
    These all have long double precision versions.

3. transcendental functions usually found in UNIX math libraries

    double erf(double x);
	error function
    double erfc(double x);
	1.0 - erf(x)
    double erfcx(double x);
	x^2 * erfc(x), this is usually not found in standard math libs.
    double gamma(double x)
    double lgamma(double x)
	gamma function and log gamma.
    double j0(double x)
    double j1(double x)
    double jn(int n, double x)
    double y0(double x)
    double y1(double x)
    double yn(int n, double x)
	Bessel functions.
    double cbrt(double x)
	cube root
    double expm1(double x)
	exp(x) - 1.0, is accurate for x close to zero
    double log1p(double x)
	log(1.0+x), is accurate for x close to zero
    and some more. The last three have long double precision versions.

There are also some more support functions (starting with an
underscore).  Those are only needed internally, but might be useful
for low level code. Note that the djgpp libm[n] function pow10 has a
different prototype than the Borland function of the same
name. (double pow10(double) in libm[n] vs. double pow10(int)). 
pow10l(x) may yield some less accuracy than powl(10.0L,x), 
but will run faster.

Known Bugs
    frexp and frexpl won't work for denormalized arguments.
    (This should only happen with extremely small arguments.)
    hypot may fail for very large arguments.
    the functions don't have support for matherr, they rather
    return NaN or Inf (or 0 in the case of underflow). This
    is easier, should be faster and I prefer it :-)

I have used this library for almost 2 years now, and did not have any
major problems. The accuracy should be at least as good as in libm.  I
have tested the library with the famous paranoia program, (Professor
W. M. Kahan, C Version by David M. Gay and Thos Sumner) with a package
called elefunt (C Version of the Code & Waite ELEFUNT Elementary
Function Test Package by Ken Stoner and Nelson H.F. Beebe, see README
in tests/elefunt) and with a consistency test for math functions by
Stephen L. Moshier (from the cephes math library).  I wish to thank
all the authors of these packages.

I have not done very much testing with a math emulator, but as I remember
the test programs did run correctly and gave almost the same accuracy
with the math emulator of Bill Metzenthen. (This was some versions ago).
The accuracy was worse, but still appropriate if you don't have
the highest demands, with the math emulator of DJ Delorie.


The test packages are included in the tests subdirectory.  To compile
them, you need dmake (available from Simtel archives), or you have to
modify the makefiles. I have adapted the test packages to work in long
double precision as well. To compile the long double versions call
"dmake LDOUBLE=1".  Paranoia will give some flaws in (normal) double
precision, but not in long double precision. These flaw are probably
harmless and due to extra precision (on the 87 stack, there are always
long double precision numbers).  The extra precision may depend on
compilation flags or on the compiler version.  Another complaint is
for 0^0 not beeing 1 (pow(0,0) = NaN). I think, this is a question of
taste or philosophy.  In the directory test/db are 2 programs to test
the functions not coverd by the other test packages. These two programs
are not self explaining, you have to read the source to understand
them.  They are of very limited use, anyway, so you might not bother.

(If you have too much spare time, you might have fun in trying to understand
the source of paranoia :-)

The package also compiles under Borland-C (and probably Turbo-C), but
not all functions may work with the Borland supplied math emulator. I
think the 387 opcodes are not supported, so you might need a real 387,
a 287 will not be enough. To recompile the library with Borland-C you
need dmake:

> dmake -f makefile.bcc MODEL=l

Note that the test programs may give peculiar output (like 0.000... !=
0.000....). This is because of bugs in the Borland printf for very
large or very small arguments.  You will have to modify
tests/cephes/mtst.c to make it work with Borland-C.

Feel free to send comments, complaints, bugs and patches to

    buers AT dg1 DOT chemie DOT uni-konstanz DOT de

I will only have extremely limited time to support this library,
but I will try to fix reported bugs. It might take some time for me
to reply to email.

I have now idea, if the following is really needed, but ...

THIS SOFTWARE IS PROVIDED ``AS IS'' AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH DAMAGE.


I hope, this library will be useful to some of you.

    Dieter Buerssner

- Raw text -


  webmaster     delorie software   privacy  
  Copyright © 2019   by DJ Delorie     Updated Jul 2019