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_log1p.S,v 1.13 2003/09/16 18:17:11 wennmach Exp $") 16 17/* 18 * The log1p() 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 * log1p() is implemented by testing the range of the argument. 27 * If it is appropriate for fyl2xp1, this instruction is used. 28 * Else, we compute log1p(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 */ 42 43.text 44 .align 4 45ENTRY(log1p) 46 XMM_ONE_ARG_DOUBLE_PROLOGUE 47 fldl ARG_DOUBLE_ONE 48 fabs 49 fld1 /* ... x 1 */ 50 fadd %st(0) /* ... x 2 */ 51 fadd %st(0) /* ... x 4 */ 52 fld1 /* ... 4 1 */ 53 fdivp /* ... x 0.25 */ 54 fcompp 55 fnstsw %ax 56 andb $69,%ah 57 jne use_fyl2x 58 jmp use_fyl2xp1 59 60 .align 4 61use_fyl2x: 62 fldln2 63 fldl ARG_DOUBLE_ONE 64 fld1 65 faddp 66 fyl2x 67 XMM_DOUBLE_EPILOGUE 68 ret 69 70 .align 4 71use_fyl2xp1: 72 fldln2 73 fldl ARG_DOUBLE_ONE 74 fyl2xp1 75 XMM_DOUBLE_EPILOGUE 76 ret 77