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]
Content-Type: text/plain; charset=ISO-8859-15; format=flowed
Content-Transfer-Encoding: 7bit
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) */