e_logf.c revision 241755
10Sduke/* e_logf.c -- float version of e_log.c. 211884Sykantser * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 30Sduke */ 40Sduke 50Sduke/* 60Sduke * ==================================================== 70Sduke * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 80Sduke * 90Sduke * Developed at SunPro, a Sun Microsystems, Inc. business. 100Sduke * Permission to use, copy, modify, and distribute this 110Sduke * software is freely granted, provided that this notice 120Sduke * is preserved. 130Sduke * ==================================================== 140Sduke */ 150Sduke 160Sduke#include <sys/cdefs.h> 170Sduke__FBSDID("$FreeBSD: head/lib/msun/src/e_logf.c 241755 2012-10-19 22:46:48Z imp $"); 180Sduke 192362Sohair#include "math.h" 202362Sohair#include "math_private.h" 212362Sohair 220Sduke/* __ieee754_log(x) 230Sduke * Return the logrithm of x 240Sduke * 250Sduke * Method : 260Sduke * 1. Argument Reduction: find k and f such that 270Sduke * x = 2^k * (1+f), 280Sduke * where sqrt(2)/2 < 1+f < sqrt(2) . 2916958Siignatyev * 300Sduke * 2. Approximation of log(1+f). 310Sduke * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 320Sduke * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 330Sduke * = 2s + s*R 340Sduke * We use a special Reme algorithm on [0,0.1716] to generate 350Sduke * a polynomial of degree 8 to approximate R The maximum error 360Sduke * of this polynomial approximation is bounded by 2**-34.24. In 370Sduke * other words, 380Sduke * 2 4 6 8 390Sduke * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s 400Sduke * (the values of Lg1 to Lg7 are listed in the program) 410Sduke * and 420Sduke * | 2 8 | -34.24 430Sduke * | Lg1*s +...+Lg4*s - R(z) | <= 2 440Sduke * | | 450Sduke * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. 460Sduke * In order to guarantee error in log below 1ulp, we compute log 470Sduke * by 480Sduke * log(1+f) = f - s*(f - R) (if f is not too large) 490Sduke * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 500Sduke * 510Sduke * 3. Finally, log(x) = k*ln2 + log(1+f). 520Sduke * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 530Sduke * Here ln2 is split into two floating point number: 540Sduke * ln2_hi + ln2_lo, 550Sduke * where n*ln2_hi is always exact for |n| < 2000. 560Sduke * 570Sduke * Special cases: 580Sduke * log(x) is NaN with signal if x < 0 (including -INF) ; 590Sduke * log(+INF) is +INF; log(0) is -INF with signal; 600Sduke * log(NaN) is that NaN with no signal. 610Sduke * 620Sduke * Accuracy: 630Sduke * according to an error analysis, the error is always less than 640Sduke * 1 ulp (unit in the last place). 650Sduke * 660Sduke * Constants: 670Sduke * The hexadecimal values are the intended ones for the following 680Sduke * constants. The decimal values may be used, provided that the 690Sduke * compiler will convert from decimal to binary accurately enough 700Sduke * to produce the hexadecimal values shown. 710Sduke */ 720Sduke 730Sdukestatic const float 740Sdukeln2_hi = 6.9313812256e-01, /* 0x3f317180 */ 750Sdukeln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ 760Sduketwo25 = 3.355443200e+07, /* 0x4c000000 */ 770Sduke/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ 780SdukeLg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ 790SdukeLg2 = 0xccce13.0p-25, /* 0.40000972152 */ 800SdukeLg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ 810SdukeLg4 = 0xf89e26.0p-26; /* 0.24279078841 */ 820Sduke 830Sdukestatic const float zero = 0.0; 840Sduke 850Sdukefloat 860Sduke__ieee754_logf(float x) 870Sduke{ 880Sduke float hfsq,f,s,z,R,w,t1,t2,dk; 890Sduke int32_t k,ix,i,j; 900Sduke 910Sduke GET_FLOAT_WORD(ix,x); 920Sduke 930Sduke k=0; 940Sduke if (ix < 0x00800000) { /* x < 2**-126 */ 950Sduke if ((ix&0x7fffffff)==0) 960Sduke return -two25/zero; /* log(+-0)=-inf */ 970Sduke if (ix<0) return (x-x)/zero; /* log(-#) = NaN */ 980Sduke k -= 25; x *= two25; /* subnormal number, scale up x */ 99 GET_FLOAT_WORD(ix,x); 100 } 101 if (ix >= 0x7f800000) return x+x; 102 k += (ix>>23)-127; 103 ix &= 0x007fffff; 104 i = (ix+(0x95f64<<3))&0x800000; 105 SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */ 106 k += (i>>23); 107 f = x-(float)1.0; 108 if((0x007fffff&(0x8000+ix))<0xc000) { /* -2**-9 <= f < 2**-9 */ 109 if(f==zero) { 110 if(k==0) { 111 return zero; 112 } else { 113 dk=(float)k; 114 return dk*ln2_hi+dk*ln2_lo; 115 } 116 } 117 R = f*f*((float)0.5-(float)0.33333333333333333*f); 118 if(k==0) return f-R; else {dk=(float)k; 119 return dk*ln2_hi-((R-dk*ln2_lo)-f);} 120 } 121 s = f/((float)2.0+f); 122 dk = (float)k; 123 z = s*s; 124 i = ix-(0x6147a<<3); 125 w = z*z; 126 j = (0x6b851<<3)-ix; 127 t1= w*(Lg2+w*Lg4); 128 t2= z*(Lg1+w*Lg3); 129 i |= j; 130 R = t2+t1; 131 if(i>0) { 132 hfsq=(float)0.5*f*f; 133 if(k==0) return f-(hfsq-s*(hfsq+R)); else 134 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); 135 } else { 136 if(k==0) return f-s*(f-R); else 137 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); 138 } 139} 140