Print this page
5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>


  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 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak log2 = __log2
  31 
  32 /* INDENT OFF */
  33 /*
  34  * log2(x) = log(x)/log2
  35  *
  36  * Base on Table look-up algorithm with product polynomial
  37  * approximation for log(x).
  38  *
  39  * By K.C. Ng, Nov 29, 2004
  40  *
  41  * (a). For x in [1-0.125, 1+0.125], from log.c we have
  42  *      log(x) =  f + ((a1*f^2) *
  43  *                 ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
  44  *                 (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
  45  *                 ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
  46  *      where f = x - 1.
  47  *      (i) modify a1 <- a1 / log2
  48  *      (ii) 1/log2 = 1.4426950408889634...
  49  *                  = 1.5 - 0.057304959... (4 bit shift)
  50  *           Let lv = 1.5 - 1/log2, then


  71  *      z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
  72  *      significant bits. Then
  73  *          log2(x) = n + log2(z).
  74  *
  75  * Special cases:
  76  *      log2(x) is NaN with signal if x < 0 (including -INF) ;
  77  *      log2(+INF) is +INF; log2(0) is -INF with signal;
  78  *      log2(NaN) is that NaN with no signal.
  79  *
  80  * Maximum error observed: less than 0.84 ulp
  81  *
  82  * Constants:
  83  * The hexadecimal values are the intended ones for the following constants.
  84  * The decimal values may be used, provided that the compiler will convert
  85  * from decimal to binary accurately enough to produce the hexadecimal values
  86  * shown.
  87  */
  88 /* INDENT ON */
  89 
  90 #include "libm.h"
  91 #include "libm_synonyms.h"
  92 #include "libm_protos.h"
  93 
  94 extern const double _TBL_log[];
  95 
  96 static const double P[] = {
  97 /* ONE   */  1.0,
  98 /* TWO52 */  4503599627370496.0,
  99 /* LN10V */  1.4426950408889634073599246810018920433347,   /* 1/log10 */
 100 /* ZERO  */  0.0,
 101 /* A1    */ -9.6809362455249638217841932228967194640116e-02,
 102 /* A2    */  1.99628461483039965074226529395673424005508422852e+0000,
 103 /* A3    */  2.26812367662950720159642514772713184356689453125e+0000,
 104 /* A4    */ -9.05030639084976384900471657601883634924888610840e-0001,
 105 /* A5    */ -1.48275767132434044270894446526654064655303955078e+0000,
 106 /* A6    */  1.88158320939722756293122074566781520843505859375e+0000,
 107 /* A7    */  1.83309386046986411145098827546462416648864746094e+0000,
 108 /* A8    */  1.24847063988317086291601754055591300129890441895e+0000,
 109 /* A9    */  1.98372421445537705508854742220137268304824829102e+0000,
 110 /* A10   */ -3.94711735767898475035764249696512706577777862549e-0001,
 111 /* A11   */  3.07890395362954372160402272129431366920471191406e+0000,




  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 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak __log2 = log2
  31 
  32 /* INDENT OFF */
  33 /*
  34  * log2(x) = log(x)/log2
  35  *
  36  * Base on Table look-up algorithm with product polynomial
  37  * approximation for log(x).
  38  *
  39  * By K.C. Ng, Nov 29, 2004
  40  *
  41  * (a). For x in [1-0.125, 1+0.125], from log.c we have
  42  *      log(x) =  f + ((a1*f^2) *
  43  *                 ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
  44  *                 (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
  45  *                 ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
  46  *      where f = x - 1.
  47  *      (i) modify a1 <- a1 / log2
  48  *      (ii) 1/log2 = 1.4426950408889634...
  49  *                  = 1.5 - 0.057304959... (4 bit shift)
  50  *           Let lv = 1.5 - 1/log2, then


  71  *      z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
  72  *      significant bits. Then
  73  *          log2(x) = n + log2(z).
  74  *
  75  * Special cases:
  76  *      log2(x) is NaN with signal if x < 0 (including -INF) ;
  77  *      log2(+INF) is +INF; log2(0) is -INF with signal;
  78  *      log2(NaN) is that NaN with no signal.
  79  *
  80  * Maximum error observed: less than 0.84 ulp
  81  *
  82  * Constants:
  83  * The hexadecimal values are the intended ones for the following constants.
  84  * The decimal values may be used, provided that the compiler will convert
  85  * from decimal to binary accurately enough to produce the hexadecimal values
  86  * shown.
  87  */
  88 /* INDENT ON */
  89 
  90 #include "libm.h"

  91 #include "libm_protos.h"
  92 
  93 extern const double _TBL_log[];
  94 
  95 static const double P[] = {
  96 /* ONE   */  1.0,
  97 /* TWO52 */  4503599627370496.0,
  98 /* LN10V */  1.4426950408889634073599246810018920433347,   /* 1/log10 */
  99 /* ZERO  */  0.0,
 100 /* A1    */ -9.6809362455249638217841932228967194640116e-02,
 101 /* A2    */  1.99628461483039965074226529395673424005508422852e+0000,
 102 /* A3    */  2.26812367662950720159642514772713184356689453125e+0000,
 103 /* A4    */ -9.05030639084976384900471657601883634924888610840e-0001,
 104 /* A5    */ -1.48275767132434044270894446526654064655303955078e+0000,
 105 /* A6    */  1.88158320939722756293122074566781520843505859375e+0000,
 106 /* A7    */  1.83309386046986411145098827546462416648864746094e+0000,
 107 /* A8    */  1.24847063988317086291601754055591300129890441895e+0000,
 108 /* A9    */  1.98372421445537705508854742220137268304824829102e+0000,
 109 /* A10   */ -3.94711735767898475035764249696512706577777862549e-0001,
 110 /* A11   */  3.07890395362954372160402272129431366920471191406e+0000,