This is the mail archive of the libc-hacker@sources.redhat.com mailing list for the glibc project.
Note that libc-hacker is a closed list. You may look at the archives of this list, but subscription and posting are not open.
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Hi!
According to ISO C99, log2 has to return domain range error if argument is
negative and error may be signalled if argument is 0 (similarly to log or
log10). So I guess we need a log2 wrapper, implemented below:
2001-06-04 Jakub Jelinek <jakub@redhat.com>
* math/Makefile (libm-calls): Add e_log2, w_log2, remove s_log2.
* math/math_private.h (__ieee754_log2, __ieee754_log2f,
__ieee754_log2l): New prototypes.
* sysdeps/generic/w_log2.c: New file.
* sysdeps/generic/w_log2f.c: New file.
* sysdeps/generic/w_log2l.c: New file.
* sysdeps/generic/s_log2l.c: Move...
* sysdeps/generic/e_log2l.c: ...to here. Rename to __ieee754_log2l.
* sysdeps/ieee754/k_standard.c (__kernel_standard): Handle log2(0)
and log2(x < 0).
* sysdeps/i386/fpu/s_log2.S: Move...
* sysdeps/i386/fpu/e_log2.S: ...to here. Rename to __ieee754_log2.
* sysdeps/i386/fpu/s_log2f.S: Move...
* sysdeps/i386/fpu/e_log2f.S: ...to here. Rename to __ieee754_log2f.
* sysdeps/i386/fpu/s_log2l.S: Move...
* sysdeps/i386/fpu/e_log2l.S: ...to here. Rename to __ieee754_log2l.
* sysdeps/m68k/fpu/s_log2.S: Move...
* sysdeps/m68k/fpu/e_log2.S: ...to here. Rename to __ieee754_log2.
* sysdeps/m68k/fpu/s_log2f.S: Move...
* sysdeps/m68k/fpu/e_log2f.S: ...to here. Rename to __ieee754_log2f.
* sysdeps/m68k/fpu/s_log2l.S: Move...
* sysdeps/m68k/fpu/e_log2l.S: ...to here. Rename to __ieee754_log2l.
* sysdeps/ieee754/dbl-64/s_log2.c: Move...
* sysdeps/ieee754/dbl-64/e_log2.c: ...to here. Rename to
__ieee754_log2.
* sysdeps/ieee754/flt-32/s_log2f.c: Move...
* sysdeps/ieee754/flt-32/e_log2f.c: ...to here. Rename to
__ieee754_log2f.
--- libc/math/Makefile.jj Tue Mar 20 13:45:34 2001
+++ libc/math/Makefile Mon Jun 4 10:43:00 2001
@@ -53,11 +53,11 @@ libm-calls = e_acos e_acosh e_asin e_ata
w_tgamma w_hypot w_j0 w_j1 w_jn w_lgamma w_lgamma_r \
w_log w_log10 w_pow w_remainder w_scalb w_sinh w_sqrt \
s_signbit s_fpclassify s_fmax s_fmin s_fdim s_nan s_trunc \
- s_remquo s_log2 e_exp2 s_round s_nearbyint s_sincos \
+ s_remquo e_log2 e_exp2 s_round s_nearbyint s_sincos \
conj cimag creal cabs carg s_cexp s_csinh s_ccosh s_clog \
s_catan s_casin s_ccos s_csin s_ctan s_ctanh s_cacos \
s_casinh s_cacosh s_catanh s_csqrt s_cpow s_cproj s_clog10 \
- s_fma s_lrint s_llrint s_lround s_llround e_exp10
+ s_fma s_lrint s_llrint s_lround s_llround e_exp10 w_log2
dbl-only-routines := branred doasin dosincos halfulp mpa mpatan2 \
mpatan mpexp mplog mpsqrt mptan sincos32 slowexp \
slowpow
--- libc/math/math_private.h.jj Wed May 23 09:21:45 2001
+++ libc/math/math_private.h Mon Jun 4 10:52:07 2001
@@ -169,6 +169,7 @@ extern double __ieee754_gamma_r (double,
extern double __ieee754_lgamma (double);
extern double __ieee754_gamma (double);
extern double __ieee754_log10 (double);
+extern double __ieee754_log2 (double);
extern double __ieee754_sinh (double);
extern double __ieee754_hypot (double,double);
extern double __ieee754_j0 (double);
@@ -211,6 +212,7 @@ extern float __ieee754_gammaf_r (float,i
extern float __ieee754_lgammaf (float);
extern float __ieee754_gammaf (float);
extern float __ieee754_log10f (float);
+extern float __ieee754_log2f (float);
extern float __ieee754_sinhf (float);
extern float __ieee754_hypotf (float,float);
extern float __ieee754_j0f (float);
@@ -250,6 +252,7 @@ extern long double __ieee754_gammal_r (l
extern long double __ieee754_lgammal (long double);
extern long double __ieee754_gammal (long double);
extern long double __ieee754_log10l (long double);
+extern long double __ieee754_log2l (long double);
extern long double __ieee754_sinhl (long double);
extern long double __ieee754_hypotl (long double,long double);
extern long double __ieee754_j0l (long double);
--- libc/sysdeps/generic/w_log2.c.jj Mon Jun 4 10:17:22 2001
+++ libc/sysdeps/generic/w_log2.c Mon Jun 4 11:47:28 2001
@@ -0,0 +1,32 @@
+/*
+ * wrapper log2(X)
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+double
+__log2 (double x) /* wrapper log2 */
+{
+#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, 48); /* log2 (0) */
+ else
+ return __kernel_standard (x, x, 49); /* log2 (x < 0) */
+ }
+ else
+ return z;
+#endif
+}
+weak_alias (__log2, log2)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__log2, __log2l)
+weak_alias (__log2, log2l)
+#endif
--- libc/sysdeps/generic/w_log2f.c.jj Mon Jun 4 10:17:22 2001
+++ libc/sysdeps/generic/w_log2f.c Mon Jun 4 11:47:21 2001
@@ -0,0 +1,30 @@
+/*
+ * wrapper log2(X)
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+float
+__log2f (float x) /* wrapper log2f */
+{
+#ifdef _IEEE_LIBM
+ return __ieee754_log2f (x);
+#else
+ float z;
+ z = __ieee754_log2f (x);
+ if (_LIB_VERSION == _IEEE_ || __isnanf (x)) return z;
+ if (x <= 0.0f)
+ {
+ if (x == 0.0f)
+ /* log2f (0) */
+ return __kernel_standard ((double) x, (double) x, 148);
+ else
+ /* log2f (x < 0) */
+ return __kernel_standard ((double) x, (double) x, 149);
+ }
+ else
+ return z;
+#endif
+}
+weak_alias (__log2f, log2f)
--- libc/sysdeps/generic/w_log2l.c.jj Mon Jun 4 10:17:22 2001
+++ libc/sysdeps/generic/w_log2l.c Mon Jun 4 11:47:40 2001
@@ -0,0 +1,28 @@
+/*
+ * wrapper log2l(X)
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+long double
+__log2l (long double x) /* wrapper log2l */
+{
+#ifdef _IEEE_LIBM
+ return __ieee754_log2l (x);
+#else
+ long double z;
+ z = __ieee754_log2l (x);
+ if (_LIB_VERSION == _IEEE_ || __isnanl (x)) return z;
+ if (x <= 0.0)
+ {
+ if (x == 0.0)
+ return __kernel_standard (x, x, 248); /* log2l (0) */
+ else
+ return __kernel_standard (x, x, 249); /* log2l (x < 0) */
+ }
+ else
+ return z;
+#endif
+}
+weak_alias (__log2l, log2l)
--- libc/sysdeps/generic/s_log2l.c.jj Mon Oct 13 05:52:31 1997
+++ libc/sysdeps/generic/s_log2l.c Mon Jun 4 10:41:10 2001
@@ -1,15 +0,0 @@
-#include <math.h>
-#include <stdio.h>
-#include <errno.h>
-
-long double
-__log2l (long double x)
-{
- fputs ("__log2l not implemented\n", stderr);
- __set_errno (ENOSYS);
- return 0.0;
-}
-weak_alias (__log2l, log2l)
-
-stub_warning (log2l)
-#include <stub-tag.h>
--- libc/sysdeps/generic/e_log2l.c.jj Mon Jun 4 10:41:18 2001
+++ libc/sysdeps/generic/e_log2l.c Mon Jun 4 10:41:42 2001
@@ -0,0 +1,14 @@
+#include <math.h>
+#include <stdio.h>
+#include <errno.h>
+
+long double
+__ieee754_log2l (long double x)
+{
+ fputs ("__ieee754_log2l not implemented\n", stderr);
+ __set_errno (ENOSYS);
+ return 0.0;
+}
+
+stub_warning (log2l)
+#include <stub-tag.h>
--- libc/sysdeps/i386/fpu/s_log2.S.jj Fri Nov 12 16:30:30 1999
+++ libc/sysdeps/i386/fpu/s_log2.S Mon Jun 4 10:34:33 2001
@@ -1,68 +0,0 @@
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
- * Public domain.
- *
- * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
- */
-
-#include <machine/asm.h>
-
-#ifdef __ELF__
- .section .rodata
-#else
- .text
-#endif
- .align ALIGNARG(4)
- ASM_TYPE_DIRECTIVE(one,@object)
-one: .double 1.0
- ASM_SIZE_DIRECTIVE(one)
- /* It is not important that this constant is precise. It is only
- a value which is known to be on the safe side for using the
- fyl2xp1 instruction. */
- ASM_TYPE_DIRECTIVE(limit,@object)
-limit: .double 0.29
- ASM_SIZE_DIRECTIVE(limit)
-
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%edx)
-#else
-#define MO(op) op
-#endif
-
- .text
-ENTRY(__log2)
-#ifdef PIC
- call 1f
-1: popl %edx
- addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
-#endif
- fldl MO(one)
- fldl 4(%esp) // x : 1
- fxam
- fnstsw
- fld %st // x : x : 1
- sahf
- jc 3f // in case x is NaN or ħInf
-4: fsub %st(2), %st // x-1 : x : 1
- fld %st // x-1 : x-1 : x : 1
- fabs // |x-1| : x-1 : x : 1
- fcompl MO(limit) // x-1 : x : 1
- fnstsw // x-1 : x : 1
- andb $0x45, %ah
- jz 2f
- fstp %st(1) // x-1 : 1
- fyl2xp1 // log(x)
- ret
-
-2: fstp %st(0) // x : 1
- fyl2x // log(x)
- ret
-
-3: jp 4b // in case x is ħInf
- fstp %st(1)
- fstp %st(1)
- ret
-END (__log2)
-weak_alias (__log2, log2)
--- libc/sysdeps/i386/fpu/s_log2f.S.jj Fri Nov 12 16:30:30 1999
+++ libc/sysdeps/i386/fpu/s_log2f.S Mon Jun 4 10:34:33 2001
@@ -1,68 +0,0 @@
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
- * Public domain.
- *
- * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
- */
-
-#include <machine/asm.h>
-
-#ifdef __ELF__
- .section .rodata
-#else
- .text
-#endif
- .align ALIGNARG(4)
- ASM_TYPE_DIRECTIVE(one,@object)
-one: .double 1.0
- ASM_SIZE_DIRECTIVE(one)
- /* It is not important that this constant is precise. It is only
- a value which is known to be on the safe side for using the
- fyl2xp1 instruction. */
- ASM_TYPE_DIRECTIVE(limit,@object)
-limit: .double 0.29
- ASM_SIZE_DIRECTIVE(limit)
-
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%edx)
-#else
-#define MO(op) op
-#endif
-
- .text
-ENTRY(__log2f)
-#ifdef PIC
- call 1f
-1: popl %edx
- addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
-#endif
- fldl MO(one)
- flds 4(%esp) // x : 1
- fxam
- fnstsw
- fld %st // x : x : 1
- sahf
- jc 3f // in case x is NaN or ħInf
-4: fsub %st(2), %st // x-1 : x : 1
- fld %st // x-1 : x-1 : x : 1
- fabs // |x-1| : x-1 : x : 1
- fcompl MO(limit) // x-1 : x : 1
- fnstsw // x-1 : x : 1
- andb $0x45, %ah
- jz 2f
- fstp %st(1) // x-1 : 1
- fyl2xp1 // log(x)
- ret
-
-2: fstp %st(0) // x : 1
- fyl2x // log(x)
- ret
-
-3: jp 4b // in case x is ħInf
- fstp %st(1)
- fstp %st(1)
- ret
-END (__log2f)
-weak_alias (__log2f, log2f)
--- libc/sysdeps/i386/fpu/s_log2l.S.jj Fri Nov 12 16:30:30 1999
+++ libc/sysdeps/i386/fpu/s_log2l.S Mon Jun 4 10:34:33 2001
@@ -1,68 +0,0 @@
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
- * Public domain.
- *
- * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
- */
-
-#include <machine/asm.h>
-
-#ifdef __ELF__
- .section .rodata
-#else
- .text
-#endif
- .align ALIGNARG(4)
- ASM_TYPE_DIRECTIVE(one,@object)
-one: .double 1.0
- ASM_SIZE_DIRECTIVE(one)
- /* It is not important that this constant is precise. It is only
- a value which is known to be on the safe side for using the
- fyl2xp1 instruction. */
- ASM_TYPE_DIRECTIVE(limit,@object)
-limit: .double 0.29
- ASM_SIZE_DIRECTIVE(limit)
-
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%edx)
-#else
-#define MO(op) op
-#endif
-
- .text
-ENTRY(__log2l)
-#ifdef PIC
- call 1f
-1: popl %edx
- addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
-#endif
- fldl MO(one)
- fldt 4(%esp) // x : 1
- fxam
- fnstsw
- fld %st // x : x : 1
- sahf
- jc 3f // in case x is NaN or ħInf
-4: fsub %st(2), %st // x-1 : x : 1
- fld %st // x-1 : x-1 : x : 1
- fabs // |x-1| : x-1 : x : 1
- fcompl MO(limit) // x-1 : x : 1
- fnstsw // x-1 : x : 1
- andb $0x45, %ah
- jz 2f
- fstp %st(1) // x-1 : 1
- fyl2xp1 // log(x)
- ret
-
-2: fstp %st(0) // x : 1
- fyl2x // log(x)
- ret
-
-3: jp 4b // in case x is ħInf
- fstp %st(1)
- fstp %st(1)
- ret
-END (__log2l)
-weak_alias (__log2l, log2l)
--- libc/sysdeps/i386/fpu/e_log2f.S.jj Mon Jun 4 10:34:55 2001
+++ libc/sysdeps/i386/fpu/e_log2f.S Mon Jun 4 10:35:10 2001
@@ -0,0 +1,67 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ *
+ * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
+ */
+
+#include <machine/asm.h>
+
+#ifdef __ELF__
+ .section .rodata
+#else
+ .text
+#endif
+ .align ALIGNARG(4)
+ ASM_TYPE_DIRECTIVE(one,@object)
+one: .double 1.0
+ ASM_SIZE_DIRECTIVE(one)
+ /* It is not important that this constant is precise. It is only
+ a value which is known to be on the safe side for using the
+ fyl2xp1 instruction. */
+ ASM_TYPE_DIRECTIVE(limit,@object)
+limit: .double 0.29
+ ASM_SIZE_DIRECTIVE(limit)
+
+
+#ifdef PIC
+#define MO(op) op##@GOTOFF(%edx)
+#else
+#define MO(op) op
+#endif
+
+ .text
+ENTRY(__ieee754_log2f)
+#ifdef PIC
+ call 1f
+1: popl %edx
+ addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
+#endif
+ fldl MO(one)
+ flds 4(%esp) // x : 1
+ fxam
+ fnstsw
+ fld %st // x : x : 1
+ sahf
+ jc 3f // in case x is NaN or ħInf
+4: fsub %st(2), %st // x-1 : x : 1
+ fld %st // x-1 : x-1 : x : 1
+ fabs // |x-1| : x-1 : x : 1
+ fcompl MO(limit) // x-1 : x : 1
+ fnstsw // x-1 : x : 1
+ andb $0x45, %ah
+ jz 2f
+ fstp %st(1) // x-1 : 1
+ fyl2xp1 // log(x)
+ ret
+
+2: fstp %st(0) // x : 1
+ fyl2x // log(x)
+ ret
+
+3: jp 4b // in case x is ħInf
+ fstp %st(1)
+ fstp %st(1)
+ ret
+END (__ieee754_log2f)
--- libc/sysdeps/i386/fpu/e_log2l.S.jj Mon Jun 4 10:34:55 2001
+++ libc/sysdeps/i386/fpu/e_log2l.S Mon Jun 4 10:35:24 2001
@@ -0,0 +1,67 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ *
+ * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
+ */
+
+#include <machine/asm.h>
+
+#ifdef __ELF__
+ .section .rodata
+#else
+ .text
+#endif
+ .align ALIGNARG(4)
+ ASM_TYPE_DIRECTIVE(one,@object)
+one: .double 1.0
+ ASM_SIZE_DIRECTIVE(one)
+ /* It is not important that this constant is precise. It is only
+ a value which is known to be on the safe side for using the
+ fyl2xp1 instruction. */
+ ASM_TYPE_DIRECTIVE(limit,@object)
+limit: .double 0.29
+ ASM_SIZE_DIRECTIVE(limit)
+
+
+#ifdef PIC
+#define MO(op) op##@GOTOFF(%edx)
+#else
+#define MO(op) op
+#endif
+
+ .text
+ENTRY(__ieee754_log2l)
+#ifdef PIC
+ call 1f
+1: popl %edx
+ addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
+#endif
+ fldl MO(one)
+ fldt 4(%esp) // x : 1
+ fxam
+ fnstsw
+ fld %st // x : x : 1
+ sahf
+ jc 3f // in case x is NaN or ħInf
+4: fsub %st(2), %st // x-1 : x : 1
+ fld %st // x-1 : x-1 : x : 1
+ fabs // |x-1| : x-1 : x : 1
+ fcompl MO(limit) // x-1 : x : 1
+ fnstsw // x-1 : x : 1
+ andb $0x45, %ah
+ jz 2f
+ fstp %st(1) // x-1 : 1
+ fyl2xp1 // log(x)
+ ret
+
+2: fstp %st(0) // x : 1
+ fyl2x // log(x)
+ ret
+
+3: jp 4b // in case x is ħInf
+ fstp %st(1)
+ fstp %st(1)
+ ret
+END (__ieee754_log2l)
--- libc/sysdeps/i386/fpu/e_log2.S.jj Mon Jun 4 10:34:55 2001
+++ libc/sysdeps/i386/fpu/e_log2.S Mon Jun 4 10:35:38 2001
@@ -0,0 +1,67 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ *
+ * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
+ */
+
+#include <machine/asm.h>
+
+#ifdef __ELF__
+ .section .rodata
+#else
+ .text
+#endif
+ .align ALIGNARG(4)
+ ASM_TYPE_DIRECTIVE(one,@object)
+one: .double 1.0
+ ASM_SIZE_DIRECTIVE(one)
+ /* It is not important that this constant is precise. It is only
+ a value which is known to be on the safe side for using the
+ fyl2xp1 instruction. */
+ ASM_TYPE_DIRECTIVE(limit,@object)
+limit: .double 0.29
+ ASM_SIZE_DIRECTIVE(limit)
+
+
+#ifdef PIC
+#define MO(op) op##@GOTOFF(%edx)
+#else
+#define MO(op) op
+#endif
+
+ .text
+ENTRY(__ieee754_log2)
+#ifdef PIC
+ call 1f
+1: popl %edx
+ addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
+#endif
+ fldl MO(one)
+ fldl 4(%esp) // x : 1
+ fxam
+ fnstsw
+ fld %st // x : x : 1
+ sahf
+ jc 3f // in case x is NaN or ħInf
+4: fsub %st(2), %st // x-1 : x : 1
+ fld %st // x-1 : x-1 : x : 1
+ fabs // |x-1| : x-1 : x : 1
+ fcompl MO(limit) // x-1 : x : 1
+ fnstsw // x-1 : x : 1
+ andb $0x45, %ah
+ jz 2f
+ fstp %st(1) // x-1 : 1
+ fyl2xp1 // log(x)
+ ret
+
+2: fstp %st(0) // x : 1
+ fyl2x // log(x)
+ ret
+
+3: jp 4b // in case x is ħInf
+ fstp %st(1)
+ fstp %st(1)
+ ret
+END (__ieee754_log2)
--- libc/sysdeps/ieee754/dbl-64/s_log2.c.jj Wed Jul 14 01:52:42 1999
+++ libc/sysdeps/ieee754/dbl-64/s_log2.c Mon Jun 4 10:30:38 2001
@@ -1,136 +0,0 @@
-/* Adapted for log2 by Ulrich Drepper <drepper@cygnus.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.
- * ====================================================
- */
-
-/* __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 "math.h"
-#include "math_private.h"
-
-#ifdef __STDC__
-static const double
-#else
-static double
-#endif
-ln2 = 0.69314718055994530942,
-two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
-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 __log2(double x)
-#else
- double __log2(x)
- double x;
-#endif
-{
- double hfsq,f,s,z,R,w,t1,t2,dk;
- int32_t k,hx,i,j;
- u_int32_t lx;
-
- EXTRACT_WORDS(hx,lx,x);
-
- k=0;
- if (hx < 0x00100000) { /* x < 2**-1022 */
- if (((hx&0x7fffffff)|lx)==0)
- return -two54/(x-x); /* log(+-0)=-inf */
- if (hx<0) return (x-x)/(x-x); /* 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) {
- hfsq=0.5*f*f;
- return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
- } else {
- return dk-((s*(f-R))-f)/ln2;
- }
-}
-
-weak_alias (__log2, log2)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__log2, __log2l)
-weak_alias (__log2, log2l)
-#endif
--- libc/sysdeps/ieee754/dbl-64/e_log2.c.jj Mon Jun 4 10:30:46 2001
+++ libc/sysdeps/ieee754/dbl-64/e_log2.c Mon Jun 4 10:31:08 2001
@@ -0,0 +1,130 @@
+/* Adapted for log2 by Ulrich Drepper <drepper@cygnus.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.
+ * ====================================================
+ */
+
+/* __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 "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const double
+#else
+static double
+#endif
+ln2 = 0.69314718055994530942,
+two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
+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 hfsq,f,s,z,R,w,t1,t2,dk;
+ int32_t k,hx,i,j;
+ u_int32_t lx;
+
+ EXTRACT_WORDS(hx,lx,x);
+
+ k=0;
+ if (hx < 0x00100000) { /* x < 2**-1022 */
+ if (((hx&0x7fffffff)|lx)==0)
+ return -two54/(x-x); /* log(+-0)=-inf */
+ if (hx<0) return (x-x)/(x-x); /* 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) {
+ hfsq=0.5*f*f;
+ return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
+ } else {
+ return dk-((s*(f-R))-f)/ln2;
+ }
+}
--- libc/sysdeps/ieee754/flt-32/s_log2f.c.jj Wed Jul 14 02:01:54 1999
+++ libc/sysdeps/ieee754/flt-32/s_log2f.c Mon Jun 4 10:32:15 2001
@@ -1,91 +0,0 @@
-/* e_logf.c -- float version of e_log.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- * adapted for log2 by Ulrich Drepper <drepper@cygnus.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.
- * ====================================================
- */
-
-
-#include "math.h"
-#include "math_private.h"
-
-#ifdef __STDC__
-static const float
-#else
-static float
-#endif
-ln2 = 0.69314718055994530942,
-two25 = 3.355443200e+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 __log2f(float x)
-#else
- float __log2f(x)
- float x;
-#endif
-{
- float hfsq,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/(x-x); /* log(+-0)=-inf */
- if (ix<0) return (x-x)/(x-x); /* 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) {
- hfsq=(float)0.5*f*f;
- return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
- } else {
- return dk-((s*(f-R))-f)/ln2;
- }
-}
-weak_alias (__log2f, log2f)
--- libc/sysdeps/ieee754/flt-32/e_log2f.c.jj Mon Jun 4 10:32:21 2001
+++ libc/sysdeps/ieee754/flt-32/e_log2f.c Mon Jun 4 10:32:43 2001
@@ -0,0 +1,90 @@
+/* e_logf.c -- float version of e_log.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * adapted for log2 by Ulrich Drepper <drepper@cygnus.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.
+ * ====================================================
+ */
+
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+ln2 = 0.69314718055994530942,
+two25 = 3.355443200e+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 hfsq,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/(x-x); /* log(+-0)=-inf */
+ if (ix<0) return (x-x)/(x-x); /* 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) {
+ hfsq=(float)0.5*f*f;
+ return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
+ } else {
+ return dk-((s*(f-R))-f)/ln2;
+ }
+}
--- libc/sysdeps/ieee754/k_standard.c.jj Wed Jul 14 01:44:31 1999
+++ libc/sysdeps/ieee754/k_standard.c Mon Jun 4 10:15:50 2001
@@ -88,6 +88,8 @@ static double zero = 0.0; /* used as con
* 45-- exp2 underflow
* 46-- exp10 overflow
* 47-- exp10 underflow
+ * 48-- log2(0)
+ * 49-- log2(x<0)
*/
@@ -937,7 +939,42 @@ static double zero = 0.0; /* used as con
__set_errno (ERANGE);
}
break;
- /* #### Last used is 47/147/247 ### */
+ case 48:
+ case 148:
+ case 248:
+ /* log2(0) */
+ exc.type = SING;
+ exc.name = type < 100 ? "log2" : (type < 200
+ ? "log2f" : "log2l");
+ if (_LIB_VERSION == _SVID_)
+ exc.retval = -HUGE;
+ else
+ exc.retval = -HUGE_VAL;
+ if (_LIB_VERSION == _POSIX_)
+ __set_errno (ERANGE);
+ else if (!matherr(&exc)) {
+ __set_errno (EDOM);
+ }
+ break;
+ case 49:
+ case 149:
+ case 249:
+ /* log2(x<0) */
+ exc.type = DOMAIN;
+ exc.name = type < 100 ? "log2" : (type < 200
+ ? "log2f" : "log2l");
+ if (_LIB_VERSION == _SVID_)
+ exc.retval = -HUGE;
+ else
+ exc.retval = NAN;
+ if (_LIB_VERSION == _POSIX_)
+ __set_errno (EDOM);
+ else if (!matherr(&exc)) {
+ __set_errno (EDOM);
+ }
+ break;
+
+ /* #### Last used is 49/149/249 ### */
}
return exc.retval;
}
--- libc/sysdeps/m68k/fpu/s_log2.c.jj Tue Mar 25 02:33:29 1997
+++ libc/sysdeps/m68k/fpu/s_log2.c Mon Jun 4 10:39:50 2001
@@ -1,2 +0,0 @@
-#define FUNC log2
-#include <s_atan.c>
--- libc/sysdeps/m68k/fpu/s_log2f.c.jj Tue Mar 25 02:33:29 1997
+++ libc/sysdeps/m68k/fpu/s_log2f.c Mon Jun 4 10:39:50 2001
@@ -1,2 +0,0 @@
-#define FUNC log2f
-#include <s_atanf.c>
--- libc/sysdeps/m68k/fpu/s_log2l.c.jj Tue Mar 25 02:33:30 1997
+++ libc/sysdeps/m68k/fpu/s_log2l.c Mon Jun 4 10:39:50 2001
@@ -1,2 +0,0 @@
-#define FUNC log2l
-#include <s_atanl.c>
--- libc/sysdeps/m68k/fpu/e_log2.c.jj Mon Jun 4 10:40:09 2001
+++ libc/sysdeps/m68k/fpu/e_log2.c Mon Jun 4 10:40:22 2001
@@ -0,0 +1,2 @@
+#define FUNC __ieee754_log2
+#include <e_acos.c>
--- libc/sysdeps/m68k/fpu/e_log2f.c.jj Mon Jun 4 10:40:09 2001
+++ libc/sysdeps/m68k/fpu/e_log2f.c Mon Jun 4 10:40:35 2001
@@ -0,0 +1,2 @@
+#define FUNC __ieee754_log10f
+#include <e_acosf.c>
--- libc/sysdeps/m68k/fpu/e_log2l.c.jj Mon Jun 4 10:40:09 2001
+++ libc/sysdeps/m68k/fpu/e_log2l.c Mon Jun 4 10:40:50 2001
@@ -0,0 +1,2 @@
+#define FUNC __ieee754_log10l
+#include <e_acosl.c>
Jakub
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |