Home | History | Annotate | Line # | Download | only in i387 
      1 /*
      2  * Written by J.T. Conklin <jtc (at) NetBSD.org>.
      3  * Public domain.
      4  */
      5 
      6 /*
      7  * Modified by Lex Wennmacher <wennmach (at) NetBSD.org>
      8  * Still public domain.
      9  */
     10 
     11 #include <machine/asm.h>
     12 
     13 #include "abi.h"
     14 
     15 RCSID("$NetBSD: s_log1p.S,v 1.14 2024/07/16 14:52:49 riastradh Exp $")
     16 
     17 /*
     18  * The log1p() function is provided to compute an accurate value of
     19  * log(1 + x), even for tiny values of x. The i387 FPU provides the
     20  * fyl2xp1 instruction for this purpose. However, the range of this
     21  * instruction is limited to:
     22  * 		-(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1
     23  *                         -0.292893 <= x <= 0.414214
     24  * at least on older processor versions.
     25  *
     26  * log1p() is implemented by testing the range of the argument.
     27  * If it is appropriate for fyl2xp1, this instruction is used.
     28  * Else, we compute log1p(x) = ln(2)*ld(1 + x) the traditional way
     29  * (using fyl2x).
     30  *
     31  * The range testing costs speed, but as the rationale for the very
     32  * existence of this function is accuracy, we accept that.
     33  *
     34  * In order to reduce the cost for testing the range, we check if
     35  * the argument is in the range
     36  *                             -0.25 <= x <= 0.25
     37  * which can be done with just one conditional branch. If x is
     38  * inside this range, we use fyl2xp1. Outside of this range,
     39  * the use of fyl2x is accurate enough.
     40  */
     41 
     42 WEAK_ALIAS(log1p, _log1p)
     43 
     44 .text
     45 	.align	4
     46 ENTRY(_log1p)
     47 	XMM_ONE_ARG_DOUBLE_PROLOGUE
     48 	fldl	ARG_DOUBLE_ONE
     49 	fabs
     50 	fld1				/* ... x 1 */
     51 	fadd	%st(0)			/* ... x 2 */
     52 	fadd	%st(0)			/* ... x 4 */
     53 	fld1				/* ... 4 1 */
     54 	fdivp				/* ... x 0.25 */
     55 	fcompp
     56 	fnstsw	%ax
     57 	andb	$69,%ah
     58 	jne	use_fyl2x
     59 	jmp	use_fyl2xp1
     60 
     61 	.align	4
     62 use_fyl2x:
     63 	fldln2
     64         fldl	ARG_DOUBLE_ONE
     65         fld1
     66         faddp
     67         fyl2x
     68 	XMM_DOUBLE_EPILOGUE
     69         ret
     70 
     71 	.align	4
     72 use_fyl2xp1:
     73 	fldln2
     74 	fldl	ARG_DOUBLE_ONE
     75 	fyl2xp1
     76 	XMM_DOUBLE_EPILOGUE
     77 	ret
     78 END(_log1p)
     79