Mail Archives: djgpp-workers/2015/07/18/13:19:37
| X-Authentication-Warning: | delorie.com: mail set sender to djgpp-workers-bounces using -f
|
| X-Recipient: | djgpp-workers AT delorie DOT com
|
| Message-ID: | <55AA8AFB.7030502@gmx.de>
|
| Date: | Sat, 18 Jul 2015 19:20:59 +0200
|
| From: | "Juan Manuel Guerrero (juan DOT guerrero AT gmx DOT de) [via djgpp-workers AT delorie DOT com]" <djgpp-workers AT delorie DOT com>
|
| User-Agent: | Mozilla/5.0 (X11; U; Linux x86_64; de; rv:1.9.2.13) Gecko/20101206 SUSE/3.1.7 Thunderbird/3.1.7
|
| MIME-Version: | 1.0
|
| To: | djgpp-workers AT delorie DOT com
|
| Subject: | Implementation of log2[f]
|
| X-Provags-ID: | V03:K0:s1Y7HgVrSDHzbJ5RUH3N/hCsK9rzyfa32fGW5XMNfiWfgIhH/mt
|
| O0HuQb6VRmSm9wNSeermBjw7aNmGCnzbT/BwfFC2TAOTGPRO4ZdnYMXVG36H43jB1Pt9mTf
|
| N3KIY671MAZIEJ16KNRriVeCuf+3f87cuYVb1RMVE5weF5OQA9sZzEgSty4gm8k3QWwWyp0
|
| OyyApcCKeGJukDT2TyjEg==
|
| X-UI-Out-Filterresults: | notjunk:1;V01:K0:Wytzh2dxFYM=:NohVtx625/H9DBHihA7noj
|
| JKInphWtywndALOvcfS+iqZzPanqzsgaHyQgCqOLqPxcJKqtV/2ZMTwv+0Jq/BWkD9EhYTmn0
|
| KRS5hAm82WNNLwrNvHgundefDH+VPRdc68E2v4cHcy0VuB6/XkYcngVg3xF+G9Md7cdWXUd4a
|
| kUARvgWfpLd9T6D96Pe/k7wV0UnEoExfnJ12TSfDStwZz+3gtI7Hz6lQ5prJ1hTvYmehtGY50
|
| s/9vJk8ZBp03+MuhKQ2jLWUfpLzehqUC4AbOBwXPmV6wREZeQrqUeSKlrJ5jxOrA76Diu7k8c
|
| etL4tcutHMnfSATHtWDrCt5+ajqbrkxEKfzAsYBOMdqK3dlQ9bPXlM/oQe8ENvD5ApmqUVDFi
|
| WFuwJ1KwFfpH8z9Fv0oSfYMGhhnkU17jVqW+tFnuqEYJ5OULvgqxDUaUW1b7uxCy4tfFa6831
|
| 3cHIoV/RsFjxbCkq168ZJT/ASyS0PgxRtQv9eRwRCD9XM4KOPjkomLAZnrk1utYPTvRacKdk2
|
| DjJdxJF7nGMY7o+fOIl02FkXDrtOf6firOUFdriHJVj4pC6f+FaBcuL9xSZu17qYjNdQUjUNu
|
| vpUD/NiCn/192ycYAo2A3Mb/tvrPmy2CW+WRbAd/6SDeD0fffyvtke3DpT4cCBLdQsYKMA1ES
|
| KWhdEebx0PHjc5Innzy2g+gOG2qg50leQrWSN8F0jYQOXJhj9ae9U1rdt4Z9TzukNnYg=
|
| Reply-To: | djgpp-workers AT delorie DOT com
|
OFYI, I have commited the patch below into the HEAD and v2_05_1 branches.
Regards,
Juan M. Guerrero
Index: djgpp/include/math.h
===================================================================
RCS file: /cvs/djgpp/djgpp/include/math.h,v
retrieving revision 1.26
diff -p -U 5 -r1.26 math.h
--- djgpp/include/math.h 23 May 2015 13:10:46 -0000 1.26
+++ djgpp/include/math.h 18 Jul 2015 15:33:41 -0000
@@ -1,5 +1,6 @@
+/* Copyright (C) 2015 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 2013 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 2003 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1999 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1998 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */
@@ -33,10 +34,11 @@ double floor(double _x);
double fmod(double _x, double _y);
double frexp(double _x, int *_pexp);
double ldexp(double _x, int _exp);
double log(double _y);
double log10(double _x);
+double log2(double _x);
double modf(double _x, double *_pint);
double pow(double _x, double _y);
double sin(double _x);
double sinh(double _x);
double sqrt(double _x);
@@ -72,10 +74,11 @@ extern float tanhf(float);
extern float expf(float);
extern float frexpf(float, int *);
extern float ldexpf(float, int);
extern float logf(float);
extern float log10f(float);
+extern float log2f(float);
extern float modff(float, float *);
extern float powf(float, float);
extern float sqrtf(float);
extern float ceilf(float);
extern float fabsf(float);
@@ -172,11 +175,10 @@ double atanh(double _x);
double cbrt(double _x);
double exp2(double _x);
double expm1(double _x);
double hypot(double _x, double _y);
double log1p(double _x);
-double log2(double _x);
/* These are in libm.a (Cygnus). You must link -lm to get these */
/* See libm/math.h for comments */
extern double erf(double);
Index: djgpp/include/libm/math.h
===================================================================
RCS file: /cvs/djgpp/djgpp/include/libm/math.h,v
retrieving revision 1.11
diff -p -U 5 -r1.11 math.h
--- djgpp/include/libm/math.h 16 Nov 2013 21:48:24 -0000 1.11
+++ djgpp/include/libm/math.h 18 Jul 2015 15:33:41 -0000
@@ -1,5 +1,6 @@
+/* Copyright (C) 2015 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 2013 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1998 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */
/* Provided by Cygnus Support (jtc AT cygnus DOT com) */
@@ -270,10 +271,11 @@ extern float tanhf __P((float));
extern float expf __P((float));
extern float frexpf __P((float, int *));
extern float ldexpf __P((float, int));
extern float logf __P((float));
extern float log10f __P((float));
+extern float log2f __P((float));
extern float modff __P((float, float *));
extern float powf __P((float, float));
extern float sqrtf __P((float));
@@ -296,11 +298,10 @@ extern float jnf __P((int, float));
extern float lgammaf __P((float));
extern float nanf __P((void));
extern float y0f __P((float));
extern float y1f __P((float));
extern float ynf __P((int, float));
-#define log2f(x) (logf (x) / (float) M_LOG2_E)
extern float acoshf __P((float));
extern float asinhf __P((float));
extern float atanhf __P((float));
extern float cbrtf __P((float));
Index: djgpp/src/docs/kb/wc205.txi
===================================================================
RCS file: /cvs/djgpp/djgpp/src/docs/kb/wc205.txi,v
retrieving revision 1.12
diff -p -U 5 -r1.12 wc205.txi
--- djgpp/src/docs/kb/wc205.txi 27 Jun 2015 07:56:10 -0000 1.12
+++ djgpp/src/docs/kb/wc205.txi 18 Jul 2015 15:33:42 -0000
@@ -403,5 +403,11 @@ files will be overwritten nor the user w
name.
@cindex @file{stdbool.h} added
@cindex @acronym{C99} compliance, @file{stdbool.h}
The @acronym{C99} file @file{stdbool.h} was added.
+
+@cindex @file{math.h}, added @code{log2f} and @code{log2}
+@findex log2f AT r{, added @code{log2f} to @file{libm.a}}
+@findex log2 AT r{, added @code{log2} to @file{libm.a}}
+The base 2 logarithm function for arguments of type @code{double} and @code{float}
+have been added to the mathematical library @file{libm.a}.
Index: djgpp/src/libm/math/ef_log.c
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/ef_log.c,v
retrieving revision 1.1
diff -p -U 5 -r1.1 ef_log.c
--- djgpp/src/libm/math/ef_log.c 6 Jul 1998 16:03:26 -0000 1.1
+++ djgpp/src/libm/math/ef_log.c 18 Jul 2015 15:33:43 -0000
@@ -6,11 +6,11 @@
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "fdlibm.h"
@@ -70,18 +70,18 @@ static float zero = 0.0;
}
R = f*f*((float)0.5-(float)0.33333333333333333*f);
if(k==0) return f-R; else {dk=(float)k;
return dk*ln2_hi-((R-dk*ln2_lo)-f);}
}
- s = f/((float)2.0+f);
+ s = f/((float)2.0+f);
dk = (float)k;
z = s*s;
i = ix-(0x6147a<<3);
w = z*z;
j = (0x6b851<<3)-ix;
- t1= w*(Lg2+w*(Lg4+w*Lg6));
- t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+ t1= w*(Lg2+w*(Lg4+w*Lg6));
+ t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
i |= j;
R = t2+t1;
if(i>0) {
hfsq=(float)0.5*f*f;
if(k==0) return f-(hfsq-s*(hfsq+R)); else
Index: djgpp/src/libm/math/fdlibm.h
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/fdlibm.h,v
retrieving revision 1.1
diff -p -U 5 -r1.1 fdlibm.h
--- djgpp/src/libm/math/fdlibm.h 4 Oct 1998 10:48:42 -0000 1.1
+++ djgpp/src/libm/math/fdlibm.h 18 Jul 2015 15:33:43 -0000
@@ -65,10 +65,11 @@ extern double __ieee754_cosh __P((double
extern double __ieee754_fmod __P((double,double));
extern double __ieee754_pow __P((double,double));
extern double __ieee754_lgamma_r __P((double,int *));
extern double __ieee754_gamma_r __P((double,int *));
extern double __ieee754_log10 __P((double));
+extern double __ieee754_log2 __P((double));
extern double __ieee754_sinh __P((double));
extern double __ieee754_hypot __P((double,double));
extern double __ieee754_j0 __P((double));
extern double __ieee754_j1 __P((double));
extern double __ieee754_y0 __P((double));
@@ -112,10 +113,11 @@ extern float __ieee754_coshf __P((float)
extern float __ieee754_fmodf __P((float,float));
extern float __ieee754_powf __P((float,float));
extern float __ieee754_lgammaf_r __P((float,int *));
extern float __ieee754_gammaf_r __P((float,int *));
extern float __ieee754_log10f __P((float));
+extern float __ieee754_log2f __P((float));
extern float __ieee754_sinhf __P((float));
extern float __ieee754_hypotf __P((float,float));
extern float __ieee754_j0f __P((float));
extern float __ieee754_j1f __P((float));
extern float __ieee754_y0f __P((float));
Index: djgpp/src/libm/math/makefile
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/makefile,v
retrieving revision 1.12
diff -p -U 5 -r1.12 makefile
--- djgpp/src/libm/math/makefile 11 May 2015 07:40:40 -0000 1.12
+++ djgpp/src/libm/math/makefile 18 Jul 2015 15:33:44 -0000
@@ -21,10 +21,11 @@ SRC += e_j0.c
SRC += e_j1.c
SRC += e_jn.c
SRC += er_lgamma.c
SRC += e_log.c
SRC += e_log10.c
+SRC += e_log2.c
SRC += e_pow.c
SRC += e_rem_pio2.c
SRC += e_remainder.c
SRC += e_scalb.c
SRC += e_sinh.c
@@ -45,10 +46,11 @@ SRC += w_j1.c
SRC += w_jn.c
SRC += w_lgamma.c
SRC += wr_lgamma.c
SRC += w_log.c
SRC += w_log10.c
+SRC += w_log2.c
SRC += w_pow.c
SRC += w_remainder.c
SRC += w_scalb.c
SRC += w_sinh.c
SRC += w_sqrt.c
@@ -103,10 +105,11 @@ SRC += ef_j0.c
SRC += ef_j1.c
SRC += ef_jn.c
SRC += erf_lgamma.c
SRC += ef_log.c
SRC += ef_log10.c
+SRC += ef_log2.c
SRC += ef_pow.c
SRC += ef_rem_pio2.c
SRC += ef_remainder.c
SRC += ef_scalb.c
SRC += ef_sinh.c
@@ -127,10 +130,11 @@ SRC += wf_j1.c
SRC += wf_jn.c
SRC += wf_lgamma.c
SRC += wrf_lgamma.c
SRC += wf_log.c
SRC += wf_log10.c
+SRC += wf_log2.c
SRC += wf_pow.c
SRC += wf_remainder.c
SRC += wf_scalb.c
SRC += wf_sinh.c
SRC += wf_sqrt.c
@@ -190,13 +194,13 @@ SRC += finitel.c
chobj = w_acos.def w_acosh.def w_asin.def s_asinh.def \
s_atan.def w_atan2.def w_atanh.def w_j0.def \
s_copysign.def w_cosh.def s_erf.def w_exp.def \
s_fabs.def s_floor.def w_fmod.def s_frexp.def \
w_gamma.def w_hypot.def s_ldexp.def w_log.def \
- w_log10.def s_log1p.def s_matherr.def s_modf.def \
- w_pow.def w_remainder.def s_sin.def w_sinh.def \
- s_cbrt.def w_sqrt.def s_tan.def s_tanh.def \
+ w_log10.def w_log2.def s_log1p.def s_matherr.def \
+ s_modf.def w_pow.def w_remainder.def s_sin.def \
+ w_sinh.def s_cbrt.def w_sqrt.def s_tan.def s_tanh.def \
s_infinity.def s_isnan.def s_scalbn.def s_nextafter.def \
s_nan.def s_ilogb.def s_expm1.def trunc.def \
lrint.def llrint.def round.def lround.def llround.def
CFLAGS = -D_USE_LIBM_MATH_H
Index: djgpp/src/libm/math/math.texi
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/math.texi,v
retrieving revision 1.5
diff -p -U 5 -r1.5 math.texi
--- djgpp/src/libm/math/math.texi 27 Oct 2013 20:48:01 -0000 1.5
+++ djgpp/src/libm/math/math.texi 18 Jul 2015 15:33:44 -0000
@@ -55,10 +55,11 @@ available when you include @file{math.h}
* ldexp:: Load exponent
* llrint:: Round to integer
* llround:: Round to nearest integer
* log:: Natural logarithms
* log10:: Base 10 logarithms
+* log2:: Base 2 logarithms
* log1p:: Log of 1 + X
* lrint:: Round to integer
* lround:: Round to nearest integer
* matherr:: Modifiable math error handler
* modf:: Split fractional and integer parts
@@ -189,10 +190,13 @@ The library is set to X/Open mode by def
@page
@include w_log10.def
@page
+@include w_log2.def
+
+@page
@include s_log1p.def
@page
@include s_matherr.def
diff -aprNU5 djgpp.orig/src/libm/math/e_log2.c djgpp/src/libm/math/e_log2.c
--- djgpp.orig/src/libm/math/e_log2.c 1969-12-31 18:57:24 -0502
+++ djgpp/src/libm/math/e_log2.c 2015-07-12 13:45:52 -0502
@@ -0,0 +1,134 @@
+
+/* @(#)e_log10.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_log2(x)
+ * Return the logarithm to base 2 of x
+ *
+ * Method :
+ * 1. Argument Reduction: find k and f such that
+ * x = 2^k * (1+f),
+ * where sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ * 2. Approximation of log(1+f).
+ * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ * = 2s + s*R
+ * We use a special Reme algorithm on [0,0.1716] to generate
+ * a polynomial of degree 14 to approximate R The maximum error
+ * of this polynomial approximation is bounded by 2**-58.45. In
+ * other words,
+ * 2 4 6 8 10 12 14
+ * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
+ * (the values of Lg1 to Lg7 are listed in the program)
+ * and
+ * | 2 14 | -58.45
+ * | Lg1*s +...+Lg7*s - R(z) | <= 2
+ * | |
+ * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ * In order to guarantee error in log below 1ulp, we compute log
+ * by
+ * log(1+f) = f - s*(f - R) (if f is not too large)
+ * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
+ *
+ * 3. Finally, log(x) = k + log(1+f).
+ * = k+(f-(hfsq-(s*(hfsq+R))))
+ *
+ * Special cases:
+ * log2(x) is NaN with signal if x < 0 (including -INF) ;
+ * log2(+INF) is +INF; log(0) is -INF with signal;
+ * log2(NaN) is that NaN with no signal.
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+static const double
+#else
+static double
+#endif
+ln2 = 0.69314718055994530942,
+two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
+Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
+Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
+Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
+Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
+Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
+Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+#ifdef __STDC__
+static const double zero = 0.0;
+#else
+static double zero = 0.0;
+#endif
+
+#ifdef __STDC__
+ double __ieee754_log2(double x)
+#else
+ double __ieee754_log2(x)
+ double x;
+#endif
+{
+ double dk,f,R,s,t1,t2,w,z;
+ __int32_t i,j,k,hx;
+ __uint32_t lx;
+
+ EXTRACT_WORDS(hx,lx,x);
+
+ k=0;
+ if (hx < 0x00100000) { /* x < 2**-1022 */
+ if (((hx&0x7fffffff)|lx)==0)
+ return -two54/zero; /* log(+-0)=-inf */
+ if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
+ k -= 54; x *= two54; /* subnormal number, scale up x */
+ GET_HIGH_WORD(hx,x);
+ }
+ if (hx >= 0x7ff00000) return x+x;
+ k += (hx>>20)-1023;
+ hx &= 0x000fffff;
+ i = (hx + 0x95f64) & 0x100000;
+ SET_HIGH_WORD (x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
+ k += (i >> 20);
+ dk = (double) k;
+ f = x - 1.0;
+ if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */
+ if (f == zero) return dk;
+ R = f * f * (0.5 - 0.33333333333333333 * f);
+ return dk - (R - f) / ln2;
+ }
+ s = f / (2.0 + f);
+ z = s * s;
+ i = hx - 0x6147a;
+ w = z * z;
+ j = 0x6b851 - hx;
+ t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+ t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+ i |= j;
+ R = t2 + t1;
+ if (i > 0) {
+ double hfsq = 0.5 * f * f;
+ return dk - ((hfsq - (s * (hfsq + R))) - f) / ln2;
+ }
+ else
+ return dk - ((s * (f - R)) - f) / ln2;
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
diff -aprNU5 djgpp.orig/src/libm/math/ef_log2.c djgpp/src/libm/math/ef_log2.c
--- djgpp.orig/src/libm/math/ef_log2.c 1969-12-31 18:57:24 -0502
+++ djgpp/src/libm/math/ef_log2.c 2015-07-12 13:20:16 -0502
@@ -0,0 +1,87 @@
+/* ef_log.c -- float version of e_log.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian AT cygnus DOT com.
+ * adapted for log2 by Juan M. Guerrero <juan DOT guerrero AT gmx DOT net>
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+ln2 = 0.69314718055994530942,
+two25 = 3.3554432000e+07, /* 0x4c000000 */
+Lg1 = 6.6666668653e-01, /* 3F2AAAAB */
+Lg2 = 4.0000000596e-01, /* 3ECCCCCD */
+Lg3 = 2.8571429849e-01, /* 3E924925 */
+Lg4 = 2.2222198546e-01, /* 3E638E29 */
+Lg5 = 1.8183572590e-01, /* 3E3A3325 */
+Lg6 = 1.5313838422e-01, /* 3E1CD04F */
+Lg7 = 1.4798198640e-01; /* 3E178897 */
+
+#ifdef __STDC__
+static const float zero = 0.0;
+#else
+static float zero = 0.0;
+#endif
+
+#ifdef __STDC__
+ float __ieee754_log2f(float x)
+#else
+ float __ieee754_log2f(x)
+ float x;
+#endif
+{
+ float f,s,z,R,w,t1,t2,dk;
+ __int32_t k,ix,i,j;
+
+ GET_FLOAT_WORD(ix,x);
+
+ k=0;
+ if (ix < 0x00800000) { /* x < 2**-126 */
+ if ((ix&0x7fffffff)==0)
+ return -two25/zero; /* log(+-0)=-inf */
+ if (ix<0) return (x-x)/zero; /* log(-#) = NaN */
+ k -= 25; x *= two25; /* subnormal number, scale up x */
+ GET_FLOAT_WORD(ix,x);
+ }
+ if (ix >= 0x7f800000) return x+x;
+ k += (ix>>23)-127;
+ ix &= 0x007fffff;
+ i = (ix+(0x95f64<<3))&0x800000;
+ SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */
+ k += (i>>23);
+ dk = (float)k;
+ f = x-(float)1.0;
+ if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */
+ if(f==zero) return dk;
+ R = f*f*((float)0.5-(float)0.33333333333333333*f);
+ return dk-(R-f)/ln2;
+ }
+ s = f/((float)2.0+f);
+ z = s*s;
+ i = ix-(0x6147a<<3);
+ w = z*z;
+ j = (0x6b851<<3)-ix;
+ t1= w*(Lg2+w*(Lg4+w*Lg6));
+ t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+ i |= j;
+ R = t2+t1;
+ if(i>0) {
+ float hfsq=(float)0.5*f*f;
+ return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
+ } else
+ return dk-((s*(f-R))-f)/ln2;
+}
diff -aprNU5 djgpp.orig/src/libm/math/wf_log2.c djgpp/src/libm/math/wf_log2.c
--- djgpp.orig/src/libm/math/wf_log2.c 1969-12-31 18:57:24 -0502
+++ djgpp/src/libm/math/wf_log2.c 2015-07-12 14:19:50 -0502
@@ -0,0 +1,64 @@
+/* wf_log2.c -- float version of w_log2.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian AT cygnus DOT com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * wrapper log2f(X)
+ */
+
+#include "fdlibm.h"
+#include <libc/ieee.h>
+
+#ifdef __STDC__
+ float log2f(float x) /* wrapper log2f */
+#else
+ float log2f(x) /* wrapper log2f */
+ float x;
+#endif
+{
+#ifdef _IEEE_LIBM
+ return __ieee754_log2f(x);
+#else
+ _float_long_union ux;
+ float z;
+
+ ux.f = x;
+
+ z = __ieee754_log2f(x);
+ if(_LIB_VERSION == _IEEE_ || isnanf(ux.l)) return z;
+ if(x<=(float)0.0) {
+ if(x==(float)0.0)
+ /* log2(0) */
+ return (float)__kernel_standard((double)x,(double)x,127);
+ else
+ /* log2(x<0) */
+ return (float)__kernel_standard((double)x,(double)x,128);
+ } else
+ return z;
+#endif
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+ double log2(double x)
+#else
+ double log2(x)
+ double x;
+#endif
+{
+ return (double) log2f((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
diff -aprNU5 djgpp.orig/src/libm/math/w_log2.c djgpp/src/libm/math/w_log2.c
--- djgpp.orig/src/libm/math/w_log2.c 1969-12-31 18:57:24 -0502
+++ djgpp/src/libm/math/w_log2.c 2015-07-12 13:16:50 -0502
@@ -0,0 +1,83 @@
+
+/* @(#)w_log10.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+ <<log2>>, <<log2f>>---base 2 logarithms
+
+INDEX
+log2
+INDEX
+log2f
+
+ANSI_SYNOPSIS
+ #include <math.h>
+ double log2(double <[x]>);
+ float log2f(float <[x]>);
+
+TRAD_SYNOPSIS
+ #include <math.h>
+ double log2(<[x]>)
+ double <[x]>;
+
+ float log2f(<[x]>)
+ float <[x]>;
+
+DESCRIPTION
+<<log2>> returns the base 2 logarithm of <[x]>.
+It is implemented as <<log(<[x]>) / log(2)>>.
+
+<<log2f>> is identical, save that it takes and returns <<float>> values.
+
+RETURNS
+<<log2>> and <<log2f>> return the calculated value.
+
+See the description of <<log>> for information on errors.
+
+PORTABILITY
+<<log2>> is ANSI C. <<log2f>> is an extension.
+
+ */
+
+/*
+ * wrapper log2(X)
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+ double log2(double x) /* wrapper log2 */
+#else
+ double log2(x) /* wrapper log2 */
+ double x;
+#endif
+{
+#ifdef _IEEE_LIBM
+ return __ieee754_log2(x);
+#else
+ double z;
+ z = __ieee754_log2(x);
+ if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
+ if(x<=0.0) {
+ if(x==0.0)
+ return __kernel_standard(x,x,27); /* log2(0) */
+ else
+ return __kernel_standard(x,x,28); /* log2(x<0) */
+ } else
+ return z;
+#endif
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
- Raw text -