Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP Path: utzoo!utgpu!water!watnot!watmath!clyde!rutgers!husc6!seismo!nbires!hao!noao!mcdsun!fnf From: fnf@mcdsun.UUCP Newsgroups: net.sources Subject: Portable Math Library in C (Part 3 of 6) Message-ID: <291@mcdsun.UUCP> Date: Fri, 10-Apr-87 18:42:53 EST Article-I.D.: mcdsun.291 Posted: Fri Apr 10 18:42:53 1987 Date-Received: Sun, 12-Apr-87 21:36:28 EST Organization: Motorola Microcomputer Division, Tempe, Az. Lines: 1770 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/src/atan.c <<'END_OF_funcs/src/atan.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 * atan double precision arc tangent X * X * KEY WORDS X * X * atan X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns double precision arc tangent of double precision X * floating point argument. X * X * USAGE X * X * double atan (x) X * double x; X * X * REFERENCES X * X * Fortran 77 user's guide, Digital Equipment Corp. pp B-3 X * X * Computer Approximations, J.F. Hart et al, John Wiley & Sons, X * 1968, pp. 120-130. X * X * RESTRICTIONS X * X * The maximum relative error for the approximating polynomial X * 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 arctangent(x) from: X * X * (1) If x < 0 then negate x, perform steps X * 2, 3, and 4, and negate the returned X * result. This makes use of the identity X * atan(-x) = -atan(x). X * X * (2) If argument x > 1 at this point then X * test to be sure that x can be inverted X * without underflowing. If not, reduce X * x to largest possible number that can X * be inverted, issue warning, and continue. X * Perform steps 3 and 4 with arg = 1/x X * and subtract returned result from X * pi/2. This makes use of the identity X * atan(x) = pi/2 - atan(1/x) for x>0. X * X * (3) At this point 0 <= x <= 1. If X * x > tan(pi/12) then perform step 4 X * with arg = (x*sqrt(3)-1)/(sqrt(3)+x) X * and add pi/6 to returned result. This X * last transformation maps arguments X * greater than tan(pi/12) to arguments X * in range 0...tan(pi/12). X * X * (4) At this point the argument has been X * transformed so that it lies in the X * range -tan(pi/12)...tan(pi/12). X * Since very small arguments may cause X * underflow in the polynomial evaluation, X * a final check is performed. If the X * argument is less than the underflow X * bound, the function returns x, since X * atan(x) approaches asin(x) which X * approaches x, as x goes to zero. X * X * (5) atan(x) is now computed by one of the X * approximations given in the cited X * reference (Hart). That is: X * X * atan(x) = x * SUM [ C[i] * x**(2*i) ] X * over i = {0,1,...8}. X * X * Where: X * X * C[0] = .9999999999999999849899 X * C[1] = -.333333333333299308717 X * C[2] = .1999999999872944792 X * C[3] = -.142857141028255452 X * C[4] = .11111097898051048 X * C[5] = -.0909037114191074 X * C[6] = .0767936869066 X * C[7] = -.06483193510303 X * C[8] = .0443895157187 X * (coefficients from HART table #4945 pg 267) X * X */ X X#include X#include X#include "pml.h" X Xstatic double atan_coeffs[] = { X .9999999999999999849899, /* P0 must be first */ X -.333333333333299308717, X .1999999999872944792, X -.142857141028255452, X .11111097898051048, X -.0909037114191074, X .0767936869066, X -.06483193510303, X .0443895157187 /* Pn must be last */ X}; X Xstatic char funcname[] = "atan"; X X#define LAST_BOUND 0.2679491924311227074725 /* tan (PI/12) */ X X Xdouble atan (x) Xdouble x; X{ X register int order; X double xt2; X double t1; X double t2; X extern double poly (); X auto struct exception xcpt; X X DBUG_ENTER (funcname); X DBUG_3 ("atanin", "arg %le", x); X if (x < 0.0) { X xcpt.retval = -(atan (-x)); X } else if (x > 1.0) { X if (x < MAXDOUBLE && x > -MAXDOUBLE) { X xcpt.retval = HALFPI - atan (1.0 / x); X } else { 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 = EDOM; X xcpt.retval = 0.0; X } X } X } else if (x > LAST_BOUND) { X t1 = x * SQRT3 - 1.0; X t2 = SQRT3 + x; X xcpt.retval = SIXTHPI + atan (t1 / t2); X } else if (x < X16_UNDERFLOWS) { 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 = EDOM; X xcpt.retval = x; X } X } else { X xt2 = x * x; X order = sizeof (atan_coeffs) / sizeof(double); X order -= 1; X xcpt.retval = x * poly (order, atan_coeffs, xt2); X } X DBUG_3 ("atanout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/atan.c if test 5187 -ne `wc -c funcs/src/cos.c <<'END_OF_funcs/src/cos.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 * cos double precision cosine X * X * KEY WORDS X * X * cos X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns double precision cosine of double precision X * floating point argument. X * X * USAGE X * X * double cos (x) X * double x; X * X * REFERENCES X * X * Computer Approximations, J.F. Hart et al, John Wiley & Sons, X * 1968, pp. 112-120. X * X * RESTRICTIONS X * X * The sin and cos routines are interactive in the sense that X * in the process of reducing the argument to the range -PI/4 X * to PI/4, each may call the other. Ultimately one or the X * other uses a polynomial approximation on the reduced X * argument. The sin approximation has a maximum relative error X * of 10**(-17.59) and the cos approximation has a maximum X * relative error of 10**(-16.18). X * X * These error bounds assume 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 cos(x) from: X * X * (1) Reduce argument x to range -PI to PI. X * X * (2) If x > PI/2 then call cos recursively X * using relation cos(x) = -cos(x - PI). X * X * (3) If x < -PI/2 then call cos recursively X * using relation cos(x) = -cos(x + PI). X * X * (4) If x > PI/4 then call sin using X * relation cos(x) = sin(PI/2 - x). X * X * (5) If x < -PI/4 then call cos using X * relation cos(x) = sin(PI/2 + x). X * X * (6) If x would cause underflow in approx X * evaluation arithmetic then return X * sqrt(1.0 - x**2). X * X * (7) By now x has been reduced to range X * -PI/4 to PI/4 and the approximation X * from HART pg. 119 can be used: X * X * cos(x) = ( p(y) / q(y) ) X * Where: X * X * y = x * (4/PI) X * X * p(y) = SUM [ Pj * (y**(2*j)) ] X * over j = {0,1,2,3} X * X * q(y) = SUM [ Qj * (y**(2*j)) ] X * over j = {0,1,2,3} X * X * P0 = 0.12905394659037374438571854e+7 X * P1 = -0.3745670391572320471032359e+6 X * P2 = 0.134323009865390842853673e+5 X * P3 = -0.112314508233409330923e+3 X * Q0 = 0.12905394659037373590295914e+7 X * Q1 = 0.234677731072458350524124e+5 X * Q2 = 0.2096951819672630628621e+3 X * Q3 = 1.0000... X * (coefficients from HART table #3843 pg 244) X * X * X * **** NOTE **** The range reduction relations used in X * this routine depend on the final approximation being valid X * over the negative argument range in addition to the positive X * argument range. The particular approximation chosen from X * HART satisfies this requirement, although not explicitly X * stated in the text. This may not be true of other X * approximations given in the reference. X * X */ X X# include X# include X# include "pml.h" X Xstatic double cos_pcoeffs[] = { X 0.12905394659037374438e7, X -0.37456703915723204710e6, X 0.13432300986539084285e5, X -0.11231450823340933092e3 X}; X Xstatic double cos_qcoeffs[] = { X 0.12905394659037373590e7, X 0.23467773107245835052e5, X 0.20969518196726306286e3, X 1.0 X}; X Xstatic char funcname[] = "cos"; X X Xdouble cos (x) Xdouble x; X{ X auto double y; X auto double yt2; X auto double rtnval; X extern double mod (); X extern double sin (); X extern double sqrt (); X extern double poly (); X struct exception xcpt; X X DBUG_ENTER (funcname); X DBUG_3 ("cosin", "arg %le", x); X if (x < -PI || x > PI) { X x = mod (x, TWOPI); X if (x > PI) { X x = x - TWOPI; X } else if (x < -PI) { X x = x + TWOPI; X } X } X if (x > HALFPI) { X xcpt.retval = -(cos (x - PI)); X } else if (x < -HALFPI) { X xcpt.retval = -(cos (x + PI)); X } else if (x > FOURTHPI) { X xcpt.retval = sin (HALFPI - x); X } else if (x < -FOURTHPI) { X xcpt.retval = sin (HALFPI + x); X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) { X xcpt.retval = sqrt (1.0 - (x * x)); X } else { X y = x / FOURTHPI; X yt2 = y * y; X xcpt.retval = poly (3, cos_pcoeffs, yt2) / poly (3, cos_qcoeffs, yt2); X } X DBUG_3 ("cosout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/cos.c if test 5011 -ne `wc -c funcs/src/log.c <<'END_OF_funcs/src/log.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 * log double precision natural log X * X * KEY WORDS X * X * log X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision natural log of double precision X * floating point argument. X * X * USAGE X * X * double log (x) X * double x; X * X * REFERENCES X * X * Computer Approximations, J.F. Hart et al, John Wiley & Sons, X * 1968, pp. 105-111. X * X * RESTRICTIONS X * X * The absolute error in the approximating polynomial is X * 10**(-19.38). Note that contrary to DEC's assertion X * in the F4P user's guide, the error is absolute and not X * relative. X * X * This error bound 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 log(X) from: X * X * (1) If argument is zero then flag an error X * and return minus infinity (or rather its X * machine representation). X * X * (2) If argument is negative then flag an X * error and return minus infinity. X * X * (3) Given that x = m * 2**k then extract X * mantissa m and exponent k. X * X * (4) Transform m which is in range [0.5,1.0] X * to range [1/sqrt(2), sqrt(2)] by: X * X * s = m * sqrt(2) X * X * (5) Compute z = (s - 1) / (s + 1) X * X * (6) Now use the approximation from HART X * page 111 to find log(s): X * X * log(s) = z * ( P(z**2) / Q(z**2) ) X * X * Where: X * X * P(z**2) = SUM [ Pj * (z**(2*j)) ] X * over j = {0,1,2,3} X * X * Q(z**2) = SUM [ Qj * (z**(2*j)) ] X * over j = {0,1,2,3} X * X * P0 = -0.240139179559210509868484e2 X * P1 = 0.30957292821537650062264e2 X * P2 = -0.96376909336868659324e1 X * P3 = 0.4210873712179797145 X * Q0 = -0.120069589779605254717525e2 X * Q1 = 0.19480966070088973051623e2 X * Q2 = -0.89111090279378312337e1 X * Q3 = 1.0000 X * X * (coefficients from HART table #2705 pg 229) X * X * (7) Finally, compute log(x) from: X * X * log(x) = (k * log(2)) - log(sqrt(2)) + log(s) X * X */ X X#include X#include X#include "pml.h" X Xstatic double log_pcoeffs[] = { X -0.24013917955921050986e2, X 0.30957292821537650062e2, X -0.96376909336868659324e1, X 0.4210873712179797145 X}; X Xstatic double log_qcoeffs[] = { X -0.12006958977960525471e2, X 0.19480966070088973051e2, X -0.89111090279378312337e1, X 1.0000 X}; X Xstatic char funcname[] = "log"; X X Xdouble log (x) Xdouble x; X{ X auto int k; X auto double s; X auto double z; X auto double zt2; X auto double pqofz; X auto struct exception xcpt; X extern double frexp (); X extern double poly (); X X DBUG_ENTER (funcname); X DBUG_3 ("login", "arg %le", x); X if (x == 0.0) { X xcpt.type = SING; X xcpt.name = funcname; X xcpt.arg1 = x; X if (!matherr (&xcpt)) { X fprintf (stderr, "%s: SINGULARITY error\n", funcname); X errno = EDOM; X xcpt.retval = -MAXDOUBLE; X } X } else if (x < 0.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 = -MAXDOUBLE; X } X } else { X s = SQRT2 * frexp (x, &k); X DBUG_3 ("log", "k = %d", k); X DBUG_3 ("log", "s = %le", s); X z = (s - 1.0) / (s + 1.0); X DBUG_3 ("log", "z = %le", z); X zt2 = z * z; X DBUG_3 ("log", "zt2 = %le", zt2); X pqofz = z * (poly (3, log_pcoeffs, zt2) / poly (3, log_qcoeffs, zt2)); X DBUG_3 ("pqofz", "pqofz = %le", pqofz); X DBUG_3 ("log", "k = %d", k); X DBUG_3 ("log", "LN2 = %le", LN2); X x = k * LN2; X DBUG_3 ("log", "x = %le", x); X x -= LNSQRT2; X DBUG_3 ("log", "x = %le", x); X x += pqofz; X DBUG_3 ("log", "x = %le", x); X xcpt.retval = x; X } X DBUG_3 ("logout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/log.c if test 4649 -ne `wc -c funcs/src/pml.h <<'END_OF_funcs/src/pml.h' 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 * This file gets included with all of the floating point math X * library routines when they are compiled. Note that X * this is the proper place to put machine dependencies X * whenever possible. X * X * It should be pointed out that for simplicity's sake, the X * environment parameters are defined as floating point constants, X * rather than octal or hexadecimal initializations of allocated X * storage areas. This means that the range of allowed numbers X * may not exactly match the hardware's capabilities. For example, X * if the maximum positive double precision floating point number X * is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is X * defined to be 1.11E100 then the numbers between 1.11E100 and X * 1.11...E100 are considered to be undefined. For most X * applications, this will cause no problems. X * X * An alternate method is to allocate a global static "double" variable, X * say "maxdouble", and use a union declaration and initialization X * to initialize it with the proper bits for the EXACT maximum value. X * This was not done because the only compilers available to the X * author did not fully support union initialization features. X * X */ X X#ifndef NO_DBUG X# include X#else X# define DBUG_ENTER(a1) X# define DBUG_RETURN(a1) return(a1) X# define DBUG_VOID_RETURN return X# define DBUG_EXECUTE(keyword,a1) X# define DBUG_2(keyword,format) X# define DBUG_3(keyword,format,a1) X# define DBUG_4(keyword,format,a1,a2) X# define DBUG_5(keyword,format,a1,a2,a3) X# define DBUG_PUSH(a1) X# define DBUG_POP() X# define DBUG_PROCESS(a1) X# define DBUG_FILE (stderr) X#endif X X#include X#include X#include X X/* X * MAXDOUBLE => Maximum double precision number X * MINDOUBLE => Minimum double precision number X * DMAXEXP => Maximum exponent of a double X * DMINEXP => Minimum exponent of a double X */ X X#define MINDOUBLE (1.0/MAXDOUBLE) X#define LOG2_MAXDOUBLE (DMAXEXP + 1) X#define LOG2_MINDOUBLE (DMINEXP - 1) X#define LOGE_MAXDOUBLE (LOG2_MAXDOUBLE / LOG2E) X#define LOGE_MINDOUBLE (LOG2_MINDOUBLE / LOG2E) X X/* X * The following are hacks which should be fixed when I understand all X * the issues a little better. |tanh(TANH_MAXARG)| = 1.0 X */ X#define TANH_MAXARG 16 X#define SQRT_MAXDOUBLE 1.304380e19 X X#define PI 3.14159265358979323846 X#define TWOPI (2.0 * PI) X#define HALFPI (PI / 2.0) X#define FOURTHPI (PI / 4.0) X#define SIXTHPI (PI / 6.0) X#define LOG2E 1.4426950408889634074 /* Log to base 2 of e */ X#define LOG10E 0.4342944819032518276 X#define SQRT2 1.4142135623730950488 X#define SQRT3 1.7320508075688772935 X#define LN2 0.6931471805599453094 X#define LNSQRT2 0.3465735902799726547 X X X/* X * MC68000 HARDWARE DEPENDENCIES X * X * cc -DIEEE => uses IEEE floating point format X * X */ X X#ifdef IEEE X#define X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */ X#define X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows */ X#endif X X#define TRUE 1 /* This really should be in stdio.h */ X#define FALSE 0 /* This too */ X#define VOID void END_OF_funcs/src/pml.h if test 3819 -ne `wc -c funcs/src/sin.c <<'END_OF_funcs/src/sin.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 * sin double precision sine X * X * KEY WORDS X * X * sin X * machine independent routines X * trigonometric functions X * math libraries X * X * DESCRIPTION X * X * Returns double precision sine of double precision X * floating point argument. X * X * USAGE X * X * double sin (x) X * double x; X * X * REFERENCES X * X * Computer Approximations, J.F. Hart et al, John Wiley & Sons, X * 1968, pp. 112-120. X * X * RESTRICTIONS X * X * The sin and cos routines are interactive in the sense that X * in the process of reducing the argument to the range -PI/4 X * to PI/4, each may call the other. Ultimately one or the X * other uses a polynomial approximation on the reduced X * argument. The sin approximation has a maximum relative error X * of 10**(-17.59) and the cos approximation has a maximum X * relative error of 10**(-16.18). X * X * These error bounds assume 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 sin(x) from: X * X * (1) Reduce argument x to range -PI to PI. X * X * (2) If x > PI/2 then call sin recursively X * using relation sin(x) = -sin(x - PI). X * X * (3) If x < -PI/2 then call sin recursively X * using relation sin(x) = -sin(x + PI). X * X * (4) If x > PI/4 then call cos using X * relation sin(x) = cos(PI/2 - x). X * X * (5) If x < -PI/4 then call cos using X * relation sin(x) = -cos(PI/2 + x). X * X * (6) If x is small enough that polynomial X * evaluation would cause underflow X * then return x, since sin(x) X * approaches x as x approaches zero. X * X * (7) By now x has been reduced to range X * -PI/4 to PI/4 and the approximation X * from HART pg. 118 can be used: X * X * sin(x) = y * ( p(y) / q(y) ) X * Where: X * X * y = x * (4/PI) X * X * p(y) = SUM [ Pj * (y**(2*j)) ] X * over j = {0,1,2,3} X * X * q(y) = SUM [ Qj * (y**(2*j)) ] X * over j = {0,1,2,3} X * X * P0 = 0.206643433369958582409167054e+7 X * P1 = -0.18160398797407332550219213e+6 X * P2 = 0.359993069496361883172836e+4 X * P3 = -0.2010748329458861571949e+2 X * Q0 = 0.263106591026476989637710307e+7 X * Q1 = 0.3927024277464900030883986e+5 X * Q2 = 0.27811919481083844087953e+3 X * Q3 = 1.0000... X * (coefficients from HART table #3063 pg 234) X * X * X * **** NOTE **** The range reduction relations used in X * this routine depend on the final approximation being valid X * over the negative argument range in addition to the positive X * argument range. The particular approximation chosen from X * HART satisfies this requirement, although not explicitly X * stated in the text. This may not be true of other X * approximations given in the reference. X * X */ X X#include X#include X#include "pml.h" X Xstatic double sin_pcoeffs[] = { X 0.20664343336995858240e7, X -0.18160398797407332550e6, X 0.35999306949636188317e4, X -0.20107483294588615719e2 X}; X Xstatic double sin_qcoeffs[] = { X 0.26310659102647698963e7, X 0.39270242774649000308e5, X 0.27811919481083844087e3, X 1.0 X}; X Xstatic char funcname[] = "sin"; X X Xdouble sin (x) Xdouble x; X{ X double y; X double yt2; X double rtnval; X extern double mod (); X extern double cos (); X extern double poly (); X auto struct exception xcpt; X X DBUG_ENTER (funcname); X DBUG_3 ("sinin", "arg %le", x); X if (x < -PI || x > PI) { X x = mod (x, TWOPI); X if (x > PI) { X x = x - TWOPI; X } else if (x < -PI) { X x = x + TWOPI; X } X } X if (x > HALFPI) { X xcpt.retval = -(sin (x - PI)); X } else if (x < -HALFPI) { X xcpt.retval = -(sin (x + PI)); X } else if (x > FOURTHPI) { X xcpt.retval = cos (HALFPI - x); X } else if (x < -FOURTHPI) { X xcpt.retval = -(cos (HALFPI + x)); X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) { X xcpt.retval = x; X } else { X y = x / FOURTHPI; X yt2 = y * y; X xcpt.retval = y * (poly (3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2)); X } X DBUG_3 ("sinout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/sin.c if test 4983 -ne `wc -c funcs/src/sqrt.c <<'END_OF_funcs/src/sqrt.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 * sqrt double precision square root X * X * KEY WORDS X * X * sqrt X * machine independent routines X * math libraries X * X * DESCRIPTION X * X * Returns double precision square root of double precision X * floating point argument. X * X * USAGE X * X * double sqrt (x) X * double x; X * X * REFERENCES X * X * Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7. X * X * Computer Approximations, J.F. Hart et al, John Wiley & Sons, X * 1968, pp. 89-96. X * X * RESTRICTIONS X * X * The relative error 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 * X * INTERNALS X * X * Computes square root by: X * X * (1) Range reduction of argument to [0.5,1.0] X * by application of identity: X * X * sqrt(x) = 2**(k/2) * sqrt(x * 2**(-k)) X * X * k is the exponent when x is written as X * a mantissa times a power of 2 (m * 2**k). X * It is assumed that the mantissa is X * already normalized (0.5 =< m < 1.0). X * X * (2) An approximation to sqrt(m) is obtained X * from: X * X * u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m) X * X * P0 = 0.594604482 X * P1 = 2.54164041 X * Q0 = 2.13725758 X * Q1 = 1.0 X * X * (coefficients from HART table #350 pg 193) X * X * (3) Three applications of Heron's iteration are X * performed using: X * X * y[n+1] = 0.5 * (y[n] + (m/y[n])) X * X * where y[0] = u = approx sqrt(m) X * X * (4) If the value of k was odd then y is either X * multiplied by the square root of two or X * divided by the square root of two for k positive X * or negative respectively. This rescales y X * by multiplying by 2**frac(k/2). X * X * (5) Finally, y is rescaled by int(k/2) which X * is equivalent to multiplication by 2**int(k/2). X * X * The result of steps 4 and 5 is that the value X * of y between 0.5 and 1.0 has been rescaled by X * 2**(k/2) which removes the original rescaling X * done prior to finding the mantissa square root. X * X * NOTES X * X * The Convergent Technologies compiler optimizes division X * by powers of two to "arithmetic shift right" instructions. X * This is ok for positive integers but gives different X * results than other C compilers for negative integers. X * For example, (-1)/2 is -1 on the CT box, 0 on every other X * machine I tried. X * X */ X X#include X#include X#include "pml.h" X X#define P0 0.594604482 /* Approximation coeff */ X#define P1 2.54164041 /* Approximation coeff */ X#define Q0 2.13725758 /* Approximation coeff */ X#define Q1 1.0 /* Approximation coeff */ X X#define ITERATIONS 3 /* Number of iterations */ X Xstatic char funcname[] = "sqrt"; X X Xdouble sqrt (x) Xdouble x; X{ X auto int k; X register int bugfix; X register int kmod2; X register int count; X auto int exponent; X auto double m; X auto double u; X auto double y; X auto double rtnval; X auto struct exception xcpt; X extern double frexp (); X extern double ldexp (); X X DBUG_ENTER ("sqrt"); X DBUG_3 ("sqrtin", "arg %le", x); X if (x == 0.0) { X rtnval = 0.0; X } else if (x < 0.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 { X m = frexp (x, &k); X u = (P0 + (P1 * m)) / (Q0 + (Q1 * m)); X for (count = 0, y = u; count < ITERATIONS; count++) { X y = 0.5 * (y + (m / y)); X } X if ((kmod2 = (k % 2)) < 0) { X y /= SQRT2; X } else if (kmod2 > 0) { X y *= SQRT2; X } X bugfix = 2; X xcpt.retval = ldexp (y, k/bugfix); X } X DBUG_3 ("sqrtout", "result %le", xcpt.retval); X DBUG_RETURN (xcpt.retval); X} END_OF_funcs/src/sqrt.c if test 4591 -ne `wc -c funcs/unused/scale.c <<'END_OF_funcs/unused/scale.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 * scale scale a double precision number by power of 2 X * X * KEY WORDS X * X * scale X * math libraries X * machine dependent routines X * X * SYNOPSIS X * X * double scale (value, scale) X * double value; X * integer scale; X * X * DESCRIPTION X * X * Adds a specified integer to a double precision number's X * exponent, effectively multiplying by a power of 2 for positive X * scale values and divided by a power of 2 for negative X * scale values. X * X * AUTHOR X * X * Fred Fish X * 1346 W 10th Pl X * Tempe, Az. 85281 X * X * INTERNALS X * X * This routine is highly machine dependent. As such, no X * attempt was made to make it general, hence it may have X * to be completely rewritten when transportation of the X * floating point library is attempted. X * X * For the DECSYSTEM-20 the double precision word format is: X * X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM X * X * For the PDP-11 and the 68000 the double precision word X * format is: X * X * WORD N => SEEEEEEEEMMMMMMM X * WORD N+1 => MMMMMMMMMMMMMMMM X * WORD N+2 => MMMMMMMMMMMMMMMM X * WORD N+3 => MMMMMMMMMMMMMMMM X * X * Where: S => Sign bit X * E => Exponent X * X => Ignored (set to 0) X * M => Mantissa bit X * X * Scale extracts the exponent, shifts it into the lower X * order bits, adds the scale value to it, checks for X * exponent overflow or underflow, shifts it back into the X * high bits, and inserts it back into the value. X * X * NOTE: On the DECSYSTEM-20, the exponent for negative X * numbers is complemented. This is not true for the PDP-11. X * X * NOTE: On the 68000, the mantissa is the same for both X * positive and negative numbers, only the sign is different. X * On the PDP-11, the mantissa is two's complement. Both X * use hidden bit normalization. X * X */ X X X#include X#include X#include "pml.h" X X#ifdef PDP10 X#define EXP_MASK 0377000000000 /* Mask for exponent */ X#define MANT_MASK 0400777777777 /* Mask for mantissa */ X#define LEXP_MASK 0377 /* Mask for shifted exponent */ X#define EXP_SHIFTS 27 /* Shifts for exp in LSBs */ X#endif X X#ifdef pdp11 X#define EXP_MASK 0x7F800000 /* Mask for exponent */ X#define MANT_MASK 0x807FFFFF /* Mask for mantissa */ X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */ X#define LEXP_MASK 0377 /* Mask for shifted exponent */ X#endif X X#ifdef mc68000 X#define EXP_MASK 0x7F800000 /* Mask for exponent */ X#define MANT_MASK 0x807FFFFF /* Mask for mantissa */ X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */ X#define LEXP_MASK 0377 /* Mask for shifted exponent */ X#endif X Xstatic union { X double dval; X long lval[2]; X} share; X X X Xdouble scale(value,scale) Xdouble value; Xregister int scale; X{ X register long temp1, temp2, *lpntr; X X lpntr = &share.lval[0]; X share.dval = value; X temp1 = *lpntr; X temp2 = (temp1 & EXP_MASK) >> EXP_SHIFTS; X#ifdef PDP10 X if (value < 0.0) { X temp2 = ~temp2 & LEXP_MASK; X } X#endif X temp2 += scale; X if (temp2 > MAX_EXPONENT+128) { X pmlerr(SCALE_OVERFLOW); X } else if (temp2 < 0) { X pmlerr(SCALE_UNDERFLOW); X } else { X#ifdef PDP10 X if (value < 0.0) { X temp2 = ~temp2 & LEXP_MASK; X } X#endif X temp1 &= MANT_MASK; X temp2 = temp2 << EXP_SHIFTS; X *lpntr = temp1 | temp2; X } X return (share.dval); X} END_OF_funcs/unused/scale.c if test 4092 -ne `wc -c funcs/unused/xexp.c <<'END_OF_funcs/unused/xexp.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 * xexp extract double precision number's exponent X * X * KEY WORDS X * X * xexp X * math libraries X * machine dependent routines X * X * SYNOPSIS X * X * int xexp(value) X * double value; X * X * DESCRIPTION X * X * Extracts exponent from a double precision number and X * returns it as a normal length integer. X * X * AUTHOR X * X * Fred Fish X * 1346 W 10th Pl. X * Tempe, Az 85281 X * X * INTERNALS X * X * This routine is highly machine dependent. As such, no X * attempt was made to make it general, hence it may have X * to be completely rewritten when transportation of the X * floating point library is attempted. X * X * For the DECSYSTEM-20 the double precision word format is: X * X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM X * X * For the PDP-11 and mc68000 (non IEEE format) the double X * precision word format is: X * X * WORD N => SEEEEEEEEMMMMMMM X * WORD N+1 => MMMMMMMMMMMMMMMM X * WORD N+2 => MMMMMMMMMMMMMMMM X * WORD N+3 => MMMMMMMMMMMMMMMM X * X * For the mc68000 using IEEE format the double precision word X * format is: X * X * WORD N => SEEEEEEEEEEEMMMM X * WORD N+1 => MMMMMMMMMMMMMMMM X * WORD N+2 => MMMMMMMMMMMMMMMM X * WORD N+3 => MMMMMMMMMMMMMMMM X * X * Where: S => Sign bit X * E => Exponent X * X => Ignored (set to 0) X * M => Mantissa bit X * X * X */ X X#include X#include X#include "pml.h" X X X#ifdef PDP10 X#define EXP_MASK 0377000000000 /* Mask for exponent */ X#define EXP_SHIFTS 27 /* Shifts to get into LSB's */ X#define EXP_BIAS 128 /* Exponent bias */ X#endif X X#ifdef mc68000 X#ifdef IEEE X#define EXP_MASK 0x7FF00000 /* Mask for exponent */ X#define EXP_SHIFTS 20 /* Shifts to get into LSB's */ X#define EXP_BIAS 1023 /* Exponent bias */ X#else X#define EXP_MASK 0x7F800000 /* Mask for exponent */ X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */ X#define EXP_BIAS 128 /* Exponent bias */ X#endif X#endif X X#ifdef pdp11 X#define EXP_MASK 0x7F800000 /* Mask for exponent */ X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */ X#define EXP_BIAS 128 /* Exponent bias */ X#endif X Xunion dtol { X double dval; X int ival[2]; X}; X X Xint xexp (value) Xunion dtol value; X{ X register int *ipntr; X X if (value.dval == 0.0) { X return (0); X } else { X ipntr = &value.ival[0]; X#ifdef PDP10 X if (value.dval < 0.0) { /* Exponent is complemented */ X *ipntr = ~*ipntr; X } X#endif X *ipntr &= EXP_MASK; X *ipntr >>= EXP_SHIFTS; X *ipntr -= EXP_BIAS; X return (*ipntr); X } X} END_OF_funcs/unused/xexp.c if test 3294 -ne `wc -c funcs/unused/xmant.c <<'END_OF_funcs/unused/xmant.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 * xmant extract double precision number's mantissa X * X * KEY WORDS X * X * xmant X * math libraries X * machine dependent routines X * X * SYNOPSIS X * X * double xmant (value) X * double value; X * X * DESCRIPTION X * X * Extracts a double precision number's mantissa and returns it X * as a double precision number with a normalized mantissa and X * a zero exponent. X * X * AUTHOR X * X * Fred Fish X * 1346 W 10th Pl X * Tempe, Az 85281 X * X * INTERNALS X * X * This routine is highly machine dependent. As such, no X * attempt was made to make it general, hence it may have X * to be completely rewritten when transportation of the X * floating point library is attempted. X * X * For the DECSYSTEM-20 the double precision word format is: X * X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM X * X * For the PDP-11 and the mc68000 the double precision word X * format is: X * X * WORD N => SEEEEEEEEMMMMMMM X * WORD N+1 => MMMMMMMMMMMMMMMM X * WORD N+2 => MMMMMMMMMMMMMMMM X * WORD N+3 => MMMMMMMMMMMMMMMM X * X * For the mc68000 using IEEE format the double precision word X * format is: X * X * WORD N => SEEEEEEEEEEEMMMM X * WORD N+1 => MMMMMMMMMMMMMMMM X * WORD N+2 => MMMMMMMMMMMMMMMM X * WORD N+3 => MMMMMMMMMMMMMMMM X * X * Where: S => Sign bit X * E => Exponent X * X => Ignored (set to 0) X * M => Mantissa bit X * X * NOTE: Beware of 0.0; on some machines which use excess 128 X * notation for the exponent, if the mantissa is zero the exponent X * is also. X * X */ X X#include X#include X#include "pml.h" X X#ifdef PDP10 X#define MANT_MASK 0400777777777 /* Mantissa extraction mask */ X#define ZPOS_MASK 0200000000000 /* Positive # mask for exp = 0 */ X#define ZNEG_MASK 0177000000000 /* Negative # mask for exp = 0 */ X#endif X X#ifdef pdp11 X#define MANT_MASK 0x807FFFFF /* Mantissa extraction mask */ X#define ZPOS_MASK 0x40000000 /* Positive # mask for exp = 0 */ X#define ZNEG_MASK 0x40000000 /* Negative # mask for exp = 0 */ X#endif X X#ifdef mc68000 X#ifdef IEEE X#define MANT_MASK 0x800FFFFF /* Mantissa extraction mask */ X#define ZPOS_MASK 0x3FF00000 /* Positive # mask for exp = 0 */ X#define ZNEG_MASK 0x3FF00000 /* Negative # mask for exp = 0 */ X#else X#define MANT_MASK 0x807FFFFF /* Mantissa extraction mask */ X#define ZPOS_MASK 0x40000000 /* Positive # mask for exp = 0 */ X#define ZNEG_MASK 0x40000000 /* Negative # mask for exp = 0 */ X#endif X#endif X Xunion dtol { X double dval; X int ival[2]; X}; X X Xdouble xmant (value) Xunion dtol value; X{ X register int *ipntr; X X ipntr = &value.ival[0]; X *ipntr &= MANT_MASK; X *ipntr |= ZPOS_MASK; X return (value.dval); X} X X X/****************** OLD WAY *********** Xdouble xmant (value) Xunion dtol value; X{ X register long *ipntr; X register int neg; X X if (value == 0.0) { X return (value); X } else { X if (value < 0.0) { X neg = TRUE; X } else { X neg = FALSE; X } X ipntr = &value.ival[0]; X *ipntr &= MANT_MASK; X if (neg) { X *ipntr |= ZNEG_MASK; X } else { X *ipntr |= ZPOS_MASK; X } X return (value.dval); X } X} X**********************/ END_OF_funcs/unused/xmant.c if test 3950 -ne `wc -c tests/c2c.dat <<'END_OF_tests/c2c.dat' Xcsin 1.0e00 2.0e00 3.165778547525405883e00 1.959601044654846191e00 Xcsin 0.0e00 2.0e00 0.0e00 3.626860409975051879e00 Xcsin 1.0e00 0.0e00 8.414709866046905517e-01 0.0e00 Xcsin -2.0e00 2.0e00 -3.420954853296279907e00 -1.509306490421295166e00 Xcsin -1.0e00 -2.0e00 -3.165778547525405883e00 -1.959601044654846191e00 Xcsin 1.0e00 -2.0e00 3.165778547525405883e00 -1.959601044654846191e00 Xccos 1.0e00 2.0e00 2.032723009586334228e00 -3.051897794008255004e00 Xccos 0.0e00 2.0e00 3.762195706367492675e00 0.0e00 Xccos 1.0e00 0.0e00 5.403023064136505127e-01 0.0e00 Xccos -2.0e00 2.0e00 -1.565625846385955810e00 3.297894805669784545e00 Xccos -1.0e00 -2.0e00 2.032723009586334228e00 -3.051897794008255004e00 Xccos 1.0e00 -2.0e00 2.032723009586334228e00 3.051897794008255004e00 Xcexp 1.0e00 2.0e00 -1.131204381585121154e00 2.471726655960083007e00 Xcexp 0.0e00 2.0e00 -4.161468371748924255e-01 9.092974215745925903e-01 Xcexp 1.0e00 0.0e00 2.718281835317611694e00 0.0e00 Xcexp -2.0e00 2.0e00 -5.631934991106390953e-02 1.230600243434309959e-01 Xcexp -1.0e00 -2.0e00 -1.530918665230274200e-01 -3.345118276774883270e-01 Xcexp 1.0e00 -2.0e00 -1.131204381585121154e00 -2.471726655960083007e00 Xctanh 1.0e00 2.0e00 1.166736245155334472e00 -2.434581927955150604e-01 Xctanh 0.0e00 2.0e00 0.0e00 -2.185039848089218139e00 Xctanh 1.0e00 0.0e00 7.615941539406776428e-01 0.0e00 Xctanh -2.0e00 2.0e00 -1.023835599422454834e00 -2.839295566082000732e-02 Xctanh -1.0e00 -2.0e00 -1.166736245155334472e00 2.434581927955150604e-01 Xctanh 1.0e00 -2.0e00 1.166736245155334472e00 2.434581927955150604e-01 Xctan 1.0e00 2.0e00 3.381283208727836608e-02 1.014793619513511657e00 Xctan 0.0e00 2.0e00 0.0e00 9.640275761485099792e-01 Xctan 1.0e00 0.0e00 1.557407721877098083e00 0.0e00 Xctan -2.0e00 2.0e00 2.839295566082000732e-02 1.023835599422454834e00 Xctan -1.0e00 -2.0e00 -3.381283208727836608e-02 -1.014793619513511657e00 Xctan 1.0e00 -2.0e00 3.381283208727836608e-02 -1.014793619513511657e00 Xcsqrt 0.0 0.0 0.000000000000000000e00 0.000000000000000000e00 Xcsqrt 1.0 2.0 1.272019654512405395e00 7.861513718962669372e-01 Xcsqrt 0.0 2.0 1.000000000000000000e00 1.000000000000000000e00 Xcsqrt 1.0 0.0 1.000000000000000000e00 0.000000000000000000e00 Xcsqrt -2.0 2.0 6.435942575335502624e-01 1.553773969411849975e00 Xcsqrt -1.0 2.0 7.861513718962669372e-01 1.272019654512405395e00 Xcsqrt 1.0 2.0 1.272019654512405395e00 7.861513718962669372e-01 Xcsinh 1.0e00 2.0e00 -4.890562593936920166e-01 1.403119236230850219e00 Xcsinh 0.0e00 2.0e00 0.0e00 9.092974215745925903e-01 Xcsinh 1.0e00 0.0e00 1.175201192498207092e00 0.0e00 Xcsinh -2.0e00 2.0e00 1.509306475520133972e00 3.420954853296279907e00 Xcsinh -1.0e00 -2.0e00 4.890562593936920166e-01 -1.403119236230850219e00 Xcsinh 1.0e00 -2.0e00 -4.890562593936920166e-01 -1.403119236230850219e00 Xcrcp 1.0e00 2.0e00 1.999999992549419403e-01 -3.999999985098838806e-01 Xcrcp 0.0e00 2.0e00 0.0e00 -5.0e-01 Xcrcp 1.0e00 0.0e00 1.0e00 0.0e00 Xcrcp -2.0e00 2.0e00 -2.5e-01 -2.5e-01 Xcrcp -1.0e00 -2.0e00 -1.999999992549419403e-01 3.999999985098838806e-01 Xcrcp 1.0e00 -2.0e00 1.999999992549419403e-01 3.999999985098838806e-01 Xclog 1.0e00 2.0e00 8.047189563512802124e-01 1.107148721814155578e00 Xclog 0.0e00 2.0e00 6.931471824645996093e-01 1.570796325802803039e00 Xclog 1.0e00 0.0e00 0.0e00 0.0e00 Xclog -2.0e00 2.0e00 1.039720773696899414e00 2.356194496154785156e00 Xclog -1.0e00 -2.0e00 8.047189563512802124e-01 -2.034443944692611694e00 Xclog 1.0e00 -2.0e00 8.047189563512802124e-01 -1.107148721814155578e00 Xccosh 1.0e00 2.0e0 -6.421481221914291381e-01 1.068607419729232788e00 Xccosh 0.0e00 2.0e00 -4.161468371748924255e-01 0.0e00 Xccosh 1.0e00 0.0e00 1.543080642819404602e00 0.0e00 Xccosh -2.0e00 2.0e00 -1.565625831484794616e00 -3.297894805669784545e00 Xccosh -1.0e00 -2.0e00 -6.421481221914291381e-01 1.068607419729232788e00 Xccosh 1.0e00 -2.0e00 -6.421481221914291381e-01 -1.068607419729232788e00 Xcatan 1.0e00 2.0e00 1.338972523808479309e00 4.023594781756401062e-01 Xcatan 0.0e00 2.0e00 1.570796325802803039e00 5.493061468005180358e-01 /* -real part */ Xcatan 1.0e00 0.0e00 7.853981629014015197e-01 0.0e00 Xcatan -2.0e00 2.0e00 -1.311223268508911132e00 2.388778626918792724e-01 Xcatan -1.0e00 -2.0e00 -1.338972523808479309e00 -4.023594819009304046e-01 Xcatan 1.0e00 -2.0e00 1.338972523808479309e00 -4.023594819009304046e-01 Xcasin 1.0 2.0 4.270785786211490631e-01 1.528570890426635742e00 Xcasin 0.0 2.0 0.000000000000000000e00 1.443635478615760803e00 Xcasin 1.0 0.0 1.570796325802803039e00 0.000000000000000000e00 Xcasin -2.0 2.0 -7.542491331696510314e-01 1.734324455261230468e00 Xcasin -1.0 -2.0 -4.270785823464393615e-01 -1.528570920228958129e00 Xcasin 1.0 -2.0 4.270785823464393615e-01 -1.528570920228958129e00 Xcacos 1.0 2.0 1.143717750906944274e00 -1.528570920228958129e00 Xcacos 0.0 2.0 1.570796325802803039e00 -1.443635478615760803e00 Xcacos 1.0 0.0 0.000000000000000000e00 0.000000000000000000e00 Xcacos -2.0 2.0 2.325045466423034668e00 -1.734324529767036438e00 Xcacos -1.0 -2.0 1.997874900698661804e00 1.528570890426635742e00 Xcacos 1.0 -2.0 1.143717750906944274e00 1.528570890426635742e00 END_OF_tests/c2c.dat if test 5354 -ne `wc -c tests/cc2c.dat <<'END_OF_tests/cc2c.dat' Xcadd 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 4.683000028133392334e00 6.440000057220458984e00 Xcadd -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 2.003005027770996093e00 -1.887652993202209472e00 Xcadd 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 3.719695717096328735e00 4.449983000755310058e00 Xcadd 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 4.954332172870635986e00 0.000000000000000000e00 Xcadd -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 -4.234000027179718017e00 -6.334455549716949462e00 Xcadd 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 1.227000000000000000e04 1.300034515380859375e03 Xcdiv 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 4.609979800879955291e-01 1.263787318021059036e-01 Xcdiv -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 -4.863302707672119140e-01 6.862613558769226074e-02 Xcdiv 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 5.729839205741882324e-02 -4.476501746103167533e-02 Xcdiv 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 4.663971476256847381e-01 0.000000000000000000e00 Xcdiv -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 5.215928852558135986e-01 8.269466646015644073e-02 Xcdiv 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 3.765932763671875000e04 -3.554461171875000000e05 Xcmult 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 -5.596120119094848632e00 1.293470001220703125e01 Xcmult -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 4.705165147781372070e00 1.302998638153076171e01 Xcmult 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 8.510770574212074279e-01 2.011436134576797485e00 Xcmult 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 5.323814511299133300e00 0.000000000000000000e00 Xcmult -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 -5.635822236537933349e00 1.193936669826507568e01 Xcmult 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 -4.487600040435791015e01 4.235604019165039062e02 Xcsubt 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 -2.437000006437301635e00 -1.760000020265579223e00 Xcsubt -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 -4.694344997406005859e00 6.356786966323852539e00 Xcsubt 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 -2.948183119297027587e00 -4.249976992607116699e00 Xcsubt 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 -1.802817046642303466e00 0.000000000000000000e00 Xcsubt -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 1.766000002622604370e00 1.665544450283050537e00 Xcsubt 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 1.227000000000000000e04 1.299965484619140625e03 END_OF_tests/cc2c.dat if test 3825 -ne `wc -c