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 76