contrib/openbsd_libm/arch/amd64/s_log1p.S

   1 /*      $OpenBSD: s_log1p.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 log1p() 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  * log1p() is implemented by testing the range of the argument.
  26  * If it is appropriate for fyl2xp1, this instruction is used.
  27  * Else, we compute log1p(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
  44 ENTRY(log1p)
  45         XMM_ONE_ARG_DOUBLE_PROLOGUE
  46         fldl    ARG_DOUBLE_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
  60 use_fyl2x:
  61         fldln2
  62         fldl    ARG_DOUBLE_ONE
  63         fld1
  64         faddp
  65         fyl2x
  66         XMM_DOUBLE_EPILOGUE
  67         ret
  68
  69         .align  4
  70 use_fyl2xp1:
  71         fldln2
  72         fldl    ARG_DOUBLE_ONE
  73         fyl2xp1
  74         XMM_DOUBLE_EPILOGUE
  75         ret