lib/libm/arch/i386/s_log1pf.S

   1 /*
   2  * Written by J.T. Conklin <jtc@NetBSD.org>.
   3  * Public domain.
   4  *
   5  * $NetBSD: s_log1pf.S,v 1.10 2003/09/16 18:17:11 wennmach Exp $
   6  * $DragonFly: src/lib/libm/arch/i386/s_log1pf.S,v 1.1 2005/07/26 21:15:19 joerg Exp $
   7  */
   8
   9 /*
  10  * Modified by Lex Wennmacher <wennmach@NetBSD.org>
  11  * Still public domain.
  12  */
  13
  14 #include <machine/asm.h>
  15
  16 #include "abi.h"
  17
  18 /*
  19  * The log1pf() function is provided to compute an accurate value of
  20  * log(1 + x), even for tiny values of x. The i387 FPU provides the
  21  * fyl2xp1 instruction for this purpose. However, the range of this
  22  * instruction is limited to:
  23  *              -(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1
  24  *                         -0.292893 <= x <= 0.414214
  25  * at least on older processor versions.
  26  *
  27  * log1pf() is implemented by testing the range of the argument.
  28  * If it is appropriate for fyl2xp1, this instruction is used.
  29  * Else, we compute log1pf(x) = ln(2)*ld(1 + x) the traditional way
  30  * (using fyl2x).
  31  *
  32  * The range testing costs speed, but as the rationale for the very
  33  * existence of this function is accuracy, we accept that.
  34  *
  35  * In order to reduce the cost for testing the range, we check if
  36  * the argument is in the range
  37  *                             -0.25 <= x <= 0.25
  38  * which can be done with just one conditional branch. If x is
  39  * inside this range, we use fyl2xp1. Outside of this range,
  40  * the use of fyl2x is accurate enough.
  41  *
  42  */
  43
  44 .text
  45         .align  4
  46 ENTRY(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
  62 use_fyl2x:
  63         fldln2
  64         flds    ARG_FLOAT_ONE
  65         fld1
  66         faddp
  67         fyl2x
  68         XMM_FLOAT_EPILOGUE
  69         ret
  70
  71         .align  4
  72 use_fyl2xp1:
  73         fldln2
  74         flds    ARG_FLOAT_ONE
  75         fyl2xp1
  76         XMM_FLOAT_EPILOGUE
  77         ret