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 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