Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP Path: utzoo!utgpu!water!watnot!watmath!clyde!rutgers!seismo!nbires!hao!noao!mcdsun!fnf From: fnf@mcdsun.UUCP Newsgroups: net.sources Subject: Portable Math Library in C (Part 2 of 6) Message-ID: <290@mcdsun.UUCP> Date: Fri, 10-Apr-87 18:41:36 EST Article-I.D.: mcdsun.290 Posted: Fri Apr 10 18:41:36 1987 Date-Received: Sun, 12-Apr-87 06:12:02 EST Organization: Motorola Microcomputer Division, Tempe, Az. Lines: 2186 Keywords: math library This is a portable math library written entirely in C. Since it has been several years since I had any interest in doing any more work on it, and people may find it useful, I have decided to post it. There should be a lead-in posting (part 0 of 6?) containing a README file and commands to make the directories for the regular parts 1-6. Be sure to read the README file in 'part 0'. -Fred Fish =============== Cut here and feed to the shell ======================== #! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh funcs/environ/testfrexp.c <<'END_OF_funcs/environ/testfrexp.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FILE X * X * testfrexp.c test the runtime environment function frexp X * X * DESCRIPTION X * X * This simple minded program is provided to aid in testing X * the "frexp" function assumed to be provided in the runtime X * environment. If not provided, a suitable substitute can X * be coded in C, however the necessary code is very machine X * dependent and is generally almost trivial to code in assembly X * language for the specific host machine. X * X * The frexp() function takes two arguments, the first is a double X * and the second is a pointer to an int. Frexp() returns the X * mantissa of the first arg, and stores the exponent indirectly X * in the location pointed to by the second arg. X * X * See "frexp(3C)" in the Unix System V User's Manual for more X * information. X * X * This program is typically used as: X * X * testfrexp junkfile X * diff testfrexp.out junkfile X * X */ X X#include X Xextern double frexp (); X Xmain () X{ X double input; /* Input value */ X double dmant; /* Mantissa, 0.5 LE |dmant| LT 1.0 */ X int intexp; /* Exponent as power of 2 */ X X while (scanf ("%le", &input) == 1) { X dmant = frexp (input, &intexp); X printf ("%le %d\n", dmant, intexp); X } X} END_OF_funcs/environ/testfrexp.c if test 2017 -ne `wc -c funcs/environ/testldexp.c <<'END_OF_funcs/environ/testldexp.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FILE X * X * testldexp.c test the runtime environment function ldexp X * X * DESCRIPTION X * X * This simple minded program is provided to aid in testing X * the "ldexp" function assumed to be provided in the runtime X * environment. If not provided, a suitable substitute can X * be coded in C, however the necessary code is very machine X * dependent and is generally almost trivial to code in assembly X * language for the specific host machine. X * X * The ldexp() function takes two arguments, the first is a double X * which is the mantissa of the returned double value. The second X * is an int which is the exponent of the returned double. The X * result is (mantissa * (2 ** exp)). X * X * See "frexp(3C)" in the Unix System V User's Manual for more X * information. X * X * This program is typically used as: X * X * testldexp junkfile X * diff testldexp.out junkfile X * X */ X X#include X Xextern double ldexp (); X Xmain () X{ X double dmant; /* Mantissa, 0.5 LE |dmant| LT 1.0 */ X int intexp; /* Exponent as power of 2 */ X double result; /* dmant times 2 to the intexp */ X X while (scanf ("%le %d", &dmant, &intexp) == 2) { X result = ldexp (dmant, intexp); X printf ("%le\n", result); X } X} END_OF_funcs/environ/testldexp.c if test 2027 -ne `wc -c funcs/environ/testmodf.c <<'END_OF_funcs/environ/testmodf.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FILE X * X * testmodf.c test the runtime environment function modf X * X * This simple minded program is provided to aid in testing X * the "modf" function assumed to be provided in the runtime X * environment. If not provided, a suitable substitute can X * be coded in C, however the necessary code is very machine X * dependent and is generally almost trivial to code in assembly X * language for the specific host machine. X * X * The modf() function takes two arguments. The first is a double value X * and the second is a pointer to a double. Modf() returns the X * signed fraction part of the first argument, and stores the integral X * part (as a double) indirectly in the location pointed to by the X * second argument. Note that both the direct and indirect result will X * have the same sign, and: X * X * + = X * X * See "frexp(3C)" in the Unix System V User's Manual for more X * information. X * X */ X Xextern double modf (); X Xmain () X{ X double input; /* Input value */ X double frac; X double ipart; X X while (scanf ("%le", &input) == 1) { X frac = modf (input, &ipart); X printf ("%le %le\n", frac, ipart); X } X} END_OF_funcs/environ/testmodf.c if test 1977 -ne `wc -c funcs/src/acos.c <<'END_OF_funcs/src/acos.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * acos double precision arc cosine X * X * KEY WORDS X * X * acos X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns double precision arc cosine of double precision X * floating point argument. X * X * USAGE X * X * double acos (x) X * double x; X * X * REFERENCES X * X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-1. X * X * RESTRICTIONS X * X * The maximum relative error for the approximating polynomial X * in atan is 10**(-16.82). However, this assumes exact arithmetic X * in the polynomial evaluation. Additional rounding and X * truncation errors may occur as the argument is reduced X * to the range over which the polynomial approximation X * is valid, and as the polynomial is evaluated using X * finite-precision arithmetic. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes arccosine(x) from: X * X * (1) If x < -1.0 or x > +1.0 then call X * matherr and return 0.0 by default. X * X * (2) If x = 0.0 then acos(x) = PI/2. X * X * (3) If x = 1.0 then acos(x) = 0.0 X * X * (4) If x = -1.0 then acos(x) = PI. X * X * (5) If 0.0 < x < 1.0 then X * acos(x) = atan(Y) X * Y = sqrt[1-(x**2)] / x X * X * (4) If -1.0 < x < 0.0 then X * acos(x) = atan(Y) + PI X * Y = sqrt[1-(x**2)] / x X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "acos"; X X Xdouble acos (x) Xdouble x; X{ X double y; X extern double atan(); X extern double sqrt(); X auto struct exception xcpt; X X DBUG_ENTER (funcname); X DBUG_3 ("acosin", "arg %le", x); X if ( x > 1.0 || x < -1.0) { X xcpt.type = DOMAIN; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: DOMAIN error\n", funcname); X errno = EDOM; X xcpt.retval = 0.0; X } X } else if (x == 0.0) { X xcpt.retval = HALFPI; X } else if (x == 1.0) { X xcpt.retval = 0.0; X } else if (x == -1.0) { X xcpt.retval = PI; X } else { X y = atan ( sqrt (1.0 - (x * x)) / x ); X if (x > 0.0) { X xcpt.retval = y; X } else { X xcpt.retval = y + PI; X } X } X DBUG_3 ("acosout", "result %le", x); X DBUG_RETURN (x); X} END_OF_funcs/src/acos.c if test 2942 -ne `wc -c funcs/src/acosh.c <<'END_OF_funcs/src/acosh.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * acosh double precision hyperbolic arc cosine X * X * KEY WORDS X * X * acosh X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision hyperbolic arc cosine of double precision X * floating point number. X * X * USAGE X * X * double acosh (x) X * double x; X * X * RESTRICTIONS X * X * The range of the ACOSH function is all real numbers greater X * than or equal to 1.0 however large arguments may cause X * overflow in the x squared portion of the function evaluation. X * X * For precision information refer to documentation of the X * floating point library primatives called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes acosh(x) from: X * X * 1. If x < 1.0 then report illegal X * argument and return zero. X * X * 2. If x > sqrt(MAXDOUBLE) then X * set x = sqrt(MAXDOUBLE and X * continue after reporting overflow. X * X * 3. acosh(x) = log [x+sqrt(x**2 - 1)] X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "acosh"; X X Xdouble acosh (x) Xdouble x; X{ X auto struct exception xcpt; X extern double log (); X extern double sqrt (); X X DBUG_ENTER (funcname); X DBUG_3 ("acoshin", "arg %le", x); X if (x < 1.0) { X xcpt.type = DOMAIN; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: DOMAIN error\n", funcname); X errno = ERANGE; X xcpt.retval = 0.0; X } X } else if (x > SQRT_MAXDOUBLE) { X xcpt.type = OVERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: OVERFLOW error\n", funcname); X errno = ERANGE; X x = SQRT_MAXDOUBLE; X xcpt.retval = log (2* SQRT_MAXDOUBLE); X } X } else { X xcpt.retval = log (x + sqrt (x*x - 1.0)); X } X DBUG_3 ("acoshout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} X END_OF_funcs/src/acosh.c if test 2642 -ne `wc -c funcs/src/asin.c <<'END_OF_funcs/src/asin.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * asin double precision arc sine X * X * KEY WORDS X * X * asin X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns double precision arc sine of double precision X * floating point argument. X * X * If argument is less than -1.0 or greater than +1.0, calls X * matherr with a DOMAIN error. If matherr does not handle X * the error then prints error message and returns 0. X * X * USAGE X * X * double asin (x) X * double x; X * X * REFERENCES X * X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-2. X * X * RESTRICTIONS X * X * For precision information refer to documentation of the floating X * point library primatives called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes arcsine(x) from: X * X * (1) If x < -1.0 or x > +1.0 then X * call matherr and return 0.0 by default. X * X * (2) If x = 0.0 then asin(x) = 0.0 X * X * (3) If x = 1.0 then asin(x) = PI/2. X * X * (4) If x = -1.0 then asin(x) = -PI/2 X * X * (5) If -1.0 < x < 1.0 then X * asin(x) = atan(y) X * y = x / sqrt[1-(x**2)] X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "asin"; X X Xdouble asin (x) Xdouble x; X{ X extern double atan (); X extern double sqrt (); X struct exception xcpt; X X DBUG_ENTER (funcname); X DBUG_3 ("asinin", "arg %le", x); X if ( x > 1.0 || x < -1.0) { X xcpt.type = DOMAIN; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: DOMAIN error\n", funcname); X errno = EDOM; X xcpt.retval = 0.0; X } X } else if (x == 0.0) { X xcpt.retval = 0.0; X } else if (x == 1.0) { X xcpt.retval = HALFPI; X } else if (x == -1.0) { X xcpt.retval = -HALFPI; X } else { X xcpt.retval = atan ( x / sqrt (1.0 - (x * x)) ); X } X DBUG_3 ("asinout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/asin.c if test 2675 -ne `wc -c funcs/src/asinh.c <<'END_OF_funcs/src/asinh.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * asinh double precision hyperbolic arc sine X * X * KEY WORDS X * X * asinh X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision hyperbolic arc sine of double precision X * floating point number. X * X * USAGE X * X * double asinh (x) X * double x; X * X * RESTRICTIONS X * X * The domain of the ASINH function is the entire real axis X * however the evaluation of x squared may cause overflow X * for large magnitudes. X * X * For precision information refer to documentation of the X * floating point library routines called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes asinh(x) from: X * X * 1. Let xmax = sqrt(MAXDOUBLE - 1) X * X * 2. If x < -xmax or xmax < x then X * let x = xmax and flag overflow. X * X * 3. asinh(x) = log [x+sqrt(x**2 + 1)] X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "asinh"; X X Xdouble asinh (x) Xdouble x; X{ X auto struct exception xcpt; X extern double log (); X extern double sqrt (); X X DBUG_ENTER (funcname); X DBUG_3 ("asinhin", "arg %le", x); X if (x < -SQRT_MAXDOUBLE || x > SQRT_MAXDOUBLE) { X xcpt.type = OVERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: OVERFLOW error\n", funcname); X errno = ERANGE; X xcpt.retval = log (2 * SQRT_MAXDOUBLE); X } X } else { X xcpt.retval = log (x + sqrt(x*x + 1.0)); X } X DBUG_3 ("asinhout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/asinh.c if test 2315 -ne `wc -c funcs/src/atan2.c <<'END_OF_funcs/src/atan2.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * atan2 double precision arc tangent of two arguments X * X * KEY WORDS X * X * atan2 X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns double precision arc tangent of two X * double precision floating point arguments ( atan(Y/X) ). X * X * USAGE X * X * double atan2(x,y) X * double x; X * double y; X * X * REFERENCES X * X * Fortran 77 user's guide, Digital Equipment Corp. pp B-4. X * X * RESTRICTIONS X * X * Note that the argument usage is exactly the reverse of the X * common FORTRAN usage where atan2(x,y) computes atan(x/y). X * The usage here is less confusing than remembering that x is X * really y and y is really x. X * X * For precision information refer to documentation of the X * other floating point library routines called. X * X * PROGRAMMER X * X * Fred Fish X * Tempe, Az 85281 X * X * INTERNALS X * X * Computes atan(y/x) from: X * X * 1. If x = 0 then X * atan(x,y) = PI/2 * (sign(y)) X * X * 2. If x > 0 then X * atan(x,y) = atan(y/x) X * X * 3. If x < 0 then atan2(x,y) = X * PI*(sign(y)) + atan(y/x) X * X */ X X#include X#include X#include "pml.h" X X Xdouble atan2 (x, y) Xdouble x; Xdouble y; X{ X double result; X extern double sign(); X extern double atan(); X X ENTER ("atan2"); X DEBUG4 ("atan2in", "x = %le y = %le", x, y); X if (x == 0.0) { X result = sign (HALFPI, y); X } else if (x > 0.0) { X result = atan (y/x); X } else { X result = atan (y/x) + sign (PI, y); X } X DEBUG3 ("atan2out", "result %le", result); X LEAVE (); X return (result); X} END_OF_funcs/src/atan2.c if test 2363 -ne `wc -c funcs/src/atanh.c <<'END_OF_funcs/src/atanh.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * atanh double precision hyperbolic arc tangent X * X * KEY WORDS X * X * atanh X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision hyperbolic arc tangent of double precision X * floating point number. X * X * USAGE X * X * double atanh (x) X * double x; X * X * RESTRICTIONS X * X * The range of the atanh function is -1.0 to +1.0 exclusive. X * Certain pathological cases near 1.0 may cause overflow X * in evaluation of 1+x / 1-x, depending upon machine exponent X * range and mantissa precision. X * X * For precision information refer to documentation of the X * other floating point library routines called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes atanh(x) from: X * X * 1. If x <= -1.0 or x >= 1.0 X * then report argument out of X * range and return 0.0 X * X * 2. atanh(x) = 0.5 * log((1+x)/(1-x)) X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "atanh"; X X Xdouble atanh (x) Xdouble x; X{ X auto struct exception xcpt; X extern double log (); X X DBUG_ENTER (funcname); X DBUG_3 ("atanhin", "arg %le", x); X if (x <= -1.0 || x >= 1.0) { X xcpt.type = DOMAIN; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: DOMAIN error\n", funcname); X errno = ERANGE; X xcpt.retval = 0.0; X } X } else { X xcpt.retval = 0.5 * log ((1+x) / (1-x)); X } X DBUG_3 ("atanhout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/atanh.c if test 2301 -ne `wc -c funcs/src/cacos.c <<'END_OF_funcs/src/cacos.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * cacos complex double precision arc cosine X * X * KEY WORDS X * X * cacos X * complex functions X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Computes double precision complex arc cosine of X * a double precision complex argument. X * X * USAGE X * X * COMPLEX cacos (z) X * COMPLEX z; X * X * PROGRAMMER X * X * Fred Fish X * Tempe,Az 85281 X * (602) 966-8871 X * X * INTERNALS X * X * Computes complex arc cosine of z = x + jy from: X * X * cacos(z) = -j * clog(z + j * csqrt(1-z*z)) X * X */ X X#include X#include X#include "pml.h" X X XCOMPLEX cacos (z) XCOMPLEX z; X{ X COMPLEX temp; X double swaptemp; X extern COMPLEX cmult (), csqrt (), clog (); X X ENTER ("cacos"); X DEBUG4 ("cacosin", "arg %le %le", z.real, z.imag); X temp = cmult(z, z); X temp.real = 1.0 - temp.real; X temp.imag = -temp.imag; X temp = csqrt (temp); X swaptemp = temp.real; X temp.real = -temp.imag; X temp.imag = swaptemp; X temp.real += z.real; X temp.imag += z.imag; X temp = clog (temp); X z.real = temp.imag; X z.imag = -temp.real; X DEBUG4 ("cacosout", "result %le %le", z.real, z.imag); X LEAVE (); X return (z); X} X END_OF_funcs/src/cacos.c if test 1983 -ne `wc -c funcs/src/cdiv.c <<'END_OF_funcs/src/cdiv.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * cdiv double precision complex division X * X * KEY WORDS X * X * cdiv X * complex functions X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Computes double precision complex result of division of X * first double precision complex argument by second double X * precision complex argument. X * X * USAGE X * X * COMPLEX cdiv (numerator, denominator) X * COMPLEX numerator; X * COMPLEX denominator; X * X * PROGRAMMER X * X * Fred Fish X * Tempe, Az 85281 X * (602) 966-8871 X * X * INTERNALS X * X * Computes cdiv(znum,zden) from: X * X * 1. Let znum = a + j b X * Let zden = c + j d X * X * 2. denom = c*c + d*d X * X * 3. If denom is zero then log error, X * set r_cdiv = maximum floating value, X * i_cdiv = 0, and go to step 5. X * X * 4. r_cdiv = (a*c + b*d) / denom X * i_cdiv = (c*b - a*d) / denom X * X * 5. Then cdiv(znum,zden) = r_cdiv + j i_cdiv X * X */ X X#include X#include X#include "pml.h" X X XCOMPLEX cdiv (znum, zden) XCOMPLEX znum; XCOMPLEX zden; X{ X COMPLEX result; X double denom; X X ENTER ("cdiv"); X DEBUG4 ("cdivin", "arg1 %le %le", znum.real, znum.imag); X DEBUG4 ("cdivin", "arg2 %le %le", zden.real, zden.imag); X denom = (zden.real * zden.real) + (zden.imag * zden.imag); X if (denom == 0.0) { X/**** X pmlerr(C_DIV_ZERO); X result.real = MAX_POS_DBLF; X******/ X result.real = 0.0; /* TERRIBLY WRONG! */ X result.imag = 0.0; X } else { X result.real = ((znum.real*zden.real)+(znum.imag*zden.imag)) / denom; X result.imag = ((zden.real*znum.imag)-(znum.real*zden.imag)) / denom; X } X DEBUG4 ("cdivout", "result %le %le", result.real, result.imag); X LEAVE (); X return (result); X} END_OF_funcs/src/cdiv.c if test 2460 -ne `wc -c funcs/src/cmult.c <<'END_OF_funcs/src/cmult.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * cmult double precision complex multiplication X * X * KEY WORDS X * X * cmult X * complex functions X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Computes double precision complex result of multiplication of X * first double precision complex argument by second double X * precision complex argument. X * X * USAGE X * X * COMPLEX cmult (z1, z2) X * COMPLEX z1; X * COMPLEX z2; X * X * PROGRAMMER X * X * Fred Fish X * Tempe, Az 85281 X * (602) 966-8871 X * X * INTERNALS X * X * Computes cmult(z1,z2) from: X * X * 1. Let z1 = a + j b X * Let z2 = c + j d X * X * 2. r_cmult = (a*c - b*d) X * i_cmult = (a*d + c*b) X * X * 3. Then cmult(z1,z2) = r_cmult + j i_cmult X * X */ X X#include X#include X#include "pml.h" X X XCOMPLEX cmult (z1, z2) XCOMPLEX z1; XCOMPLEX z2; X{ X COMPLEX result; X X ENTER ("cmult"); X DEBUG4 ("cmultin", "arg1 %le %le", z1.real, z1.imag); X DEBUG4 ("cmultin", "arg2 %le %le", z2.real, z2.imag); X result.real = (z1.real * z2.real) - (z1.imag * z2.imag); X result.imag = (z1.real * z2.imag) + (z2.real * z1.imag); X DEBUG4 ("cmultout", "result %le %le", result.real, result.imag); X LEAVE (); X return (result); X} END_OF_funcs/src/cmult.c if test 1994 -ne `wc -c funcs/src/cosh.c <<'END_OF_funcs/src/cosh.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * cosh double precision hyperbolic cosine X * X * KEY WORDS X * X * cosh X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision hyperbolic cosine of double precision X * floating point number. X * X * USAGE X * X * double cosh (x) X * double x; X * X * REFERENCES X * X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-4 X * X * RESTRICTIONS X * X * Inputs greater than log(MAXDOUBLE) result in overflow. X * Inputs less than log(MINDOUBLE) result in underflow. X * X * For precision information refer to documentation of the X * floating point library routines called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes hyperbolic cosine from: X * X * cosh(X) = 0.5 * (exp(X) + exp(-X)) X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "cosh"; X X Xdouble cosh (x) Xdouble x; X{ X auto struct exception xcpt; X extern double exp (); X X DBUG_ENTER (funcname); X DBUG_3 ("coshin", "arg %le", x); X if (x > LOGE_MAXDOUBLE) { X xcpt.type = OVERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: OVERFLOW error\n", funcname); X errno = ERANGE; X xcpt.retval = MAXDOUBLE; X } X } else if (x < LOGE_MINDOUBLE) { X xcpt.type = UNDERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: UNDERFLOW error\n", funcname); X errno = ERANGE; X xcpt.retval = MINDOUBLE; X } X } else { X x = exp (x); X xcpt.retval = 0.5 * (x + 1.0/x); X } X DBUG_3 ("coshout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/cosh.c if test 2424 -ne `wc -c funcs/src/csqrt.c <<'END_OF_funcs/src/csqrt.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * csqrt complex double precision square root X * X * KEY WORDS X * X * csqrt X * machine independent routines X * complex functions X * math libraries X * X * DESCRIPTION X * X * Computes double precision complex square root of X * a double precision complex argument. X * X * USAGE X * X * COMPLEX csqrt (z) X * COMPLEX z; X * X * REFERENCES X * X * Fortran 77 user's guide, Digital Equipment Corp. pp B-12 X * X * RESTRICTIONS X * X * The relative error in the double precision square root X * computation is 10**(-30.1) after three applications X * of Heron's iteration for the square root. X * X * However, this assumes exact arithmetic in the iterations X * and initial approximation. Additional errors may occur X * due to truncation, rounding, or machine precision limits. X * X * PROGRAMMER X * X * Fred Fish X * Tempe, Az 85281 X * (602) 966-8871 X * X * INTERNALS X * X * Computes complex square root of z = x + j y from: X * X * 1. If z = 0 + j 0 then return z. X * X * 2. root = sqrt((dabs(x) + cabs(z)) / 2) X * X * 3. q = y / (2 * root) X * X * 4. If x >= 0 then X * csqrt(z) = (root,q) X * X * 5. If x < 0 and y >= 0 then X * csqrt(z) = (q,root) X * X * 6. If x < 0 and y < 0 then X * csqrt(z) = (-q,root) X * X */ X X#include X#include X#include "pml.h" X X XCOMPLEX csqrt (z) XCOMPLEX z; X{ X double root, q; X extern double dabs(), sqrt(), cabs (); X X ENTER ("csqrt"); X DEBUG4 ("csqrtin", "arg %le %le", z.real, z.imag); X if (z.real != 0.0 || z.imag != 0.0) { X root = sqrt (0.5 * (dabs (z.real) + cabs (z))); X q = z.imag / (2.0 * root); X if (z.real >= 0.0) { X z.real = root; X z.imag = q; X } else if (z.imag < 0.0) { X z.real = -q; X z.imag = -root; X } else { X z.real = q; X z.imag = root; X } X } X DEBUG4 ("csqrtout", "result %le %le", z.real, z.imag); X LEAVE (); X return (z); X} END_OF_funcs/src/csqrt.c if test 2679 -ne `wc -c funcs/src/csubt.c <<'END_OF_funcs/src/csubt.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * csubt double precision complex subtraction X * X * KEY WORDS X * X * csubt X * machine independent routines X * complex functions X * math libraries X * X * DESCRIPTION X * X * Computes double precision complex result of subtraction of X * second double precision complex argument from first double X * precision complex argument. X * X * Note that the complex subtraction function is so simple that X * it would not normally be called as a function but simply X * done "inline". It is supplied mostly for completeness. X * X * USAGE X * X * COMPLEX csubt (z1, z2) X * COMPLEX z1; X * COMPLEX z2; X * X * PROGRAMMER X * X * Fred Fish X * Tempe, Az X * (602) 966-8871 X * X * INTERNALS X * X * Computes csubt (z1, z2) from: X * X * 1. Let z1 = a + j b X * Let z2 = c + j d X * X * 2. Then csubt(z1,z2) = (a - c) + j (b - d) X * X */ X X#include X#include X#include "pml.h" X X XCOMPLEX csubt (z1, z2) XCOMPLEX z1; XCOMPLEX z2; X{ X ENTER ("csubt"); X DEBUG4 ("csubtin", "arg %le %le", z1.real, z1.imag); X DEBUG4 ("csubtin", "arg2 %le %le", z2.real, z2.imag); X z1.real -= z2.real; X z1.imag -= z2.imag; X DEBUG4 ("csubtout", "result %le %le", z1.real, z1.imag); X LEAVE (); X return (z1); X} END_OF_funcs/src/csubt.c if test 2005 -ne `wc -c funcs/src/poly.c <<'END_OF_funcs/src/poly.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * poly double precision polynomial evaluation X * X * KEY WORDS X * X * poly X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Evaluates a polynomial and returns double precision X * result. Is passed a the order of the polynomial, X * a pointer to an array of double precision polynomial X * coefficients (in ascending order), and the independent X * variable. X * X * USAGE X * X * double poly (order, coeffs, x) X * int order; X * double *coeffs; X * double x; X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Evalates the polynomial using recursion and the form: X * X * P(x) = P0 + x(P1 + x(P2 +...x(Pn))) X * X */ X X#include X#include X#include "pml.h" X X Xdouble poly (order, coeffs, x) Xregister int order; Xdouble *coeffs; Xdouble x; X{ X auto double curr_coeff; X auto double rtn_value; X X DBUG_ENTER ("poly"); X DBUG_5 ("polyin", "args %d %#x %le", order, coeffs, x); X if (order <= 0) { X rtn_value = *coeffs; X } else { X curr_coeff = *coeffs; /* Bug in Unisoft's compiler. Does not */ X coeffs++; /* generate good code for *coeffs++ */ X rtn_value = curr_coeff + x * poly (--order, coeffs, x); X } X DBUG_3 ("polyout", "result %le", rtn_value); X DBUG_RETURN (rtn_value); X} END_OF_funcs/src/poly.c if test 2049 -ne `wc -c funcs/src/sinh.c <<'END_OF_funcs/src/sinh.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * sinh double precision hyperbolic sine X * X * KEY WORDS X * X * sinh X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision hyperbolic sine of double precision X * floating point number. X * X * USAGE X * X * double sinh (x) X * double x; X * X * REFERENCES X * X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-5 X * X * RESTRICTIONS X * X * Inputs greater than ln(MAXDOUBLE) result in overflow. X * Inputs less than ln(MINDOUBLE) result in underflow. X * X * For precision information refer to documentation of the X * floating point library routines called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes hyperbolic sine from: X * X * sinh(x) = 0.5 * (exp(x) - 1.0/exp(x)) X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "sinh"; X X Xdouble sinh (x) Xdouble x; X{ X extern double exp (); X auto struct exception xcpt; X X DBUG_ENTER (funcname); X DBUG_3 ("sinhin", "arg %le", x); X if (x > LOGE_MAXDOUBLE) { X xcpt.type = OVERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: OVERFLOW error\n", funcname); X errno = ERANGE; X xcpt.retval = MAXDOUBLE; X } X } else if (x < LOGE_MINDOUBLE) { X xcpt.type = UNDERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: UNDERFLOW error\n", funcname); X errno = ERANGE; X xcpt.retval = MINDOUBLE; X } X } else { X x = exp (x); X xcpt.retval = 0.5 * (x - (1.0 / x)); X } X DBUG_3 ("sinhout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/sinh.c if test 2423 -ne `wc -c funcs/src/tan.c <<'END_OF_funcs/src/tan.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * tan Double precision tangent X * X * KEY WORDS X * X * tan X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns tangent of double precision floating point number. X * X * USAGE X * X * double tan (x) X * double x; X * X * INTERNALS X * X * Computes the tangent from tan(x) = sin(x) / cos(x). X * X * If cos(x) = 0 and sin(x) >= 0, then returns largest X * floating point number on host machine. X * X * If cos(x) = 0 and sin(x) < 0, then returns smallest X * floating point number on host machine. X * X * REFERENCES X * X * Fortran IV plus user's guide, Digital Equipment Corp. pp. B-8 X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "tan"; X Xdouble tan (x) Xdouble x; X{ X double sinx; X double cosx; X auto struct exception xcpt; X extern double sin (); X extern double cos(); X X DBUG_ENTER (funcname); X DBUG_3 ("tanin", "arg %le", x); X sinx = sin (x); X cosx = cos (x); X if (cosx == 0.0) { X xcpt.type = OVERFLOW; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: OVERFLOW error\n", funcname); X errno = ERANGE; X if (sinx >= 0.0) { X xcpt.retval = MAXDOUBLE; X } else { X xcpt.retval = -MAXDOUBLE; X } X } X } else { X xcpt.retval = sinx / cosx; X } X DBUG_3 ("tanout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/tan.c if test 2208 -ne `wc -c funcs/src/tanh.c <<'END_OF_funcs/src/tanh.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X X/* X * FUNCTION X * X * tanh double precision hyperbolic tangent X * X * KEY WORDS X * X * tanh X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision hyperbolic tangent of double precision X * floating point number. X * X * USAGE X * X * double tanh (x) X * double x; X * X * REFERENCES X * X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-5 X * X * RESTRICTIONS X * X * For precision information refer to documentation of the X * floating point library routines called. X * X * PROGRAMMER X * X * Fred Fish X * X * INTERNALS X * X * Computes hyperbolic tangent from: X * X * tanh(x) = 1.0 for x > TANH_MAXARG X * = -1.0 for x < -TANH_MAXARG X * = sinh(x) / cosh(x) otherwise X * X */ X X#include X#include X#include "pml.h" X Xstatic char funcname[] = "tanh"; X X Xdouble tanh (x) Xdouble x; X{ X auto struct exception xcpt; X register int positive; X extern double sinh (); X extern double cosh (); X X DBUG_ENTER (funcname); X DBUG_3 ("tanhin", "arg %le", x); X if (x > TANH_MAXARG || x < -TANH_MAXARG) { X if (x > 0.0) { X positive = 1; X } else { X positive = 0; X } X xcpt.type = PLOSS; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: PLOSS error\n", funcname); X errno = ERANGE; X if (positive) { X xcpt.retval = 1.0; X } else { X xcpt.retval = -1.0; X } X } X } else { X xcpt.retval = sinh (x) / cosh (x); X } X DBUG_3 ("tanhout", "result %le", xcpt.retval); X return (xcpt.retval); X} END_OF_funcs/src/tanh.c if test 2312 -ne `wc -c funcs/unused/dint.c <<'END_OF_funcs/unused/dint.c' X/************************************************************************ X * * X * N O T I C E * X * * X * Copyright Abandoned, 1987, Fred Fish * X * * X * This previously copyrighted work has been placed into the * X * public domain by the author (Fred Fish) and may be freely used * X * for any purpose, private or commercial. I would appreciate * X * it, as a courtesy, if this notice is left in all copies and * X * derivative works. Thank you, and enjoy... * X * * X * The author makes no warranty of any kind with respect to this * X * product and explicitly disclaims any implied warranties of * X * merchantability or fitness for any particular purpose. * X * * X ************************************************************************ X */ X X/* X * FUNCTION X * X * dint double precision integer portion X * X * KEY WORDS X * X * dint X * machine dependent routines X * math libraries X * X * DESCRIPTION X * X * Returns integer portion of double precision number as double X * precision number. X * X * USAGE X * X * double dint(x) X * double x; X * X * PROGRAMMER X * X * Fred Fish X * Tempe, Az 85281 X * (602) 966-8871 X * X * RESTRICTIONS X * X * The current DEC-20 C system treats double as float. This X * routine will need to be modified when true double precision X * is implemented. X * X */ X X#include X#include X#include "pml.h" X X X#ifdef PDP10 X X#define W1_FBITS 27 /* Number of fractional bits in word 1 */ X#define WORD_MASK 0777777777777 /* Mask for all 36 bits of first word */ X Xdouble dint(x) Xdouble x; X{ X register int *vpntr, exponent; X int xexp(); X X vpntr = &x; X if ((exponent=xexp(x)) <= 0) { X x = 0.0; X } else if (exponent <= W1_FBITS) { X *vpntr &= (WORD_MASK << (W1_FBITS - exponent)); X } else { X pmlerr(DINT_2BIG); X x = 0.0; X } X if (x < 0.0) { X x += 1.0; X } X return (x); X} X X#else X X#define W1_FBITS 24 /* (NOTE HIDDEN BIT NORMALIZATION) */ X#define W2_FBITS 32 /* Number of fractional bits in word 2 */ X#define WORD_MASK 0xFFFFFFFF /* Mask for each long word of double */ X Xstatic union { X double dval; X long lval[2]; X} share; X Xdouble dint(x) Xdouble x; X{ X int exponent, xexp(), fbitdown; X register long *lpntr; X X lpntr = &share.lval[0]; X share.dval = x; X if ((exponent=xexp(x)) <= 0) { X share.dval = 0.0; X } else if (exponent <= W1_FBITS) { X *lpntr++ &= (WORD_MASK << (W1_FBITS - exponent)); X *lpntr++ = 0; X } else if (exponent <= (fbitdown = W1_FBITS+W2_FBITS)) { X lpntr++; X *lpntr++ &= (WORD_MASK << (fbitdown - exponent)); X } else { X pmlerr(DINT_2BIG); X } X return (share.dval); X} X#endif END_OF_funcs/unused/dint.c if test 2630 -ne `wc -c