>Number:         20399
>Category:       port-i386
>Synopsis:       libm387: inaccuracies in s_log1p.S and s_log1pf.S
>Confidential:   no
>Severity:       non-critical
>Priority:       medium
>Responsible:    port-i386-maintainer
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Mon Feb 17 15:44:00 PST 2003
>Closed-Date:
>Last-Modified:
>Originator:     
>Release:        NetBSD 1.6
>Organization:
	
>Environment:
	
	
System: NetBSD m3.ods.org 1.6 NetBSD 1.6 (Monostatos) #0: Tue Feb 11 17:09:36 GMT 2003 jtb@m3.ods.org:/usr/src/sys/arch/i386/compile/Monostatos i386
Architecture: i386
Machine: i386
>Description:
        libm387 gives bad results in s_log1p.S and s_log1pf.S for small
        values.  The fyl2xp1 instruction needs to be used for accuracy.

>How-To-Repeat:
        Run the tests for pkgsrc/math/gsl.
>Fix:
        I made the following changes:

--- abi.h~	Mon Feb 17 23:20:29 2003
+++ abi.h	Mon Feb 17 23:09:34 2003
@@ -8,7 +8,7 @@
  * The x86-64 ABI specifies that float, double and long double
  * arguments are passed in SSE2 (xmm) registers. Unfortunately,
  * there is no way to push those on to the FP stack, which is
- * where he fancier instructions get their arguments from.
+ * where the fancier instructions get their arguments from.
  *
  * Define some prologues and epilogues to store and retrieve
  * xmm regs to local variables.
--- s_log1p.S~	Mon Feb 17 23:20:59 2003
+++ s_log1p.S	Mon Feb 17 23:07:14 2003
@@ -10,17 +10,30 @@
 RCSID("$NetBSD: s_log1p.S,v 1.9 2001/06/19 00:26:30 fvdl Exp $")
 
 /*
- * Since the fyl2xp1 instruction has such a limited range:
- *	-(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1
- * it's not worth trying to use it.
+ * The fyl2xp1 instruction has a very limited range:
+ *	-(1 - (sqrt(2) / 2)) <= x <= (1 - (sqrt(2) / 2))
+ * but is much more accurate than fyl2x for values near zero.
  */
 
+	.section	.rodata
+	.align	4
+X:	.double	0.293            /* approximately 1 - sqrt(2)/2 */
+
 ENTRY(log1p)
-	XMM_ONE_ARG_DOUBLE_PROLOGUE
+        XMM_ONE_ARG_DOUBLE_PROLOGUE
 	fldln2
 	fldl	ARG_DOUBLE_ONE
+	fld	%st(0)
+	fabs
+	fcompl	X
+	fnstsw	%ax
+	sahf
+	jb	1f
 	fld1
 	faddp
 	fyl2x
-	XMM_DOUBLE_EPILOGUE
+       	XMM_DOUBLE_EPILOGUE
+	ret
+1:	fyl2xp1
+       	XMM_DOUBLE_EPILOGUE
 	ret
--- s_log1pf.S~	Mon Feb 17 23:21:03 2003
+++ s_log1pf.S	Mon Feb 17 23:07:12 2003
@@ -10,17 +10,30 @@
 RCSID("$NetBSD: s_log1pf.S,v 1.6 2001/06/19 00:26:30 fvdl Exp $")
 
 /*
- * Since the fyl2xp1 instruction has such a limited range:
- *	-(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1
- * it's not worth trying to use it.
+ * The fyl2xp1 instruction has a very limited range:
+ *	-(1 - (sqrt(2) / 2)) <= x <= (1 - (sqrt(2) / 2))
+ * but is much more accurate than fyl2x for values near zero.
  */
 
+	.section	.rodata
+	.align	4
+X:	.float	0.293
+
 ENTRY(log1pf)
-	XMM_ONE_ARG_FLOAT_PROLOGUE
+       	XMM_ONE_ARG_FLOAT_PROLOGUE
 	fldln2
 	flds	ARG_FLOAT_ONE
+	fld	%st(0)
+	fabs
+	fcompl	X
+	fnstsw	%ax
+	sahf
+	jb	1f
 	fld1
 	faddp
 	fyl2x
-	XMM_FLOAT_EPILOGUE
+	ret
+        XMM_FLOAT_EPILOGUE
+1:	fyl2xp1
+        XMM_FLOAT_EPILOGUE
 	ret

>Release-Note:
>Audit-Trail:
>Unformatted: