1/*
2 * Written by J.T. Conklin <jtc@NetBSD.org>.
3 * Public domain.
4 */
5
6/*
7 * Modified by Lex Wennmacher <wennmach@NetBSD.org>
8 * Still public domain.
9 */
10
11#include <machine/asm.h>
12
13#include "abi.h"
14
15RCSID("$NetBSD: s_log1pf.S,v 1.11 2024/07/16 14:52:49 riastradh Exp $")
16
17/*
18 * The log1pf() 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 * log1pf() is implemented by testing the range of the argument.
27 * If it is appropriate for fyl2xp1, this instruction is used.
28 * Else, we compute log1pf(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
42WEAK_ALIAS(log1pf, _log1pf)
43
44.text
45          .align    4
46ENTRY(_log1pf)
47          XMM_ONE_ARG_FLOAT_PROLOGUE
48          flds      ARG_FLOAT_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
62use_fyl2x:
63          fldln2
64        flds        ARG_FLOAT_ONE
65        fld1
66        faddp
67        fyl2x
68          XMM_FLOAT_EPILOGUE
69        ret
70
71          .align    4
72use_fyl2xp1:
73          fldln2
74          flds      ARG_FLOAT_ONE
75          fyl2xp1
76          XMM_FLOAT_EPILOGUE
77          ret
78END(_log1pf)
79