1/* $OpenBSD: s_log1p.S,v 1.6 2017/08/19 18:27:19 deraadt 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 "DEFS.h" 13 14/* 15 * The log1p() function is provided to compute an accurate value of 16 * log(1 + x), even for tiny values of x. The i387 FPU provides the 17 * fyl2xp1 instruction for this purpose. However, the range of this 18 * instruction is limited to: 19 * -(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1 20 * -0.292893 <= x <= 0.414214 21 * at least on older processor versions. 22 * 23 * log1p() is implemented by testing the range of the argument. 24 * If it is appropriate for fyl2xp1, this instruction is used. 25 * Else, we compute log1p(x) = ln(2)*ld(1 + x) the traditional way 26 * (using fyl2x). 27 * 28 * The range testing costs speed, but as the rationale for the very 29 * existence of this function is accuracy, we accept that. 30 * 31 * In order to reduce the cost for testing the range, we check if 32 * the argument is in the range 33 * -0.25 <= x <= 0.25 34 * which can be done with just one conditional branch. If x is 35 * inside this range, we use fyl2xp1. Outside of this range, 36 * the use of fyl2x is accurate enough. 37 * 38 */ 39 40ENTRY(log1p) 41 fldl 4(%esp) 42 fabs 43 fld1 /* ... x 1 */ 44 fadd %st(0) /* ... x 2 */ 45 fadd %st(0) /* ... x 4 */ 46 fld1 /* ... 4 1 */ 47 fdivp /* ... x 0.25 */ 48 fcompp 49 fnstsw %ax 50 andb $69,%ah 51 jne use_fyl2x 52 jmp use_fyl2xp1 53 54 .align 4,0xcc 55use_fyl2x: 56 fldln2 57 fldl 4(%esp) 58 fld1 59 faddp 60 fyl2x 61 ret 62 63 .align 4,0xcc 64use_fyl2xp1: 65 fldln2 66 fldl 4(%esp) 67 fyl2xp1 68 ret 69END_STD(log1p) 70