1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 /* 22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak logf = __logf 30 31 /* 32 * Algorithm: 33 * 34 * Let y = x rounded to six significant bits. Then for any choice 35 * of e and z such that y = 2^e z, we have 36 * 37 * log(x) = e log(2) + log(z) + log(1+(x-y)/y) 38 * 39 * Note that (x-y)/y = (x'-y')/y' for any scaled x' = sx, y' = sy; 40 * in particular, we can take s to be the power of two that makes 41 * ulp(x') = 1. 42 * 43 * From a table, obtain l = log(z) and r = 1/y'. For |s| <= 2^-6, 44 * approximate log(1+s) by a polynomial p(s) where p(s) := s+s*s* 45 * (K1+s*(K2+s*K3)). Then we compute the expression above as 46 * e*ln2 + l + p(r*(x'-y')) all evaluated in double precision. 47 * 48 * When x is subnormal, we first scale it to the normal range, 49 * adjusting e accordingly. 50 * 51 * Accuracy: 52 * 53 * The largest error is less than 0.6 ulps. 54 */ 55 56 #include "libm.h" 57 58 /* 59 * For i = 0, ..., 12, 60 * TBL[2i] = log(1 + i/32) and TBL[2i+1] = 2^-23 / (1 + i/32) 61 * 62 * For i = 13, ..., 32, 63 * TBL[2i] = log(1/2 + i/64) and TBL[2i+1] = 2^-23 / (1 + i/32) 64 */ 65 static const double TBL[] = { 66 0.000000000000000000e+00, 1.192092895507812500e-07, 67 3.077165866675368733e-02, 1.155968868371212153e-07, 68 6.062462181643483994e-02, 1.121969784007352926e-07, 69 8.961215868968713805e-02, 1.089913504464285680e-07, 70 1.177830356563834557e-01, 1.059638129340277719e-07, 71 1.451820098444978890e-01, 1.030999260979729787e-07, 72 1.718502569266592284e-01, 1.003867701480263102e-07, 73 1.978257433299198675e-01, 9.781275040064102225e-08, 74 2.231435513142097649e-01, 9.536743164062500529e-08, 75 2.478361639045812692e-01, 9.304139672256097884e-08, 76 2.719337154836417580e-01, 9.082612537202380448e-08, 77 2.954642128938358980e-01, 8.871388989825581272e-08, 78 3.184537311185345887e-01, 8.669766512784091150e-08, 79 -3.522205935893520934e-01, 8.477105034722222546e-08, 80 -3.302416868705768671e-01, 8.292820142663043248e-08, 81 -3.087354816496132859e-01, 8.116377160904255122e-08, 82 -2.876820724517809014e-01, 7.947285970052082892e-08, 83 -2.670627852490452536e-01, 7.785096460459183052e-08, 84 -2.468600779315257843e-01, 7.629394531250000159e-08, 85 -2.270574506353460753e-01, 7.479798560049019504e-08, 86 -2.076393647782444896e-01, 7.335956280048077330e-08, 87 -1.885911698075500298e-01, 7.197542010613207272e-08, 88 -1.698990367953974734e-01, 7.064254195601851460e-08, 89 -1.515498981272009327e-01, 6.935813210227272390e-08, 90 -1.335313926245226268e-01, 6.811959402901785336e-08, 91 -1.158318155251217008e-01, 6.692451343201754014e-08, 92 -9.844007281325252434e-02, 6.577064251077586116e-08, 93 -8.134563945395240081e-02, 6.465588585805084723e-08, 94 -6.453852113757117814e-02, 6.357828776041666578e-08, 95 -4.800921918636060631e-02, 6.253602074795082293e-08, 96 -3.174869831458029812e-02, 6.152737525201612732e-08, 97 -1.574835696813916761e-02, 6.055075024801586965e-08, 98 0.000000000000000000e+00, 5.960464477539062500e-08, 99 }; 100 101 static const double C[] = { 102 6.931471805599452862e-01, 103 -2.49887584306188944706e-01, 104 3.33368809981254554946e-01, 105 -5.00000008402474976565e-01 106 }; 107 108 #define ln2 C[0] 109 #define K3 C[1] 110 #define K2 C[2] 111 #define K1 C[3] 112 113 float 114 logf(float x) 115 { 116 double v, t; 117 float f; 118 int hx, ix, i, exp, iy; 119 120 hx = *(int *)&x; 121 ix = hx & ~0x80000000; 122 123 if (ix >= 0x7f800000) /* nan or inf */ 124 return ((hx < 0)? x * 0.0f : x * x); 125 126 exp = 0; 127 if (hx < 0x00800000) { /* negative, zero, or subnormal */ 128 if (hx <= 0) { 129 f = 0.0f; 130 return ((ix == 0)? -1.0f / f : f / f); 131 } 132 133 /* subnormal; scale by 2^149 */ 134 f = (float)ix; 135 ix = *(int *)&f; 136 exp = -149; 137 } 138 139 exp += (ix - 0x3f320000) >> 23; 140 ix &= 0x007fffff; 141 iy = (ix + 0x20000) & 0xfffc0000; 142 i = iy >> 17; 143 t = ln2 * (double)exp + TBL[i]; 144 v = (double)(ix - iy) * TBL[i + 1]; 145 v += (v * v) * (K1 + v * (K2 + v * K3)); 146 f = (float)(t + v); 147 return (f); 148 } 149