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