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