From ae89720834d70256a18af17ca64ff49892431c63 Mon Sep 17 00:00:00 2001 From: Keith Marshall Date: Mon, 21 Nov 2016 19:24:19 +0000 Subject: [PATCH] Replace defective powf() and powl() function implementations. --- mingwrt/ChangeLog | 14 + mingwrt/Makefile.in | 2 +- mingwrt/include/math.h | 4 - mingwrt/mingwex/math/pow_generic.sx | 423 +++++++++++++++++++ mingwrt/mingwex/math/powf.c | 3 - mingwrt/mingwex/math/powl.c | 803 ------------------------------------ 6 files changed, 438 insertions(+), 811 deletions(-) create mode 100644 mingwrt/mingwex/math/pow_generic.sx delete mode 100644 mingwrt/mingwex/math/powf.c delete mode 100644 mingwrt/mingwex/math/powl.c diff --git a/mingwrt/ChangeLog b/mingwrt/ChangeLog index b5ea2ae..3a71f3a 100644 --- a/mingwrt/ChangeLog +++ b/mingwrt/ChangeLog @@ -1,3 +1,17 @@ +2016-11-21 Keith Marshall + + Replace defective powf() and powl() function implementations. + + * mingwex/math/powf.c mingwex/math/powl.c: Delete; replaced by... + * mingwex/math/pow_generic.sx: ...this new file; it implements... + (__x87pow): ...this generic power function back-end, serving... + (pow, powf, powl): ...each of these front-end entry points. + + * Makefile.in (libmingwex.a): Add x87pow.$OBJEXT + + * include/math.h (powf): Remove inline implementation; it no longer + offers any inline advantage. + 2016-11-18 Keith Marshall Implement unit tests for power functions. diff --git a/mingwrt/Makefile.in b/mingwrt/Makefile.in index bf5a702..c95d788 100644 --- a/mingwrt/Makefile.in +++ b/mingwrt/Makefile.in @@ -462,7 +462,7 @@ libmingwex.a: $(addsuffix .$(OBJEXT), cosf cosl acosf acosl sinf sinl asinf \ powf powl powi powif powil remainder remainderf remainderl remquo remquof \ remquol rint rintf rintl round roundf roundl scalbn scalbnf scalbnl signbit \ signbitf signbitl sqrtf sqrtl tgamma tgammaf tgammal trunc truncf truncl \ - x87cvt x87cvtf x87log x87log1p) + x87cvt x87cvtf x87log x87log1p x87pow) # Replacement I/O functions in libmingwex.a, providing better POSIX # compatibility than their Microsoft equivalents. diff --git a/mingwrt/include/math.h b/mingwrt/include/math.h index af94513..2ac331c 100644 --- a/mingwrt/include/math.h +++ b/mingwrt/include/math.h @@ -650,10 +650,6 @@ extern long double __cdecl hypotl (long double, long double); /* 7.12.7.4 The pow functions. Double in C89 */ extern float __cdecl powf (float, float); -#ifndef __NO_INLINE__ -__CRT_INLINE float __cdecl powf (float x, float y) - {return (float) pow (x, y);} -#endif extern long double __cdecl powl (long double, long double); /* 7.12.7.5 The sqrt functions. Double in C89. */ diff --git a/mingwrt/mingwex/math/pow_generic.sx b/mingwrt/mingwex/math/pow_generic.sx new file mode 100644 index 0000000..65436b4 --- /dev/null +++ b/mingwrt/mingwex/math/pow_generic.sx @@ -0,0 +1,423 @@ +/* + * pow_generic.sx + * + * Generic implementation for the pow(), powl(), and powf() functions. + * + * $Id$ + * + * Written by Keith Marshall + * Copyright (C) 2016, MinGW.org Project + * + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice (including the next + * paragraph) shall be included in all copies or substantial portions of the + * Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + * + */ +#undef __function + +#if defined _powl_source +# define __function _powl +# define __yoffset __xoffset+12 +# define __fldx fldt __xoffset(%esp) +# define __fldy fldt __yoffset(%esp) + +#elif defined _powf_source +# define __function _powf +# define __yoffset __xoffset+4 +# define __fldx flds __xoffset(%esp) +# define __fldy flds __yoffset(%esp) +# define ___x87cvt ___x87cvtf + +#elif defined _pow_source +# define __function _pow +# define __yoffset __xoffset+8 +# define __fldx fldl __xoffset(%esp) +# define __fldy fldl __yoffset(%esp) +#endif + +#ifdef __MINGW64__ +# define __xoffset 8 +# define esp rsp +# define eax rax +# define edx rdx +#else +# define __xoffset 4 +#endif + +.text +.align 4 +#ifdef __function +/* A specific front-end entry point name has been identified; thus, + * we are assembling the front-end stub implementation for just one + * of the three supported functions, with C language prototypes: + * + * double pow (double x, double y); + * long double powl (long double x, long double y); + * float powf (float x, float y); + */ +.globl __function +.def __function; .scl 2; .type 32; .endef + +__function: +/* First, load x and y into the FPU, using the appropriate operand + * size specification for the specified front-end entry point, then + * hand off control to the generic back-end function. + */ + __fldx /* x */ + __fldy /* y ; x */ + +#ifdef _powl_source +/* For the long double powl (long double, long double) form of the + * primary function, we can simply delegate computation of the REAL10 + * result to the ___x87pow() handler, with that returning directly + * to powl()'s own caller... + */ + jmp ___x87pow + +#else +/* ...whereas for each of double pow (double, double) form, and its + * float powf (float, float) sibling, we must call the backend handler + * to compute an intermediate REAL10 result... + */ + call ___x87pow + +/* ...then return that via the appropriate type conversion/validation + * handler, to obtain the ultimately required REAL8 or REAL4 result. + */ + jmp ___x87cvt + +#endif +#else +/* No specific function entry point identified; implement the generic + * back-end, which is common to all supported front-end entry points; + * it also provides the error reporting API. + */ +#include "errno.sx" + +.align 4 +.globl ___x87pow +.def ___x87pow; .scl 2; .type 32; .endef + +___x87pow: + fxam /* classify y input value */ + fnstsw %ax /* copy FPU flags to CPU flags */ + fld1 /* +1.0 ; y ; x */ + sahf /* examine ZF = C3 and PF = C2 */ + jnz 30f /* y is non-zero */ + jp 30f /* y is non-zero, denormalized */ + +/* In the case where y is zero, then POSIX says that the value of x is + * irrelevant, (even if it is indefinite); the return value is +1.0 + */ + fstp %st(2) /* y ; x^y = 1.0 */ + fstp %st(0) /* x^y */ + ret + +/* When y is non-zero, proceed to consideration of the x argument value; + * (this is necessary, even if the value of y is indeterminate). + */ +30: movb %ah, %dl /* save y classification flags */ + fucomp %st(2) /* y ; x */ + fnstsw %ax /* copy FPU flags to CPU flags */ + fxch /* x ; y */ + sahf /* examine ZF = C3 and PF = C2 */ + jp 32f /* return x^y = indeterminate x */ + jnz 40f /* x != +1.0 */ + +/* For this specific case, where x == +1.0, POSIX says that the return + * value shall be +1.0, (even if the value of y is indeterminate). + */ + fld1 /* x^y = 1.0 ; y ; x */ +31: fstp %st(1) /* x^y ; x */ +32: fstp %st(1) /* x^y */ + ret + +/* In any other case, if either x or y is NaN, then POSIX requires that + * NaN shall be returned; first check for x being NaN, or infinite. + */ +40: fxam /* classify x input value */ + fnstsw %ax /* copy FPU flags to CPU flags */ + sahf /* test for infinity or NaN */ + jnc 50f /* x is finite, so pass it on */ + jnp 32b /* return x^y = x as NaN */ + +/* We've identified that x is infinite; how we handle this boundary + * condition depends on whether it's a +ve infinity, or a -ve. + */ + testb $0x02, %ah /* x ; y */ + jnz 42f /* x is -ve infinity */ + +/* In the case where x is +ve infinity, POSIX stipulates that the return + * value should be +ve infinity when y > 0.0, or +0.0 when y < 0.0; first + * deal with the y > 0.0 case. + */ + testb $0x02, %dl /* check if y is +ve, or -ve? */ + jz 32b /* +ve: return x^y = x = +ve infinity */ + +/* Alternatively, in the case when x is +ve infinity, and y < 0.0, we + * substitute 0.0 for x, then return it as the value for x^y. + */ +41: fldz /* 0.0 ; x ; y */ + jmp 31b /* return x^y = 0.0 */ + +/* Similarly, in the case where x is -ve infinity, we must again return + * infinity for y > 0.0, or 0.0 for y < 0.0; however, in this case, the + * sign of the returned value must be -ve if y is an odd valued integer, + * or +ve for any other value of y. + */ +42: testb $0x02, %dl /* check if y is +ve, or -ve? */ + jz 43f /* when +ve, return signed infinity */ + +/* Fall through when x is -ve infinity and y < 0.0; substitute -0.0 for + * the infinite value of x, then adjust the sign depending on whether y + * is an odd valued integer, or any other value. + */ + fldz /* 0.0 ; x ; y */ + fstp %st(1) /* 0.0 ; y */ + fchs /* -0.0 ; y */ + +/* Determine if y is non-integral, or an even valued integer, in either + * of which cases we force a +ve return value, or an odd valued integer, + * in whiich case we leave the sign of the return value as it is; begin + * by checking if y/2 is an integer, which asserts that y itself is an + * even valued integer. + */ +43: fld1 /* 1.0 ; x ; y */ + fchs /* -1.0 ; x ; y */ + fld %st(2) /* y ; -1.0 ; x ; y */ + fscale /* y/2 ; -1.0 ; x ; y */ + fst %st(1) /* y/2 ; y/2 ; x ; y */ + frndint /* int(y/2) ; y/2 ; x ; y */ + fucompp /* x ; y */ + fnstsw %ax /* check if int(y/2) == y/2 ? */ + sahf /* hence y is an even valued integer */ + je 44f /* so go force +ve x^y return value */ + +/* When we've established that y is not an even valued integer, we must + * still confirm the possibility that it is an odd valued integer; i.e. + * if it is an integer, it must be odd valued. + */ + fld %st(1) /* y ; x ; y */ + fld %st(0) /* y ; y ; x ; y */ + frndint /* int(y) ; y ; x ; y */ + fucompp /* x ; y */ + fnstsw %ax /* check if int(y) == y ? */ + sahf /* hence y is an odd valued integer */ + je 32b /* so return x^y as is */ + +/* When y is either an even valued integer, or not an integer at all: + */ +44: fabs /* make x^y value +ve */ + jmp 32b /* then return it */ + +/* When x is finite, we still need to check the possibility that y may + * be NaN, or may be infinite. + */ +50: xchgb %dl, %ah /* reload y classification flags */ + movb %ah, %dh /* save a copy */ + sahf /* check for y finite, infinite, or NaN */ + jnc 60f /* y is also finite */ + jp 52f /* y is infinite, but not NaN */ + +/* y is NaN; pop x off FPU stack, and return x^y as NaN value of y. + */ +51: fstp %st(0) /* y */ + ret /* return x^y = y */ + +/* We've identified x as finite, but y as infinite; POSIX defines + * boundary conditions about the range -1.0 < x < +1.0, which may be + * differentiated by comparison between +1.0 and |x|. We've already + * that x != +1.0, so if we now identify that |x| == +1.0, then this + * must represent x == -1.0, a boundary condition for which POSIX + * prescribes a return value of +1.0 + */ +52: fabs /* |x| ; y */ + fld1 /* 1.0 ; |x| ; y */ + fucom %st(1) /* check for |x| == 1.0 */ + fnstsw %ax /* copy FPU flags to CPU flags */ + sahf /* if ZF == 1 -> |x| == 1.0 */ + je 31b /* return x^y = |x| = +1.0 */ + +/* When |x| != 1.0, we have no further use for the comparative values + * of 1.0 and |x|, on the FPU stack; discard them, then check the flag + * state to establish whether x lies within the boundary range. + */ + fstp %st(0) /* |x| ; y */ + fstp %st(0) /* y */ + jnb 53f /* -1.0 < x < +1.0 */ + +/* This represents the POSIX boundary condition where y is infinite, + * and |x| > 1.0; for this condition, POSIX specifies a return value + * of x^y = 0.0 if y is -ve infinity, otherwise x^y = y. + */ + test $0x02, %dh /* if y is -ve infinity */ + jnz 41b /* then go return +0.0 */ + ret /* else return -ve infinity */ + +/* Here, we have -1.0 < x < +1.0, and y is infinite; for this case, + * POSIX prescribes a return value of +0.0 when y is +ve infinity, or + * +ve infinity when y itself -ve infinity. + */ +53: test $0x02, %dh /* if y is +ve infinity */ + jz 41b /* then go return +0.0 */ + fabs /* else force to +ve infinity */ + ret /* and return it */ + +/* We've now established that both x and y are finite, but we must + * still consider the special restrictions which apply when x == 0.0 + * or x < 0.0 + */ +60: movb %dl, %ah /* review x value classification */ + sahf /* examining ZF = C3 and PF = C2 */ + jnz 70f /* x is non-zero */ + jp 70f /* x is non-zero, denormalized */ + +/* When x is zero, the return value is (possibly signed) zero for + * all y > 0.0, but infinite, and reported as a pole error, for any + * y < 0.0 + */ + testb $0x02, %dh /* if y is -ve? */ + jnz 61f /* then go process the pole error */ + +/* For the case where y > 0.0, the sign of the original x value is + * preserved when y is an odd valued integer, or forced to +0.0 for + * any other +ve value of y; (obviously, if x is already +0.0, the + * sign preservation condition becomes irrelevant). + */ + testb $0x02, %ah /* if x == +0.0 */ + jz 32b /* then just go return it */ + jmp 43b /* else go ajust sign */ + +/* For the pole error case, we must substitute an infinity for the + * original value of x; a convenient way to achieve this is to take + * the logarithm of x, (which is -ve infinity by definition). + */ +61: fld1 /* 1.0 ; 0.0 ; y */ + fxch %st(1) /* 0.0 ; 1.0 ; y */ + fyl2x /* 1.0 * log2(0.0) = -inf ; y */ + +/* To diagnose the pole error, we will set errno = ERANGE, (which + * is compliant with POSIX); on Win32, we call __errno() to get a + * pointer to errno itself, but note that we haven't done with EDX + * yet, so we must guard against possible modification during the + * execution of __errno(). + */ + pushl %edx /* we must save this */ + errno ERANGE /* because this may change it */ + popl %edx /* restore saved value */ + +/* The returned infinity must preserve the sign of the original x, + * when y is an odd valued integer, otherwise it is forced to +inf; + * (obviously, if x is +0.0, we may just force +inf anyway). + */ + testb $0x02, %dl /* if x is -0.0 */ + jnz 43b /* then go do signed return */ + fabs /* else force to +inf */ + jmp 32b /* and return it */ + +/* When both x and y are finite and non-zero, then we must check + * for a possible domain error condition, which occurs when x < 0 + * and y has any value which is not an integer. + */ +70: testb $0x02, %ah /* if x > 0.0 */ + jz 80f /* then result is computable */ + +/* Here, x < 0.0; the result may still be computable, if (and only + * if) the value of y is an integer. + */ + fld %st(1) /* y ; x ; y */ + fld %st(0) /* y ; y ; x ; y */ + frndint /* int(y) ; y ; x ; y */ + fucompp /* x ; y */ + fnstsw %ax /* copy FPU flags to CPU flags */ + sahf /* to test if y == int(y) ? */ + je 71f /* then result is computable */ + +/* Fall through when x < 0.0 and y is not an integer; in this case + * we must set errno to report a domain error, and return NaN. + */ + fsqrt /* NaN ; y */ + errno EDOM /* set errno = EDOM */ + jmp 32b /* return NaN */ + +/* When x < 0.0 and y is an integer, we may still compute x^y + * according to the relationship x^y = -1^y * 2^(y * log2(|x|)) + */ +71: fabs /* |x| ; y */ + fld %st(1) /* y ; |x| ; y */ + fxch %st(1) /* |x| ; y ; y */ + call 80f /* |x|^y ; y */ + fchs /* assume y is odd valued */ + jmp 43b /* adjust if even valued */ + +/* When x > 0.0, and y is finite, we may proceed to compute x^y, + * according to the relationship: x^y = 2^(y * log2(x)); first we + * compute log2(x), preferring the FYL2XP1 method for values of x + * close to zero, but falling back on FYL2X for x > 1.29 + */ +80: call ___x87log /* y*log2(x) */ + +/* Having computed the value of y * log2(x), we may now compute + * the final result as 2^(y * log2(x)). We must compute this in + * stages, combining 2^frac(y * log2(x)) * 2^int(y * log2(x)) to + * yielding the final result for x^y; first separate y * log2(x) + * into fractional and integer parts: + */ + fld %st /* y*log2(x) ; y*log2(x) */ + frndint /* int(y*log2(x)) ; y*log2(x) */ + fxch %st(1) /* y*log2(x) ; int(y*log2(x)) */ + fsub %st(1), %st /* frac(y*log2(x)) ; int(y*log2(x)) */ + +/* Now compute the intermediate 2^frac(y * log(x)) - 1.0 result: + */ + f2xm1 /* 2^frac(y*log2(x))-1 ; int(y*log2(x)) */ + +/* Add the 1.0 deficit, to yield the 2^frac(y * log2(x)) result: + */ + fld1 /* 1 ; 2^frac(y*log2(x))-1 ; int(y*log2(x)) */ + faddp /* 2^frac(y*log2(x)) ; int(y*log2(x)) */ + +/* Finally, multiply by 2^int(y * log2(x)), to yield the x^y result: + */ + fscale /* x^y ; int(y*log2(x)) */ + fstp %st(1) /* x^y */ + +/* At this point, the value of x^y should not be zero; if it is, then + * the computation has underflowed, in which case POSIX recommends that + * errno should be set to ERANGE. Alternatively, if the result becomes + * infinite then the computation has overflowed, in which case POSIX + * requires that errno be so set. Check if either is appropriate. + */ + fxam /* classify x^y result */ + fnstsw %ax /* copy FPU status flags */ + sahf /* to test via CPU flags, hence report */ + jbe 81f /* ZF -> underflow; CF -> overflow */ + ret /* else return x^y, errno unchanged */ + +/* Here, we provide an alternative function return, for use when either + * overflow or underflow is detected during the computation of x^y; it + * returns whatever x^y value has been computed, after having set errno + * to indicate the ERANGE condition. + */ +81: errno ERANGE /* set errno = ERANGE */ + ret /* return x^y, errno = ERANGE */ +#endif + +/* vim: set autoindent filetype=asm formatoptions=croql: */ +/* $RCSfile$: end of file */ diff --git a/mingwrt/mingwex/math/powf.c b/mingwrt/mingwex/math/powf.c deleted file mode 100644 index 1af4d2d..0000000 --- a/mingwrt/mingwex/math/powf.c +++ /dev/null @@ -1,3 +0,0 @@ -#include -float powf (float x, float y) - {return (float) pow (x, y);} diff --git a/mingwrt/mingwex/math/powl.c b/mingwrt/mingwex/math/powl.c deleted file mode 100644 index 03d7122..0000000 --- a/mingwrt/mingwex/math/powl.c +++ /dev/null @@ -1,803 +0,0 @@ -/* - * powl.c - * - * Power function, long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, z, powl(); - * - * z = powl( x, y ); - * - * - * - * DESCRIPTION: - * - * Computes x raised to the yth power. Analytically, - * - * x**y = exp( y log(x) ). - * - * Following Cody and Waite, this program uses a lookup table - * of 2**-i/32 and pseudo extended precision arithmetic to - * obtain several extra bits of accuracy in both the logarithm - * and the exponential. - * - * - * - * ACCURACY: - * - * The relative error of pow(x,y) can be estimated - * by y dl ln(2), where dl is the absolute error of - * the internally computed base 2 logarithm. At the ends - * of the approximation interval the logarithm equal 1/32 - * and its relative error is about 1 lsb = 1.1e-19. Hence - * the predicted relative error in the result is 2.3e-21 y . - * - * Relative error: - * arithmetic domain # trials peak rms - * - * IEEE +-1000 40000 2.8e-18 3.7e-19 - * .001 < x < 1000, with log(x) uniformly distributed. - * -1000 < y < 1000, y uniformly distributed. - * - * IEEE 0,8700 60000 6.5e-18 1.0e-18 - * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. - * - * - * ERROR MESSAGES: - * - * message condition value returned - * pow overflow x**y > MAXNUM INFINITY - * pow underflow x**y < 1/MAXNUM 0.0 - * pow domain x<0 and y noninteger 0.0 - * - */ - -/* -Cephes Math Library Release 2.7: May, 1998 -Copyright 1984, 1991, 1998 by Stephen L. Moshier -*/ - -/* -Modified for mingw -2002-07-22 Danny Smith -*/ - -#ifdef __MINGW32__ -#include "cephes_mconf.h" -#else -#include "mconf.h" - -static char fname[] = {"powl"}; -#endif - -#ifndef _SET_ERRNO -#define _SET_ERRNO(x) -#endif - - -/* Table size */ -#define NXT 32 -/* log2(Table size) */ -#define LNXT 5 - -#ifdef UNK -/* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z) - * on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1 - */ -static long double P[] = { - 8.3319510773868690346226E-4L, - 4.9000050881978028599627E-1L, - 1.7500123722550302671919E0L, - 1.4000100839971580279335E0L, -}; -static long double Q[] = { -/* 1.0000000000000000000000E0L,*/ - 5.2500282295834889175431E0L, - 8.4000598057587009834666E0L, - 4.2000302519914740834728E0L, -}; -/* A[i] = 2^(-i/32), rounded to IEEE long double precision. - * If i is even, A[i] + B[i/2] gives additional accuracy. - */ -static long double A[33] = { - 1.0000000000000000000000E0L, - 9.7857206208770013448287E-1L, - 9.5760328069857364691013E-1L, - 9.3708381705514995065011E-1L, - 9.1700404320467123175367E-1L, - 8.9735453750155359320742E-1L, - 8.7812608018664974155474E-1L, - 8.5930964906123895780165E-1L, - 8.4089641525371454301892E-1L, - 8.2287773907698242225554E-1L, - 8.0524516597462715409607E-1L, - 7.8799042255394324325455E-1L, - 7.7110541270397041179298E-1L, - 7.5458221379671136985669E-1L, - 7.3841307296974965571198E-1L, - 7.2259040348852331001267E-1L, - 7.0710678118654752438189E-1L, - 6.9195494098191597746178E-1L, - 6.7712777346844636413344E-1L, - 6.6261832157987064729696E-1L, - 6.4841977732550483296079E-1L, - 6.3452547859586661129850E-1L, - 6.2092890603674202431705E-1L, - 6.0762367999023443907803E-1L, - 5.9460355750136053334378E-1L, - 5.8186242938878875689693E-1L, - 5.6939431737834582684856E-1L, - 5.5719337129794626814472E-1L, - 5.4525386633262882960438E-1L, - 5.3357020033841180906486E-1L, - 5.2213689121370692017331E-1L, - 5.1094857432705833910408E-1L, - 5.0000000000000000000000E-1L, -}; -static long double B[17] = { - 0.0000000000000000000000E0L, - 2.6176170809902549338711E-20L, --1.0126791927256478897086E-20L, - 1.3438228172316276937655E-21L, - 1.2207982955417546912101E-20L, --6.3084814358060867200133E-21L, - 1.3164426894366316434230E-20L, --1.8527916071632873716786E-20L, - 1.8950325588932570796551E-20L, - 1.5564775779538780478155E-20L, - 6.0859793637556860974380E-21L, --2.0208749253662532228949E-20L, - 1.4966292219224761844552E-20L, - 3.3540909728056476875639E-21L, --8.6987564101742849540743E-22L, --1.2327176863327626135542E-20L, - 0.0000000000000000000000E0L, -}; - -/* 2^x = 1 + x P(x), - * on the interval -1/32 <= x <= 0 - */ -static long double R[] = { - 1.5089970579127659901157E-5L, - 1.5402715328927013076125E-4L, - 1.3333556028915671091390E-3L, - 9.6181291046036762031786E-3L, - 5.5504108664798463044015E-2L, - 2.4022650695910062854352E-1L, - 6.9314718055994530931447E-1L, -}; - -#define douba(k) A[k] -#define doubb(k) B[k] -#define MEXP (NXT*16384.0L) -/* The following if denormal numbers are supported, else -MEXP: */ -#ifdef DENORMAL -#define MNEXP (-NXT*(16384.0L+64.0L)) -#else -#define MNEXP (-NXT*16384.0L) -#endif -/* log2(e) - 1 */ -#define LOG2EA 0.44269504088896340735992L -#endif - - -#ifdef IBMPC -static const uLD P[] = { -{ { 0xb804,0xa8b7,0xc6f4,0xda6a,0x3ff4, XPD } }, -{ { 0x7de9,0xcf02,0x58c0,0xfae1,0x3ffd, XPD } }, -{ { 0x405a,0x3722,0x67c9,0xe000,0x3fff, XPD } }, -{ { 0xcd99,0x6b43,0x87ca,0xb333,0x3fff, XPD } } -}; -static const uLD Q[] = { -{ { 0x6307,0xa469,0x3b33,0xa800,0x4001, XPD } }, -{ { 0xfec2,0x62d7,0xa51c,0x8666,0x4002, XPD } }, -{ { 0xda32,0xd072,0xa5d7,0x8666,0x4001, XPD } } -}; -static const uLD A[] = { -{ { 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD } }, -{ { 0x033a,0x722a,0xb2db,0xfa83,0x3ffe, XPD } }, -{ { 0xcc2c,0x2486,0x7d15,0xf525,0x3ffe, XPD } }, -{ { 0xf5cb,0xdcda,0xb99b,0xefe4,0x3ffe, XPD } }, -{ { 0x392f,0xdd24,0xc6e7,0xeac0,0x3ffe, XPD } }, -{ { 0x48a8,0x7c83,0x06e7,0xe5b9,0x3ffe, XPD } }, -{ { 0xe111,0x2a94,0xdeec,0xe0cc,0x3ffe, XPD } }, -{ { 0x3755,0xdaf2,0xb797,0xdbfb,0x3ffe, XPD } }, -{ { 0x6af4,0xd69d,0xfcca,0xd744,0x3ffe, XPD } }, -{ { 0xe45a,0xf12a,0x1d91,0xd2a8,0x3ffe, XPD } }, -{ { 0x80e4,0x1f84,0x8c15,0xce24,0x3ffe, XPD } }, -{ { 0x27a3,0x6e2f,0xbd86,0xc9b9,0x3ffe, XPD } }, -{ { 0xdadd,0x5506,0x2a11,0xc567,0x3ffe, XPD } }, -{ { 0x9456,0x6670,0x4cca,0xc12c,0x3ffe, XPD } }, -{ { 0x36bf,0x580c,0xa39f,0xbd08,0x3ffe, XPD } }, -{ { 0x9ee9,0x62fb,0xaf47,0xb8fb,0x3ffe, XPD } }, -{ { 0x6484,0xf9de,0xf333,0xb504,0x3ffe, XPD } }, -{ { 0x2590,0xd2ac,0xf581,0xb123,0x3ffe, XPD } }, -{ { 0x4ac6,0x42a1,0x3eea,0xad58,0x3ffe, XPD } }, -{ { 0x0ef8,0xea7c,0x5ab4,0xa9a1,0x3ffe, XPD } }, -{ { 0x38ea,0xb151,0xd6a9,0xa5fe,0x3ffe, XPD } }, -{ { 0x6819,0x0c49,0x4303,0xa270,0x3ffe, XPD } }, -{ { 0x11ae,0x91a1,0x3260,0x9ef5,0x3ffe, XPD } }, -{ { 0x5539,0xd54e,0x39b9,0x9b8d,0x3ffe, XPD } }, -{ { 0xa96f,0x8db8,0xf051,0x9837,0x3ffe, XPD } }, -{ { 0x0961,0xfef7,0xefa8,0x94f4,0x3ffe, XPD } }, -{ { 0xc336,0xab11,0xd373,0x91c3,0x3ffe, XPD } }, -{ { 0x53c0,0x45cd,0x398b,0x8ea4,0x3ffe, XPD } }, -{ { 0xd6e7,0xea8b,0xc1e3,0x8b95,0x3ffe, XPD } }, -{ { 0x8527,0x92da,0x0e80,0x8898,0x3ffe, XPD } }, -{ { 0x7b15,0xcc48,0xc367,0x85aa,0x3ffe, XPD } }, -{ { 0xa1d7,0xac2b,0x8698,0x82cd,0x3ffe, XPD } }, -{ { 0x0000,0x0000,0x0000,0x8000,0x3ffe, XPD } } -}; -static const uLD B[] = { -{ { 0x0000,0x0000,0x0000,0x0000,0x0000, XPD } }, -{ { 0x1f87,0xdb30,0x18f5,0xf73a,0x3fbd, XPD } }, -{ { 0xac15,0x3e46,0x2932,0xbf4a,0xbfbc, XPD } }, -{ { 0x7944,0xba66,0xa091,0xcb12,0x3fb9, XPD } }, -{ { 0xff78,0x40b4,0x2ee6,0xe69a,0x3fbc, XPD } }, -{ { 0xc895,0x5069,0xe383,0xee53,0xbfbb, XPD } }, -{ { 0x7cde,0x9376,0x4325,0xf8ab,0x3fbc, XPD } }, -{ { 0xa10c,0x25e0,0xc093,0xaefd,0xbfbd, XPD } }, -{ { 0x7d3e,0xea95,0x1366,0xb2fb,0x3fbd, XPD } }, -{ { 0x5d89,0xeb34,0x5191,0x9301,0x3fbd, XPD } }, -{ { 0x80d9,0xb883,0xfb10,0xe5eb,0x3fbb, XPD } }, -{ { 0x045d,0x288c,0xc1ec,0xbedd,0xbfbd, XPD } }, -{ { 0xeded,0x5c85,0x4630,0x8d5a,0x3fbd, XPD } }, -{ { 0x9d82,0xe5ac,0x8e0a,0xfd6d,0x3fba, XPD } }, -{ { 0x6dfd,0xeb58,0xaf14,0x8373,0xbfb9, XPD } }, -{ { 0xf938,0x7aac,0x91cf,0xe8da,0xbfbc, XPD } }, -{ { 0x0000,0x0000,0x0000,0x0000,0x0000, XPD } } -}; -static const uLD R[] = { -{ { 0xa69b,0x530e,0xee1d,0xfd2a,0x3fee, XPD } }, -{ { 0xc746,0x8e7e,0x5960,0xa182,0x3ff2, XPD } }, -{ { 0x63b6,0xadda,0xfd6a,0xaec3,0x3ff5, XPD } }, -{ { 0xc104,0xfd99,0x5b7c,0x9d95,0x3ff8, XPD } }, -{ { 0xe05e,0x249d,0x46b8,0xe358,0x3ffa, XPD } }, -{ { 0x5d1d,0x162c,0xeffc,0xf5fd,0x3ffc, XPD } }, -{ { 0x79aa,0xd1cf,0x17f7,0xb172,0x3ffe, XPD } } -}; - -/* 10 byte sizes versus 12 byte */ -#define douba(k) (A[(k)].ld) -#define doubb(k) (B[(k)].ld) -#define MEXP (NXT*16384.0L) -#ifdef DENORMAL -#define MNEXP (-NXT*(16384.0L+64.0L)) -#else -#define MNEXP (-NXT*16384.0L) -#endif -static const -union -{ - unsigned short L[6]; - long double ld; -} log2ea = {{0xc2ef,0x705f,0xeca5,0xe2a8,0x3ffd, XPD}}; - -#define LOG2EA (log2ea.ld) -/* -#define LOG2EA 0.44269504088896340735992L -*/ -#endif - -#ifdef MIEEE -static long P[] = { -0x3ff40000,0xda6ac6f4,0xa8b7b804, -0x3ffd0000,0xfae158c0,0xcf027de9, -0x3fff0000,0xe00067c9,0x3722405a, -0x3fff0000,0xb33387ca,0x6b43cd99, -}; -static long Q[] = { -/* 0x3fff0000,0x80000000,0x00000000, */ -0x40010000,0xa8003b33,0xa4696307, -0x40020000,0x8666a51c,0x62d7fec2, -0x40010000,0x8666a5d7,0xd072da32, -}; -static long A[] = { -0x3fff0000,0x80000000,0x00000000, -0x3ffe0000,0xfa83b2db,0x722a033a, -0x3ffe0000,0xf5257d15,0x2486cc2c, -0x3ffe0000,0xefe4b99b,0xdcdaf5cb, -0x3ffe0000,0xeac0c6e7,0xdd24392f, -0x3ffe0000,0xe5b906e7,0x7c8348a8, -0x3ffe0000,0xe0ccdeec,0x2a94e111, -0x3ffe0000,0xdbfbb797,0xdaf23755, -0x3ffe0000,0xd744fcca,0xd69d6af4, -0x3ffe0000,0xd2a81d91,0xf12ae45a, -0x3ffe0000,0xce248c15,0x1f8480e4, -0x3ffe0000,0xc9b9bd86,0x6e2f27a3, -0x3ffe0000,0xc5672a11,0x5506dadd, -0x3ffe0000,0xc12c4cca,0x66709456, -0x3ffe0000,0xbd08a39f,0x580c36bf, -0x3ffe0000,0xb8fbaf47,0x62fb9ee9, -0x3ffe0000,0xb504f333,0xf9de6484, -0x3ffe0000,0xb123f581,0xd2ac2590, -0x3ffe0000,0xad583eea,0x42a14ac6, -0x3ffe0000,0xa9a15ab4,0xea7c0ef8, -0x3ffe0000,0xa5fed6a9,0xb15138ea, -0x3ffe0000,0xa2704303,0x0c496819, -0x3ffe0000,0x9ef53260,0x91a111ae, -0x3ffe0000,0x9b8d39b9,0xd54e5539, -0x3ffe0000,0x9837f051,0x8db8a96f, -0x3ffe0000,0x94f4efa8,0xfef70961, -0x3ffe0000,0x91c3d373,0xab11c336, -0x3ffe0000,0x8ea4398b,0x45cd53c0, -0x3ffe0000,0x8b95c1e3,0xea8bd6e7, -0x3ffe0000,0x88980e80,0x92da8527, -0x3ffe0000,0x85aac367,0xcc487b15, -0x3ffe0000,0x82cd8698,0xac2ba1d7, -0x3ffe0000,0x80000000,0x00000000, -}; -static long B[51] = { -0x00000000,0x00000000,0x00000000, -0x3fbd0000,0xf73a18f5,0xdb301f87, -0xbfbc0000,0xbf4a2932,0x3e46ac15, -0x3fb90000,0xcb12a091,0xba667944, -0x3fbc0000,0xe69a2ee6,0x40b4ff78, -0xbfbb0000,0xee53e383,0x5069c895, -0x3fbc0000,0xf8ab4325,0x93767cde, -0xbfbd0000,0xaefdc093,0x25e0a10c, -0x3fbd0000,0xb2fb1366,0xea957d3e, -0x3fbd0000,0x93015191,0xeb345d89, -0x3fbb0000,0xe5ebfb10,0xb88380d9, -0xbfbd0000,0xbeddc1ec,0x288c045d, -0x3fbd0000,0x8d5a4630,0x5c85eded, -0x3fba0000,0xfd6d8e0a,0xe5ac9d82, -0xbfb90000,0x8373af14,0xeb586dfd, -0xbfbc0000,0xe8da91cf,0x7aacf938, -0x00000000,0x00000000,0x00000000, -}; -static long R[] = { -0x3fee0000,0xfd2aee1d,0x530ea69b, -0x3ff20000,0xa1825960,0x8e7ec746, -0x3ff50000,0xaec3fd6a,0xadda63b6, -0x3ff80000,0x9d955b7c,0xfd99c104, -0x3ffa0000,0xe35846b8,0x249de05e, -0x3ffc0000,0xf5fdeffc,0x162c5d1d, -0x3ffe0000,0xb17217f7,0xd1cf79aa, -}; - -#define douba(k) (*(long double *)&A[3*(k)]) -#define doubb(k) (*(long double *)&B[3*(k)]) -#define MEXP (NXT*16384.0L) -#ifdef DENORMAL -#define MNEXP (-NXT*(16384.0L+64.0L)) -#else -#define MNEXP (-NXT*16382.0L) -#endif -static long L[3] = {0x3ffd0000,0xe2a8eca5,0x705fc2ef}; -#define LOG2EA (*(long double *)(&L[0])) -#endif - - -#define F W -#define Fa Wa -#define Fb Wb -#define G W -#define Ga Wa -#define Gb u -#define H W -#define Ha Wb -#define Hb Wb - -#ifndef __MINGW32__ -extern long double MAXNUML; -#endif - -static VOLATILE long double z; -static long double w, W, Wa, Wb, ya, yb, u; - -#ifdef __MINGW32__ -static __inline__ long double reducl( long double ); -extern long double __powil ( long double, int ); -extern long double powl ( long double x, long double y); -#else -#ifdef ANSIPROT -extern long double floorl ( long double ); -extern long double fabsl ( long double ); -extern long double frexpl ( long double, int * ); -extern long double ldexpl ( long double, int ); -extern long double polevll ( long double, void *, int ); -extern long double p1evll ( long double, void *, int ); -extern long double __powil ( long double, int ); -extern int isnanl ( long double ); -extern int isfinitel ( long double ); -static long double reducl( long double ); -extern int signbitl ( long double ); -#else -long double floorl(), fabsl(), frexpl(), ldexpl(); -long double polevll(), p1evll(), __powil(); -static long double reducl(); -int isnanl(), isfinitel(), signbitl(); -#endif /* __MINGW32__ */ - -#ifdef INFINITIES -extern long double INFINITYL; -#else -#define INFINITYL MAXNUML -#endif - -#ifdef NANS -extern long double NANL; -#endif -#ifdef MINUSZERO -extern long double NEGZEROL; -#endif - -#endif /* __MINGW32__ */ - -#ifdef __MINGW32__ - -/* No error checking. We handle Infs and zeros ourselves. */ -static __inline__ long double -__fast_ldexpl (long double x, int expn) -{ - long double res; - __asm__ ("fscale" - : "=t" (res) - : "0" (x), "u" ((long double) expn)); - return res; -} - -#define ldexpl __fast_ldexpl - -#endif - - -long double powl( long double x, long double y ) -{ - /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */ - int i, nflg, iyflg, yoddint; - long e; - - if( y == 0.0L ) - return( 1.0L ); - -#ifdef NANS - if( isnanl(x) ) - { - _SET_ERRNO (EDOM); - return( x ); - } - if( isnanl(y) ) - { - _SET_ERRNO (EDOM); - return( y ); - } -#endif - - if( y == 1.0L ) - return( x ); - - if( isinfl(y) && (x == -1.0L || x == 1.0L) ) - return( y ); - - if( x == 1.0L ) - return( 1.0L ); - - if( y >= MAXNUML ) - { - _SET_ERRNO (ERANGE); -#ifdef INFINITIES - if( x > 1.0L ) - return( INFINITYL ); -#else - if( x > 1.0L ) - return( MAXNUML ); -#endif - if( x > 0.0L && x < 1.0L ) - return( 0.0L ); -#ifdef INFINITIES - if( x < -1.0L ) - return( INFINITYL ); -#else - if( x < -1.0L ) - return( MAXNUML ); -#endif - if( x > -1.0L && x < 0.0L ) - return( 0.0L ); - } - if( y <= -MAXNUML ) - { - _SET_ERRNO (ERANGE); - if( x > 1.0L ) - return( 0.0L ); -#ifdef INFINITIES - if( x > 0.0L && x < 1.0L ) - return( INFINITYL ); -#else - if( x > 0.0L && x < 1.0L ) - return( MAXNUML ); -#endif - if( x < -1.0L ) - return( 0.0L ); -#ifdef INFINITIES - if( x > -1.0L && x < 0.0L ) - return( INFINITYL ); -#else - if( x > -1.0L && x < 0.0L ) - return( MAXNUML ); -#endif - } - if( x >= MAXNUML ) - { -#if INFINITIES - if( y > 0.0L ) - return( INFINITYL ); -#else - if( y > 0.0L ) - return( MAXNUML ); -#endif - return( 0.0L ); - } - - w = floorl(y); - /* Set iyflg to 1 if y is an integer. */ - iyflg = 0; - if( w == y ) - iyflg = 1; - - /* Test for odd integer y. */ - yoddint = 0; - if( iyflg ) - { - ya = fabsl(y); - ya = floorl(0.5L * ya); - yb = 0.5L * fabsl(w); - if( ya != yb ) - yoddint = 1; - } - - if( x <= -MAXNUML ) - { - if( y > 0.0L ) - { -#ifdef INFINITIES - if( yoddint ) - return( -INFINITYL ); - return( INFINITYL ); -#else - if( yoddint ) - return( -MAXNUML ); - return( MAXNUML ); -#endif - } - if( y < 0.0L ) - { -#ifdef MINUSZERO - if( yoddint ) - return( NEGZEROL ); -#endif - return( 0.0 ); - } - } - - - nflg = 0; /* flag = 1 if x<0 raised to integer power */ - if( x <= 0.0L ) - { - if( x == 0.0L ) - { - if( y < 0.0 ) - { -#ifdef MINUSZERO - if( signbitl(x) && yoddint ) - return( -INFINITYL ); -#endif -#ifdef INFINITIES - return( INFINITYL ); -#else - return( MAXNUML ); -#endif - } - if( y > 0.0 ) - { -#ifdef MINUSZERO - if( signbitl(x) && yoddint ) - return( NEGZEROL ); -#endif - return( 0.0 ); - } - if( y == 0.0L ) - return( 1.0L ); /* 0**0 */ - else - return( 0.0L ); /* 0**y */ - } - else - { - if( iyflg == 0 ) - { /* noninteger power of negative number */ - mtherr( fname, DOMAIN ); - _SET_ERRNO (EDOM); -#ifdef NANS - return(NANL); -#else - return(0.0L); -#endif - } - nflg = 1; - } - } - - /* Integer power of an integer. */ - - if( iyflg ) - { - i = w; - w = floorl(x); - if( (w == x) && (fabsl(y) < 32768.0) ) - { - w = __powil( x, (int) y ); - return( w ); - } - } - - - if( nflg ) - x = fabsl(x); - - /* separate significand from exponent */ - x = frexpl( x, &i ); - e = i; - - /* find significand in antilog table A[] */ - i = 1; - if( x <= douba(17) ) - i = 17; - if( x <= douba(i+8) ) - i += 8; - if( x <= douba(i+4) ) - i += 4; - if( x <= douba(i+2) ) - i += 2; - if( x >= douba(1) ) - i = -1; - i += 1; - - - /* Find (x - A[i])/A[i] - * in order to compute log(x/A[i]): - * - * log(x) = log( a x/a ) = log(a) + log(x/a) - * - * log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a - */ - x -= douba(i); - x -= doubb(i/2); - x /= douba(i); - - - /* rational approximation for log(1+v): - * - * log(1+v) = v - v**2/2 + v**3 P(v) / Q(v) - */ - z = x*x; - w = x * ( z * polevll( x, P, 3 ) / p1evll( x, Q, 3 ) ); - w = w - ldexpl( z, -1 ); /* w - 0.5 * z */ - - /* Convert to base 2 logarithm: - * multiply by log2(e) = 1 + LOG2EA - */ - z = LOG2EA * w; - z += w; - z += LOG2EA * x; - z += x; - - /* Compute exponent term of the base 2 logarithm. */ - w = -i; - w = ldexpl( w, -LNXT ); /* divide by NXT */ - w += e; - /* Now base 2 log of x is w + z. */ - - /* Multiply base 2 log by y, in extended precision. */ - - /* separate y into large part ya - * and small part yb less than 1/NXT - */ - ya = reducl(y); - yb = y - ya; - - /* (w+z)(ya+yb) - * = w*ya + w*yb + z*y - */ - F = z * y + w * yb; - Fa = reducl(F); - Fb = F - Fa; - - G = Fa + w * ya; - Ga = reducl(G); - Gb = G - Ga; - - H = Fb + Gb; - Ha = reducl(H); - w = ldexpl( Ga + Ha, LNXT ); - - /* Test the power of 2 for overflow */ - if( w > MEXP ) - { - _SET_ERRNO (ERANGE); - mtherr( fname, OVERFLOW ); - return( INFINITYL ); - } - - if( w < MNEXP ) - { - _SET_ERRNO (ERANGE); - mtherr( fname, UNDERFLOW ); - return( 0.0L ); - } - - e = w; - Hb = H - Ha; - - if( Hb > 0.0L ) - { - e += 1; - Hb -= (1.0L/NXT); /*0.0625L;*/ - } - - /* Now the product y * log2(x) = Hb + e/NXT. - * - * Compute base 2 exponential of Hb, - * where -0.0625 <= Hb <= 0. - */ - z = Hb * polevll( Hb, R, 6 ); /* z = 2**Hb - 1 */ - - /* Express e/NXT as an integer plus a negative number of (1/NXT)ths. - * Find lookup table entry for the fractional power of 2. - */ - if( e < 0 ) - i = 0; - else - i = 1; - i = e/NXT + i; - e = NXT*i - e; - w = douba( e ); - z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */ - z = z + w; - z = ldexpl( z, i ); /* multiply by integer power of 2 */ - - if( nflg ) - { - /* For negative x, - * find out if the integer exponent - * is odd or even. - */ - w = ldexpl( y, -1 ); - w = floorl(w); - w = ldexpl( w, 1 ); - if( w != y ) - z = -z; /* odd exponent */ - } - - return( z ); -} - -static __inline__ long double -__convert_inf_to_maxnum(long double x) -{ - if (isinf(x)) - return (x > 0.0L ? MAXNUML : -MAXNUML); - else - return x; -} - - -/* Find a multiple of 1/NXT that is within 1/NXT of x. */ -static __inline__ long double reducl( long double x ) -{ - long double t; - - /* If the call to ldexpl overflows, set it to MAXNUML. - * This avoids Inf - Inf = Nan result when calculating the 'small' - * part of a reduction. Instead, the small part becomes Inf, - * causing under/overflow when adding it to the 'large' part. - * There must be a cleaner way of doing this. - */ - t = __convert_inf_to_maxnum (ldexpl( x, LNXT )); - t = floorl( t ); - t = ldexpl( t, -LNXT ); - return(t); -} -- 2.11.0